orientation of DMPLEX for tetrehedral meshes
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? 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
DMPlex internal ordering is not the same as VTK ordering. The code you supply will have an inward-facing normal (as evaluated by DMPlex geometry routines). The IO routines translate to the internal ordering convention. If you are creating a mesh manually or writing a new file reader, you should use (or convert to) the DMPlex convention. Matteo Semplice via petsc-users <[email protected]> writes:
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?
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
Thanks! is there a document explaining the DMPLex internal ordering? (and by the way what's the reason to have the factors 0.5 in the jacobians? I would not expect them if the reference tet is (0,0,0), (1,0,0), etc 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])); } } Thanks Matteo Il 02/07/26 17:22, Jed Brown ha scritto:
DMPlex internal ordering is not the same as VTK ordering. The code you supply will have an inward-facing normal (as evaluated by DMPlex geometry routines). The IO routines translate to the internal ordering convention. If you are creating a mesh manually or writing a new file reader, you should use (or convert to) the DMPlex convention.
Matteo Semplice via petsc-users<[email protected]> writes:
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?
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
-- 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
See DMPlexCreateReferenceCell_Internal (in plexcreate.c), and look at the coordinates. One common convention for the reference line segment is (-1, 1), which is the natural domain for Legendre polynomials and quadrature. Some packages use this, and some use (0, 1). There is a similar choice of reference element for triangles and tetrahedra. Matteo Semplice <[email protected]> writes:
Thanks!
is there a document explaining the DMPLex internal ordering?
(and by the way what's the reason to have the factors 0.5 in the jacobians? I would not expect them if the reference tet is (0,0,0), (1,0,0), etc
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])); } }
Thanks
Matteo
Il 02/07/26 17:22, Jed Brown ha scritto:
DMPlex internal ordering is not the same as VTK ordering. The code you supply will have an inward-facing normal (as evaluated by DMPlex geometry routines). The IO routines translate to the internal ordering convention. If you are creating a mesh manually or writing a new file reader, you should use (or convert to) the DMPlex convention.
Matteo Semplice via petsc-users<[email protected]> writes:
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?
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
-- 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
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... >
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
participants (3)
-
Jed Brown -
Matteo Semplice -
Matthew Knepley