[petsc-users] [C/Fortran] Calling DMPlexComputeCellGeometryFEM in a loop

[Noam T.]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250512/62b630b5/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mesh_quad.msh
Type: model/mesh
Size: 1441 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250512/62b630b5/attachment.msh>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: testF.F90
Type: text/x-fortran
Size: 2290 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250512/62b630b5/attachment.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: testC.c
Type: text/x-csrc
Size: 2455 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250512/62b630b5/attachment-0001.bin>