[petsc-users] orientation of DMPLEX for tetrehedral meshes

[Jed Brown]

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 <matteo.semplice at uninsubria.it> 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<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
>
> -- 
> 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