Finite Element Programming

This article is devoted to the own implementation (Solker Joker FEM ) of the finite element method for systems of diffusion-reaction equations.

It is usually preferable to use ready-made solutions, but if there are specific features in the task, then the task is easier to solve using a simple library.

Introduction

The finite element method [1-4] is a group of effective methods for solving problems of mathematical physics, continuum mechanics. FEM programming issues are devoted to sources [5–9].

Formulation of the problem

Robin’s boundary problem for the system $N$ diffusion-reaction equations in a bounded domain $\ Omega \ subset \ mathbb R ^ 2$ with border $\ Gamma$ :

- a_{i} D e l t a u_{i} (x) + s u m_{k = 1}^{N} q_{i k} (x) u_{k} (x) = g_{i} (x), x i n O m e g a,

$-a_i \ Delta u_i (x) + \ sum_ {k = 1} ^ N q_ {ik} (x) u_k (x) = g_i (x), \; \; x \ in \ Omega,$

a_{i} f r a c p a r t i a l u_{i} (x) p a r t i a l m a t h b f n + b_{i} (x) u_{i} (x) = b_{i} (x) w_{i} (x), x i n G a m m a .

$a_i \ frac {\ partial u_i (x)} {\ partial \ mathbf n} + b_i (x) u_i (x) = b_i (x) w_i (x), \; \; x \ in \ Gamma.$

Here $i = 1, 2, \ ldots, N$ , $\ mathbf n$ - vector of external normal to $\ Gamma$ .

Weak language

Multiply each of $N$ equations for an arbitrary test function $v$ , we integrate the obtained identities by $\ Omega$ , apply the Green formula and take into account the boundary conditions. Will get

a_{i} i n t l i m i t s_{} O m e g a n a b l a u_{i} n a b l a v d x + i n t l i m i t s_{} G a m m a b_{i} u_{i} v d G a m m a + s u m_{k = 1}^{N} i n t l i m i t s_{} O m e g a q_{i k} u_{k} v d x = i n t l i m i t s_{} O m e g a g_{i} v d x + i n t l i m i t s_{} G a m m a b_{i} w_{i} v d G a m m a, i = 1, l d o t s, N .

$a_i \ int \ limits_ \ Omega \ nabla u_i \ nabla v \, dx + \ int \ limits_ \ Gamma b_i u_i v \, d \ Gamma + \ sum_ {k = 1} ^ N \ int \ limits_ \ Omega q_ { ik} u_k v \, dx = \ int \ limits_ \ Omega g_i v \, dx + \ int \ limits_ \ Gamma b_i w_i v \, d \ Gamma, \; \; i = 1, \ ldots, N.$

We introduce the Sobolev functional space $V = H ^ 1 (\ Omega)$ . Functions $u_1, \ ldots, u_N \ in V$ are called a weak solution of the original boundary value problem if the indicated identities are satisfied for an arbitrary test function $v \ in V$ .

Weak wording can be written as follows: find functions $u_1, \ ldots, u_N \ in V$ such that

A_{i} (u_{i}, v) + s u m_{k = 1}^{N} B_{i k} (u_{k}, v) = F_{i} (v) f o r a l l v i n V, i = 1, l d o t s, N,

$A_i (u_i, v) + \ sum_ {k = 1} ^ N B_ {ik} (u_k, v) = F_i (v) \; \; \ forall v \ in V, \; i = 1, \ ldots, N,$

Where $A_i, B_ {ik}$ - bilinear forms, $F_i$ - linear form:

A_{i} (u, v) = a_{i} i n t l i m i t s_{} O m e g a n a b l a u n a b l a v d x + i n t l i m i t s_{} G a m m a b_{i} u v d G a m m a, q u a d B_{i k} (u, v) = s u m_{k = 1}^{N} i n t l i m i t s_{} O m e g a q_{i k} u v d x,

$A_i (u, v) = a_i \ int \ limits_ \ Omega \ nabla u \ nabla v \, dx + \ int \ limits_ \ Gamma b_i uv \, d \ Gamma, \ quad B_ {ik} (u, v) = \ sum_ {k = 1} ^ N \ int \ limits_ \ Omega q_ {ik} uv \, dx,$

F_{i} (v) = i n t l i m i t s_{} O m e g a g_{i} v d x + i n t l i m i t s_{} G a m m a b_{i} w_{i} v d G a m m a .

$F_i (v) = \ int \ limits_ \ Omega g_i v \, dx + \ int \ limits_ \ Gamma b_i w_i v \, d \ Gamma.$

FEM algorithm

We introduce in space $V$ finite-dimensional subspace $V_m$ . According to the Galerkin method, we arrive at a discrete problem: find functions $u_1, \ ldots, u_N \ in V_m$ such that

A_{i} (u_{i}, v) + s u m_{k = 1}^{N} B_{i k} (u_{k}, v) = F_{i} (v) f o r a l l v i n V_{m}, i = 1, l d o t s, N .

$A_i (u_i, v) + \ sum_ {k = 1} ^ N B_ {ik} (u_k, v) = F_i (v) \; \; \ forall v \ in V_m, \; i = 1, \ ldots, N.$

Let be $\ phi_1, \ ldots, \ phi_m$ - basis $V_m$ . If we put $u_i (x) = \ displaystyle \ sum_ {s = 1} ^ m c_ {is} \ phi_s (x)$ , $v (x) = \ phi_j (x)$ , then the task is reduced to SLAE with respect to $Nm$ unknown $c_ {ks}$ :

s u m_{s = 1}^{m} A_{i} (p h i_{s}, p h i_{j}) c_{i s} + s u m_{k = 1}^{N} s u m_{s = 1}^{m} B_{i k} (p h i_{s}, p h i_{j}) c_{k s} = F_{i} (p h i_{j}), i = 1, l d o t s, N, j = 1, l d o t s, m,

$\ sum_ {s = 1} ^ m A_i (\ phi_s, \ phi_j) c_ {is} + \ sum_ {k = 1} ^ N \ sum_ {s = 1} ^ m B_ {ik} (\ phi_s, \ phi_j) c_ {ks} = F_i (\ phi_j), \; \; i = 1, \ ldots, N, \; j = 1, \ ldots, m,$

s u m_{s = 1}^{m} l e f t (a_{i} i n t l i m i t s_{O m e g a} n a b l a p h i_{s} c d o t n a b l a p h i_{j} d x + i n t l i m i t s_{} G a m m a b_{i} p h i_{s} p h i_{j} r i g h t) c_{i s} + s u m_{k = 1}^{N} s u m_{s = 1}^{m} l e f t (i n t_{} O m e g a q_{i k} p h i_{s} p h i_{j} d x r i g h t) c_{k s} =

$\ sum_ {s = 1} ^ m \ left (a_i \ int \ limits _ {\ Omega} \ nabla \ phi_s \ cdot \ nabla \ phi_j \ dx + \ int \ limits_ \ Gamma b_i \ phi_s \ phi_j \ right) c_ {is} + \ sum_ {k = 1} ^ N \ sum_ {s = 1} ^ m \ left (\ int_ \ Omega q_ {ik} \ phi_s \ phi_j \, dx \ right) c_ {ks} =$

= i n t l i m i t s_{} O m e g a g_{i} p h i_{j} d x + i n t l i m i t s_{} G a m m a b_{i} w_{i} p h i_{j} d G a m m a, i = 1, l d o t s, N, j = 1, l d o t s, m .

$= \ int \ limits_ \ Omega g_i \ phi_j \, dx + \ int \ limits_ \ Gamma b_i w_i \ phi_j \, d \ Gamma, \; \; i = 1, \ ldots, N, \; j = 1, \ ldots, m.$

The simplest example of a subspace $V_m$ - the space of continuous piecewise linear functions. Region $\ Omega$ is divided into finite elements (usually triangles), and on each triangle the functions of $V_m$ are linear. In the FEM, the basic functions are chosen so that in most parts of the domain they vanish, then there will be many zeros in the SLAE matrix. As basic functions $\ phi_j$ take "pyramids", which are equal to 1 in one grid node and 0 in the other nodes. Each node of the grid has its own basic function, and the dimension $m$ subspace is equal to the number of nodes in the grid.

Note that if $s$ th and $j$ nodes are not connected by an edge, then $\ varphi_s \ varphi_j \ equiv 0$ , $\ nabla \ varphi_s \ cdot \ nabla \ varphi_j \ equiv 0$ so when summing over $s$ it is enough to go through only the indices $s$ for which between $s$ m and $j$ nodes are edge either $s = j$ .

Thus, the solution of the boundary value problem by the finite element method is reduced to the construction and solution of the corresponding SLAE.

Piecewise linear functions

To set the function $u \ in V_m$ , you need to specify its values in the grid nodes. We denote the coordinates $i$ -th node through $(X_i, Y_i)$ and the value of the function $u$ at $i$ node - through $U_i = u (X_i, Y_i)$ .

To find the value of a function $u$ at an arbitrary point $(x, y)$ belonging to the triangle with $i$ th $j$ th and $k$ vertices, we introduce the barycentric coordinates $L_1, L_2, L_3$ that take values from 0 to 1 and are related to Cartesian coordinates $x, y$ the following relations (see [3, Section 4.7.1], [4, p. 35]):

b e g i n g a t h e r e d x = L_{1} X_{i} + L_{2} X_{j} + L_{3} X_{k}, y = L_{1} Y_{i} + L_{2} Y_{j} + L_{3} Y_{k}, 1 = L_{1} + L_{2} + L_{3} . e n d g a t h e r e d

$\ begin {gathered} x = L_1 X_i + L_2 X_j + L_3 X_k, \\ y = L_1 Y_i + L_2 Y_j + L_3 Y_k, \\ 1 = L_1 + L_2 + L_3. \ end {gathered}$

From here

b e g i n g a t h e r e d L_{1} (x, y) = f r a c 1 2 A (a_{i} + b_{i} x + c_{i} y), a_{i} = X_{j} Y_{k} - X_{k} Y_{j}, b_{i} = Y_{j} - Y_{k}, c_{i} = X_{k} - X_{j}, L_{2} (x, y) = f r a c 1 2 A (a_{j} + b_{j} x + c_{j} y), a_{j} = X_{k} Y_{i} - X_{i} Y_{k}, b_{j} = Y_{k} - Y_{i}, c_{j} = X_{i} - X_{k}, L_{3} (x, y) = f r a c 1 2 A (a_{k} + b_{k} x + c_{k} y), a_{k} = X_{i} Y_{j} - X_{j} Y_{i}, b_{k} = Y_{i} - Y_{j}, c_{k} = X_{j} - X_{i}, e n d g a t h e r e d

$\ begin {gathered} L_1 (x, y) = \ frac {1} {2A} (a_i + b_i x + c_i y), \; a_i = X_j Y_k - X_k Y_j, \; b_i = Y_j - Y_k, \; c_i = X_k - X_j, \\ L_2 (x, y) = \ frac {1} {2A} (a_j + b_j x + c_j y), \; a_j = X_k Y_i - X_i Y_k, \; b_j = Y_k - Y_i, \; c_j = X_i - X_k, \\ L_3 (x, y) = \ frac {1} {2A} (a_k + b_k x + c_k y), \; a_k = X_i Y_j - X_j Y_i, \; b_k = Y_i - Y_j, \; c_k = X_j - X_i, \\ \ end {gathered}$

Where $A = \ dfrac {1} {2} \ begin {vmatrix} 1 & X_i & Y_i \\ 1 & X_j & Y_j \\ 1 & X_k & Y_k \ end {vmatrix}$ - area of a triangle with a sign.

The value of the function can be found by the formula $u (x, y) = U_i L_1 (x, y) + U_j L_2 (x, y) + U_k L_3 (x, y)$ .

Numerical integration

Let us reduce the integral over the boundary edge with $i$ m and $j$ ends of the integral over a segment $[- 1, 1]$ :

i n t l i m i t s_{i}^{j} f (x, y) d l = f r a c | i j | 2 i n t l i m i t s_{- 1}^{1} f (x (t), y (t)) d t .

$\ int \ limits_i ^ jf (x, y) \, dl = \ frac {| i \, j |} {2} \ int \ limits _ {- 1} ^ 1 f (x (t), y (t) ) \, dt.$

Here $| i \, j |$ - edge length, $x (t) = \ dfrac {X_i + X_j} {2} + \ dfrac {X_j - X_i} {2} t$ ,

$y (t) = \ dfrac {Y_i + Y_j} {2} + \ dfrac {Y_j - Y_i} {2} t$ .

To calculate the last integral, apply the Gauss quadrature formula:

i n t l i m i t s_{- 1}^{1} v a r p h i (s) d s a p p r o x s u m_{j = 1}^{n} w_{j} v a r p h i (x i_{j}) .

$\ int \ limits _ {- 1} ^ 1 \ varphi (s) \, ds \ approx \ sum_ {j = 1} ^ n w_j \ varphi (\ xi_j).$

Points $\ xi_j \ in [-1, 1]$ and weight factors $w_j$ taken from the table [3, sec. 5.9.2].

We will reduce the integral of the triangle with $i$ th $j$ th and $k$ vertices of a triangle integral $D = \ {(L_1, L_2) \ in [0,1] \ times [0,1] \ colon L_1 + L_2 \ leq 1 \}$ :

i n t l i m i t s_{i j k} f (x, y) d x d y = 2 | i j k | i n t l i m i t s_{D} f (x (L_{1}, L_{2}), y (L_{1}, L_{2})) d L_{1} d L_{2} .

$\ int \ limits_ {i \, j \, k} f (x, y) \, dxdy = 2 | i \, j \, k | \ int \ limits_D f (x (L_1, L_2), y (L_1, L_2)) \, dL_1dL_2.$

Here $| i \, j \, k | = | A |$ - area of a triangle, $x (L_1, L_2) = L_1X_i + L_2X_j + (1 - L_1 - L_2) X_k$ ,

$y (L_1, L_2) = L_1Y_i + L_2Y_j + (1 - L_1 - L_2) Y_k$ .

To calculate the last integral, we apply the quadrature formula [3, Sec. 5.11].

Software implementation

General considerations

Testing is difficult in this subject domain [10], so it is important that the code be written very clearly. OOP allows you to increase the readability of the code; to use OOP, you need to select a set of concepts and formulate an algorithm in these terms. This makes it possible in the first place to pay attention to the algorithm, and only then to the implementation details.

Triangulation

To build the grid, the FreeFem ++ package is used, the triangulation is saved in the .mesh format.

The library introduces the Mesh class, which contains mesh information that is read from the file.

Build SLAU

To write the algorithm, classes of basic functions and integrands of functions (products of functions) are introduced, which are inherited from the base classes Integrand and BoundaryIntegrand . In these classes there are support (carrier) and boundary_support (boundary carrier) fields, respectively. The carrier is a list of triangles on which the function is not equal to 0. The boundary carrier is a list of boundary edges on which the function is not equal to 0. Also, the Value method is defined for classes of integrand functions, which calculates the value of the function at a given point.

The Integrate and BoundaryIntegrate functions take as input an integrand function, iterate through the carrier list and call the integration function along a given triangle (or boundary edge), which is defined in the Integrand base class ( BoundaryIntegrand ) and calls the virtual function Value .

Thus, the SLAU build code has a simple form:

 for (int i = 0; i < data.N; ++i) { for (int j = 0; j < data.mesh.nodes_num; ++j) { for (auto s : data.mesh.adjacent_nodes[j]) { sys.AddCoeff(i, j, i, s, data.a[i] * Integrate(grad(phi[s]) * grad(phi[j])) + BoundaryIntegrate(data.b[i] * phi[s] * phi[j]) ); } for (int k = 0; k < data.N; ++k) { for (auto s : data.mesh.adjacent_nodes[j]) { sys.AddCoeff(i, j, k, s, Integrate(data.q[i][k] * phi[s] * phi[j]) ); } } sys.AddRhs(i, j, Integrate(data.g[i] * phi[j]) + BoundaryIntegrate(data.b[i] * data.w[i] * phi[j]) ); } }

Decision SLAU

To solve the SLAE, an iterative method is used, implemented in the MTL4 library.

Access to the result

To find the solution value at a given point, it is necessary in a short time to find the triangle containing the given point. For this you can use a hash table.

Divide the region enclosing rectangle. $\ Omega$ on $nm$ blocks where $n = [W / a]$ , $m = [H / b]$ , $W, H$ - the length and width of the enclosing rectangle, $a, b$ - the average size of the enclosing triangles of triangles triangles.

For each of $nm$ blocks fill the list of triangles that intersect it. For this, let's iterate through all the triangles and find their enclosing rectangles,

for which it is easy to find the blocks crossing them.

When executing a query for a given point, we find the block to which it belongs, and we look at the list of triangles corresponding to this block, checking which of them the given point belongs to.

Testing

Testing the program was carried out on tasks with a known exact solution. With double grinding of the grid, the error decreases by about 4 times or 2 times, which indicates that the method has the 2nd or 1st order of accuracy. The reduction in the order of convergence to the 1st may be due to the inconsistency of the boundary conditions at the joints of the parts of the boundary on which different values of the function are given. $w_i$ .

Conclusion

The advantage of the developed library over large solvers is its simplicity, as well as the ability to take into account the specific features of the problem (piecewise constant coefficient in the boundary condition). In the future, we plan to implement a similar library for 3-dimensional tasks, which will be used to solve optimal control problems for the boundary coefficient in the complex heat exchange model, where this coefficient is piecewise constant [11].

Sources

Alekseev G.V. Numerical methods of mathematical physics. Introduction to the finite element method. Vladivostok: Dal'nauka, 1999.
Dautov R.Z., Karchevsky M.M. Introduction to the theory of the finite element method. Kazan: KSU, 2004.
Zienkiewicz OC, Taylor RL, Zhu JZ The finite element method. Elsevier, 2005.
Segerlind L. Application of the finite element method. M .: Mir, 1979.
The finite element method on the example of the Poisson equation // Habrahabr, 07.27.2015.
Writing the FEM of the calculator in less than 180 lines of code //

Habrahabr, 12/01/2015.
Practice of using Freefem ++ // Habrahabr, 06.06.2013.
Gockenbach MS Understanding and finite element method. SIAM, 2006.
Smith IM, Griffiths DV, Margetts L. Programming the finite element method. Wiley, 2014.
Why unit tests do not work in scientific applications // Habrahabr, 04/26/2010.
Grenkin G.V. Algorithm for Solving the Boundary Optimal Control Problem in the Model of Complex Heat Exchange, Dal'nevost. Math journals 2016. V. 16. No. 1. P. 24-38.

Source: https://habr.com/ru/post/344564/

All Articles