Hello, I ran into a seemingly non-consistent behavior of the function
[DMPlexComputeCellGeometryFEM](https://urldefense.us/v3/__https://petsc.org/release/manualpages/DMPlex/DMPlexComputeCellGeometryFEM/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_rJABo0_$ ) ( [DM](https://urldefense.us/v3/__https://petsc.org/release/manualpages/DM/DM/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_kyYWLFi$ ) dm , [PetscInt](https://urldefense.us/v3/__https://petsc.org/release/manualpages/Sys/PetscInt/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_jAVrwAd$ ) cell , [PetscQuadrature](https://urldefense.us/v3/__https://petsc.org/release/manualpages/DT/PetscQuadrature/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_jVjnRvX$ ) quad , [PetscReal](https://urldefense.us/v3/__https://petsc.org/release/manualpages/Sys/PetscReal/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_heQHCRx$ ) * v , [PetscReal](https://urldefense.us/v3/__https://petsc.org/release/manualpages/Sys/PetscReal/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_heQHCRx$ ) J [], [PetscReal](https://urldefense.us/v3/__https://petsc.org/release/manualpages/Sys/PetscReal/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_heQHCRx$ ) invJ [], [PetscReal](https://urldefense.us/v3/__https://petsc.org/release/manualpages/Sys/PetscReal/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_heQHCRx$ ) * detJ ) Calling it in a loop happens to crash, only some times, when the input npoints is not equal to one. Attached is a mesh file (msh format), consisting of four quadrilaterals. Small examples in both C and Fortran also attached; they read the mesh, create a Gauss tensor quadrature rule and then loop through all cells in order to get the mapped quadrature points. When the quadrature created has a single point per dimension (npoints = 1), there seems to be no issue. The program runs and the output is the expected. On the C side, For npoints = 2 the program hangs, usually no errors; sometimes *** stack smashing detected ***: terminated, [:1345101] *** Process received signal *** [:1345101] Signal: Aborted (6) [:1345101] Signal code: (-6) or LeakSanitizer: detected memory leaks (when compiled with -fsanitize=addres). On the Fortran side, similar behavior for npoints = 1 or 2. Either no error at all (no error meaning the one that only mentions "probably memory access out of bounds"), or sometimes "corrupted size vs. prev_size", followed by a stack trace. Question: Do the arguments J, invJ need to be allocated (e.g. allocate(J(4))) before calling DMPlexComputeCellGeometryFEM? Allocation seems to be the cause for triggering an error. --- NOTE: The interface to DMPlexComputeCellGeometryFEM, under $PETSC_ARCH$/ftn/dm/petscdmplex.h90, read: interface DMPlexComputeCellGeometryFEM subroutine DMPlexComputeCellGeometryFEM(a,b,c,d,e,f,g, z) import tDM,tPetscQuadrature DM :: a PetscInt :: b PetscQuadrature :: c PetscReal :: d #### should it be d(*) ? PetscReal :: e(*) PetscReal :: f(*) PetscReal :: g PetscErrorCode z end subroutine end interface Should argument d (referring to the array of the mapped integration points v) be defined differently? --- NOTE: Tested with PETSc 3.23.1 (from the tarball) and the latest git. Configured with ./configure --with-fortran-bindings=1 Using gcc/gfortran 11.4.0, OpenMPI 4.1.2 Thanks, Noam
mesh_quad.msh
Description: Mesh model
program test
#include <petsc/finclude/petscdm.h>
#include <petsc/finclude/petscdmplex.h>
use PETScDMplex
implicit none
DM :: dm
PetscInt :: c
PetscInt :: plexdim, Nc, npoints
PetscInt :: cell, cStart, cEnd
PetscInt :: ierr
PetscReal :: vv(8), detJ ! change size of v accordingly
PetscQuadrature :: quad
PetscReal, dimension(:), pointer :: J, invJ
PetscCallA(PetscInitialize(ierr))
PetscCallA(DMPlexCreateFromFile(PETSC_COMM_WORLD, "mesh_quad.msh", "test", PETSC_TRUE, dm, ierr))
plexdim = 2
Nc = 2
npoints = 2 ! npoints = 1 works all the time
PetscCallA(PetscDTGaussTensorQuadrature(plexdim, Nc, npoints, -1.0d0, +1.0d0, quad, ierr))
PetscCallA(DMPlexGetHeightStratum(dm, 0, cStart, cEnd, ierr))
! Allocating J, invJ before calling DMPlexComputeCellGeometryFEM
!
! allocate(J(4))
! allocate(invJ(4))
!
! sometimes works, but others does not, with a message:
!
! "corrupted size vs. prev_size" , followed by a stack trace
!
do c = cStart, cEnd-1
PetscCallA(DMPlexComputeCellGeometryFEM(dm, c, quad, vv, J, invJ, detJ, ierr))
CHKERRA(ierr)
print *, "Cell ", c
print "(a, *(f12.6))", "v0 =", vv
end do
PetscCallA(PetscQuadratureDestroy(quad, ierr))
PetscCallA(PetscFinalize(ierr))
end program test
!
! EXPECTED OUTPUT
! ---------------
!
! npoints = 1
!
! Cell 0
! v0 = 2.750000 2.500000
! Cell 1
! v0 = 7.750000 2.750000
! Cell 2
! v0 = 7.750000 7.750000
! Cell 3
! v0 = 2.750000 7.500000
!
! ---------------------------------
!
! npoints = 2
!
! Cell 0
! v0 = 4.110042 1.178633 1.101283 0.934616 4.565384 4.398717 1.223291 3.488034
! Cell 1
! v0 = 8.988034 1.101283 6.223291 1.223291 9.110042 4.110042 6.678633 4.565384
! Cell 2
! v0 = 8.988034 8.988034 9.110042 6.223291 6.223291 9.110042 6.678633 6.678633
! Cell 3
! v0 = 1.101283 8.821367 4.110042 9.065384 1.223291 5.601283 4.565384 6.511966
!
! ---------------------------------
#include <stdio.h>
#include <petscdmplex.h>
#include "petscdt.h"
int main(int argc, char *argv[]) {
DM dm;
PetscInt dim, Nc, npoints;
PetscQuadrature quad;
PetscInt cStart, cEnd;
PetscCall(PetscInitialize(&argc, &argv, NULL, NULL));
PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD, "mesh_quad.msh", "test", PETSC_TRUE, &dm));
dim = 2;
Nc = 2;
npoints = 2; // npoints = 1 works
PetscCall(PetscDTGaussTensorQuadrature(dim, Nc, npoints, -1.0, +1.0, &quad));
PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
for (PetscInt s = cStart; s < cEnd; s++) {
PetscReal vv[8], J[4], invJ[4], detJ; // vv[2] if npoints==1
PetscCall(DMPlexComputeCellGeometryFEM(dm, s, quad, vv, J, invJ, &detJ));
printf("Cell [%d] -- v =", s);
for (int k = 0; k < (npoints == 1 ? 2 : 8); k++) {
printf("%12.6f ", vv[k]);
}
printf("\n");
}
PetscCall(PetscQuadratureDestroy(&quad));
PetscCall(PetscFinalize());
}
/*
* EXPECTED OUTPUT
* ---------------
*
* npoints = 1
*
* Cell [0] -- v = 2.750000 2.500000
* Cell [1] -- v = 7.750000 2.750000
* Cell [2] -- v = 7.750000 7.750000
* Cell [3] -- v = 2.750000 7.500000
*
* ---------------------------------
*
* npoints = 2
*
* Cell [0] -- v = 4.110042 1.178633 1.101283 0.934616 4.565384 4.398717 1.223291 3.488034
* Cell [1] -- v = 8.988034 1.101283 6.223291 1.223291 9.110042 4.110042 6.678633 4.565384
* Cell [2] -- v = 8.988034 8.988034 9.110042 6.223291 6.223291 9.110042 6.678633 6.678633
* Cell [3] -- v = 1.101283 8.821367 4.110042 9.065384 1.223291 5.601283 4.565384 6.511966
*
* ---------------------------------
*
* Frequent output for npoints = 2 (wrong values)
* ----------------------------------------------
*
* Cell [0] -- v = 0.394338 -2.894338 2.788675 -0.788675 0.105662 -2.894338 2.211325 -0.788675
* Cell [1] -- v = 0.105662 -2.105662 2.605662 0.394338 0.394338 -2.105662 2.894338 0.394338
* Cell [2] -- v = -2.394338 0.394338 0.105662 -2.105662 -2.105662 0.394338 0.394338 -2.105662
* Cell [3] -- v = 0.105662 2.894338 -2.788675 0.788675 0.394338 2.894338 -2.211325 0.788675
*
*/
