On Thu, Jul 2, 2026 at 10:55 AM Matteo Semplice via petsc-users <petsc-users@mcs.anl.gov> 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://gitlab.com/petsc/petsc/-/blob/main/src/dm/impls/plex/plexinterpolate.c?ref_type=heads#L194

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://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