The Mesh Topology

The MeshTopology class provides topological information about the mesh, i.e. it gives access to its entities such as edges, faces or cells and neighborly relation, e.g. the neighboring cells to each edge.

All the information is stored using incidence matrices. So, for every pair \((d,d')\) of topological dimensions \(0 \leq d,d' \leq D\) where \(D \in \{ 1,2,3 \}\) is the spatial dimension of the mesh, there exists a matrix \(M_{d,d'}\) that for the \(d\)-dimensional entities stores their incident entities of topological dimension \(d'\). The entries \(m_{ij}\) are either \(1\) or \(0\) depending on whether the \(j\)-th \(d'\)-dimensional mesh entity is incident to the \(i\)-th mesh entity of topological dimension \(d\).

To explain how to get to these matrices we will consider the following simple sample mesh in 2D.

which consists of 4 vertices, 5 edges and 2 triangular cells. For this mesh the incidence relation \(D \to 0\), i.e. the vertices incident to each cell \((D=2)\), would be stored in a matrix

\[M_{2,0} = \begin{pmatrix} 1 & 1 & 1 & 0 \\ 0 & 1 & 1 & 1 \end{pmatrix}\]

This matrix shows that e.g. the first cell, which corresponds to the first row in the matrix, is incident to the vertices with indices \(1,2,3\), as those are the non-zero column indices which represent the mesh vertices.

Note that this mesh is very simple and coarse. If the meshes get finer the number of non-zero entries will become very small compared to the total number of entries, i.e. the incidence matrices will be sparse. To make use of this fact we will utilize the sparse matrix capabilities of the scipy package which provides sparse matrix classes that only store the non-zero entries, hence saving lots of memory.

Similarly to [Log09] the construction of these matrices relies on the three algorithms build, transpose and intersection, which will be explained in the following sections.

Initialize

A prerequisite for the algorithms to work is, that we already have the incidence matrix for the relation \(D \to 0\). This matrix is constructed using the following auxiliary intialization method.

First, we need the connectivity array that defines the mesh cells (entities of topological dimension \(D\)) row-wise via their node indices. For the sample mesh this array would be

\[\begin{bmatrix} 1 & 2 & 3 \\ 2 & 3 & 4 \end{bmatrix},\]

i.e. the first cell is define by the node indices \(1,2,3\) and the second cell is defined by the nodes \(2,3,4\). This array has to be supplied by the user and passed as an argument when creating an instance of the MeshTopology class.

In the end we will represent the incidence relation using a sparse matrix in coordinate format (COO) which stores its non-zero entries by their ‘cordinates’ \((i,j)\), i.e. their row and column indices. To do so, three 1-dimensional arrays are used. The first one is the data array which contains the values of the non-zero entries. In our case this will be an array of \(1\)‘s. The second array will contain the row indices and the third one the respective column indices.

To create the row index array we need to know the number of vertices that define the cells of the mesh. For a mesh of triangles this number is equal to \(nV = 3\) (in 3D for tetrahedra it would be \(nV = 4\)). Remember that we want to construct the incidence relation between the mesh cells and the mesh vertices, so if each cell is defined by \(nV\) vertices there will be that many non-zero entries in every row of the matrix. Hence, for each row index \(i\) (which corresponds to one cell in the mesh) there will be \(nV\) coordinate tuples \((i,j)\) that have this index as a first component. Therefore, the number of vertices tells us how often each cell index appears in the row array.

For our sample mesh the row index array would be:

\[rows = \begin{bmatrix} 1 & 1 & 1 & 2 & 2 & 2 \end{bmatrix},\]

containing each row index \(nV = 3\) times.

The columns of the incidence matrix we create correspond to the vertices of the mesh. So, each index in the column array is the column index of the non-zero entry corresponding to the respective index in the row array. As those non-zero entries represent the incident vertices for each cell the column indices are given by the connectivity array that has been passed as an argument.

Again, for our sample mesh this column index array would be:

\[cols = \begin{bmatrix} 1 & 2 & 3 & 2 & 3 & 4 \end{bmatrix}.\]

As stated before, the remaining data array contains the values of the non-zero entries in the matrix. For us those values are all equal to \(1\), so the data array would be

\[data = \begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 \end{bmatrix}.\]

and we can now create the incidence matrix for the relation \(D \to 0\) by passing the three array to the scipy class for a sparse matrix in COO format.

Build

The build algorithm computes the incidence relation \(d \to 0\) for a topological dimension \(0 < d < D\), i.e. it determines the incident vertices for every \(d\)-dimensional mesh entity.

Starting from the incidence relation \(D \to 0\) that was constructed using the initialization method above, we first have to determine the local vertex sets for every \(d\)-dimensional subentity of the mesh cells, i.e. we need for each subentity of topological dimension \(d\) of every mesh cell to vertex indices that define it.

For example, assume we want to build the incidence relation \(1 \to 0\), i.e. the incident vertices of every edge (topological dimension \(d = 1\)). As every triangular cell is incident to three edges and each edge is defined by two vertex indices, the local vertex set for the edges of every cell will consist of three 2-tuples of indices.

For our sample mesh the local vertex set for the edges of the first cell would be

\[\begin{bmatrix} (1,2) & (1,3) & (2,3) \end{bmatrix}\]

because the first edge of the cell is defined by the vertex indices \(1\) and \(2\), the second edge is defined by the indices \(1,3\) and the third edge by the indices \(2,3\). Analogously, the local vertex for the edges of the second cell would be

\[\begin{bmatrix} (2,3) & (2,4) & (3,4) \end{bmatrix}.\]

Obviously, there are duplicates in the local vertex sets because edges will be incident to two neighboring triangular cells if they are not part of the boundary of the domain. So, the next step is to determine the unique vertex tuples in all the local vertex sets. For us this would result in the following array, where the tuples are now written as the rows.

\[\begin{bmatrix} 1 & 2 \\ 1 & 3 \\ 2 & 3 \\ 2 & 4 \\ 3 & 4 \end{bmatrix}.\]

This array is now a connectivity array that defines the edges of our sample mesh just like the one we had in the initialization method for the cells. That means we can use the same procedure to construct the sparse matrix for the incidence relation \(d \to 0\).

As there are five edges in the mesh and each edge is incident to two vertices the row index array for the edges will contain each index from \(1\) to \(5\) twice.

\[rows = \begin{bmatrix} 1 & 1 & 2 & 2 & 3 & 3 & 4 & 4 & 5 & 5 \end{bmatrix}\]

The column index array just consists of the indices given by the connectivity array flattened row-wise

\[cols = \begin{bmatrix} 1 & 2 & 1 & 3 & 2 & 3 & 2 & 4 & 3 & 4 \end{bmatrix}\]

and the data array is an array of \(1\)‘s with the same length as the row and column index arrays.

\[data = \begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \end{bmatrix}\]

Transpose

The transpose algorithm computes the incidence relation \(d \to d'\) for \(d < d'\) from \(d' \to d\).

This is the simplest and most straight forward of the algorithms. Assuming we already have the incidence matrix \(M_{d',d}\) for the relation \(d' \to d\) we get the incidence matrix \(M_{d,d'}\) simply by transposing

\[M_{d,d'} = M_{d',d}^{T}\]

For example, we already know the incident vertices to every mesh edge

\[M_{1,0} = \begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 1 \end{pmatrix}.\]

If we now want to know the incident edges for every mesh vertex we just have to transpose this matrix to get

\[M_{1,0}^{T} = M_{0,1} = \begin{pmatrix} 1 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 \\ 0 & 1 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 & 1 \\ \end{pmatrix}.\]

where we can see for each of the four mesh vertices the edges that are incident to them.

Intersection

The intersection algorithm computes the incidence relation \(d \to d'\) for \(d \geq d'\) from the two relations \(d \to d''\) and \(d'' \to d'\).

The auxiliary topological dimension \(d''\) is set to \(0\), so we assume that we already have the incidence matrices \(M_{d,0}, M_{0,d'}\) for the incidence relations \(d \to 0\) and \(0 \to d'\).

Then, the intersection is computed as the dot product of these two matrices. The resulting intersection matrix \(I_{d,d'}\) describes for every \(d\)-dimensional mesh entity (represented by the rows) how many vertices (\(d'' = 0\)-dimensional entities) it shares with each of the mesh entities of topological dimension \(d\) (represented by the columns).

Next, we have to distinguish between two cases to define when the incidence relation \(d \to d'\) holds for two entities.

The case \(d = d'\)

In case \(d = d'\), e.g. if we want the neighboring (incident) cells for every mesh cell (\(2 \to 2\)), we define two entities of topological dimension \(d\) incident to each other if they share a \(d - 1 = \tilde{d}\)-dimensional mesh entity. So, for two triangles (\(d = 2\)) to be incident they must share an edge (\(\tilde{d} = 1\)).

Simplicial entities of dimension \(d\), i.e. triangles for \(d = 2\) or tetrahedra for \(d = 3\), which are the only type of entities that currently are supported by Py\(\small\mathbb{SOFE}\), are defined by \(d + 1\) vertices. So, as the computed intersection matrix \(I_{d,d'}\) tells us the number of vertices shared between each of the \(d\)-dimensional entities and every entity of topological dimension \(d'\) the incidence matrix \(M_{d,d'}\) for \(d = d'\) is defined by

\[(M_{d,d})_{ij} = m_{ij} = \begin{cases} 1 & (I_{d,d})_{ij} = d \\ 0 & \text{else} \end{cases}\]

Suppose, we want to know the neighboring edges for every edge of our sample mesh, i.e. we are interested in the incidence relation \(1 \to 1\). The necessary incidence matrices for computing the intersection are

\[M_{1,0} = \begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 1 \end{pmatrix} \quad M_{0,1} = \begin{pmatrix} 1 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 \\ 0 & 1 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 & 1 \\ \end{pmatrix}.\]

Then the intersection matrix is

\[M_{1,0} \cdot M_{0,1} = I_{1,1} = \begin{pmatrix} 2 & 1 & 1 & 1 & 0 \\ 1 & 2 & 1 & 0 & 1 \\ 1 & 1 & 2 & 1 & 1 \\ 1 & 0 & 1 & 2 & 1 \\ 0 & 1 & 1 & 1 & 2 \end{pmatrix}\]

and therefore the resulting incidence matrix would be

\[M_{1,1} = \begin{pmatrix} 0 & 1 & 1 & 1 & 0 \\ 1 & 0 & 1 & 0 & 1 \\ 1 & 1 & 0 & 1 & 1 \\ 1 & 0 & 1 & 0 & 1 \\ 0 & 1 & 1 & 1 & 0 \end{pmatrix}.\]

A special case is the relation \(0 \to 0\) for which we define that vertices are incident to themselves only. So the incidence matrix \(M_{0,0}\) would be the identity matrix.

The case \(d > d'\)

In case \(d > d'\) we define a entity of topological dimension \(d'\) to be incident to a \(d\)-dimensional one if all its vertices are incident to this entity of dimension \(d\).

Since the \(d'\)-dimensional mesh entities are defined by \(d' + 1\) vertices the incidence matrix \(M_{d,d'}\) is defined similarly to the previous section via the intersection matrix

\[(M_{d,d'})_{ij} = m_{ij} = \begin{cases} 1 & (I_{d,d'})_{ij} = d' + 1 \\ 0 & \text{else} \end{cases}\]

Suppose, we want to know the incident edges for every cell of our sample mesh, i.e. we are interested in the incidence relation \(2 \to 1\). The necessary incidence matrices for computing the intersection are

\[M_{2,0} = \begin{pmatrix} 1 & 1 & 1 & 0 \\ 0 & 1 & 1 & 1 \end{pmatrix} \quad M_{0,1} = \begin{pmatrix} 1 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 \\ 0 & 1 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 & 1 \\ \end{pmatrix}.\]

Then the intersection matrix is

\[M_{2,0} \cdot M_{0,1} = I_{2,1} = \begin{pmatrix} 2 & 2 & 2 & 1 & 1 \\ 1 & 1 & 2 & 2 & 2 \end{pmatrix}\]

and therefore the resulting incidence matrix would be

\[M_{2,1} = \begin{pmatrix} 1 & 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 & 1 \end{pmatrix}.\]

Computing Incidence Relations

Now that we have the algorithms build, transpose and intersection we can construct the incidence matrix for any relation \(d \to d'\) by a combination of those algorithms as follows

Compute relation (d -> d')
--------------------------

if relation (d -> 0) does not exist
   build relation (d -> 0)

if relation (d' -> 0) does not exist
   build relation (d' -> 0)

if d < d'
   compute relation (d' -> d)
   transpose relation (d' -> d)
else
   intersect relation (d -> dd)