Decoding the DMPlex

The purpose of this notebook is to provide a tutorial for the DMPlex features in PETSc. DMPlex and its sub-objects are an attempt to properly abstract out the concept of grids and the assignment of degree of freedom information to entities in that grid. The hope is that this will allow for easy implementation of different discretization ideas and subsequently lead to their fair comparison. To be able to use DMPlex you need to speak its language and understand how to get the information your methods need. In this tutorial I will explain and demonstrate the functionality of the DMPlex API.



In [1]:

    
import sys,petsc4py
petsc4py.init(sys.argv)
from petsc4py import PETSc
import numpy as np

One way to create a DMPlex is to specify coordinates of vertices and cell connectivities. Here we encode a simple 2 by 2 element mesh of quads.



In [2]:

    
dim = 2
coords = np.asarray([[0.0, 0.0],
                     [0.5, 0.0],
                     [1.0, 0.0],
                     [0.0, 0.5],
                     [0.5, 0.5],
                     [1.0, 0.5],
                     [0.0, 1.0],
                     [0.5, 1.0],
                     [1.0, 1.0]])
cells = np.asarray([[0,1,4,3],
                    [1,2,5,4],
                    [3,4,7,6],
                    [4,5,8,7]],dtype='int32')

Now we initialize the DMPlex using this mesh information.



In [3]:

    
plex = PETSc.DMPlex().createFromCellList(dim,cells,coords)

PETSc converts the mesh into their abstraction of a topology which they store as a Hasse Diagram. Essentially this is a list of integers which encodes all the entities of each dimenion. We can use the view method to see what the DMPlex has encoded (broken, prints to the terminal but not here, capture magic doesn't work).



In [4]:

    
plex.view()

The view reveals that we have a mesh in 2 dimensions and that we have 9 0-cells (vertices), 12 1-cells (edges), and 4 2-cells (quads). All mesh entities are stored as integers in a single array called a chart. Each entity in the chart is called a point. (At this point it would be good to make some kind of plot with all points numbered, I suggest sketching one.)

Cones and Supports

The cone of a point consists of the points of a dimension lower which make up that entity. So we can loop through the chart and print out each points' cone.



In [17]:

    
pStart,pEnd = plex.getChart()
for i in range(pStart,pEnd):
    print "point =", i, "\tcone =", plex.getCone(i)









    



point = 0 	cone = [13 14 15 16]
point = 1 	cone = [17 18 19 14]
point = 2 	cone = [15 20 21 22]
point = 3 	cone = [19 23 24 20]
point = 4 	cone = []
point = 5 	cone = []
point = 6 	cone = []
point = 7 	cone = []
point = 8 	cone = []
point = 9 	cone = []
point = 10 	cone = []
point = 11 	cone = []
point = 12 	cone = []
point = 13 	cone = [4 5]
point = 14 	cone = [5 8]
point = 15 	cone = [8 7]
point = 16 	cone = [7 4]
point = 17 	cone = [5 6]
point = 18 	cone = [6 9]
point = 19 	cone = [9 8]
point = 20 	cone = [ 8 11]
point = 21 	cone = [11 10]
point = 22 	cone = [10  7]
point = 23 	cone = [ 9 12]
point = 24 	cone = [12 11]

Note that the numbering is completely different from our original mesh encoding. Here we summarize what we observe from the cones of the chart entities.

Points 0 through 3 correspond to quad cells. So their cones are made up of lists of 4 integers which refer to the lower dimensional entities which make up that cell--the edges.
Points 4 through 12 correspond to vertices. These are the lowest dimensional object we have and thus they are empty.
Points 13 through 24 correspond to edges. Each edge is made up of two vertices.

Similarly, each point has a support. The support of a point is the list which consists of points of a higher dimension which contain the point in its cone. So we can now repeat the above exercise but for the support.



In [18]:

    
for i in range(pStart,pEnd):
    print "point =", i, "\tsupport =", plex.getSupport(i)









    



point = 0 	support = []
point = 1 	support = []
point = 2 	support = []
point = 3 	support = []
point = 4 	support = [13 16]
point = 5 	support = [13 14 17]
point = 6 	support = [17 18]
point = 7 	support = [15 16 22]
point = 8 	support = [14 15 19 20]
point = 9 	support = [18 19 23]
point = 10 	support = [21 22]
point = 11 	support = [20 21 24]
point = 12 	support = [23 24]
point = 13 	support = [0]
point = 14 	support = [0 1]
point = 15 	support = [0 2]
point = 16 	support = [0]
point = 17 	support = [1]
point = 18 	support = [1]
point = 19 	support = [1 3]
point = 20 	support = [2 3]
point = 21 	support = [2]
point = 22 	support = [2]
point = 23 	support = [3]
point = 24 	support = [3]

Points 0 through 3 (quads) have no support, there is nothing of higher dimension in this mesh
Points 4 through 12 (vertices) have at least 2 edges in their support and the middle (8) has 4 edges
Points 13 through 24 (edges) have at least 1 cell in their support as many as 2

So the DMPlex is a dimension-independent, low memory, abstract representation of a topology which we use to represent grid objects. The DMPlex knows nothing about elements, basis functions, fluxes, degrees of freedom, etc. It is just the topology itself and is completely general. DMPlex can be used to construct any kind of topological relation. Here we created one from a cell list and then accessed its cone/support information. A DMPlex can also be built by hand using the appropriate set routines or with other kinds of constructors available in the API.

Labeling

DMPlex provides support for the labeling of points. This can be helpful if you would like to flag certain entities for some reason. By default, the DMPlex comes with a label called 'depth'. This labels each entity based on how deep it is in the chart. You could also think of it as the dimensionality of the objects. Here we can check that the label does exist.



In [22]:

    
for i in range(plex.getNumLabels()):
    name = plex.getLabelName(i)
    print "label name = %s" % name, "\tlabel size = %d" % plex.getLabelSize(name)









    



label name = boundary 	label size = 1
label name = depth 	label size = 3

So the label 'depth' does exist and we see that there are 3 different entries. Now we will loop over each item in the DMPlex and print the value of the label.



In [20]:

    
for i in range(pStart,pEnd):
    print "point =",i, "\tlabel(depth) = %d" % plex.getLabelValue("depth",i)









    



point = 0 	label(depth) = 2
point = 1 	label(depth) = 2
point = 2 	label(depth) = 2
point = 3 	label(depth) = 2
point = 4 	label(depth) = 0
point = 5 	label(depth) = 0
point = 6 	label(depth) = 0
point = 7 	label(depth) = 0
point = 8 	label(depth) = 0
point = 9 	label(depth) = 0
point = 10 	label(depth) = 0
point = 11 	label(depth) = 0
point = 12 	label(depth) = 0
point = 13 	label(depth) = 1
point = 14 	label(depth) = 1
point = 15 	label(depth) = 1
point = 16 	label(depth) = 1
point = 17 	label(depth) = 1
point = 18 	label(depth) = 1
point = 19 	label(depth) = 1
point = 20 	label(depth) = 1
point = 21 	label(depth) = 1
point = 22 	label(depth) = 1
point = 23 	label(depth) = 1
point = 24 	label(depth) = 1

The depths listed match our intuition of which entities were which when we looked at the cones and support. The DMPlex has support for identifying the range of indices in the chart which correspond to each value of the depth, the so-called depth stratum. (I do not understand what height stratum is for yet)



In [9]:

    
for i in range(plex.getDepth()+1):
    print "depth = %d" % i,"\tdepth stratum = ",plex.getDepthStratum(i),"\theight stratum = ",plex.getHeightStratum(i)









    



depth = 0 	depth stratum =  (4, 13) 	height stratum =  (0, 4)
depth = 1 	depth stratum =  (13, 25) 	height stratum =  (13, 25)
depth = 2 	depth stratum =  (0, 4) 	height stratum =  (4, 13)

We can also use labels to mark, say, boundary edges. These are the edges with only 1 entry in the support.



In [23]:

    
plex.createLabel("boundary")
for i in range(pStart,pEnd):
    if plex.getLabelValue("depth",i) == 1: # this is an edge
        if plex.getSupportSize(i) == 1:    # only one cell has it as an edge
            plex.setLabelValue("boundary",i,1)
    print "point =", i, "\tlabel(boundary) = %d" % plex.getLabelValue("boundary",i)









    



point = 0 	label(boundary) = -1
point = 1 	label(boundary) = -1
point = 2 	label(boundary) = -1
point = 3 	label(boundary) = -1
point = 4 	label(boundary) = -1
point = 5 	label(boundary) = -1
point = 6 	label(boundary) = -1
point = 7 	label(boundary) = -1
point = 8 	label(boundary) = -1
point = 9 	label(boundary) = -1
point = 10 	label(boundary) = -1
point = 11 	label(boundary) = -1
point = 12 	label(boundary) = -1
point = 13 	label(boundary) = 1
point = 14 	label(boundary) = -1
point = 15 	label(boundary) = -1
point = 16 	label(boundary) = 1
point = 17 	label(boundary) = 1
point = 18 	label(boundary) = 1
point = 19 	label(boundary) = -1
point = 20 	label(boundary) = -1
point = 21 	label(boundary) = 1
point = 22 	label(boundary) = 1
point = 23 	label(boundary) = 1
point = 24 	label(boundary) = 1

The default values are set to -1 and the boundary edges were marked with a 1. Labels aren't so useful on their own--the useful part about labeling things is that you can also get index sets of all entities with the same value of label. This means that if we wanted an index set which maps vertex numbers to the chartID we can get the PETSc IS for the 'depth' label for a value of 0 (again view doesn't print here).



In [11]:

    
vis = plex.getStratumIS("depth",0)
vis.view()

Note that there are also routines which would appear to do the same as the above. However, their index sets return a local to global mapping of vertices (smallest depth in chart) and cells (largest depth in chart) useful in parallel.



In [12]:

    
vis = plex.getVertexNumbering()
vis.view()
cis = plex.getCellNumbering()
cis.view()

Meets and Joins

Many times in performing grid operations, you need to know how lower and/or higher dimensional items are connected to each other. In the PETSc parlance, these are called meets and joins. A meet of a set of points is the intersection of the points' cones and a join is the intersection of the points' support. Here we demonstrate the concept with a few examples.



In [13]:

    
# Two cells, meet is the common edge, no join
pnts = [0,1]
print "meet =",plex.getMeet(pnts),"\tjoin =",plex.getJoin(pnts)









    



meet = [14] 	join = []



In [14]:

    
# Two edges, meet is the common vertex, join is the cell to which they are both connected
pnts = [14,19]
print "meet =",plex.getMeet(pnts),"\tjoin =",plex.getJoin(pnts)









    



meet = [8] 	join = [1]



In [15]:

    
# Two vertices, no meet, join is the common edge to which they are both connected
pnts = [8,11]
print "meet =",plex.getMeet(pnts),"\tjoin =",plex.getJoin(pnts)









    



meet = [] 	join = [20]

Transitive Closure

The transitive closure of a point in the DMPlex is a list of all reachable points from the given point, by default in the 'in-edge' direction. The transitive closure is then a set created by recursively taking the union on all points in the cone and its cones. In other words, it is all points of lower or equal dimension that this point can "reach". The routine also returns an array parallel to the closure array which defines how the points are oriented relative to the give point (e.g. for a cell, you might need to flip some edges to follow a right-handed convention). There is a flag in the routine, useCone, which if set to False will perform the same operation but in the 'out-edge' direction (that is, instead of recursively operating on cones, it will use the supports).



In [16]:

    
for i in range(pStart,pEnd):
    coneclose,orient = plex.getTransitiveClosure(i)
    suppclose,orient = plex.getTransitiveClosure(i,useCone=False)
    print "point =",i,"\tclosure(cone) =",coneclose,"\tclosure(supp) =",suppclose









    



point = 0 	closure(cone) = [ 0 13 14 15 16  4  5  8  7] 	closure(supp) = [0]
point = 1 	closure(cone) = [ 1 17 18 19 14  5  6  9  8] 	closure(supp) = [1]
point = 2 	closure(cone) = [ 2 15 20 21 22  7  8 11 10] 	closure(supp) = [2]
point = 3 	closure(cone) = [ 3 19 23 24 20  8  9 12 11] 	closure(supp) = [3]
point = 4 	closure(cone) = [4] 	closure(supp) = [ 4 13 16  0]
point = 5 	closure(cone) = [5] 	closure(supp) = [ 5 13 14 17  0  1]
point = 6 	closure(cone) = [6] 	closure(supp) = [ 6 17 18  1]
point = 7 	closure(cone) = [7] 	closure(supp) = [ 7 15 16 22  0  2]
point = 8 	closure(cone) = [8] 	closure(supp) = [ 8 14 15 19 20  0  1  2  3]
point = 9 	closure(cone) = [9] 	closure(supp) = [ 9 18 19 23  1  3]
point = 10 	closure(cone) = [10] 	closure(supp) = [10 21 22  2]
point = 11 	closure(cone) = [11] 	closure(supp) = [11 20 21 24  2  3]
point = 12 	closure(cone) = [12] 	closure(supp) = [12 23 24  3]
point = 13 	closure(cone) = [13  4  5] 	closure(supp) = [13  0]
point = 14 	closure(cone) = [14  5  8] 	closure(supp) = [14  0  1]
point = 15 	closure(cone) = [15  8  7] 	closure(supp) = [15  0  2]
point = 16 	closure(cone) = [16  7  4] 	closure(supp) = [16  0]
point = 17 	closure(cone) = [17  5  6] 	closure(supp) = [17  1]
point = 18 	closure(cone) = [18  6  9] 	closure(supp) = [18  1]
point = 19 	closure(cone) = [19  9  8] 	closure(supp) = [19  1  3]
point = 20 	closure(cone) = [20  8 11] 	closure(supp) = [20  2  3]
point = 21 	closure(cone) = [21 11 10] 	closure(supp) = [21  2]
point = 22 	closure(cone) = [22 10  7] 	closure(supp) = [22  2]
point = 23 	closure(cone) = [23  9 12] 	closure(supp) = [23  3]
point = 24 	closure(cone) = [24 12 11] 	closure(supp) = [24  3]



In [16]: