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

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@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