On Mon, May 12, 2025 at 5:29 PM Noam T. via petsc-users < petsc-users@mcs.anl.gov> wrote:
> 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. > Yes, this is what happens when interfaces evolve. When I originally wrote the function, I had only implemented affine geometry. When we implemented other coordinate spaces, the meaning of the arrays (x, J, JInv, detJ) changed. Now they held the evaluations at each quadrature point since they were not constant anymore. You need to allocate room for npoints copies of all of them. I need to fix the documentation. Thanks, Matt > 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 > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!eAWWlW5LS8g4psK5Lef5MqtppJYejYI7g1LadD-JXWhiZVrx_jhz5HQ2T-VWxPh-ZdIPhN9KiG6PbssvMu8E$ <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!eAWWlW5LS8g4psK5Lef5MqtppJYejYI7g1LadD-JXWhiZVrx_jhz5HQ2T-VWxPh-ZdIPhN9KiG6Pbq7RNIV5$ >
