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