Simulation und wissenschaftliches Rechnen II (SiwiR II) Christoph - PowerPoint PPT Presentation

Example 1: Poisson’s Equation in 1D − u ′′ = on ]0 , 1[ , f u (0) = 0 , u (1) = 0 . Discretize this equation by A h x h = ˜ f h where �� 1 � �� 1 � ˜ v ′ q v ′ A h = f h = I h ( f ) v p dx p dx , . p,q ∈ ◦ p ∈ ◦ 0 0 Ω h Ω h This means     2 − 1 4 1     f ( h ) x 1     − 1 2 − 1 1 4 1         . .         1 = h  ... ... ...     ... ... ...    .    .     .      6     h         . .         − 1 2 − 1 1 4 1             x m − 1 f (1 − h ) − 1 2 1 4 – p. 27/123

Finite Elements - Central Difference The discretization of u �→ u ′ by finite elements leads to a discretization similar to the central difference discretization   0 1     − 1 0 1     1 ... ... ...   = A h .   2     − 1 0 1   − 1 0 How do we get something similar to FD upwind? – p. 28/123

Streamline Diffusion Discretization Consider the convection diffusion equation in 1D: − u ′′ − bu ′ = f, u (0) = u (1) = 0 ◦ Multiply this equation by v = v h − ρhv ′ h sgn b , where v h ∈ V h and integrate. Assuming b > 0 , this yields � 1 � 1 � 1 � � u ′ v ′ h + hρbu ′ v ′ h − bu ′ v h u ′′ v ′ f ( v h − ρhv ′ dx + ρh h dx = h ) dx. 0 0 0 In the streamline diffusion discretization, we neglect the term of third order and replace u by u h : � 1 � 1 � � (1 + hρb ) u ′ h v ′ h − bu ′ f ( v h − ρhv ′ dx = h ) dx. h v h 0 0 – p. 29/123

Streamline Diffusion Discretization Let ρ = 1 2 . Then, the stencil corresponding to the term � 1 � 1 � � � 1 � hρbu ′ h v ′ h − bu ′ 2 hu ′ h v ′ h − u ′ dx = b h v h h v h dx 0 0 is   1 − 1     0 1 − 1     ... ... ...   b .       0 1 − 1   0 1 This shows that the finite element streamline diffusion discretization is similar to the FD upwind discretization. – p. 30/123

Finite Elements in 2D/3D Definition 4. T = { T 1 , . . . , T M } is a conform triangulation of Ω if Ω = � M i =1 T i , T i is a triangle or quadrangle (in 2D) or tetrahedron, hexahedron, prism, or pyramid (in 3D) T i ∩ T j is either empty or one common corner or one common edge. – p. 31/123

Finite Elements in 2D/3D Remark. Let us write T h , if the diameter h T of every element T ∈ T h is less or equal h : h T ≤ h. A family of triangulations {T h } is called quasi-uniform, if there exists a constant ρ > 0 such that the radius ρ T of the largest inner ball of every triangle T ∈ T h satisfies ρ T > ρh. – p. 32/123

Good and Bad Triangles bad approximation good bad for Gauss − Seidel – p. 33/123

Linear Elements in 2D Definition 5. Let T h be a triangulation of Ω . Then, let V h be the space of linear finite elements defined as follows: � � � � � � � v ∈ C 0 (Ω) T is linear for every T ∈ T h V h = � v ◦ V h ∩ H 1 = 0 (Ω) V h � � � � T ( x, y ) = a + bx + cy . v T is linear means that v – p. 34/123

Bilinear Elements in 2D Definition 6 (Bilinear elements on a Cartesian 2D grid) . Let Ω =]0 , 1[ 2 , h = 1 m and � � � � � T h = [ ih, ( i + 1) h ] × [ jh, ( j + 1) h ] � i, j = 0 , . . . , m − 1 . The space of bilinear finite elements on Ω is defined as follows � � � � � v ∈ C 0 (Ω) � � T is bilinear for every T ∈ T H V h = � v . � � � � v T is bilinear means that v T ( x, y ) = a + bx + cy + dxy . – p. 35/123

FE Discretization of Poisson’s equation − ∆ u = f � � = 0 . u δ Ω Multiplication with v h and integration leads to: ◦ FE Discretization: Find u h ∈ V h such that � � ◦ ∇ u h ∇ v h d ( x, y ) = ∀ v h ∈ f v h d ( x, y ) V h . (12) Ω Ω – p. 36/123

Nodal Basis Functions Definition 7. Let V h be the space of linear or bilinear finite elements on T h and N h the set of corners of T h . Then, define the nodal basis function v q ∈ V h at the point q by: � 1 if x = q for x ∈ N h v q ( x ) = if x � = q 0 Observe that � � � � � V h = span � q ∈ N h v q This means that every function u h ∈ V h can be represented as � u h = λ q v q q ∈N h – p. 37/123

Stiffness matrix � � ∇ v q ∇ v p d ( x, y ) , a p.q := f p := f v p d ( x, y ) Ω Ω 0 N h := N h ∩ Ω A h := ( a p,q ) , 0 p,q ∈ N h � = u h λ q v q . 0 q ∈ N h Then, (12) implies U h = ( λ q ) 0 q ∈ N h = where A h U h F h F h = ( f p ) 0 p ∈ N h The matrix A h is called the stiffness matrix of the FE discretization. – p. 38/123

Bilinear Elements on a Structured Grid Consider the structured grid on Ω =]0 , 1[ 2 : � � � � � T h = [ ih, ( i + 1) h ] × [ jh, ( j + 1) h ] � i, j = 0 , . . . , m − 1 . N h is the set of corresponding nodal points (corner points). Observe that the nodal basis functions can be decomposed as v p x p y ( x, y ) = v p x ( x ) · v p y ( y ) . Thus, � � 1 � 1 ∂v p x ∂v q x ∇ v q x q y ∇ v p x p y d ( x, y ) = ∂x d x v p y v q y d y ∂x Ω 0 0 � 1 � 1 ∂v p y ∂v q y + d y v p x v q x d x. ∂y ∂y 0 0 – p. 39/123

Bilinear Elements on a Structured Grid This shows that the discretization stencil for Poisson’s equation is:     − 1 1 � � � �     1 1 · h  · h     + − 1 − 1 2 4 2 1 4 1    6 6 h h − 1 1   − 1 − 1 − 1   1   = − 1 − 1 8   3 − 1 − 1 − 1 and for the right hand side the stencil is:     1 1 4 1 � �   h 2   h · h     =  . 1 4 1 4 4 16 4    6 6 36 1 1 4 1 – p. 40/123

Local Stiffness Matrix Since Ω = � M i =1 T i , we obtain � � M � ∇ v q ∇ v p d ( x, y ) = ∇ v q ∇ v p d ( x, y ) . Ω T i i =1 For linear or bilinear elements, we obtain � ∇ v q ∇ v p d ( x, y ) � = 0 ⇔ p, q ∈ T i . T i Definition 8. The matrix �� � ∇ v q ∇ v p d ( x, y ) T i p,q ∈ T i is called local stiffness matrix at T i . – p. 41/123

Reference Element To calculate the local stiffness matrices we need a reference element ˆ T and a mapping Ψ i : ˆ T → T i for every i . Example 1. A reference element for triangles is: ˆ T = { ( ξ, η ) | ξ + η ≤ 1 ξ, η ≥ 0 } . and If T i consists of the corners P 1 , P 2 , P 3 , then Ψ i ( ξ, η ) = P 1 + ( P 2 − P 1 ) ξ + ( P 3 − P 1 ) η. Example 2. A reference element for quadrangles is: ˆ T = { ( ξ, η ) | 0 ≤ ξ, η ≤ 1 } . – p. 42/123

Calculation of Local Stiffness Matrices Now, the local stiffness element can be calculated by � ∇ v T q ∇ v p d ( x, y ) = T i � � � T � � ( D Ψ i ) − T ∇ ˆ ( D Ψ i ) − T ∇ ˆ | det D Ψ i | d ( ξ, η ) . = v q v p ˆ T Example 3. Consider triangles. Then, describe the mapping Ψ i by      a  c  ξ +  η. Ψ i ( ξ, η ) = P 1 + b d Then,    a c  . D Ψ i = b d – p. 43/123

Numerical Integration Calculate the integral � � � T � � ( D Ψ i ) − T ∇ ˆ ( D Ψ i ) − T ∇ ˆ | det D Ψ i | d ( ξ, η ) . v q v p ˆ T by Gauss quadrature rule. Example 4. Consider triangles and choose the above reference element. Then, the first Gauss quadrature rule implies: � � � T � � ( D Ψ i ) − T ∇ ˆ ( D Ψ i ) − T ∇ ˆ | det D Ψ i | d ( ξ, η ) v q v p ˆ T � 1 � � � T � � 1 3 , 1 ( D Ψ i ) − T ∇ ˆ ( D Ψ i ) − T ∇ ˆ ≈ | det D Ψ i | v q v p . 2 3 – p. 44/123

Calculation of Stiffness Matrix P 2 The stiffness matrix of this triangulation is: A P 3 A h = B P 1 C   l A 11 + l B l A l A 13 + l B l B 0 11 12 13 14 P 4 P 5     l A l A l A 0 0   21 22 23     l A 31 + l B l A l A 33 + l B l B 34 + l C l C   31 32 33 34 35   { 1 , 2 , 3 } N ( A ) =   l B l B 43 + l C l B 44 + l C l C 0   41 43 44 45 { 1 , 3 , 4 } = N ( B ) l C l C l C 0 0 53 54 55 { 4 , 3 , 5 } N ( C ) = – p. 45/123

Algorithmic Calculation of Stiffness Matrix The calculation of stiffness matrix has to be performed in two steps: 1. Step: Calculate the local stiffness matrix. 2. Step: Calculate the stiffness matrix. But there exist two approaches: 1. First compute and store the whole local stiffness matrix. Then, calculate the stiffness matrix. Advantage: The local stiffness matrices can be used for coarsening the local stiffness matrices in a multigrid algorithm. Faster code for some non-linear problems. 2. After the calculation of the local stiffness matrix of one element T , add these integrals to the whole stiffness matrix. Advantage: Less storage requirement. – p. 46/123

Data Structure for (Local) Stiffness Matrix 1. Local stiffness matrix: Let n T be the number of degrees of freedom for an element T ∈ T h (3 for triangle). Then, for every T ∈ T h a n T × n T matrix has to be stored. Data structure: list or array for storing T ∈ T h . Each entry must contain n T and a pointer to a n T × n T matrix. 2. Stiffness matrix: For every unknown (grid point) the discretization stencil has no fixed size. Data structure: Sparse matrix. – p. 47/123

Sparse Matrix Format 1.row 2.row k.row n k − 1 n 2 + n 1 1 + n k − 2 n 1 + ... + n 1 n 1 n 2 n k i = 1 i = 2 i = k n k − 1 + ... n 2 + ... n 1 1 – p. 48/123

Algorithm (Stiffness Matrix Calculation) Let N ( T ) be the corner points of the triangle T . 1. Calculate local stiffness matrix ( l T ij ) i,j ∈N ( T ) for every finite element T ∈ T h . 2. Calculate the number of neighbour points m i for every point i . This gives the value n i = � i s =1 m s + 1 in the sparse matrix of the the stiffness matrix. 3. Go to every grid point i and iterate over the neighbour element T ∈ T h , i ∈ N ( T ) . Add the l T ij to a ij for every j ∈ N ( T ) . ⇒ Later we explain how to obtain a suitable data structure for the discretization grid. – p. 49/123

Data Structure of the Discretization Grid Array (or list) of objects of type Triangle . Every triangle has an id . class Triangulation_grid { int number_triangles; Triangle* triangles; // id is number in list int number_points; Points* points; // id is number in list } class Triangle { int id_point_1, id_point_2, id_point_3; } class Points { double x,y; // coordinate of point int number_neighbour_points; int* id_neighbour_points; } – p. 50/123

Weak Formulation of a PDE Let V be a vector space and a : V × V → R a symmetric positive definite bilinear form. a induces the “energy” norm � � u � E = a ( u, u ) . Furthermore, let f : V → R be a � · � E continuous linear functional and let V be complete with respect to � · � E . Problem 1. Find u ∈ V such that a ( u, v ) = f ( v ) ∀ v ∈ V. Theorem 1. The above problem has a unique solution u . – p. 51/123

Examples of PDEs with Weak Formulation Example 5 (Poisson’s Equation with Reaction Term) . Let V = H 1 0 (Ω) := { u ∈ L 2 (Ω) | ∇ u ∈ L 2 (Ω) } and c ≥ 0 . � � � ∇ u ∇ v + cuv a ( u, v ) = d ( x, y ) . Ω Example 6 (Linear Elasticity) . Let V = ( H 1 0 (Ω)) 3 , u ∈ V , C a suitable 6 × 6 matrix and Du the vector of symmetric derivatives (see section ?? ). � ( Du ) T C Dv d ( x, y, z ) . a ( u, v ) = Ω Example 7 (Maxwell’s Equations) . Let V be a suitable vector space similar to ( H 1 0 (Ω)) 3 and c ≥ 0 . � ( ∇ × u ) T ( ∇ × v ) + cuv d ( x, y, z ) . a ( u, v ) = Ω – p. 52/123

General Convergence Theory Let V h be a subspace of V . Problem 2. Find u h ∈ V h such that a ( u h , v h ) = f ( v h ) ∀ v h ∈ V h . Theorem 2. � u − u h � E ≤ v h ∈ V h � u − v h � E inf Example 8. Consider Poisson’s equation. Let V h be the space of linear elements corresponding to a a familiy of quasi-uniform triangulations. Furthermore, assume that u ∈ C 2 (Ω) ( H 2 (Ω) ) is the weak solution of Poisson’s equation. Then, there is a constant C such that � u − u h � E ≤ hC for every h , where u h ∈ V h is the finite element solution. Remark: In case of H 2 (Ω) -regularity, one can prove � u − u h � L 2 ≤ h 2 C . – p. 53/123

Operator Formulation Let V h be a finite element space and ( v q ) q ∈N h the corresponding nodal basis. Let u,f be vectors of length |N h | . Then, f = Laplace_FE(u); means   � �   f = ( f p ) p ∈N h = ∇ ( u q v q ) ∇ v p d ( x, y ) Ω q ∈N h p ∈N h and f = Helm_FE(u); means   � �   f = ( f p ) p ∈N h = ( u q v q ) v p d ( x, y ) . Ω q ∈N h p ∈N h – p. 54/123

Operator Formulation Thus, Laplace_FE( ), Helm_FE( ) are operators. Let Diag_Laplace_FE( ), Diag_Helm_FE( ) be the corresponding diagonal operators. Let interior, boundary represent the interior and boundary points of the domain and product(u,v) the scalar product of u and v . – p. 55/123

Test 1: Volume Calculation Now let us implement the above operators by expression templates. Then, the code u = 1.0; f = Helm_FE(u); cout << product(u,f) << endl; calculates the volume of the domain. – p. 56/123

Test 2: Volume Calculation Now let us implement the above operators by expression templates. Then, the code u = X; f = Poisson_FE(u); cout << product(u,f) << endl; calculates the volume of the domain. Here, let X be the x-coordinate. – p. 57/123

Dirichlet Boundary Conditions The problem Problem 3. Find u ∈ H 1 (Ω) such that −△ u = on Ω , f u | Γ = g can approximatively be solved by finite elements and the Gauss-Seidel iteration as follows: u = Helm_FE(f); f = u; u = g | boundary; for(i=0;i<i_max;++i) u = u - (Laplace_FE(u)-f)/Diag_Laplace_FE() | interior_points; – p. 58/123

Using a Direct Solver Let us assume that there is a good direct solver ui=Inverse(S,fi) which calculates ui = S − 1 fi. Then, apply the code u = Helm_FE(f); f = u; u = 0.0 | interior; u = g | boundary; f = f - Laplace_FE(u) | boundary; u = 0.0 | boundary; S = Sparse_matrix(Laplace_FE, interior); fi = vector(f,interior); ui = vector(u,interior); ui = Inverse(S,fi); u = g | boundary; – p. 59/123

Boundary Conditions Let us consider the equation − div µ grad u = on Ω , f = on Γ D , u g ∂u = 0 on Γ N , ∂� n ∂u n + β ( u − u ref ) = 0 on Γ third , ∂� here Ω is a domain and ∂ Ω = Γ D ∪ Γ N ∪ Γ third is a disjunct subdivision of the boundary, where Γ D � = ∅ . Furthermore, let µ : Ω → R be a piecewise constant parameter and 0 < β ∈ R . – p. 60/123

Boundary Conditions Define the finite element space � ¯ � V h := { v h ∈ V h | v h Γ D = 0 } . Then, we obtain � � � � ∇ uµ ∇ v h d ( x, y ) + βµ uv h dσ = βµ u ref v h dσ + fv h d ( x, y ) Ω Γ third Γ third Ω for every v h ∈ ¯ V h . FE Discretization Find u h ∈ V h such that � � � � ∇ u h µ ∇ v h d ( x, y ) + βµ = u ref v h dσ + fv h d ( x, y u h v h dσ βµ Ω Γ third Γ third Ω ∀ v h ∈ ¯ V h , u h ( z ) = g ( z ) ∀ z ∈ Ω h ∩ Γ D . – p. 61/123

Boundary Conditions Observe that the bilinear form � � a ( u, v ) := ∇ uµ ∇ v d ( x, y ) + βµ uv dσ Ω Γ third is symmetric positive definite on the space � � { v ∈ H 1 (Ω) | v Γ D = 0 } . – p. 62/123

Operator Formulation Let us implement the operators in section ?? by expression templates. Let A FE be the operator corresponding to the bilinear form a ( u, v ) . The previous problem can be solved by the Gauss-Seidel iteration as follows: u = Helm_FE(f); f = u; u = g | boundary_D; for(i=0;i<i_max;++i) u = u - (A_FE(u)-f)/Diag_A_FE() | grid_space; Here grid space represents the interior points and the boundary points which are no Dirichlet boundary points. – p. 63/123

Neumann Boundary Conditions Let us consider the equation −△ u = on Ω , f (13) ∂u = 0 on Γ . (14) ∂� n A short calculation shows � f d ( x, y ) = 0 . (15) Ω Thus, we assume (15). – p. 64/123

Neumann Boundary Conditions A natural way to obtain a well-defined problem is: Problem 4. Find u ∈ H 1 (Ω) such that −△ u = f on Ω , � ∂u = 0 on Γ and u d ( x, y ) = 0 , ∂� n Ω where we assume � f d ( x, y ) = 0 . Ω – p. 65/123

Operator Formulation Let us implement the operators in section ?? by expression templates. The previous problem can be solved by Gauss-Seidel iteration as follows: Eins = 1.0; // set up for normalization IntE = Helm_FE(Eins); Eins = Eins / product(Eins,IntE); f = f - Eins * product(f,IntE); u = Helm_FE(f); f = u; for(int i=0; i<N;++i) { u = u - (A_FE(u) -f) / Diag_A_FE(); u = u - Eins * product(u,IntE); } – p. 66/123

Streamline Diffusion Discretization Consider the convection diffusion equation in 1D: −△ u − � b grad u = f � � = 0 . u ∂ Ω b : Ω → R 2 is a vector field. where � ◦ Discretization: Find u h ∈ V h such that � � � ∇ u h ◦ ∇ v h + hρ� b ◦ ∇ u h � b ◦ ∇ v h � � b � − 1 − � b ◦ ∇ u h v h d ( x, y ) 2 Ω � f ( v h − ρh� b ◦ ∇ v h � � b � − 1 = 2 ) d ( x, y ) Ω ◦ for every v h ∈ V h . – p. 67/123

Types of Grids There exist Cartesian grids block structured grids unstructured grids ... to discretize a domain Ω . – p. 68/123

Cartesian Grid Example of an Cartesian grid: { ( ih, jk ) + ( x 0 , y 0 ) | i = 0 , ..., m, j = 0 , ..., n } , Ω h,k = where h, k > 0 . Data structure: Array! – p. 69/123

Block Structured Grid Let Ω 0 = { ( ih, jh ) | i, j = 0 , ..., n } , h where h = 1 n . Furthermore, let T = { T 1 , . . . , T M } be a subdivision by quadrangles (2D) and Ψ k : [0 , 1] 2 → T k smooth bijections such that Ω = � M k =1 Ψ k ([0 , 1] 2 ) . Then, M � Ψ k (Ω 0 Ω h = h ) k =1 is a block structured grid. (Generalizations in 3D and for different mesh sizes are possible.) Data structure: Unstructured grid of quadrangles, each block by array. – p. 70/123

Simple Interpolation A simple construction of the mapping Ψ k : [0 , 1] 2 → T k is P SW + ( P SE − P SW ) η + ( P NW − P SW ) ξ + Ψ k ( η, ξ ) = ( P NE − P SE − P NW + P SW ) ξη. P NW P NE 1 P SW 0 0 1 PSE – p. 71/123

Transfinite Interpolation Let β N , β S , β W , β E : [0 , 1] → R 2 be parameterizations of the north, south, west and east boundary. Then, the transfinite interpolation is: β S ( η ) + ( β N ( η ) − β S ( η )) ξ Ψ k ( η, ξ ) = + β W ( ξ ) + ( β E ( ξ ) − β W ( ξ )) η − P SW − ( P SE − P SW ) η − ( P NW − P SW ) ξ − ( P NE − P SE − P NW + P SW ) ξη. – p. 72/123

Unstructured Grid An example of an unstructured grid is: Data structure for an unstructured grid: list or array of corners (information of coordinates) list or array of triangles, quadrangles, ... with pointers to corners and number of corners. By this information one can construct: list or array of edges and faces (in 3D). – p. 73/123

Visualization of an Unstructured Grid Examples of powerful 3D visualization programs are: AVS (commercial) OpenDx (public domain) ParaView (uses vtk Toolkit, public domain) AVS supports structured grids and unstructured grids. An unstructured grid may consist of geometric elements point (0D) line (1D) triangle, quadrangle (2D) tetrahedron, hexahedron, prism, pyramid (3D). – p. 74/123

Example of file in UCD format for AVS # UCD file format for AVS 22 44 1 0 0 0 0.292053 0.292053 0.292053 1 0.292053 0.292053 0.892053 ... 21 0.292053 0.292053 1.00005 1 1 tet 12 2 7 0 ... 43 1 tet 19 21 17 1 44 1 tet 21 19 18 1 1 1 variable 0 0.433753 1 -0.296865 ... 21 -0.369419 – p. 75/123

Example of file in dx format for OpenDx # dx file format for OpenDx unstructured grid object 1 class array type float rank 1 shape 3 items 22 data follows 0.292053 0.292053 0.292053 ... 0.292053 0.292053 1.00005 object 2 class array type int rank 1 shape 4 items 44 data follows 12 2 7 0 20 0 17 1 ... 21 19 18 1 attribute "element type" string "tetrahedra" attribute "ref" string "positions" object 3 class array type float rank 0 items 22 data follows 0.433753 ... -0.369419 attribute "dep" string "positions" object "irregular positions irregular connections" class field component "positions" value 1 component "connections" value 2 component "data" value 3 – p. 76/123 end

Example of file in dx format for OpenDx # dx file format for OpenDx structured grid object 1 class gridpositions counts 10 10 10 origin 0.005 0.000 0.005 delta 0.010 0 0 delta 0 0.010 0 delta 0 0 0.010 object 2 class gridconnections counts 100 101 99 attribute "element type" string "cubes" attribute "ref" string "positions" object 3 class array type float rank 0 items 1000 data follows -0.200 ... -0.369419 attribute "dep" string "positions" object 4 class field component "positions" value 1 component "connections" value 2 component "data" value 3 end – p. 77/123

Example of file in vtk format # vtk file format Population_Inversion ASCII DATASET UNSTRUCTURED_GRID POINTS 6655 float 0.853816 0.0 529.0 0.768435 0.0853816 530.0 ... CELLS 5000 45000 8 132 121 122 133 11 0 1 12 8 5927 5916 5917 5928 5806 5795 5796 5807 ... CELL_TYPES 5000 12 12 ... POINT_DATA 6655 SCALARS Population_Inversion float 1 LOOKUP_TABLE default 5.72747e+11 4.47913e+11 – p. 78/123

Visualization Using VTK vtk has 14 different cell types. Some of them are: type [1]: VTK VERTEX type [12]: VTK HEXAHEDRON type [14]: VTK PYRAMID In order to visualize data one can use a vtk-file viewer like paraview or one applies vtk library to visualize data directly. – p. 79/123

Interpolation between Grids Assume that a finite element function u is given on a triangulation T h 1 . How to find the values of u on a grid Ω h 2 ? How to find the triangles T p for every p ∈ Ω h 2 ? – p. 80/123

Test for one Triangle Let P 1 , P 2 , P 3 be the corners of one triangle. Is a certain point P contained in the triangle P 1 P 2 P 3 ? Let ( ξ, η ) be such that P = P 1 + ( P 2 − P 1 ) ξ + ( P 3 − P 1 ) η. Then P is contained in the triangle P 1 P 2 P 3 if and only if ξ + η ≤ 1 ξ, η ≥ 0 . and Such a test for all points P ∈ Ω h 2 and triangles T h 1 is very time consuming. – p. 81/123

From Structured to Unstructured Grid Let the structured grid be Ω h 1 = { ( x 0 + h 1 i, y 0 + h 1 j ) | i, j = 0 , ..., N } Then, by a simple indices calculation one obtains the index i ′ , j ′ such that p ∈ ( x 0 , y 0 ) + h 1 [ i ′ , i ′ + 1] × h 1 [ j ′ , j ′ + 1] . – p. 82/123

From Unstructured to Structured Grid Let the structured grid be Ω h 2 = { ( x 0 + h 2 i, y 0 + h 2 j ) | i, j = 0 , ..., N. } Now perform the following steps: 1. For every triangle T = T (( x 1 , y 1 ) , ( x 2 , y 2 ) , ( x 3 , y 3 )) ∈ T h 1 consider the quadrangle Q = [(min( x 1 , x 2 , x 3 ) , min( y 1 , y 2 , y 3 )) , (max( x 1 , x 2 , x 3 ) , max( y 1 , y 2 , y 3 ))] . 2. For every p ∈ Q ∩ Ω h 2 , set T ( p ) = T , if p ∈ T . This means store the index of T at p , if p ∈ T . 3. Test if T ( p ) is set for every p ∈ Ω h 2 . If not, then calculate the next point q ∈ Ω h 2 from p such that T ( q ) is set and T ( p ) = T ( q ) . 4. Interpolate data from Ω h 1 to Ω h 2 by using T ( p ) . – p. 83/123

From Unstructured to Unstructured Grid Construct an auxiliary structured grid such that the domain of this grid contains the domain of the two unstructured grids. The meshsize of the auxiliary structured grid should roughly be the meshsize of the two unstructured grids. Then, for every triangle T ∈ T h 1 , put T in the cell c of the structured grid, if c intersects with T . (see “From Unstructured to Structured Grid”). This means let T ∈ T c , if T ∩ c � = ∅ . For every p ∈ Ω h 2 , find the structured cell c such that p ∈ c . Then, find triangle T ∈ T c such that p ∈ T . – p. 84/123

Heating of a Body u z u y u x original body heating of the body leads to deformed body – p. 85/123

Linear Elasticity Let Ω ⊂ R 3 be the domain of the body. Let T 0 ∈ R be the original temperature of the body. Let T : Ω → R be the temperature of the body after heating.   u x    : Ω → R 3 be the deformation vector of Let � u = u y  u z the body after heating. Problem: Let Ω , T 0 , T be given. Then, calculate � u . – p. 86/123

Symmetric Derivative u : Ω → R 3 . The symmetric derivative is defined by: Definition 9. Let �   ∂u x ∂x   ∂u y   ∂y   ∂u z   � �   ∂z : Ω → R 6 u := D�   ∂y + ∂u y 1 ∂u x   2 ∂x � �     ∂u y 1 ∂z + ∂u z   2 ∂y � ∂u x � 1 ∂z + ∂u z . 2 ∂x – p. 87/123

Symmetric Derivative Another notation is � ∂u i � 1 + ∂u j := ǫ ij 2 ∂x j ∂x i ( i,j ) ∈ Φ   ǫ 11     ǫ 22     ǫ 33   := D� u ,   ǫ 12     ǫ 13   ǫ 23 where Φ = { (1 , 1) , (2 , 2) , (3 , 3) , (1 , 2) , (2 , 3) , (3 , 1) } . – p. 88/123

Symmetric Divergence � � � � 1 ∂σ ij ∂σ ij div ( σ ij ) ( i,j ) ∈ Φ := Σ j e i + Σ i e j . 2 ∂x j ∂x i � � div ( σ ij ) ( i,j ) ∈ Φ is the adjoint operator of D� u in the following sense: � � � � v d ( x, y, z ) = − div ( σ ij ) ( i,j ) ∈ Φ ( σ ij ) ( i,j ) ∈ Φ Dv d ( x, y, z ) Ω Ω Observe, that for a symmetric matrix ( σ ij ) ( i,j ) ∈ Φ : � � ∂σ ij div ( σ ij ) ( i,j ) ∈ Φ = Σ j e i . ∂x j – p. 89/123

Linear Elasticity Definition 10. Let E > 0 and 0 < ν < 1 2 be the physical constants E-Modul and Poisson ratio. Then, define the matrix   1 − ν − ν    − ν 1 − ν 0      C − 1 = 1 − ν − ν 1   ,   E 1 + ν     0 1 + ν   1 + ν – p. 90/123

Linear Elasticity The deformation vector of the body satisfies the equations:   α     α       α   ( T − T 0 ) + C − 1 σ, D� u =     0       0   0 div ( σ ) = 0 , where α is a physical constant and σ is called stress vector. – p. 91/123

Weak Formulation of Linear Elasticity Define the bilinear form a : ( H 1 (Ω)) 3 × ( H 1 (Ω)) 3 → R � ( Du ) T C Dv d ( x, y, z ) . �→ ( u, v ) Ω Let v ∈ ( H 1 0 (Ω)) 3 . Then, we obtain   α     α     �   α   a ( u, v ) = div C ( T − T 0 ) v d ( x, y, z ) .     0 Ω       0   0 – p. 92/123

Matrix of Linear Elasticity A short calculation shows that a ( � v ) = u,� � � ∂x + ν ∂u y 2 (1 − 2 ν ) ∂u y (1 − ν ) ∂u x ∂v x ∂v x ∂x + ν ∂u z ∂v x ∂x + 1 ∂v x = ∂y + F ∂x ∂y ∂z ∂x Ω 1 ∂y + 1 ∂z + 1 2 (1 − 2 ν ) ∂u x ∂v x 2 (1 − 2 ν ) ∂u z ∂v x 2 (1 − 2 ν ) ∂u x ∂v x ∂z + ∂y ∂x ∂z ∂v y 2 (1 − 2 ν ) ∂u y ∂v y ∂v y 1 2 (1 − 2 ν ) ∂u x ∂x + 1 ∂x + ν ∂u x ∂y + ∂y ∂x ∂x (1 − ν ) ∂u y ∂v y ∂v y ∂v y 2 (1 − 2 ν ) ∂u y ∂v y ∂y + ν ∂u z ∂y + 1 2 (1 − 2 ν ) ∂u z ∂z + 1 ∂z + ∂y ∂z ∂y ∂z 2 (1 − 2 ν ) ∂u y 1 2 (1 − 2 ν ) ∂u x ∂v z ∂x + 1 2 (1 − 2 ν ) ∂u z ∂v z ∂x + 1 ∂v z ∂y + ∂z ∂x ∂z � ∂z + ν ∂u y 2 (1 − 2 ν ) ∂u z 1 ∂v z ∂y + ν ∂u x ∂v z ∂v z ∂z + (1 − ν ) ∂u z ∂v z d ( x, y, z ) . ∂y ∂x ∂y ∂z ∂z – p. 93/123

RHS of Linear Elasticity The right hand side can be written as � � � (1 + ν ) α T ∆ T ∂v x ∂x + (1 + ν ) α T ∆ T ∂v y ∂y + (1 + ν ) α T ∆ T ∂v z − d ( x, y, z ) . F ∂z Ω (16) For the implementation of the right hand side, it is helpful to sort the above terms:   F 1      . F 2 (17)  F 3 – p. 94/123

Boundary Conditions for Linear Elasticity Let Γ D ⊂ ∂ Ω be the fixed boundary of the deformation process. Let the rest of the boundary be free. Then, define V = { v ∈ H 1 (Ω) | v | Γ D = 0 } . Problem 5 (Weak formulation with boundary condition) . Find u ∈ V such that   α     α     �   α   a ( u, v ) = div C ( T − T 0 ) v d ( x, y, z )     0 Ω       0   0 for every v ∈ V . – p. 95/123

FE Discretization of Linear Elasticity Let V h be the space of trilinear finite elements. Then, define v ∈ ( V h ) 3 | � � V h = { � v | Γ D = 0 } . Problem 6 (Weak formulation with boundary condition) . Find u h ∈ � V h such that   α     α     �   α   a ( u h , v h ) = div C ( T − T 0 ) v h d ( x, y, z )     0 Ω       0   0 for every v h ∈ � V h . – p. 96/123

Rigid Body Modes Consider the vector space functions                 1 0 0 y 0 z                             M := span  ,  ,  , − x  ,  , 0 1 0 z 0              − y − x 0 0 1 0 M is the kernel of the bilinear form a . This means ∀ � m ∈ M , ∀ � v ∈ V. a ( � v ) = 0 m,� Therefore, in case of Neumann boundary conditions, we have to construct V such that V ∩ M = { � 0 } . In case of pure boundary conditions, define V as follows: V = { v ∈ H 1 (Ω) | � v, � m � = 0 ∀ � m ∈ M} . – p. 97/123

Superconvergence of the Gradient The stress σ has to be calculated by D� u . The finite element theory for linear and trilinear finite elements shows � D� u − D� u h � L 2 = O ( h ) . This is a slow convergence. But one can prove the following superconvergence of the gradient in case of structured grids: Let Σ h be the cell points of the structured grid. Then, for a sufficient smooth solution � u and a not complicated boundary Γ , we obtain: u h )( p ) � = O ( h 2 ) . p ∈ Σ h � ( D� max u − D� Therefore, in case of linear elasticity, apply trilinear elements on a block-structured grid or quadratic finite elements. – p. 98/123

Fluid Dynamics Let us describe a two dimensional flow by: u x-component of the velocity vector of the flow, v y-component of the velocity vector of the flow, p pressure of the flow. (u,v) v u – p. 99/123

Navier-Stokes-Equations ∂x + ∂ ( u 2 ) + ∂ ( uv ) 1 ∂u ∂t + ∂p = Re ∆ u ∂x ∂y + ∂ ( v 2 ) ∂y + ∂ ( uv ) 1 ∂v ∂t + ∂p = Re ∆ v ∂x ∂y ∂u ∂x + ∂v = 0 ∂y There exist different kind of boundary conditions: input, output, slip, and no-slip boundary conditions. – p. 100/123