DMPlex: Unstructured Grids#

Representing Unstructured Grids#

The main advantage of DMPLEX in representing topology is that it treats all the different pieces of a mesh, e.g. cells, faces, edges, and vertices, in the same way. This allows the interface to be small and simple, while remaining flexible and general. This also allows “dimension independent programming”, which means that the same algorithm can be used unchanged for meshes of different shapes and dimensions.

All pieces of the mesh (vertices, edges, faces, and cells) are treated as points (or mesh entities), which are each identified by a PetscInt. A mesh is built by relating points to other points, in particular specifying a “covering” relation among the points. For example, an edge is defined by being covered by two vertices, and a triangle can be defined by being covered by three edges (or even by three vertices). This structure is known as a Hasse Diagram, which is a Directed Acyclic Graph (DAG) representing a cell complex using the covering relation. The graph edges represent the relation, which also encodes a partially ordered set (poset).

For example, we can encode the doublet mesh as in Fig. 8,

which can also be represented as the DAG in Fig. 9.

To use the PETSc API, we consecutively number the mesh pieces. The PETSc convention in 3 dimensions is to number first cells, then vertices, then faces, and then edges. In 2 dimensions the convention is to number faces, vertices, and then edges. In terms of the labels in Fig. 8, these numberings are

First, we declare the set of points present in a mesh,

DMPlexSetChart(dm, 0, 11);

Note that a chart here corresponds to a semi-closed interval (e.g $[0,11) = \{0,1,\ldots,10\}$) specifying the range of indices we’d like to use to define points on the current rank. We then define the covering relation, which we call the cone, which are also the in-edges in the DAG. In order to preallocate correctly, we first provide sizes,

/* DMPlexSetConeSize(dm, point, number of points that cover the point); */
DMPlexSetConeSize(dm, 0, 3);
DMPlexSetConeSize(dm, 1, 3);
DMPlexSetConeSize(dm, 6, 2);
DMPlexSetConeSize(dm, 7, 2);
DMPlexSetConeSize(dm, 8, 2);
DMPlexSetConeSize(dm, 9, 2);
DMPlexSetConeSize(dm, 10, 2);
DMSetUp(dm);

and then point values (recall each point is an integer that represents a single geometric entity, a cell, face, edge, or vertex),

/* DMPlexSetCone(dm, point, [points that cover the point]); */
DMPlexSetCone(dm, 0, [6, 7, 8]);
DMPlexSetCone(dm, 1, [7, 9, 10]);
DMPlexSetCone(dm, 6, [2, 3]);
DMPlexSetCone(dm, 7, [3, 4]);
DMPlexSetCone(dm, 8, [4, 2]);
DMPlexSetCone(dm, 9, [4, 5]);
DMPlexSetCone(dm, 10, [5, 3]);

There is also an API for providing the dual relation, using DMPlexSetSupportSize() and DMPlexSetSupport(), but this can be calculated automatically using the provided DMPlexSetConeSize() and DMPlexSetCone() information and then calling

DMPlexSymmetrize(dm);

The “symmetrization” is in the sense of the DAG. Each point knows its covering (cone) and each point knows what it covers (support). Note that when using automatic symmetrization, cones will be ordered but supports will not. The user can enforce an ordering on supports by rewriting them after symmetrization using DMPlexSetSupport().

In order to support efficient queries, we construct fast search structures and indices for the different types of points using

DMPlexStratify(dm);

Grid Point Orientations#

TODO: fill this out with regard to orientations.

Should probably reference Hapla and summarize what’s there.

Dealing with Periodicity#

Plex allows you to represent periodic domains is two ways. Using the default scheme, periodic topology can be represented directly. This ensures that all topological queries can be satisfied, but then care must be taken in representing functions over the mesh, such as the coordinates. The second method is to use a non-periodic topology, but connect certain mesh points using the local-to-global map for that DM. This allows a more general set of mappings to be implemented, such as partial twists, but topological queries on the periodic boundary cease to function.

For the default scheme, a call to DMLocalizeCoordinates() (which usually happens automatically on mesh creation) creates a second, discontinuous coordinate field. These values can be accessed using DMGetCellCoordinates() and DMGetCellCoordinatesLocal(). Plex provides a convenience method, DMPlexGetCellCoordinates(), that extracts cell coordinates correctly, depending on the periodicity of the mesh. An example of its use is shown below:

const PetscScalar *array;
PetscScalar       *coords = NULL;
PetscInt           numCoords;
PetscBool          isDG;

PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &numCoords, &array, &coords));
for (PetscInt cc = 0; cc < numCoords/dim; ++cc) {
  if (cc > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " -- "));
  PetscCall(PetscPrintf(PETSC_COMM_SELF, "("));
  for (PetscInt d = 0; d < dim; ++d) {
    if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
    PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(coords[cc * dim + d])));
  }
  PetscCall(PetscPrintf(PETSC_COMM_SELF, ")"));
}
PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &numCoords, &array, &coords));

Connecting Grids to Data Using PetscSection:#

A PetscSection is used to describe the connection between the grid and data associated with the grid. Specifically, it assigns a number of dofs to each mesh entity on the DAG. See PetscSection: Connecting Grids to Data for more details. Using the mesh from Fig. 8, we provide an example of creating a PetscSection for a single field. We can lay out data for a continuous Galerkin $P_3$ finite element method,

PetscInt pStart, pEnd, cStart, cEnd, c, vStart, vEnd, v, eStart, eEnd, e;

DMPlexGetChart(dm, &pStart, &pEnd);
DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);   // cells
DMPlexGetHeightStratum(dm, 1, &eStart, &eEnd);   // edges
DMPlexGetHeightStratum(dm, 2, &vStart, &vEnd);   // vertices, equivalent to DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
PetscSectionSetChart(s, pStart, pEnd);
for(c = cStart; c < cEnd; ++c)
    PetscSectionSetDof(s, c, 1);
for(v = vStart; v < vEnd; ++v)
    PetscSectionSetDof(s, v, 1);
for(e = eStart; e < eEnd; ++e)
    PetscSectionSetDof(s, e, 2); // two dof on each edge
PetscSectionSetUp(s);

DMPlexGetHeightStratum() returns all the points of the requested height in the DAG. Since this problem is in two dimensions the edges are at height 1 and the vertices at height 2 (the cells are always at height 0). One can also use DMPlexGetDepthStratum() to use the depth in the DAG to select the points. DMPlexGetDepth(dm,&depth) returns the depth of the DAG, hence DMPlexGetDepthStratum(dm,depth-1-h,) returns the same values as DMPlexGetHeightStratum(dm,h,).

For $P_3$ elements there is one degree of freedom at each vertex, 2 along each edge (resulting in a total of 4 degrees of freedom along each edge including the vertices, thus being able to reproduce a cubic function) and 1 degree of freedom within the cell (the bubble function which is zero along all edges).

Now a PETSc local vector can be created manually using this layout,

PetscSectionGetStorageSize(s, &n);
VecSetSizes(localVec, n, PETSC_DETERMINE);
VecSetFromOptions(localVec);

When working with DMPLEX and PetscFE (see below) one can simply get the sections (and related vectors) with

DMSetLocalSection(dm, s);
DMGetLocalVector(dm, &localVec);
DMGetGlobalVector(dm, &globalVec);

Data Layout using DMPLEX and PetscFE#

A DM can automatically create the local section if given a description of the discretization, for example using a PetscFE object. We demonstrate this by creating a PetscFE that can be configured from the command line. It is a single, scalar field, and is added to the DM using DMSetField(). When a local or global vector is requested, the DM builds the local and global sections automatically.

DMPlexIsSimplex(dm, &simplex);
PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, -1, &fe);
DMSetField(dm, 0, NULL, (PetscObject) fe);
DMCreateDS(dm);

Here the call to DMSetField() declares the discretization will have one field with the integer label 0 that has one degree of freedom at each point on the DMPlex. To get the $P_3$ section above, we can either give the option -petscspace_degree 3, or call PetscFECreateLagrange() and set the degree directly.

Partitioning and Ordering#

In the same way as MatPartitioning or MatGetOrdering(), give the results of a partitioning or ordering of a graph defined by a sparse matrix, PetscPartitionerDMPlexPartition or DMPlexPermute are encoded in an IS. However, the graph is not the adjacency graph of the matrix but the mesh itself. Once the mesh is partitioned and reordered, the data layout from a PetscSection can be used to automatically derive a problem partitioning/ordering.

Evaluating Residuals#

The evaluation of a residual or Jacobian, for most discretizations has the following general form:

• Traverse the mesh, picking out pieces (which in general overlap),

• Extract some values from the current solution vector, associated with this piece,

• Calculate some values for the piece, and

• Insert these values into the residual vector

DMPlex separates these different concerns by passing sets of points from mesh traversal routines to data extraction routines and back. In this way, the PetscSection which structures the data inside a Vec does not need to know anything about the mesh inside a DMPLEX.

The most common mesh traversal is the transitive closure of a point, which is exactly the transitive closure of a point in the DAG using the covering relation. In other words, the transitive closure consists of all points that cover the given point (generally a cell) plus all points that cover those points, etc. So in 2d the transitive closure for a cell consists of edges and vertices while in 3d it consists of faces, edges, and vertices. Note that this closure can be calculated orienting the arrows in either direction. For example, in a finite element calculation, we calculate an integral over each element, and then sum up the contributions to the basis function coefficients. The closure of the element can be expressed discretely as the transitive closure of the element point in the mesh DAG, where each point also has an orientation. Then we can retrieve the data using PetscSection methods,

PetscScalar *a;
PetscInt     numPoints, *points = NULL, p;

VecGetArrayRead(u,&a);
DMPlexGetTransitiveClosure(dm,cell,PETSC_TRUE,&numPoints,&points);
for (p = 0; p <= numPoints*2; p += 2) {
  PetscInt dof, off, d;

  PetscSectionGetDof(section, points[p], &dof);
  PetscSectionGetOffset(section, points[p], &off);
  for (d = 0; d <= dof; ++d) {
    myfunc(a[off+d]);
  }
}
DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &numPoints, &points);
VecRestoreArrayRead(u, &a);

This operation is so common that we have built a convenience method around it which returns the values in a contiguous array, correctly taking into account the orientations of various mesh points:

const PetscScalar *values;
PetscInt           csize;

DMPlexVecGetClosure(dm, section, u, cell, &csize, &values);
// Do integral in quadrature loop putting the result into r[]
DMPlexVecRestoreClosure(dm, section, u, cell, &csize, &values);
DMPlexVecSetClosure(dm, section, residual, cell, &r, ADD_VALUES);

A simple example of this kind of calculation is in DMPlexComputeL2Diff_Plex() (source). Note that there is no restriction on the type of cell or dimension of the mesh in the code above, so it will work for polyhedral cells, hybrid meshes, and meshes of any dimension, without change. We can also reverse the covering relation, so that the code works for finite volume methods where we want the data from neighboring cells for each face:

PetscScalar *a;
PetscInt     points[2*2], numPoints, p, dofA, offA, dofB, offB;

VecGetArray(u,  &a);
DMPlexGetTransitiveClosure(dm, cell, PETSC_FALSE, &numPoints, &points);
assert(numPoints == 2);
PetscSectionGetDof(section, points[0*2], &dofA);
PetscSectionGetDof(section, points[1*2], &dofB);
assert(dofA == dofB);
PetscSectionGetOffset(section, points[0*2], &offA);
PetscSectionGetOffset(section, points[1*2], &offB);
myfunc(a[offA], a[offB]);
VecRestoreArray(u, &a);

This kind of calculation is used in TS Tutorial ex11.

Saving and Loading DMPlex Data with HDF5#

PETSc allows users to save/load DMPLEXs representing meshes, PetscSections representing data layouts on the meshes, and Vecs defined on the data layouts to/from an HDF5 file in parallel, where one can use different number of processes for saving and for loading.

DMViewFromOptions(dm, NULL, "-dm_view");
VecViewFromOptions(vec, NULL, "-vec_view")

$ ./myprog -dm_view hdf5:myprog.h5 -vec_view hdf5:myprog.h5::append

PetscViewer  viewer;

PetscViewerHDF5Open(PETSC_COMM_WORLD, "example.h5", FILE_MODE_WRITE, &viewer);

PetscObjectSetName((PetscObject)dm, "plexA");
PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_PETSC);
DMView(dm, viewer);
PetscViewerPopFormat(viewer);

DMPlexTopologyView(dm, viewer);
DMPlexCoordinatesView(dm, viewer);
DMPlexLabelsView(dm, viewer);

DMClone(dm, &sdm);
PetscObjectSetName((PetscObject)sdm, "dmA");
DMSetLocalSection(sdm, s);
DMPlexSectionView(dm, viewer, sdm);

Vec  vec;

DMGetLocalVector(sdm, &vec);

PetscObjectSetName((PetscObject)vec, "vecA");
DMPlexLocalVectorView(dm, viewer, sdm, vec);

PetscViewerDestroy(&viewer);

$ h5dump --contents example.h5
HDF5 "example.h5" {
FILE_CONTENTS {
 group      /
 group      /topologies
 group      /topologies/plexA
 group      /topologies/plexA/dms
 group      /topologies/plexA/dms/dmA
 dataset    /topologies/plexA/dms/dmA/order
 group      /topologies/plexA/dms/dmA/section
 dataset    /topologies/plexA/dms/dmA/section/atlasDof
 dataset    /topologies/plexA/dms/dmA/section/atlasOff
 group      /topologies/plexA/dms/dmA/vecs
 group      /topologies/plexA/dms/dmA/vecs/vecA
 dataset    /topologies/plexA/dms/dmA/vecs/vecA/vecA
 group      /topologies/plexA/labels
 group      /topologies/plexA/topology
 dataset    /topologies/plexA/topology/cells
 dataset    /topologies/plexA/topology/cones
 dataset    /topologies/plexA/topology/order
 dataset    /topologies/plexA/topology/orientation
 }
}

-dm_plex_view_hdf5_storage_version 3.0.0

$ h5dump -a /dmplex_storage_version example.h5

$ h5dump --contents example.h5
HDF5 "example.h5" {
FILE_CONTENTS {
 ...
 group      /topologies/plexA/topology
 dataset    /topologies/plexA/topology/permutation
 group      /topologies/plexA/topology/strata
 group      /topologies/plexA/topology/strata/0
 dataset    /topologies/plexA/topology/strata/0/cone_sizes
 dataset    /topologies/plexA/topology/strata/0/cones
 dataset    /topologies/plexA/topology/strata/0/orientations
 group      /topologies/plexA/topology/strata/1
 dataset    /topologies/plexA/topology/strata/1/cone_sizes
 dataset    /topologies/plexA/topology/strata/1/cones
 dataset    /topologies/plexA/topology/strata/1/orientations
 group      /topologies/plexA/topology/strata/2
 dataset    /topologies/plexA/topology/strata/2/cone_sizes
 dataset    /topologies/plexA/topology/strata/2/cones
 dataset    /topologies/plexA/topology/strata/2/orientations
 group      /topologies/plexA/topology/strata/3
 dataset    /topologies/plexA/topology/strata/3/cone_sizes
 dataset    /topologies/plexA/topology/strata/3/cones
 dataset    /topologies/plexA/topology/strata/3/orientations
 }
}

PetscViewerHDF5Open(PETSC_COMM_WORLD, "example.h5", FILE_MODE_READ, &viewer);

DMCreate(PETSC_COMM_WORLD, &dm);
DMSetType(dm, DMPLEX);
PetscObjectSetName((PetscObject)dm, "plexA");
PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_PETSC);
DMLoad(dm, viewer);
PetscViewerPopFormat(viewer);

PetscSF  sfO;

DMCreate(PETSC_COMM_WORLD, &dm);
DMSetType(dm, DMPLEX);
PetscObjectSetName((PetscObject)dm, "plexA");
PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_PETSC);
DMPlexTopologyLoad(dm, viewer, &sfO);
DMPlexCoordinatesLoad(dm, viewer, sfO);
PetscViewerPopFormat(viewer);

DM        distributedDM;
PetscInt  overlap = 1;
PetscSF   sfDist, sf;

DMPlexDistribute(dm, overlap, &sfDist, &distributedDM);
if (distributedDM) {
  DMDestroy(&dm);
  dm = distributedDM;
  PetscObjectSetName((PetscObject)dm, "plexA");
}
PetscSFCompose(sfO, sfDist, &sf);
PetscSFDestroy(&sfO);
PetscSFDestroy(&sfDist);

PetscSF  globalDataSF, localDataSF;

DMClone(dm, &sdm);
PetscObjectSetName((PetscObject)sdm, "dmA");
DMPlexSectionLoad(dm, viewer, sdm, sf, &globalDataSF, &localDataSF);

Vec  vec;

DMGetLocalVector(sdm, &vec);

PetscObjectSetName((PetscObject)vec, "vecA");
DMPlexLocalVectorLoad(dm, viewer, sdm, localDataSF, localVec);

PetscViewerDestroy(&viewer);

Metric-based mesh adaptation#

DMPLEX supports mesh adaptation using the Riemannian metric framework. The idea is to use a Riemannian metric space within the mesher. The metric space dictates how mesh resolution should be distributed across the domain. Using this information, the remesher transforms the mesh such that it is a unit mesh when viewed in the metric space. That is, the image of each of its elements under the mapping from Euclidean space into the metric space has edges of unit length.

One of the main advantages of metric-based mesh adaptation is that it allows for fully anisotropic remeshing. That is, it provides a means of controlling the shape and orientation of elements in the adapted mesh, as well as their size. This can be particularly useful for advection-dominated and directionally-dependent problems.

See [Ala10] for further details on metric-based anisotropic mesh adaptation.

The two main ingredients for metric-based mesh adaptation are an input mesh (i.e. the DMPLEX) and a Riemannian metric. The implementation in PETSc assumes that the metric is piecewise linear and continuous across elemental boundaries. Such an object can be created using the routine

DMPlexMetricCreate(DM dm, PetscInt field, Vec *metric);

A metric must be symmetric positive-definite, so that distances may be properly defined. This can be checked using

DMPlexMetricEnforceSPD(DM dm, Vec metricIn, PetscBool restrictSizes, PetscBool restrictAnisotropy, Vec metricOut, Vec determinant);

This routine may also be used to enforce minimum and maximum tolerated metric magnitudes (i.e. cell sizes), as well as maximum anisotropy. These quantities can be specified using

DMPlexMetricSetMinimumMagnitude(DM dm, PetscReal h_min);
DMPlexMetricSetMaximumMagnitude(DM dm, PetscReal h_max);
DMPlexMetricSetMaximumAnisotropy(DM dm, PetscReal a_max);

or the command line arguments

-dm_plex_metric_h_min <h_min>
-dm_plex_metric_h_max <h_max>
-dm_plex_metric_a_max <a_max>

One simple way to combine two metrics is by simply averaging them entry-by-entry. Another is to intersect them, which amounts to choosing the greatest level of refinement in each direction. These operations are available in PETSc through the routines

DMPlexMetricAverage(DM dm, PetscInt numMetrics, PetscReal weights[], Vec metrics[], Vec metricAvg);
DMPlexMetricIntersection(DM dm, PetscInt numMetrics, Vec metrics[], Vec metricInt);

However, before combining metrics, it is important that they are scaled in the same way. Scaling also allows the user to control the number of vertices in the adapted mesh (in an approximate sense). This is achieved using the $L^p$ normalization framework, with the routine

DMPlexMetricNormalize(DM dm, Vec metricIn, PetscBool restrictSizes, PetscBool restrictAnisotropy, Vec metricOut, Vec determinant);

There are two important parameters for the normalization: the normalization order $p$ and the target metric complexity, which is analogous to the vertex count. They are controlled using

DMPlexMetricSetNormalizationOrder(DM dm, PetscReal p);
DMPlexMetricSetTargetComplexity(DM dm, PetscReal target);

or the command line arguments

-dm_plex_metric_p <p>
-dm_plex_metric_target_complexity <target>

Two different metric-based mesh adaptation tools are available in PETSc:

• Pragmatic;

• Mmg/ParMmg.

Mmg is a serial package, whereas ParMmg is the MPI version. Note that surface meshing is not currently supported and that ParMmg works only in three dimensions. Mmg can be used for both two and three dimensional problems. Pragmatic, Mmg and ParMmg may be specified by the command line arguments

-dm_adaptor pragmatic
-dm_adaptor mmg
-dm_adaptor parmmg

Once a metric has been constructed, it can be used to perform metric-based mesh adaptation using the routine

DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DMLabel rgLabel, DM dmAdapt);

where bdLabel and rgLabel are boundary and interior tags to be preserved under adaptation, respectively.

Footnotes