[petsc-users] orientation of DMPLEX for tetrehedral meshes

[Jed Brown]

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 <petsc-users at mcs.anl.gov> 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