Lecture 7(a)

Finite elements in two dimensions

Unstructured grids

• Most problems have complex domains.
• Use unstructured grids.
• Use Finite Elements.
• Higher dimensions increases practical complexity.

We have spent a lot of time developing a (seemingly) quite complex method in only one dimension. The point of all the effort put in is that the conceptual difficulty of moving to unstructured grids in higher dimensions is small. The practical difficulties can be somewhat larger, however.

Reference triangle

• Split unit square into two triangles.
• Each element maps to reference triangle.
• Compute local \(K^{(e)}, F^{(e)}\) on reference element.
• Assemble.
• Solve.

In the one dimensional case, we took each element and mapped it to the reference element \(x \in [-1, 1]\). All integrals were done on the reference element.

In two dimensions, we map to a reference triangle or a reference square. Here we will use triangles. This mapping involves a translation and a scaling (the same as in one dimension), but also potentially a rotation. This increases practical complexity.

Note that the picture includes the mapping of the two local nodes with local node number 1. These are mapped to the same point in the reference triangle, starting from different points. This is the standard case.

Local calculations

Follow 1d case, need local arrays.

Stiffness matrix

\[ k^e_{ab} = \int_{\triangle^e} \text{d}\triangle^e \, \left( \partial_{x} N_a (\mathbf{x}) \partial_{x} N_b(\mathbf{x}) + \partial_{y} N_a(\mathbf{x}) \partial_{y} N_b(\mathbf{x}) \right). \]

Force vector

\[ f^e_b = \int_{\triangle^e} \text{d}\triangle^e \, N_b(\mathbf{x}) S(\mathbf{x}). \]

In the one dimensional case the computation over the reference element simplified the general code structure. That is going to happen massively again here. However, we need to know how to map to the reference element and compute the multi-dimensional integral.

Integrate over elements

\[ \begin{aligned} \int_\triangle \mathrm{d} \triangle \phi(x, y) &= \int_{0}^1 \mathrm{d}\xi_2 \int_0^{1-\xi_2} \mathrm{d}\xi_1 \phi(\xi_1, \xi_2) j(\xi_1, \xi_2) \, , \\ j = \det \left[ \frac{\partial \mathbf{x}}{\partial \symbf{\xi}} \right] &= \det \begin{pmatrix} \partial_{\xi_1} x & \partial_{\xi_2} x \\ \partial_{\xi_1} y & \partial_{\xi_2} y \end{pmatrix} \, . \end{aligned} \]

Gauss quadrature on reference triangle:

\[ \begin{split} \int_{0}^1 \mathrm{d}\xi_2 \int_0^{1-\xi_2} \mathrm{d}\xi_1 \psi(\symbf{\xi}) \simeq \tfrac{1}{6} \sum_{j=1}^3 \psi(\mathbf{x}(\symbf{\xi}^{(j)})) \, , \\ \symbf{\xi}^{(1)} = \tfrac{1}{6} \begin{pmatrix} 1 & 1 \end{pmatrix} \, , \quad \symbf{\xi}^{(2)} = \tfrac{1}{6} \begin{pmatrix} 1 & 4 \end{pmatrix} \, , \quad \symbf{\xi}^{(3)} = \tfrac{1}{6} \begin{pmatrix} 4 & 1 \end{pmatrix} \, . \end{split} \]

The key thing we need to do in the finite element method is to integrate over the reference element. Here, the reference element is the triangle. We need to map the quantity to the triangle. This change of coordinates means that the area integral over the triangle picks up a term related to the determinant of the coordinate transformation. The coordinate transformation is directly linked to the map. However, this formula will hold in more general cases.

Having performed the mapping, we can then integrate over the reference triangle using Gauss quadrature again. In this case we will not use the Gauss-Lobatto formula that evaluates at the nodes, because that is not very accurate. Instead we will use Gauss-Legendre quadrature, which gives high enough accuracy using only three evaluation points. However, those points are within the reference triangle. To evaluate an explicit function of space, such as the source function, we will need to be able to map the reference coordinates back to the original space to get these values.

Maps

• Coordinates are scalar fields.
• So use function basis:

\[ x = x^e_0 N_0(\symbf{\xi}) + x^e_1 N_1(\symbf{\xi}) + x^e_2 N_2(\symbf{\xi}) \, . \]

Linear basis on reference triangle

\[ \begin{split} N_0 = 1 - \xi_1 - \xi_2, \quad N_1 = \xi_1, \quad N_2 = \xi_2 \\ \implies \qquad \partial_{\xi_1} x = -x^e_0 + x^e_1 \, . \end{split} \]

To get the maps, we represent the coordinates in terms of the function basis. We can always do this by thinking of the coordinates themselves as scalar fields. Therefore the nodal coordinates of the element in the original coordinates can be directly written in terms of the basis functions, and the basis functions are simple linear functions of the reference coordinates. We can then compute the derivatives in order to get the explicit coordinate transformation. We can also evaluate the basis functions at reference coordinates inside the reference element to find the coordinate locations for integration.

Gridding

Now need four things:

1. List of node locations nodes (\(\mathtt{i} \to \mathbf{x}\));
2. List of nodes in elements IEN (\(\mathtt{e} \to \mathtt{i}_0, \mathtt{i}_1, \mathtt{i}_2\));
3. List of boundary nodes with values, type;
4. Location matrix LM linking local/global numbering.

Node numbers arbitrary: add all nodes not on Dirichlet boundaries to ID.

\[ \mathtt{LM}(\mathtt{a}, \mathtt{e}) = \mathtt{ID}(\mathtt{IEN}(\mathtt{e}, \mathtt{a})) \, . \]

We now have nearly everything that we need. We know that the local stiffness matrix and force vectors can be computed by integrals over the reference element. Given the global nodal coordinates of an element in the original coordinates we know how to map this to reference space. However, we still need to know how to link the local stiffness matrix, for example, to the global entries. This requires a bit more machinery than in the one dimensional case.

We need a global node number. We saw that in the original picture. It can be chosen arbitrarily.

We need a global element number, which is also arbitrary, and we need to know the three global nodes that are in this element. This information defines the grid. If building a grid by hand you need to do this consistently; if using gridding software, which is recommended, you need it to supply this information.

We need to know which of the global nodes lie on the boundary. Again, the gridding software should tell us this. From that, you can check which nodes on the boundary have which type.

Finally, we need to link the global node number to the appropriate row in the matrix, through the location matrix, or LM array, as in the one dimensional case.

Fortunately, the location matrix can be given from the rest of the information. We loop over the global node numbers. If it is a Dirichlet node then its LM entry is -1 as it is already fixed; otherwise we increment the number by one. This gives the ID array, for “integer destination”. We can then set the location matrix using the ID array and the array linking element numbers to node numbers. This can be set up in advance.

Diffusion with sources

nodes = np.loadtxt('data/ews_nodes_s5k_lc50k.txt')
IEN = np.loadtxt('data/ews_IEN_s5k_lc50k.txt', 
                 dtype=np.int64)
boundary_nodes = np.loadtxt('data/ews_bdry_s5k_lc50k.txt', 
                            dtype=np.int64)
# Make all boundary points Dirichlet
ID = np.zeros(len(nodes), dtype=np.int64)
n_eq = 0
for i in range(len(nodes[:, 1])):
    if i in boundary_nodes:
        ID[i] = -1
    else:
        ID[i] = n_eq
        n_eq += 1
N_equations = np.max(ID)+1
N_elements = IEN.shape[0]
N_nodes = nodes.shape[0]
N_dim = nodes.shape[1]
# Location matrix
LM = np.zeros_like(IEN.T)
for e in range(N_elements):
    for a in range(3):
        LM[a,e] = ID[IEN[e,a]]
# Global stiffness matrix and force vector
K = np.zeros((N_equations, N_equations))
F = np.zeros((N_equations,))
# Loop over elements
for e in range(N_elements):
    k_e = stiffness_2d(nodes[IEN[e,:],:])
    f_e = force_2d(nodes[IEN[e,:],:], S)
    for a in range(3):
        A = LM[a, e]
        for b in range(3):
            B = LM[b, e]
            if (A >= 0) and (B >= 0):
                K[A, B] += k_e[a, b]
        if (A >= 0):
            F[A] += f_e[a]
# Solve
Psi_interior = np.linalg.solve(K, F)
Psi_A = np.zeros(N_nodes)
for n in range(N_nodes):
    if ID[n] >= 0: # Otherwise Psi should be zero, and we've initialized that already.
        Psi_A[n] = Psi_interior[ID[n]]

Here is a simplified two dimensional finite element code. It uses the complex unstructured grid, but the boundary conditions are homogeneous Dirichlet boundaries. The source term is piecewise constant and concentrated over the major cities.

The first part of the code sets up the grid. The three key pieces of information - node locations, element to node number mapping, and which nodes lie on the boundary - have all been produced by gridding software and saved to files. Here I used gmsh to construct the grid. Having loaded this data the ID array works out which nodes will have rows in the final matrix, which is all those not having Dirichlet boundary conditions. In this case that is all the interior nodes. From the ID array we construct the location matrix for use in the next step.

The second part of the code should look nearly identical to the one dimensional case. It loops over all elements, computing the local stiffness matrix and force vector, and using the location matrix to update the global arrays. The key difference will be that the local arrays are now size 3, not 2, so the loop ranges are modified. In this case the boundary conditions are trivial so no modifications are done.

Finally the solution is computed by solving the linear system, and the boundary terms updated. Again, this is identical to the one dimensional case.

The major change to the one dimensional case would appear to be the setup of the grid, and indeed that is substantial. However, the real major change is in the functions that compute the local stiffness matrix and force vector. In one dimension these are simple one line functions. In higher dimensions they are more complex, involving 5-10 steps. Each step is individually simple enough, but they all need careful testing.

Summary

• Finite elements use function basis.
• Extension to higher dimensions requires
  • mapping to reference triangle/tetrahedron;
  • multi-point Gauss quadrature;
  • more complex gridding.
• Use gridding software.
• Core of algorithm conceptually identical.

The key step that needs updating for multiple dimensions is the computation of the local arrays. When implementing this, test by integrating known functions over given triangles, starting with the simplest cases and adding complexity. Once this can be done, the rest follows the one dimensional case nearly identically.

The main exception is the construction of the grid. Building a good grid is hard. Where possible use existing gridding software to do this for you.