we can therefore propose the reconstruction function g˜Δ si115_e as a function of the Lagrangian coordinates ξA,ξB si116_e, and ξC si117_e using Eq.1 as

g˜ΔξAξBξC=gΔAξA+gΔBξB+gΔCξC

si118_e  (Eq. 32.52)

where A, B, and C are the vertices of the triangle. We can then calculate the derivatives as

∂g˜Δ∂x=gΔA∂ξA∂x+gΔB∂ξB∂x+gΔC∂ξC∂x

si119_e

∂g˜Δ∂y=gΔA∂ξA∂y+gΔB∂ξB∂y+gΔC∂ξC∂y

si120_e

Using Eq. 32.5, Eq. 32.6, and Eq. 32.7 we can note down the partial derivatives of the Lagrangian coordinates as

∂ξA∂x=yB−yCxAyB−yC+xByC−yA+xCyA−yB∂ξB∂x=yC−yAxAyB−yC+xByC−yA+xCyA−yB∂ξC∂x=yA−yBxAyB−yC+xByC−yA+xCyA−yB∂ξA∂y=xC−xBxAyB−yC+xByC−yA+xCyA−yB∂ξB∂y=xA−xCxAyB−yC+xByC−yA+xCyA−yB∂ξC∂y=xB−xAxAyB−yC+xByC−yA+xCyA−yB

si121_e

Basis Functions. The basis functions Ψj are the pyramid function which is 1 at the center node j and zero on the neighboring nodes. Within each of the triangles that make up the pyramid, it is described by Eq. 32.1

ψΔξA=ξA

si122_e  (Eq. 32.53)

and is only dependent on one Lagrangian coordinate ξA. The partial derivatives required for Eq. 32.51 are thus given by

∂ψΔ∂x=∂ξA∂x

si123_e  (Eq. 32.54)

∂ψΔ∂y=∂ξA∂y

si124_e  (Eq. 32.55)

Integral Contribution of Each Triangle. Having derived Eq. 32.52, Eq. 32.54, and Eq. 32.55 we can now rewrite Eq. 32.51 for each triangle Δ as

−∫ΩΔ∂ψΔ∂xgΔA∂ξA∂x+gΔB∂ξB∂x+gΔC∂ξC∂x+∂ψΔ∂ygΔA∂ξA∂y+gΔB∂ξB∂y+gΔC∂ξC∂ydΩ=∫ΩΔρxyψΔdΩ

si125_e

which we rewrite to

−1xAyB−yC+xByC−yA+xCyA−yB2∫ΩΔgΔAyB−yC2+xC−xB2+gΔByB−yCyC−yA+xA−xCxC−xB+gΔCyB−yCyA−yB+xB−xAxC−xBdΩ=∫ΩΔρxyξAdΩ

si126_e  (Eq. 32.56)

We are now still left with the problem that the integration over the domain ΩΔ is somewhat difficult to carry out since the domain is a triangle. However, if you look very closely at Eq. 32.56 you will note one thing: All the terms within the integral are constants. We can therefore rewrite the left-hand side as

−1xAyB−yC+xByC−yA+xCyA−yB2gΔAyB−yC2+xC−xB2+gΔByB−yCyC−yA+xA−xCxC−xB+gΔCyB−yCyA−yB+xB−xAxC−xB∫ΩΔdΩ=∫ΩΔρxyξAdΩ

si127_e  (Eq. 32.57)

The integral is therefore nothing else than the constant multiplied with the area of the triangle. If you recall from Eq. 32.5, Eq. 32.6, and Eq. 32.7, the denominators of ξA, ξB, ξC, and are simply twice the area of the triangle; therefore the area of the triangle is

AΔ=xAyB−yC+xByC−yA+xCyA−yB2

si128_e

Which allows us to evaluate the integral of the left-hand side of Eq. 32.57:

−12xAyB−yC+xByC−yA+xCyA−yBgΔAyB−yC2+xC−xB2+gΔByB−yCyC−yA+xA−xCxC−xB+gΔCyB−yCyA−yB+xB−xAxC−xB=∫ΩΔρxyξAdΩ

si129_e  (Eq. 32.58)

In the next step, we use the same trick to evaluate the integral on the right-hand side by simply assuming that ρ(x, y) is constant within the triangle and given by ρ(Δ). This allows us to rewrite the right-hand side of Eq. 32.58 to

∫ΩΔρxyξAdΩ=ρΔ∫ΩΔξAdΩ

si130_e

where ∫ΩΔξAdΩ si131_e is nothing else than the volume of the pyramid that is given by

VΔ=13AΔh=13AΔ=xAyB−yC+xByC−yA+xCyA−yB6

si132_e

which we rewrite to

−1xAyB−yC+xByC−yA+xCyA−yB2∫ΩΔgΔAyB−yC2+xC−xB2+gΔByB−yCyC−yA+xA−xCxC−xB+gΔCyB−yCyA−yB+xB−xAxC−xBdΩ=∫ΩΔρxyξAdΩ

si126_e  (Eq. 32.56)

We are now still left with the problem that the integration over the domain ΩΔ is somewhat difficult to carry out since the domain is a triangle. However, if you look very closely at Eq. 32.56 you will note one thing: All the terms within the integral are constants. We can therefore rewrite the left-hand side as

−1xAyB−yC+xByC−yA+xCyA−yB2gΔAyB−yC2+xC−xB2+gΔByB−yCyC−yA+xA−xCxC−xB+gΔCyB−yCyA−yB+xB−xAxC−xB∫ΩΔdΩ=∫ΩΔρxyξAdΩ

si127_e  (Eq. 32.57)

The integral is therefore nothing else than the constant multiplied with the area of the triangle. If you recall from Eq. 32.5, Eq. 32.6, and Eq. 32.7, the denominators of ξA, ξB, ξC, and are simply twice the area of the triangle; therefore the area of the triangle is

AΔ=xAyB−yC+xByC−yA+xCyA−yB2

si128_e

Which allows us to evaluate the integral of the left-hand side of Eq. 32.57:

−12xAyB−yC+xByC−yA+xCyA−yBgΔAyB−yC2+xC−xB2+gΔByB−yCyC−yA+xA−xCxC−xB+gΔCyB−yCyA−yB+xB−xAxC−xB=∫ΩΔρxyξAdΩ

si129_e  (Eq. 32.58)

In the next step, we use the same trick to evaluate the integral on the right-hand side by simply assuming that ρ(x, y) is constant within the triangle and given by ρ(Δ). This allows us to rewrite the right-hand side of Eq. 32.58 to

∫ΩΔρxyξAdΩ=ρΔ∫ΩΔξAdΩ

si130_e

where ∫ΩΔξAdΩ si131_e is nothing else than the volume of the pyramid that is given by

VΔ=13AΔh=13AΔ=xAyB−yC+xByC−yA+xCyA−yB6

si132_e

where we set h = 1 because the height of the pyramid will be scaled by a prefactor (that we still need to determine). We can now simplify Eq. 32.58 to

12xAyB−yC+xByC−yA+xCyA−yBgΔAyB−yC2+xC−xB2+gΔByB−yCyC−yA+xA−xCxC−xB+gΔCyB−yCyA−yB+xB−xAxC−xB=−ρΔxAyB−yC+xByC−yA+xCyA−yB6

si140_e  (Eq. 32.59)

Eq. 32.59 is often referred to as the contribution equation because it describes the contribution of each triangle Δ to the overall integral in the current element j, i.e., the current pyramid.

Integral Form for Each Element. Having derived the integral contribution for each triangle we can now assemble the integral form for each element. For each triangle around the node j the contribution of all triangles is evaluated using Eq. 32.59. If you look closely you will see that Eq. 32.59 is a linear equation of the unknown values g(Δ, A), g(Δ, B), and g(Δ, C). Each triangle will thus contribute to the stiffness matrix of the specific node as well as to the load vector via the right-hand side of Eq. 32.59. It merely remains to cleverly sum up these contributions to obtain the overall stiffness matrix and the load vector upon which a linear system of equation is obtained, which can be simply solved for the unknown values of the nodes.

Returning to our mesh of Fig. 32.6 we can see that node 1 is surrounded by a total of 12 triangles. For each of these triangles, we evaluate Eq. 32.59 and thus obtain values to add to the stiffness matrix and the load vector. We proceed by doing this for all nodes for which we do not know the values.

32.6.4 Maple Worksheet

Obviously, FEM is not a method one would like to perform manually although it is possible. In the following we will implement a fully function two-dimensional FEM solver for Poisson equations of the type given in Eq. 32.50. As you will see, this implementation is generic and it can be used to solve, e.g., the heat equation, Eq. 8.2, or the Poisson equation for Poisson flow, which is the example we will use here.

We will split the implementation across several listings starting with the mesh definition.

Mesh Definition. Obviously the first thing we require is the implementation of the mesh. At this point, it is not really important which type of differential equation we are actually solving because FEM requires a mesh in all cases. Listing 32.2 shows the way we input the mesh information taken from Fig. 32.6. If a different domain is chosen or the mesh is changed, this will be the only input to the solver that needs to be adapted. All subsequent implementation steps will be identical.

Listing 32.2

[Maple] Mesh definition for two-dimensional FEM. This listing is part of the listing FEM2D.mw, which can be found in the supplementary material.

1 restart:read "Core.txt":read "CoreFunctions.txt":interface(rtablesize = infinity):
2
3 nodes:=Matrix ([[0,0],[0.125,−0.25],[0.25,−0.125],[0.25,0],[0.25,0.125],[0.125,0.25],[0,0.25],[− 0.125,0.25],[− 0.25,0.125],[− 0.25,0],[− 0.25,−0.125],[− 0.125,−0.25],[0,−0.25],[0.25,−0.25],[0.25,0.25],[− 0.25,0.25],[− 0.25,−0.25],[0,−0.5],[0.25,−0.43],[0.43,−0.25],[0.5,0],[0.43,0.25],[0.25,0.43],[0,0.5],[− 0.25,0.43],[− 0.43,0.25],[− 0.5,0],[− 0.43,−0.25], [− 0.25,−0.43]]);
4 neighbors:=Matrix([[2,3,4,5,6,7,8,9,10,11,12,13],[3,1,13,18,19,14],[4,1,2,14,20,21],[5,1,3,21],[15,6,1,4,21,22],[7,1,5,15,23,24],[8,1,6,24],[9,1,7,24,25,16],[10,1,8,16,26,27],[11,1,9,27],[12,1,10,27,28,17],[13,1,11,17,29,18],[2,1,12,18],[3,2,19,20],[23,6,5,22],[25,26,9,8],[11,28,29,12],[29,12,13,2,19],[18,2,14,20],[19,14,3,21],[20,3,4,5,22],[21,5,15,23],[22,15,6,24],[23,6,7,8,25],[24,8,16,26],[25,16,9,27],[26,9,10,11,28],[27,11,17,29],[28,17,12,18]]);
5 boundaryvalues:=Matrix ([[18,0],[19,0],[20,0],[21,0],[22,0],[23,0],[24,0],[25,0],[26,0],[27,0],[28,0],[29,0]]);

We first start by restarting the Maple server and including Core.txt and CoreFunctions.txt. We also change the interface setting to output matrices completely instead of only displaying a shortened notation. This will allow us to check more easily if we made a mistake during input.

We then proceed defining three matrices. The first matrix is termed nodes and contains the x-coordinates and the y-coordinates of the grid we are using. For our example, we can directly deduce these values from Fig. 32.6. We input the nodes in the numbering order we chose for Fig. 32.6. The order is important since the algorithm expects the nth node to be found in the nth row of the matrix nodes.

Next we define the matrix neighbors. In this matrix, we enter the neighbors of node n as a list. Similarly to nodes, the algorithm expects the neighbors of the nth node to be found at the nth position of the matrix neighbors. Therefore, the first line of the matrix contains the neighbors of node 1, which are nodes 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, and 13. Due to the way Maple handles matrices, the matrix neighbors will be expanded to a width suitable to hold the maximum number of neighbors. In this case, this is the first node, i.e., the matrix will have a total of 12 columns. For each node that has fewer neighbors, the nonrequired entries are left at 0. As the numbering in our mesh starts with 1, a zero in this matrix indicates that the there are no more neighbors for this node. We will use this fact to later check the total number of neighbors programmatically. Please note that a cleaner implementation could use a reference pointer to arrays of different lengths instead of storing the data in one matrix. However, this is simply a question of convenience and does not alter the performance of the code.

One important thing to keep in mind when assigning the neighbors is that the neighbors should be given such that they can be “walked along.” This means that if we moved along the mesh visiting the neighbors one-by-one we would circle around the node for which we query the neighbors. As an example, if we follow the order 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, and 13 around the first node, we actually make a complete tour around it. This is important since we will see later derive the subtriangles associated with this mesh point. If the neighbors are given in a way that the associated subtriangles can be derived simply by always taking two subsequent points (which when combined with the respective point will form the subtriangle) this will significantly simplify the process. As an example, the subtriangles associated with the first mesh points are the triangles (1,2,3), (1,3,4), and all the way up to (1,13, 2). Here we always take two subsequent points moving along the list of neighbors until we reach the end. If we give the neighboring points in a different order, it is more difficult to derive the subtriangles.

The last matrix we define is boundaryvalues, which is simply a list of values for which we know the actual values. These values are the boundary values, and obviously, we do not need to solve the mesh for them since these values are known. In our example, we know the values for all the mesh points 18 to 29. At each of these points, the value of the dependent function is zero due to the no-slip boundary condition. Obviously, the fact that the boundary values are bundled is simply due to the fact that the mesh was hand-designed, and thus the outmost points received the highest node numbers. However, this numbering could just as well be different and the (known) boundary values could be somewhere in the middle.

These three matrices are required for defining the mesh. We now need a couple of helper functions before implementing the FEM scheme.

Helper Functions.Listing 32.3 implements two helper functions that we require for FEM. The first is the function GetNumberOfNeighbors (see line 1). It expects the matrix neighbors defined in listing 32.2 and the index j of the node for which the number of neighbors must be queried. It first checks for the number of columns in the matrix neighbors, which is the maximum number of neighbors possible. It then checks the column entries of row j for the first value which is zero. If a zero is found the last column count is returned (see line 5). If no zero is found, the maximum number of neighbors is returned (see line 8).

Listing 32.3

[Maple] Helper functions for two-dimensional FEM. This listing is part of the listing FEM2D.mw, which can be found in the supplementary material.

 1 GetNumberOfNeighbors:=proc(j,neighbors) local nmax,i:
 2  nmax:=Size(neighbors,2):
 3  for i from 1 to nmax do
 4   if neighbors[j,i] = 0 then
 5    return i − 1:
 6   end if
 7  end do:
 8  return nmax:
 9 end proc:
10
11 GetSubtriangles:=proc(j,neighbors) local neighborscount,triangles,i:
12  neighborscount:=GetNumberOfNeighbors(j,neighbors):
13  triangles:=Matrix(neighborscount,3):
14  for i to neighborscount − 1 do
15   triangles[i,1]:=j:
16   triangles[i,2]:=neighbors[j,i]:
17   triangles[i,3]:=neighbors[j,i + 1]:
18  end do:
19  triangles[neighborscount,1]:=j:
20  triangles[neighborscount,2]:=neighbors[j,neighborscount]:
21  triangles[neighborscount,3]:=neighbors[j,1]:
22  return triangles;
23 end proc:

The second function we will implement is the function GetSubtriangles (see line 11). It accepts the index j of the node for which the subtriangles are to be returned. It also requires the matrix neighbors that we defined in listing 32.2. The number of subtriangles is the same as the number of neighbors. We then create a matrix with a line for each subtriangle. The matrix has three columns holding the number of each of the three nodes for the subtriangle. Obviously, the first entry always is the index j. As explained earlier, the subtriangles are then found by always taking two subsequent entries of the jth row in neighbors (see line 14). We can do this for all subtriangles but the very last one. For this subtriangle we have to take the last and the first entry in the array neighbors (see line 19). The function returns the subtriangle matrix triangles.

Initialize. As the next step, we implement a function that will prepare the stiffness matrix, the load vector, and the mesh information for the actual FEM implementation. This function is called Initialize and is shown in listing 32.4. This function requires the mesh definition in the form of the matrices nodes, neighbors and boundaryvalues. It then creates the stiffness matrix stiffm, the load vector loadv, and an array containing the condensed mesh information in nodesinfo.

Listing 32.4

[Maple] Function Initialize required for two-dimensional FEM. This listing is part of the listing FEM2D.mw, which can be found in the supplementary material.

 1 Initialize:=proc(nodes,neighbors,boundaryvalues,stiffm,loadv,nodesinfo)local n,nboundaries,currentindex,i,j:
 2  n:=Size(nodes,1);
 3  nboundaries:=Size(boundaryvalues,1):
 4  nodesinfo:=Matrix(n,6);
 5  stiffm:=Matrix(n–nboundaries,n–nboundaries):
 6  loadv:=Vector(1..n–nboundaries):
 7  currentindex:=1:
 8  for i to n do
 9   for j to nboundaries do
10    if boundaryvalues[j,1] = i then
11     nodesinfo[i,1]:=1:
12     nodesinfo[i,2]:=boundaryvalues[j,2]:
13    end if:
14   end do:
15   if nodesinfo[i,1] = 0 then
16    nodesinfo[i,2]:=currentindex:
17    currentindex:=currentindex + 1:
18   end if:
19   nodesinfo[i,3]:=nodes[i,1]:
20   nodesinfo[i,4]:=nodes[i,2]:
21   nodesinfo[i,5]:=GetNumberOfNeighbors(i,neighbors):
22   nodesinfo[i,6]:=GetSubtriangles(i,neighbors):
23  end do:
24  return;
25 end proc:

The function begins by checking the number of nodes in the mesh (see line 2). This is simply the number of rows in the matrix nodes. It then checks the number of boundary values given by checking the number of rows in boundaryvalues. Obviously, the number of unknowns in our mesh is simply the difference of the total number of nodes and the number of nodes for which we know the actual values, i.e., the boundary values. We then create the stiffness matrix stiffm and the load vector loadv suitable to hold the required number of unknowns (see line 5).

We also initialize the matrix nodesinfo. This matrix will then be filled with the condensed information of the three matrices nodes, neighbors, and boundaryvalues. The matrix is designed such that the information for node i is found in the ith entry of the matrix. The matrix has a total of 6 columns that contain the following information.

Column 1 will be set to 0 if the value for the respective node is not known and must be solved for. In this case, this node must be linked to a row of the stiffness matrix stiffm and the load vector loadv. This reference is set in the second column. If, on the other hand, the value of this node is known (because this known is a boundary value) this entry is set to 1.

In case node j is a boundary value, column 2 will hold the actual value. If node i is not known, this entry is a reference to the stiffness matrix stiffm and the load vector loadv. Recall that both the stiffness matrix and the load vector have fewer entries than the actual number of nodes in the mesh if some of the nodes are boundary values. As the boundary values can occur at any location of the mesh, the information for node j in the stiffness matrix may be located at a different row than row i. As an example, if we query this entry for node i = 3 and find the value 5 in the second column of nodesinfo we find the information for node i = 3 at the row 5 of the stiffness matrix and the load vector.

The way that this assignment is made is shown in line 8. The array boundaryvalues is checked for any entry for node i in the first column. If an entry is found, the value in the first column is set to 1 (see line 11) and the value of the second column is set to the value taken from boundaryvalues (see line 12). If the node i is not found as an entry in boundaryvalues, we leave the first column at 0 indicating that this node will be solved for. We then assign the value currentindex to the second column (see line 16). currentindex holds the next free reference in the stiffness matrix stiffm and the load vector loadv. This value is initialized to 1 (see line 7). After each use the value is incremented (see line 17).

The third and fourth columns in the matrix nodesinfo hold the x- and the y-coordinates, respectively, which are simply copied from nodes (see line 19 and line 20). The fifth column holds the number of neighbors, which is returned by the helper function GetNumberOfNeighbors (see line 21). The sixth column holds the subtriangles associated with node i, which are returned by the helper function GetSubtriangles (see line 22).

The function Initialize will assemble all required information from our mesh into nodesinfo and create the stiffness matrix stiffm and the load vector loadv for the subsequent calculation. We now proceed by implementing the actual FEM scheme.

Subtriangle Contribution.listing 32.5 implements the function FEM2DTriangle that calculates the individual triangle contribution using Eq. 32.59. For a given node j we need to sum up all contributions for all subtriangles. For each subtriangle we will call DoTriangleContribution, indicating the current central node j and the indices of the mesh points indexA, indexB, and indexC, which make up the subtriangle. We also require the value rho that forms the right-hand side of the differential equation. The function DoTriangleContribution will then update the stiffness matrix stiffm and the load vector loadv.

Listing 32.5

[Maple] Single triangle contribution for two-dimensional FEM. This listing is part of the listing FEM2D.mw, which can be found in the supplementary material.

 1 DoTriangleContribution:=proc(j,indexA,indexB,indexC,nodesinfo,rho,stiffm,loadv) local xA,yA,xB,yB, 2xC,yC,gA,gB,gC,denominator:
 2  xA:=nodesinfo[indexA,3]:yA:=nodesinfo[indexA,4]:
 3  xB:=nodesinfo[indexB,3]:yB:=nodesinfo[indexB,4]:
 4  xC:=nodesinfo[indexC,3]:yC:=nodesinfo[indexC,4]:
 5  denominator:=2*(xA*(yB − yC) + xB*(yC − yA) + xC*(yA − yB));
 6  gA:=((yB − yC)^2 + (xC − xB)^2)/denominator:
 7  gB:=((yB − yC)*(yC − yA) + (xA − xC)*(xC − xB))/denominator:
 8  gC:=((yB − yC)*(yA − yB) + (xB − xA)*(xC − xB))/denominator:
 9  loadv[nodesinfo[j,2]]:=loadv[nodesinfo[j,2]] − rho/6*(xA*(yB − yC) + xB*(yC − yA) + xC*(yA − yB)):
10  if nodesinfo[indexA,1] = 1 then
11   loadv[nodesinfo[j,2]]:=loadv[nodesinfo[j,2]] − nodesinfo[indexA,2]*gA:
12  else
13   stiffm[nodesinfo[j,2],nodesinfo[indexA,2]]:=stiffm[nodesinfo[j,2],nodesinfo[indexA,2]] + gA:
14  end if:
15  if nodesinfo[indexB,1] = 1 then
16   loadv[nodesinfo[j,2]]:=loadv[nodesinfo[j,2]] − nodesinfo[indexB,2]*gB:
17  else
18   stiffm[nodesinfo[j,2],nodesinfo[indexB,2]]:=stiffm[nodesinfo[j,2],nodesinfo[indexB,2]] + gB:
19  end if:
20  if nodesinfo[indexC,1] = 1 then
21   loadv[nodesinfo[j,2]]:=loadv[nodesinfo[j,2]]-nodesinfo[indexC,2]*gC:
22  else
23   stiffm[nodesinfo[j,2],nodesinfo[indexC,2]]:=stiffm[nodesinfo[j,2],nodesinfo[indexC,2]] + gC:
24  end if:
25 end proc:

We require the x- and y-coordinates of the three points of the subtriangle. These are obtained from the respective rows in nodesinfo. Column 3 holds the x-coordinate and column 4 the y-coordinate (see line 2, line 3, and line 4). We then calculate the denominator denominator and the values gA, gB, and gC according to Eq. 32.59. Using these values we update the load vector also according to Eq. 32.59 in line 9.

We then proceed by updating the stiffness matrix stiffm for each of the three points (see line 10, line 15, and line 20). However, we have to account for two cases here. If the respective node point is a boundary value (in which case the first column entry in nodesinfo is 1) we do not update the respective entry in the stiffness matrix and the value can simply be evaluated to a constant. On the other hand, if the value of the respective node is not known (in which case the first column entry in nodesinfo is 0), we update the respective value in the stiffness matrix stiffm and the load vector loadvvia the reference stored in the second column. Please recall that due to the fact that there are more nodes in the mesh than entries in the stiffness matrix and the load vector, we have to use a reference rather than accessing row i in stiffm and loadv, which may not hold the entry for node i.

Having implemented Eq. 32.59 in the function FEM2DTriangle the FEM scheme is almost complete.

Processing the Mesh. The next function we require simply scans over the mesh and processes each node sequentially. This function is ProcessMesh shown in listing 32.6. It accepts the matrix nodesinfo and the function rho, which is assumed to be a function of the node index. This function could then be used to calculate the right-hand side of Eq. 32.59 depending on the position of the current node. In addition the function requires the stiffness matrix stiffm and the load vector loadv. It starts by determining the total number of nodes in the mesh (see line 2). It then processes the nodes sequentially checking if the respective node is a boundary value (see line 4). If the value is a boundary value, it is skipped. For all nonboundary values the value of rho currentrho is updated by calling the function rho with the current node index as an argument. For each nonboundary value node we then call the function DoTriangleContribution for each of the subtriangles stored in the sixth column of nodesinfo. Please note that the stiffness matrix stiffm and the load vector loadv must be passed by reference in order to allow DoTriangleContribution to update the matrices. As you can see, the function ProcessMesh simply calls DoTriangleContribution for each subtriangle of node i.

Listing 32.6

[Maple] Mesh processing for two-dimensional FEM. This listing is part of the listing FEM2D.mw, which can be found in the supplementary material.

 1 ProcessMesh:=proc(nodesinfo,rho,stiffm,loadv) local nodescount,currentrho,i,j:
 2  nodescount:=Size(nodesinfo,1):
 3  for i to nodescount do
 4   if nodesinfo[i,1] = 0 then
 5    currentrho:=rho(i):
 6    for j to nodesinfo[i,5] do:
 7     DoTriangleContribution(i,nodesinfo[i,6][j,1],nodesinfo[i,6][j,2],nodesinfo[i,6][j,3], nodesinfo,currentrho,'stiffm','loadv'):
 8    end do:
 9   end if:
10  end do:
11  return:
12 end proc:

Solving the Mesh. The last function we require is the function SolveMesh shown in listing 32.7. This function requires the matrix nodesinfo as well as the updated stiffness matrix stiffm and the updated load vector loadv. It then calls the functions DoLUSubstitutionDensePivot as well as the second required function DoLUDecompositionSequentialDensePivot from CoreFunctions.txt, which solve the linear system of equations using sequential dense LU-decomposition (see section 3.6). This is done in line 2. The solution is then simply pasted into a matrix that contains the x-coordinates (first column), the y-coordinates (second column), and the value of the dependent functions (third column, see line 6 and line 7). Here again, we have to determine whether the respective mesh value is a boundary value or not (see line 8). In the case of a boundary value (which is indicated by a 1 in the first column of nodesinfo) the dependent function value is simply copied from the second column of nodesinfo. If the value is not a boundary value (which is indicated by a 0 in the first column of nodesinfo) the value of the dependent function is contained in the solution vector sol, which can be accessed by using the index stored in the second column of nodesinfo.

Listing 32.7

[Maple] Solving the mesh in two-dimensional FEM. This listing is part of the listing FEM2D.mw, which can be found in the supplementary material.

 1 SolveMesh:=proc(nodesinfo,stiffm,loadv) local nodescount,sol,XYZ,i:
 2  sol:=DoLUSubstitutionDensePivot(DoLUDecompositionSequentialDensePivot(stiffm),loadv):
 3  nodescount:=Size(nodesinfo,1):
 4  XYZ:=Matrix(nodescount,3):
 5  for i from 1 to nodescount do:
 6   XYZ[i,1]:=nodesinfo[i,3]:
 7   XYZ[i,2]:=nodesinfo[i,4]:
 8   if nodesinfo[i,1] = 1 then
 9     XYZ[i,3]:=nodesinfo[i,2]:
10   else
11     XYZ[i,3]:=sol[nodesinfo[i,2]]:
12   end if:
13  end do:
14  return XYZ:
15 end proc:

32.6.5 Using the Worksheet to Solve the Poisson Equation on a Circular Cross-Sections

Returning to our example given by Eq. 32.50 we now use the implemented functions to solve the Poisson equation on the circular cross-section. As you can see, listing 32.8 is rather short. We call the functions Initialize, ProcessMesh, and SolveMesh and simply plot the result which is shown in Fig. 32.8. The values used for the calculation were η = 1 mPas and dpdx=−0.001 si4_e mbarmm−1. The mesh used is shown in Fig. 32.6 with a radius of 0.5 mm. The result is surprisingly close to the analytical solution that we have already derived in Eq. 16.19, although we have used a rather coarse mesh. The analytical solution expects a center velocity of 6.25 mms−1.

Listing 32.8

[Maple] Applying two-dimensional FEM using the implemented functions. This listing is part of the listing FEM2D.mw, which can be found in the supplementary material.

1 Initialize(nodes,neighbors,boundaryvalues,'stiffm','loadv','nodesinfo'):
2 ProcessMesh(nodesinfo,i − > − 100,'stiffm','loadv'):
3 res:=SolveMesh(nodesinfo,stiffm,loadv):
4 surfdata(res,labels = ["y [mm]", "z [mm]",typeset(v[x]^empty,[mm/s])]);
f32-08-9782294740381
Fig. 32.8 Two-dimensional FEM applied to solving the Poisson equation on a circular cross-section using the implemented functions. The values used were η = 1 mPas, dpdx=−0.001 si4_e mbar mm−1 and R = 0.5 mm. The obtained numerical result is very close to the analytical solution given by Eq. 16.19.

Using our implemented two-dimensional FEM we find a center velocity of 6.47 mms−1. Even more exact solutions can be obtained by choosing a finer mesh, which becomes somewhat tedious to do manually. This is why FEM solvers always implement algorithms for mesh generation. In many applications, a solver may first approximate the solution on a coarse mesh and then refine the mesh for patches in which large changes are observed. By this iterative meshes can be obtained that gradually refine the solution where necessary.

32.7 Summary

In this section, we have studied FEM, which is one of the most commonly used techniques in CFD. As we have seen, the advantage of FEM is the freedom of choice for the mesh and its generation. The solutions we developed for one- and two-dimensional applications were very close to the analytical solutions. In particular, we used a circular cross-section in our example for two-dimensional FEM, highlighting how conveniently this method can be applied to cross-sections with curvatures. The mathematics underlying FEM are surprisingly simple, and, as we have seen, it can be conveniently implemented in any programming language. Please note that although we chose to implement our functions in Maple the listings do not use mechanisms that are specific to Maple. In fact, the listings can be directly translated to all programming languages, allowing for fast and convenient means of solving differential equations.

..................Content has been hidden....................

You can't read the all page of ebook, please click here login for view all page.
Reset