Thanks both for your replies! Matteo Il 03/07/26 01:56, Matthew Knepley ha scritto:
On Thu, Jul 2, 2026 at 10:55 AM Matteo Semplice via petsc-users <[email protected]> wrote:
Dear PETSc,
I am a bit puzzled by the 3d tet mesh representations as DMPlex. I boiled it down to a mesh of a single tet. More specifically,
* if a I create the mesh as
const PetscInt dim=3; const PetscInt numCells=1; // The number of cells, only on process 0 const PetscInt numVertices=4; // The number of vertices owned by this process, or PETSC_DECIDE, only on process 0 const PetscInt numCorners=4; // The number of vertices for each cell, only on process 0 const PetscBool interpolate=PETSC_TRUE; // should intermediate mesh entities (faces, edges) be created? // An array of numCells×numCornersnumCells×numCorners numbers, the vertices for each cell, only on process 0 const PetscInt cells[] = {0,1,2,3 }; const PetscInt spaceDim=3; // The spatial dimension used for coordinates // An array of numVertices×spaceDimnumVertices×spaceDim numbers, the coordinates of each vertex, only on process 0 const PetscReal vertexCoords[] = { 0,0,0, 1,0,0, 0,1,0, 0,0,1 }; PetscCall( DMPlexCreateFromCellListPetsc(MPI_COMM_WORLD, dim, numCells, numVertices, numCorners, interpolate, cells, spaceDim, vertexCoords, dmMesh) );
I get the representation that I would expect:
DM Object: Mesh 1 MPI process type: plex Mesh in 3 dimensions: Number of 0-cells per rank: 4 Number of 1-cells per rank: 6 Number of 2-cells per rank: 4 Number of 3-cells per rank: 1 Labels: celltype: 4 strata with value/size (0 (4), 1 (6), 3 (4), 6 (1)) depth: 4 strata with value/size (0 (4), 1 (6), 2 (4), 3 (1)) Mesh with 1 cells
with DMPlexGetCellCoordinates giving me the coordinates
CELL 0 with vertices at ( 0.000, 0, 0.000 ) ( 1.000, 0, 0.000 ) ( 0.000, 1, 0.000 ) ( 0.000, 0, 1.000 )
* on the other hand, if I load the MSH file
$MeshFormat 2.2 0 8 $EndMeshFormat $Nodes 4 1 0.0 0.0 0.0 2 1.0 0.0 0.0 3 0.0 1.0 0.0 4 0.0 0.0 1.0 $EndNodes $Elements 1 1 4 2 0 1 1 2 3 4 $EndElements
with DMPlexCreateGmshFromFile, the first and second vertices are swapped and I get
DM Object: Mesh 1 MPI process type: plex Mesh in 3 dimensions: Number of 0-cells per rank: 4 Number of 1-cells per rank: 6 Number of 2-cells per rank: 4 Number of 3-cells per rank: 1 Labels: celltype: 4 strata with value/size (0 (4), 1 (6), 3 (4), 6 (1)) depth: 4 strata with value/size (0 (4), 1 (6), 2 (4), 3 (1)) Cell Sets: 1 strata with value/size (0 (1)) Mesh with 1 cells
with DMPlexGetCellCoordinates giving me
CELL 0 with vertices at ( 1.000, 0, 0.000 ) ( 0.000, 0, 0.000 ) ( 0.000, 1, 0.000 ) ( 0.000, 0, 1.000 )
I understand that this behaviour is somehow consistent with the formulas in DMPlexComputeTetrahedronGeometry_Internal, where you compute the Jacobians swappnig the order of vertices (and with 0.5 factors that I do not understand at all... it looks as if your reference element was not the unit simplex)
if (J) { for (d = 0; d < dim; d++) { /* I orient with outward face normals */ J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d])); J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d])); J[d * dim + 2] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d])); } }
but still I am worried by this fact that different ways of creating a mesh gives a DMPlex with different orientations.
Can you explain the rationale behind this or refer me to some documentation?
1. DMPlexCreateFromCellListPetsc() is very low-level, so I let you do whatever you want, including tell me to use inverted cells. DMPlexCreateGmshFromFile() is guaranteed to get a properly oriented simplex from GMsh, so I guarantee that it produces a properly oriented Plex simplex.
2. Here is the order of vertices, and face vertices, in a Plex tetrahedron.
https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/blob/main/src/dm...
I truly believe the only sensible way to describe a simple is to have every normal outward facing. The convention in Gmsh is to have some normals inward and some outward, which I think is ridiculous. It breaks the inherent symmetry of the simplex.
3. The reference simplex has coordinates 1 and -1, and thus volume 4/3. This is the convention of FIAT (https://urldefense.us/v3/__https://www.firedrakeproject.org/fiat/index.html_... )
4. You can always get my orientation for things by calling DMPlexOrient().
5. Orientations are representors for elements of the dihedral group of the k-cell. I am almost done with a description of this for every geometric cell supported (of course you can define any cell you want in Plex). I can show you this if you want, although it is unfinished.
Thanks,
Matt
Best
Matteo
-- Prof. Matteo Semplice Università degli Studi dell’Insubria Dipartimento di Scienza e Alta Tecnologia – DiSAT Professore Associato Via Valleggio, 11 – 22100 Como (CO) – Italia tel.: +39 031 2386316
-- 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_uCf... <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfs... >
-- Prof. Matteo Semplice Università degli Studi dell’Insubria Dipartimento di Scienza e Alta Tecnologia – DiSAT Professore Associato Via Valleggio, 11 – 22100 Como (CO) – Italia tel.: +39 031 2386316