Spatial Discretisation V. Selmin

1. Spatial Discretisation V. Selmin Multidisciplinary Computation and Numerical Simulation

2. Outline • Conservation principle & differential equations • Grid topology • Finite Difference Discretisation • Finite Volume Discretisation • Finite Element Discretisation • - Approximation by finite elements • - Integral methods • - Numerical integration Outline

3. Conservation principle & differential equations Conservation principle & differential equations

4. Grid Topology Structured Grids • Mesh structure: • Domain divided into a structured assembly of • quadrilateral cells • Each interior nodal points is surrounded by • exactly the same number of mesh cells (or • elements) • Directions within the mesh can be immediately • identify by associating a curvilinear co-ordinates • system • it is possible to immediately identify the nearest • neigthbours of any node j on the mesh • Advantages: • Large number of algorithms for discretisation are • available • The algorithms can be normally implemented in • a computationallly efficient manner • Disadvantages: • Difficulty to generate grids of regions of general • shapes  multi-block grids • Very high elapsed time necessary to produce a • grid for domains of extremely complex shape Finite difference, finite volume and finite element discretisations Grid topology – structured gris

5. Grid Topology Unstructured Grids • Mesh structure: • Computational domain divided into an unstructured • assembly of computational cells • The number of cells surrounding a typical interior • node is not necessarily constant • The nodes and the elements has to be numbered • To get the necessary information on the neightbours • the numerotation of the nodes wich belong to each • element has to be stored • The concept of directionality does not exist anymore • Advantages: • Powerful tool for discretising domains of complex • shapes • Unstructured mesh methods naturally offer the • possibility of incorporating adaptivity • Disadvantages: • Alternative solution algorithms are more limited • Computational implementation places large demands • on both computer memory and CPU Finite volume and finite element discretisations Grid topology – unstructured gris

6. Finite Difference Discretisation Model differential equation Strictly structured grids Unknowns at grid nodes Divided differences & Taylor-series expansion First derivative definition Computational stencil Taylor-series expansion Examples Grid of constant mesh size Finite difference discretisation

7. Finite Volume Discretisation Model differential equation Structured/unstructuredgrids Unknowns at grid nodes or at cell centre Integral formulation & divergence theorem Integral formulation Divergence theorem Cell-centred control volume Examples Node-centred control volume Control volume types Cell-centred, node-centred, cell-vertex Cell-vertex control volume Finite volume discretisation

8. error approximate function exact function Nodal approximation if Finite Element Approximation Approximation error Approximation function • Select a finite set of functions which depend of n parameters • Compute the parameters in order to minimise the error, for example by taking • in n points , i.e. by cancelling the e(x) at this n points • These functions are select in such a way to be easily estimated on computers, • to be easily integrated or differentiated explicitely • In general, the approximated function is linear in the parameters Finite element discretisation Finite element discretisation

9. The subdomain is the element The coordinates are the nodal coordinates. The quantities are the nodal variables Finite Element Approximation Sub-domain approximation • The sub-domain nodal approximation method simplifies the construction • of u(x) and is adapted to numerical computation. • It consists in • identify a set of sub-domain of the domain • define an approximated function different on each sub-domain • by using the nodal approximation method • Each function may be defined with respect to nodal variables belonging • to other sub-domains (spline method) Finite element approximation • The finite element approximation method is a particular case of the previous • method with the following characteristics: • The nodal approximation on each sub-domains only uses nodal • variables related to nodes on and on its boundary • The approximated functionson each sub-domain are built in • order to be continuous on and they satisfy continuity conditions • between the different sub-domains Finite element discretisation

10. : shape funtion at node i of element e • The finite element approximation presents the following features: • The geometry of the elements has to be defined analytically, • that is more or less complex depending of their shapes • The interpolation functions related to each element • have to be built Finite Element Approximation Finite element approximation Properties of the shape functions For each element e Compact support Finite element discretisation

11. Finite Element Approximation Element partitioning rules • The partitioning of domain in elements has to satisfy the following rules: • Two different elements can have in common only nodes that belongs to their common boundary. • In particular, overlap between elements must be avoided. • The set of all the elements must form a domain as much as possible close to the given domain • In particular holes between elements must be avoided • The two previous rules are satisfied if the elements are built as follows: • Each element is univocally defined from the co-ordinates of the geometrical nodes of the element. • In most cases, those nodes lay on the element boundary and are common to several elements. • The boundary of a 2D or 3D elements is formed by a set of curves or surfaces. Each boundary • portion has to be univocally defined from the co-ordinates of the only geometrical nodes located • on the portion of boundary. So, the boundary portions that are common to two elements are • identically defined for the two elements • Moreover, the geometrical discretisation error at boundaries of the domain may be reduced by • decreasing the element size or by using elements that allow more complex boundaries. Finite element discretisation

12. Classical Elements Shapes 1-D Elements linear (2) quadratic (3) cubic (4) 2-D Elements Triangular elements linear (3) quadratic (6) cubic (9) Quadrilateral elements linear (4) quadratic (8) cubic (12) Finite element discretisation

13. Tetrahedral elements linear (4) quadratic (10) cubic (16) Hexahedral elements linear (8) quadratic (20) cubic (32) Prismatic elements linear (6) quadratic (15) cubic (24) Classical Elements Shapes 3-D Elements Finite element discretisation

14. In order to simplify the analytical definition of elements of complex shape, the notion of reference elements has been introduced: A reference element is an element of very simple shape that can be transformed in each physical element by using a geometrical transformation Through this transformation, the co-ordinates of each node of the physical element may be defined from the co-ordinates of the corresponding node of the reference element Reference element Physical element The transformation depends on the shape and the location of the physical element, and in particular from the co-ordinates of its geometrical nodes. The transformation is different for each physical element where are the co-ordinates of the geometrical nodes that belong to the element e Reference Elements Finite element discretisation

15. Reference element Physical element Reference Elements • The transformations must have the following properties: • It is a bijection in all the points of the reference element or of its • boundary: then to one point of corresponds one and only one • point of and inversely • The geometrical nodes of the reference element corresponds to • the geometrical nodes of the physical element • Each portion of the boundary of the reference element , defined • by the geometrical nodes, corresponds to the portion of the • boundary of the physical element defined though the related • nodes • A linear transformation with respect to the coordinates of the • geometrical nodes of the physical element is in general used where are in general polynomial function in and are named geometrical transformation functions Note that those functions are identical for all the coordinates components Finite element discretisation

16. Reference Elements Three nodes triangular element The reference element is defined analytically by: By considering the transformation linear in : • The three following properties are verified: • The geometrical nodes of of co-ordinates (0,0), (1,0) and (0,1) transform themself in the geometrical nodes • of of co-ordinates . • Each boundary of is transformed in the corresponding boundary of . • For example, the boundary going through the nodes (1,0) and (0,1), whose equation is • transformed into the boundary of going through and whose parametric equation is Finite element discretisation

17. Reference Elements • The transformation is bijective if the matrix is not singular The determinant is equal to two times the triangle area and is zero only if the nodes are aligned. In general the internal angles of a triangle or of a quadrilateral must be less that 180 degrees in order to avoid • Remarks: • The geometrical transformation may be interpreted as a simple change of variables • may also be considered as a local co-ordinates system related to each element Finite element discretisation

18. Classical Reference Elements 1-D Elements 2-D Elements Finite element discretisation

19. Classical Reference Elements 3-D Elements Finite element discretisation

20. interpolation functions nodal variables are part of the polynomial basis used to build the are independent monomial functions for a linear triangle Reference Elements Approximation on the reference element On each element it is possible to replace the approximation on the physical element by an approximation on the reference element • Properties of the shape functions: • Remarks: • In general, the interpolation functions are used only for very simple elements. • They are replaced by the functions where and are linked through the • transformation • The are independent of the geometry of the physical element. • The same functions may be used for all the elements that own the same reference element • characterised by itsform, itsgeometrical nodes, itsinterpolation nodes. Finite element discretisation

21. The equations of a physical problem to be studied are described on the physical domain and make use of the unknown functions and of its derivatives in x: , etc. As the approximation on the physical elelement is often complex, the approximation on the reference element will be systematically used which is associated to the transformation The transformation being bijective, Although the inverse transformation always exists, its is not easy to built it explicitely. (Only for linear transformations) In the case its is known explicitely, the expression of can be used in order to get the approximation on the physical element: A more simple approach is to transform the derivatives of with respect to x to derivatives with respect to through the Jacobian of the transformation Transformation of the differential operators Finite element discretisation

22. Transformation of the differential operators First derivatives We use the chain rule for the computation of derivatives of a function with respect to from the knowledge of its derivatives with respect to x : where J is the Jacobian of the geometrical transformation. In a similar manner, the derivatives of a function with respect to x may be obtained from the derivatives in It is the matrix j which is used in practice since we have to express the derivatives of u with respect to x,y,z from the derivatives of u with respect to The transformation of being assumed bijective, the inverse of J exists in any point of the reference element Finite element discretisation

23. Transformation of an integral The change of variables allows to move from the integration of a function on the physical element to a simpler integration on the reference element : where is the determinant of the Jacobian matrix . As a matter of fact, the element of volume is the product In a cartesian system, In the curvilinear system, and the mixed product may be written as Finite element discretisation

24. Singularity of the Jacobian matrix The singularity of J on a point of the reference element implies that the transformation is no more bijective This singularity appears when the reference element is strongly deformed. It is advisable to verify that maintains a constant sign in each point of the reference element Example: Singularity within a quadratic transformation in one dimension By choosing . Then the determinant of the Jacobian may be expressed according to This determinant becomes zero at the point The point is contained into the reference element only if Finite element discretisation

25. Approximation Error on an Element • The approximation error is defined as • in any point x of and in any point of • For a point x that corresponds to a point through the transformation the relation holds. • In order to characterise the value of the maximum error, we will use the maximum norm of the function . • Let define the error on each of the derivatives of order by: • The corresponding norm may be written according to • Maximum on of • for any for which • The following semi-norm is often used in the finite element method: • Strang gives the following expression for the norms and : • where • c and C depend only of of the element type and of the approximation used • the polynomial basis is complete up to the order • l is related to the maximum dimension of the element • the derivatives of the approximation up to order s are bounded • is the norm of with Finite element discretisation

26. Approximation Error on an Element • For s=0 , the norms take the form • In the case of a one-dimensional linear element ( ) • Improvement of the approximation accuracy • In order to improve the accuracy of the approximation, it is then requested to • decrease l , that means decrease the element dimensions • increase , that means to use an approximation whose polynomial basis is complete up to a higher order • The following techniques may be used • Decrease the dimension of the element and consequently increase the number of elements needed to cover • Increase the order of the approximation polynomials, which leads to an increase of the number of nodal variables • or of the degrees of freedom of each element. That can be done by: • increase the number of interpolation nodes of each element, maintaining one nodal variable per node, • that leads to the Lagrangian family of element type • increase of the number of nodal variables at each node, maintaining the same number of nodes. The additional • nodal variables are the values at nodes of , that leads to the Hermitian family of element type. Finite element discretisation

27. Physical system Physics laws, engineer science Formulation of the equations Partial Differential Equations Weighting residuals methods Transformation of the equations Integral Formulation Approximation of the unknown functions & matricial organisation System of Algebraic Equations Numerical solution Numerical solution of the system Approximated Solution Integral Formulation Transformation of the equations of a physical system Finite element discretisation

28. Integral Formulation Weighted Residual Method Let consider a stationary continuous physical models, whose behaviour is represented by a system of partial differential equations, linear or nonlinear of order m: on the domain V The boundary conditions being written according to on the boundary S The unknow variables depends from the co-ordinates x. Functions u are solution of the equilibrium problem if they satified both the equations in the field and on the boundary. We named residual the quantity defined by: which becomes zero when u is solution of the problem. The weighted residual method consists to search for the functions u that annul the integral form for any weighting functions that belongs to a set of functions . The functions u belong to the set of the admissible solutions that satisfy the boundary conditions and that are differentiable up to order m . Finite element discretisation

29. Integral Formulation • Transformation of the integral form • The integration by part allows to transform the integrals in such a way to decrease the conditions imposed • to the admissible functions u. • Let recall the integration by part formula: (use of Gauss theorem) • The integration by parts leads to integral forms named weak form that leads to the following advantages • the maximum order of the derivatives of u that appears in the integral form decreases. The conditions of • differentiability on u are less strong. • some of the boundary conditions that appear in the weak form may be taken into account into the integral • formulation, instead to be identically satisfied by u . • On the other hand, derivatives on appear in the integration by parts. That means that the conditions of • differentiability on increase. In addition, could have to satisfy conditions on a part of the boundary • in such a way that some boundary terms may disappear. • The solution of the partial differential problem is approximated by the solution of the weak integral form, even • if this solution does not satisfy the differentiability conditions of the original problem. Finite element discretisation

30. Integral Formulation • Discretisation of the integral forms • The solution of the partial differential equations has been replaced by the search for the functions u that annul • the integral form • for any function . • In order to build an approximated function u , the integral form is discretised in two steps: • Choose a n-parameters approximation of the unknown functions u , that may be expressed as • The integral equation becomes • for any . • Choose a set of n independent weighting functions . The number of weighting functions must be • equal to the number of parameters of the approximation. The choice of the function leads to different methods. Finite element discretisation

31. Integral Formulation Choice of the weighting functions 1 – Collocation per points The function is the Dirac distribution at the point . The integral form reduces to The accuracy of the solution depends of the choice of the points . The number of collocation points must be equal to the number of parameters. In practice, this method is not very used since it is difficult to implement it in a finite element approximation. In addition, it leads to non symmetrical system of equations. In another hand, it has the advantage to avoid integration on the volume, which may be interesting for some non linear problems. 2 – Collocation by sub domain We choose n sub-domains and take as functions : The accuracy of the solution depends of the choice of the sub-domains . The number sub-domains must be equal to the number of parameters. In practice, this method is not very used since the choice of the sub-domains is difficult. In another hand, as integrations on the volume are needed, it is better to use the Galerkin method Finite element discretisation

32. Integral Formulation 3 – Galerkin Method The function is built by considering the set of variations of the functions u . where are the variations of the approximation parameters The integral formulation reduces to Since W has to becomes zero for any , the previous relation is equivalent to n algebraic equations The system is symmetric if the operator L is self-adjoint. 4 – Least square method The least square method consists in the minimisation of the expression with respect the n approximation parameters. R is the residual. This method is not very used because it does not all the integration by parts, and needs more strict conditions on the approximation function that the Galerkin method. On another hand, it leads to a definite-positive symmetrical system for any operator L. Finite element discretisation

33. Finite Element Method The finite element method consists in using a finite element approximation of the unknown functions u in order to discretise an integral form W ; afterwards, to solve the resulting system of algebraic equations. We will use integral forms of Galerkin type for which the weighting functions are : Let replace this integral by a sum of integrals on each element : In order to compute each term , named elementary integral form, we use a finite element approximation of u and on each element : As the are zero in any point outside to , and as are only nodal variables of the element , each term is computed from the only variables associated to the element e . This property has contributed to the success of the finite element method. Then, becomes In general we take advantage of the integration by parts in order to decrease the maximum order of the derivatives. Finite element discretisation

34. Integral Formulation on the Reference Space Transformation of the derivatives The derivatives in x are expressed with respect to derivatives in and to the inverse of the Jacobian matrix of the geometrical transformation. Transformation of the integration domain The integration on the volume is replaced by the integration of the reference element : The integration limits in for the classical reference elements are 1-D 2-D Triangle Quadrilateral 3-D Tetrahedra Hexahedra Prism Finite element discretisation

35. Numerical Integration The functions to be integrated are in general complex polynomials or rational fractions . Their explicit integration is in general not easy. The numerical integration consists in replacing the evaluation of an integral by a linear weighted combination of the value of the function to be integrated at selected points: where: are the co-ordinates of the r integration points are the associated weighting coefficients Numerical integration in 1-D 1 – Gauss Method The method of Gauss is a method of numerical integration for which the r coefficients and the r abscisses are determined in such a way to integrate exactly polynomials of degree Replace the integral of a polynomial function by a linear combination of its values at the integration points Compute the 2r coefficients in such a way that the previous equation is verified exactly for the following polynomial That leads to the following equation Finite element discretisation

36. Numerical Integration This equation is verified for any if This system of 2r equations is linear in and non linear in ; it allows to compute the 2r parameters under the conditions: The abscisses are also the roots of the Legendre polynomials of order r: defined by the recurrence formula: The weights are written as: The integration error takes the form: Finite element discretisation

37. Numerical Integration Finite element discretisation

38. Numerical Integration 2 – Newton Cotes Method If the abscisses are a priori fixed, it remains r coefficients to evaluate in such a manner that a polynomial of degree r-1 is integrated exactly. In the method of Newton-Cotes the points are regularly spaced and symmetric with respect to : In order to compute the coefficient , let represent by a Lagrange polynomial of degree r-1 that takes the values at r integration points : Then : and the weights are the integrals of the Lagrange interpolation functions : The integration error is of the form: Therefore, it advisable to use an odd number of integration points Finite element discretisation

39. Numerical Integration Finite element discretisation

40. Numerical Integration • Remarks: • For a given number of integration points, the maximum degree of the polynomials integrated exactly by the • method of Newton-Cotes is less that those obtained by the method of Gauss. • Nevertheless, those method allows sometimes to have the integration points and the interpolation nodes coincident. • The integration of the terms that contains the interpolation functions is then simplified since are zero at all • the integration points other that . • The method of Gauss at r points integrates exactly a polynomial of order 2r-1. • The method of Newton-cotes integrates exactly a polynomial of order r-1. • In practice, we have to integrate polynomials of high order or non-polynomial functions (rational fractions). • By using the previous methods, the approximate integration as much accurate as the number of integration points • is high. Finite element discretisation

41. Numerical Integration Numerical integration in 2-D A - Product Method It consists to use in each direction and a 1-D numerical integration. If we use points in the direction and points in the direction, the method of Gauss integrates exactly the product of a polynomial in of order order and of a polynomial in of order . B - Direct Method It is possibile to extent directly to 2-D the methods described in the previous section: In particular, we can build methods of Gauss that integrate exactly all the monomials of order m: These methods often uses less points that the methods “product”. For quadrilateral elements, the methods “product” are the more often used, while for triangular elements the direct methods are more commonly used. Finite element discretisation

42. Numerical Integration 1 - Quadrilaterals The “product” method leads to where the are the coefficients of the Gauss or Newton-Cotes methods and are the associated integration points. The direct method leads to Finite element discretisation

43. Numerical Integration 2 - Triangles 2-1 Gauss-Radau method The “product” method consists first to transform the integral on the triangle into an integral on a square by using a change of variables The integral becomes: The numerical integration gives This method, named Gauss-Radau, is not often used because the location these integration points do not allow to enforce the symmetry of the triangle. On another hand, it can be effective when the variations of are very strong in the neighborhood of the node A of the triangle. Finite element discretisation

44. Numerical Integration 2-2 Direct method Finite element discretisation

45. Numerical Integration Numerical integration in 3-D 1 – Cube reference element The “product” method leads to where the are the coefficients of the Gauss or Newton-Cotes methods and are the associated integration points. The direct method leads to Finite element discretisation

46. Numerical Integration 2 - Tetrahedra reference element The method of Gauss-Radau is rarely used. The direct method leads to Finite element discretisation

47. Reference element Physical element Integral method Integration by parts Weighted residuals Galerkin method Physical space Reference space Physical element Reference element Numerical integration Gauss method Summary Function approximation PDE discretisation method Numerical integration Finite element discretisation