<div dir="ltr"><div dir="ltr">On Mon, May 12, 2025 at 5:29 PM Noam T. via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a>> wrote:</div><div class="gmail_quote gmail_quote_container"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div>Hello,</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">I ran into a seemingly non-consistent behavior of the function <span></span><pre><span><a href="https://urldefense.us/v3/__https://petsc.org/release/manualpages/DMPlex/DMPlexComputeCellGeometryFEM/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_rJABo0_$" target="_blank"><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">DMPlexComputeCellGeometryFEM</span></a></span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">(</span><span><a href="https://urldefense.us/v3/__https://petsc.org/release/manualpages/DM/DM/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_kyYWLFi$" target="_blank"><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">DM</span></a></span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)"> </span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">dm</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">,</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)"> </span><span><a href="https://urldefense.us/v3/__https://petsc.org/release/manualpages/Sys/PetscInt/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_jAVrwAd$" target="_blank"><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">PetscInt</span></a></span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)"> </span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">cell</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">,</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)"> </span><span><a href="https://urldefense.us/v3/__https://petsc.org/release/manualpages/DT/PetscQuadrature/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_jVjnRvX$" target="_blank"><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">PetscQuadrature</span></a></span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)"> </span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">quad</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">,</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)"> </span><span><a href="https://urldefense.us/v3/__https://petsc.org/release/manualpages/Sys/PetscReal/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_heQHCRx$" target="_blank"><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">PetscReal</span></a></span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)"> </span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">*</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">v</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">,</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)"> </span><span><a href="https://urldefense.us/v3/__https://petsc.org/release/manualpages/Sys/PetscReal/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_heQHCRx$" target="_blank"><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">PetscReal</span></a></span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)"> </span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">J</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">[],</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)"> </span><span><a href="https://urldefense.us/v3/__https://petsc.org/release/manualpages/Sys/PetscReal/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_heQHCRx$" target="_blank"><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">PetscReal</span></a></span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)"> </span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">invJ</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">[],</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)"> </span><span><a href="https://urldefense.us/v3/__https://petsc.org/release/manualpages/Sys/PetscReal/__;!!G_uCfscf7eWS!YEv7_aaCdN0cxV8sUhRBIACYZGNNcwQusNjEpfVep6yOAxRYE9f8kw4XiVENrjEvQEVwiheebMTQVhPw4eenYApT_heQHCRx$" target="_blank"><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">PetscReal</span></a></span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)"> </span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">*</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">detJ</span><span style="font-family:Menlo,Consolas,"Courier New",monospace;color:rgb(0,0,0)">)</span></pre>Calling it in a loop happens to crash, only some times, when the input <span style="font-family:Menlo,Consolas,"Courier New",monospace">npoints</span> is not equal to one.</div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Attached is a mesh file (<span style="font-family:Menlo,Consolas,"Courier New",monospace">msh</span> 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.</div></blockquote><div><br></div><div>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.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><div><div><div><div><div><span>When the quadrature created has a single point per dimension (<span style="font-family:Menlo,Consolas,"Courier New",monospace">npoints = 1</span>), there seems to be no issue. The program runs and the output is the expected. </span></div><div><span><br></span></div><div><span>On the C side, For npoints = 2 the program hangs, usually no errors; sometimes</span></div><div><span><br></span></div><div><span> <span style="font-family:Menlo,Consolas,"Courier New",monospace">*** stack smashing detected ***: terminated</span>,</span></div><div><span><span style="font-family:Menlo,Consolas,"Courier New",monospace">[:1345101] *** Process received signal ***</span><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">[:1345101] Signal: Aborted (6)</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">[:1345101] Signal code:  (-6)</span></div><span></span><br></span></div><div><span>or</span></div><div><span><br></span></div><div><span> <span style="font-family:Menlo,Consolas,"Courier New",monospace">LeakSanitizer: detected memory leaks</span> (when compiled with</span><span style="font-family:Menlo,Consolas,"Courier New",monospace"> -fsanitize=addres).</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"><br></span></div><div><span style="font-family:Arial,sans-serif">On the Fortran side, similar behavior for </span><span style="font-family:Menlo,Consolas,"Courier New",monospace">npoints = 1 or 2</span><span style="font-family:Arial,sans-serif">. </span><span style="font-family:Arial,sans-serif">Either no error at all (no error meaning the one that only mentions "</span><span style="font-family:Menlo,Consolas,"Courier New",monospace">probably memory access out of bounds</span><span style="font-family:Arial,sans-serif">"), or sometimes "</span><span style="font-family:Menlo,Consolas,"Courier New",monospace"><span style="font-family:Menlo,Consolas,"Courier New",monospace">corrupted size vs. prev_size</span></span><span style="font-family:Arial,sans-serif">", followed by a stack trace. </span></div><div><span style="font-family:Arial,sans-serif"><br></span></div><div><span style="font-family:Arial,sans-serif">Question: Do the arguments </span><span style="font-family:Menlo,Consolas,"Courier New",monospace">J, invJ </span><span style="font-family:Arial,sans-serif">need to be allocated (e.g. </span><span style="font-family:Menlo,Consolas,"Courier New",monospace">allocate(J(4))</span><span style="font-family:Arial,sans-serif">) before calling </span><span style="font-family:Menlo,Consolas,"Courier New",monospace"><span style="font-family:Menlo,Consolas,"Courier New",monospace">DMPlexComputeCellGeometryFEM</span></span><span style="font-family:Arial,sans-serif">? Allocation seems to be the cause for triggering an error.</span><span style="font-family:Menlo,Consolas,"Courier New",monospace"><br></span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"><br></span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">---<br></span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"><br></span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">NOTE:</span><span style="font-family:Arial,sans-serif"> The interface to </span><span style="font-family:Menlo,Consolas,"Courier New",monospace"><span><span style="font-family:Menlo,Consolas,"Courier New",monospace">DMPlexComputeCellGeometryFEM</span><span style="font-family:Arial,sans-serif">, under </span><span style="font-family:Menlo,Consolas,"Courier New",monospace">$PETSC_ARCH$/ftn/dm/petscdmplex.h90</span><span style="font-family:Arial,sans-serif">, read:</span></span></span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"><span><span style="font-family:Arial,sans-serif"><br></span></span></span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace"><span><span style="font-family:Arial,sans-serif"><span style="font-family:Menlo,Consolas,"Courier New",monospace">  interface DMPlexComputeCellGeometryFEM</span><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">  subroutine DMPlexComputeCellGeometryFEM(a,b,c,d,e,f,g, z)</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">  import tDM,tPetscQuadrature</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">  DM :: a</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">  PetscInt :: b</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">  PetscQuadrature :: c</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">  PetscReal :: d          ####   should it be d(*) ?<br></span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">  PetscReal :: e(*)</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">  PetscReal :: f(*)</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">  PetscReal :: g</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">  PetscErrorCode z</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">  end subroutine</span></div><div><span style="font-family:Menlo,Consolas,"Courier New",monospace">  end interface</span></div><span></span><br></span></span><span style="font-family:Arial,sans-serif"></span></span></div></div></div></div></div></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">Should argument <span style="font-family:Menlo,Consolas,"Courier New",monospace">d</span> (referring to the array of the mapped integration points <span style="font-family:Menlo,Consolas,"Courier New",monospace">v</span>) be defined differently? <br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)">---<br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><br></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span style="font-family:Arial,sans-serif">NOTE: Tested with PETSc 3.23.1 (from the tarball) and the latest git.</span></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span style="font-family:Arial,sans-serif">Configured with</span><span style="font-family:Menlo,Consolas,"Courier New",monospace"> ./configure --with-fortran-bindings=1</span></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span style="font-family:Menlo,Consolas,"Courier New",monospace">Using gcc/gfortran 11.4.0, OpenMPI 4.1.2<br></span></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span style="font-family:Menlo,Consolas,"Courier New",monospace"><br></span></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span style="font-family:Menlo,Consolas,"Courier New",monospace">Thanks,<br></span></div><div style="font-family:Arial,sans-serif;font-size:14px;color:rgb(0,0,0);background-color:rgb(255,255,255)"><span style="font-family:Menlo,Consolas,"Courier New",monospace">Noam<br></span></div></blockquote></div><div><br clear="all"></div><div><br></div><span class="gmail_signature_prefix">-- </span><br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!eAWWlW5LS8g4psK5Lef5MqtppJYejYI7g1LadD-JXWhiZVrx_jhz5HQ2T-VWxPh-ZdIPhN9KiG6Pbq7RNIV5$" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>