Data Structure of FEM in FEniCS

• Post by: Philip Luke K
• June 27, 2021
• 1 Comment

Storing the element connectivity, element nodes, nodal degrees of freedom etc.. are of significant importance in any numerical package which implements finite element analysis. And this constitutes the data structure of FEM. In this post, I will try to explain how FEniCS does this.

---

Before the construction/manufacturing stage, typical engineering involves 2 main stages; design and analysis. Design stage is concerned with developing the structural geometry through a Computer-Aided Design (CAD) environment. This geometry is then passed to the analysis stage where it is converted into an approximate analysis-suitable geometry, which is then meshed and analyzed using Finite Element Methods (FEM).

For simplicity, let us consider 2D triangular elements as shown below here onwards. Conversion of CAD design into analysis-suitable geometry may be achieved through commercial packages such as gmsh, Ansys, ABAQUS etc… depending upon the complexity of the geometry. In this stage, the geometry is discretized into elements and each element will have vertices. Nearby elements may or may not share some of these vertices.

This analysis-suitable-geometry needs to be meshed for analysing using FEM. This involves assigning nodes to the geometry based on the degree of element and assigning degrees of freedom (dof) to the nodes depending on the solution space. For example, linear element have nodes only at the vertices, whereas quadratic elements have nodes in between the vertices as well. Cubic elements have nodes inside the element. This is illustrated below.

• 

• 

• 

Let us understand this through a unit sized square mesh in FEnicS

from dolfin import *
mesh = UnitSquareMesh(1, 1)

which creates a 2-element unit-square mesh with 4 vertices as shown below. (FEniCS, by default, works with triangular and tetrahedral elements for 2D and 3D domains respectively) The first element has the vertices 0-1-3 whereas the second element will have the vertices 0-2-3. This can be seen using the command mesh.cells() . The command plot(mesh) visualizes the mesh. As natural to any coding language, the indexing starts from 0, which is used for the illustration as well. The corresponding cartesian coordinates of the vertices of a mesh can be obtained by using the command mesh.coordinates(). Some comparisons between various mesh entities and how they relate to the vertices can be found at Appendix A

mesh = UnitSquareMesh(1, 1)
plot(mesh)
coordinates = mesh.coordinates()
print("Vertex vs Coordinates:")
for i in range(len(coordinates)): 
    
    print(i, "\t", coordinates[i])

Now we will assign nodes and dof’s to this mesh and see how it changes. For the initial parametric study we will change the degree of the element in each iteration (linear, quadratic and cubic) while fixing the dof to be 1; i.e a scalar function space. This is achieved by the command V = FunctionSpace(mesh, "CG", p) where p is the degree of the element

for p in range(1, 4):
	V = FunctionSpace(mesh, "CG", p)
	dm = V.dofmap()
	print("Polynomial order: %d" % p)
	print("Vertex dofs: %s" % str(dm.entity_dofs(mesh, 0)))
	print("Facet dofs: %s" % str(dm.entity_dofs(mesh, 1)))
	print("Cell dofs: %s" % str(dm.entity_dofs(mesh, 2)))
	print("All DoFs (Vertex, Facet, and Cell) associated with cell 0: %s" % str(dm.cell_dofs(0)))
	print("******")

When the degree = 1 (linear element), FEniCS assigns nodes only to the vertices of the element (as illustrated before). The numbering of the nodes follow a standard FEniCS format and doesn’t follow the cell structure numbering as shown before. The results are shown below. There are 4 vertices and hence 4 nodes. Since dof/node = 1, there are 4 dof’s in total, all of which corresponds to the vertice nodes. There are no nodes in the edges or inside the cell, hence no dof is assigned there. For the first element (cell no: 0), the corresponding dof’s are 1, 3 and 2

Similarly when degree = 2 (quadratic element), there will 9 nodes in the mesh. In addition to 4 nodes in the vertices, an intermediate node will be present in the 5 edges, which are the facet dof’s.

When degree = 3 (cubic element), there will be 16 nodes in total. In addition to 4 nodes in the vertices, each edge will have 2 intermediate edges. Also one node will be present inside each cell. (4 + 2*4 + 2 = 16)

Next let us consider a case where the dof/node = 2, i.e each node can move in both the planar directions. In such a case, the solution space will be a vector. This may be achieved by just changing the command V = FunctionSpace(mesh, "CG", p) to V = VectorFunctionSpace(mesh, "CG", p)

So when degree = 1, there will be only 4 nodes in the mesh, all of them in the vertices. But each node will have 2 dof. Hence there will be a total of 8 dof in the mesh. FEniCS assigns this as follows

Notice that consecutive numbers are given to each dof of any node. Similarly when degree =2, there will be 9 nodes and thus 18 dof’s. Some of them will be the nodes on the edges.

Likewise, when degree=3, 16 nodes and thus 32 dof’s

In short, it is to be noted that the mesh function only creates vertices in FEniCS. The indices/numbering of nodes and the corresponding dof’s will be completely independent of the vertex numbering.

Appendix A

A triangular mesh consists of vertices, edges and cells, and that the edges may alternatively be referred to as facets and the cells as faces

A tetrahedral mesh consists of vertices, edges, faces and cells, and that the faces may alternatively be referred to as facets.

1. Edges vs Vertex

There are 5 edges in total for the unit square mesh we have defined above. To find out which vertices correspond to which edge, we use the following code snippet. The total number of edges in a mesh may be found by using mesh.num_edges()

edge = edges(mesh)
print("Edge vs Vertex")
for i in edge:
    print(i.index(), "\t", end =" ")
    for j in vertices(i):
        print(j.index(), end =" ")
    print("\n")

2. Triangle vs Vertex

For a 2D domain, FEniCS divides the mesh into triangular elements. To find the corresponding vertices for each triangular element, use :

print("Triangle vs Vertex")
cells = mesh.cells()
for i in range(len(cells)):
    print(i, "\t", cells[i])

3. Tetrahedra vs Vertex

3D domains are discretized into tetrahedra elements in FEniCS. For 3D, the mesh command itself changes

mesh2 = UnitCubeMesh(1, 1, 1)
plot(mesh2)
print("Tetrahedra vs Vertex")
cells2 = mesh2.cells()
for i in range(len(cells2)):
    print(i, "\t", cells2[i])

References

1. https://fenicsproject.discourse.group/t/problem-with-dof-to-vertex-map/1143/4
2. https://fenicsproject.org/pub/tutorial/html/._ftut1019.html