1D hp-ADAPTIVE FINITE ELEMENT PACKAGE FORTRAN 90
Contents
1. au v f bu v dx f cuv dx n fio das 2 6 0 0 0 0 At this point further derivation depends upon the kind of the boundary conditions being used In the case of Dirichlet boundary condition we eliminate the boundary term by restricting ourselves to only those test functions that vanish at that point For example if we assume Dirichlet boundary condition at x 0 u 0 7 0 2 7 we have to assume that v 0 0 In the case of Cauchy or Neumann boundary condition we build the boundary condition into the formulation by using it to calculate the derivative in terms of the boundary data and in the case of Cauchy boundary condition still unknown solution at that point Thus in the case of Dirichlet boundary condition at 0 and Cauchy boundary condition at x l we get Find u x u 0 y l l n au v bu v cuv dz Bu I w 0 n f x v dz yu 0 2 8 0 0 for every test function v x such that v 0 0 In the case of Neumann boundary condition at x l constant jj 0 and the boundary term simply vanishes Notice that we have kept all terms involving solution u r on the left hand side and the terms involving the test function v only have been moved to the right hand side of the equation The weak variational formulation can be shown to be equivalent to the classical form of the boundary value problem We integrate back by parts and use the Fourier lemma argument to recover the differential equation and Cauch2. Error tolerance tol 0 01 per cent of the energy norm of the solution Figure 9 Example 1 Final hp mesh with the corresponding indistinguishable numerical and exact solutions Figure 10 Example 1 Rates of convergence This is a very smooth solution The algorithm chooses from start to use p refinements only In this case the hp method reduces just to the p method Decreasing the error toller ance reveals also that as the number of degrees of freedom grows the order of approximation pis increased uniformly This is consistent with the hp approximation theory see e g 2 Example 2 A singular solution Ueract x99 Error tolerance tol 1 per cent of the energy norm of the solution Please use routine result to verify that the mesh is geometrically graded towards the singularity at x 0 with the order of approximation p increasing linearly from p 1 at the Figure 11 Example 2 Final hp mesh with the corresponding indistinguishable numerical and exact solutions Figure 12 Example 2 Rates of convergence singularity to p 3 away from it This is consistent with the well known result of Babuska and Gui see e g 2 Example 3 A solution with an internal layer u amp 44 r atan60 x 5 Error tolerance tol 1 per cent of the energy norm of the solution When the error tollerance is decreased the algorithm continues to choose exclusively the p refinements only Moreover th
3. 4 p 1 Note that except for the first two linear functions the remaining shape functions vanish at the element endpoints For that reason they are frequently called the bubble functions see Fig 1 Finally we introduce the notion of the hp interpolation over the master element Given a continuous function u 0 1 we define its hp interpolant tiny as the sum of the standard linear interpolant thy coinciding with function at the endpoints and a Hi projection of the difference D onto the span of the bubble functions introduced above More precisely dnp thy Uy 3 26 where inp 0 x1 G 1 R E 3 27 and p 1 a WHO 3 28 j 3 where the coefficients us are determined by solving the following system of equations pti dX dXx ie d ti thy d 2j SY f Ok 3 1 3 29 Dij aa m 3 29 3 2 2 A 1D parametric element of arbitrary order We consider now an arbitrary closed interval K x xg C 0 and assume that K is the image of the master element through some map zg K 0ljas o2 cx eK 3 30 The simplest choice is an affine map which may be convenietly defined through the master element linear shape functions rK vLX v nX amp L 1 rR 3 31 zr F zn TL zr Ehk where hy zg zr is the element length We assume that the map is invertible with inverse ax and define the element space of shape functions X K as the
4. one 5 2 Organization of the code The code is organized in the following subdirectories e blas_pack basic linear algebra routines e commons system common blocks e data pack exact solution and material data e da strs pack data structure routines e errest pack error estimation routines e clem pack element routines e files system files sample input files e frontsol_pack frontal solver routines e gcommons graphics common blocks e graph_1D actual graphics routines for the code e graph util graphics utilities e graph interf graphics interface routines e main program driver for the code e meshgen pack initial mesh generation routines e meshmods pack mesh modification routines e module pack data structure moduli e solver pack interface with frontal solver e utilities pack general utilities 5 3 Mesh generation and postprocessing routines Mesh generation datastrs pack meshgen f A sequence of input data is files input These include read from e MAXELES MAXNODS maximum number of elements and nodes in the mesh the corresponding memory is allocated dynamically e NRELIS number of uniform size elements in the initial mesh e norder uniform order of approximation for the initial mesh elements e XL length of interval 0 to set up the boundary value problem e ibcl ibc2 boundary conditions flags A sample input file input can be fo
5. 23 4 1 A priori error estimation oo soo o qe V D E a 23 4 2 A posteriori error estimation i wow be ovCE woe a wo Wo X E s Robe Xo 24 4 2 1 The element implicit residual method 25 5 User Manual 27 5 1 Data structure in 1Dhp90 uuo Eu ee phus Re eee OH n OC es 27 5 2 Organization of the code sad Pci e ed Tecie eee o xc ioa Be 27 5 3 Mesh generation and postprocessing routines llle 28 5 4 Processing algorithms a doxo3o a of eg ee ANE EUEOK A SERE URSUS 29 5 4 1 Evaluation of element matrices elem pack elem f 29 5 4 2 Modification of the element matrices due to boundary conditions 30 5 4 3 Assembling global matrices exer rene een Rer co 30 Dd SSOLIVOL tasses n ure Ios E Eun By X tme Xe o qut ri pet FU 32 54 5 Setting ap datas Files s cede duo Zee ty ak ow et per Sok of doo sex Geo 32 6 Adaptivity 34 6 1 p refinement unrefinement meshmods pack modord f 34 6 2 h refinement meshmods_pack break f llle 34 6 3 Natural order of elements datastrs_pack nelcon f ls 34 6 4 h unrefinement meshmods_pack cluster f een 34 6 5 The h refinements strategy 2s 35 6 6 Trading h refinements for p refinements a 35 6 7 Interactive refinements uos autom bah oe he Reg Oboe quom E gg s 36 6 8 The final hp adaptive algorithm lle 36 6 9 Examples of hp adaptive solutions 22e 36 A Interface with a frontal solver
6. Fig 2 The construction of the bubble basis functions is much easier As the element bubble shape functions vanish at the element endpoints we need simply to extend them only by One should really use a symbol Xap as the discretization depends upon the element size h hx and order of approximation p px the zero function elsewhere The support of a vertex node basis function extends over the two adjacent elements whereas for a bubble function it is restricted just to one element Finally the continuity of the approximation at the vertex nodes allows us to introduce the concept of the global hp interpolation Given a continuous function u x x a b we define its hp interpolant as the union of the contributing elements hp interpolants upy r up x where z e K 3 39 with D denoting the Ap interpolant over element A For linear first order elements the hp interpolation reduces to the standard Lagrange interpolation We emphasize that the interpolation is done ocally separately over each element 3 2 4 Element stiffness matrix and load vector Having selected an appropriate set of shape functions we calculate the element stiffness matrix and load vector d ue dx 55 f Ce x PL b Xiy exis dz da 3 40 Le j Pa The local master element coordinate system proves to be more convenient in the deriva tion of the element matrices The matrices are calculated in the local coordinate system by using the chain r
7. space of compositions of inverse zp and functions defined on the master element X K ozg Gc X K 3 32 u x where zx y and X K j Consequently the element shape functions are defined as Xk z X amp where zk z k 1 p 1 3 33 Note that in general the shape functions are no longer polynomials unless the map xx is an affine map In such a case we speak about an affine element For practical reasons most of the time it is convenient that map xy is specified using the master element shape functions i e it is a polynomial of order p OSI tx X E 3 34 In such a case we talk about an isoparametric element Coefficients zg will be identified as the geometry degrees of freedom g d o f Note that only the first two have the interpretation of coordinates of the endpoints of element K In our implementation we have restricted ourselves to the affine elements only Conse quently throughout the rest of this presentation we shall assume that the summation in 3 34 extends over the linear shape functions only The hp interpolation operator can now be generalized to an arbitrary element K The idea is to perform the interpolation procedure always on the master element Given a con tinuous function u x 7 K defined over the element K we compose it with the map transforming the master element into element K 8 woaK u zx 3 35 and find the corres
8. the point of view of minimizing the bandwidth Besides as a result of refinements unrefinements of the mesh nodes may no longer be numbered using consequtive integers We shall adopt the philosophy that we are always given an order of elements One such order called the natural order of elements is provided by routine datstrs pack nelcon and will be discussed in the next section Given the order of elements and an order of nodes for each element we can define the natural order of nodes Finally following the order of shape functions d o f for each of the nodes we can define the natural order of d o f see the discussion in the previous section On the practical level we may introduce an extra attribute for each node nod in the mesh say nbij nod equal to the number of the first corresponding d o f in the global natural order of d o f The pointers are determined in the following way initiate array nb 7 with zeros initiate d o f counter idof 0 for each element nel for each element node 1 get the global node number nod ELEM S nel nodes i skip if nbij nod Z 0 i e the node has already been visited set the counter for the first d o f of the node nbij nod idof 1 determine the number of d o f ndof corresponding to the node update the counter idof idof ndof end of loop through element nodes end of loop through elements Once the bijection between the local d o f and the global denumeration of d o f the connectivi
9. y Ir un v 2 ixl K v o em d This leads to the final estimate lirli i 4 78 5 User Manual 5 1 Data structure in 1Dhp90 We introduce two user defined structures module pack data structure D e type node e type element The attributes of a node include node type a character indicating whether the node is a vertex or middle one integer order of approximation integer boundary condition flag a real array coord containing geometrical degrees of freedom node coordinate and real array dofs containig the actual degrees of freedom Both the geometry and the actual d o f are allocated dynamically dependently upon the order of approximation for the node The entire information about a mesh is now stored in two allocatable arrays ELEMS and NODES as declared in the data structure module The module also includes a declaration for a number of integer attributes of the whole mesh The following parameters are relevant at the moment NRELIS number of elements in the initial mesh NRELES total number of elements in the mesh NRNODS total number of nodes in the mesh MAXEQNS maximum number of equations to be solved MAXNODS maximum number of nodes MAXELES maximum number of elements NREQNS actual number of equations to be solved Parameters MAXEQNS and NREQNS are placed into the code in anticipation of using the code for the solution of systems of equations In the code both parameters are set to
10. 1D hp ADAPTIVE FINITE ELEMENT PACKAGE FORTRAN 90 IMPLEMENTATION 1Dhp90 Leszek Demkowicz and Chang Wan Kim Texas Institute for Computational and Applied Mathematics The University of Texas at Austin Austin TX 78712 Abstract This manual provides a description of a 1D hp adaptive finite element code for the solution of a class of two point boundary value problems The implementation is based on hierarchial shape functions and allows for both h and p refinements of the mesh based on a posteriori error estimation The code includes an automatic hp adaptive strategy Numerical results for several test problems illustrate the method Contents 1 Introduction 4 2 A General Two Point Boundary Value Problem 5 2 1 Strong or classical form of the problem 5 2 2 Weak variational formulation of the problem oaoa 6 3 Finite Element Method 9 3 1 Galerkin approximation o ooa a 9 3 2 Finite Element Method aana a 10 3 2 1 1D master element of an arbitary order 10 3 2 2 A 1D parametric element of arbitrary order 12 3 2 3 1D hp finite element space xc e Rom RR a 14 3 2 4 Element stiffness matrix and load vector ooo a 15 3 2 5 Taking into account the boundary conditions Modified element matrices 16 3 2 6 Global stiffness matrix and load vector the assembling procedure 17 3 3 Structure of a classical FE codes iub RES eee pa ee ae ae 4 20 4 Error estimation
11. 42 1 Introduction The finite element method is a general technique for constructing approximate solutions to boundary value problems The method involves dividing the domain of the solution into a finite number of simple subdomains the finite elements and using variational concepts to construct an approximation of the solution over the collection of finite elements Because of the generality and richness of the idea underlying the method it has been used with a remarkable success in solving a wide range of problems in virtual all areas of engineering and mathematical physics The goals of this manual are e to give a brief introduction to fundamental ideas of variational formulation Galerkin method a posteriori error estimation in context of a one dimensional finite element method e to introduce the reader to the concept of hp adaptive Finite Element Methods e to provide a minimal user information for the package Consequently we have decided to write this document in a format of lecture notes Indeed it is our intention to use the notes in the ASE384P EM 394F CAM394F Finite Element Methods 4 class taught in the ASE EM Department and the CAM programs at the University of Texas at Austin 2 A General Two Point Boundary Value Problem 2 1 Strong or classical form of the problem We begin by considering a two point boundary value problem of finding a function u u r r 0 1 The strong form of the boundary value pr
12. aphical representation of the mesh and the current solution allows also for an interactive mesh modification using the mouse 6 8 The final hp adaptive algorithm Fig 8 presents the final automatic hp refinements algorithm The algorithm can be invoked from the main program menu by selecting the automatic hp refinements option Parameter tol acceptable error tollerance in per cent of the energy norm of the solution has to be input from the keyboard After each refinement the corresponding number of d o f in the mesh and the computed FE error are written to file files result automatically open by the program The data can later be used to visualize the corresponding convergence rates by selecting the option rates from the graphics program 6 9 Examples of hp adaptive solutions Problem u f x x 0 1 u 0 0 6 79 u 1 g x All convergence rates are represented using the log log scale in terms of the total number of degrees of freedom Input data and generate initial mesh Modify interactively the initial mesh Solve the problem on the current mesh Estimate the error Is the global residual acceptable YES Perform h refinements Solve the problem on the new mesh Trade h refinements for p refinements Figure 8 Flow chart for the hp adaptive FE code Example 1 A smooth solution Ucrace x sin x
13. e order of approximation p stays essentially uniform Thus we might say that the initial h refinements help to resolve the scale and to construct an optimal initial mesh only Once the mesh is determined the uniform p refinements again turn out to be asymptotically optimal Of course the point is that in practical computations we always work in the preasymptotic range only Figure 13 Example 3 Final hp mesh with the corresponding indistinguishable numerical and exact solutions Figure 14 Example 3 Rates of convergence Acknowledgements The authors are much indebted to Professors Ivo Babuska and J Tinsley Oden for numerous discussions regarding the subject of hp mesh optimization References 1 2 M Ainsworth J T Oden A Posteriori Error Estimation in Finite Element Analysis Comput Meth Appl Mech Engrg 142 1 88 1997 I Babu ka and B Q Guo Approximation Properties of the hp Version of the Finite Element Method Computer Methods in Applied Mechanics and Engineering Special Issue on p and hp Methods eds I Babu ka and J T Oden 133 319 346 1996 I Babuska and W C Rheinboldt A Posteriori Error Estimates for the Finite Element Method Int J Numer Meth Engng 12 1597 1615 1978 E Becker J T Oden and G Carey An Introduction to Finite Elements Prentice Hall 1985 L Demkowicz L J T Oden and T Strouboulis Adaptive Finite Elements for Flo
14. g natural order of elements The idea is to follow the leaves of the tree and the ordering of elements in the initial mesh compare Fig 7 6 4 h unrefinement meshmods pack cluster f A refined element can be back unrefined The corresponding entries for its element sons and their nodes are deleted from the data structure arrays We admit a situation in which the unrefined element will have a greater order of approximation than its sons The new d o f for the unrefined element are evaluated in routine elem pack project f by perform ing hp interpolation of the old representation of the solution using the clustered element shape functions Having done the interpolation projection we evaluate the corresponding interpolation error using the energy norm Figure 7 The family tree for a sequence of h refinements of elements 1 3 7 10 The dotted line indicates the natural order of elements 6 5 The h refinements strategy Having calculated the element residuals in routine errest_pack errest f we break h refine elements with the biggest residuals More precisely given a percentage perc of the total residual squared set to perc 60 in the code we reorder elements according to their residuals and refine the first M elements from the list that contribute with perc percentage to the global residual The operation is performed in routine meshmods_pack refine f 6 6 Trading h refinements for p refinements When refining the mesh we do not at
15. h SL S2 o B B Boo OEE E T FrEE 3 53 L 3 54 The 3 55 and 3 56 show the situation after assembling the third and the fourth element contributions Sh PCI Sos Sh Sh Sis Sh Sis S3 S32 S33 S31 S Sis S13 Sia 3 55 S5 Se St S Sa S5 Sh S31 Sb S58 Sa S35 Sj Sh Si Su Si S5 Sb S55 Ss S gt Si 855 Li L Li L L3 n 3 56 L Li L 3 3 Structure of a classical FE code We are now in a position to sketch the computer flow diagram of a classical finite element code as in Fig 6 A sequence of input data is read or generated in the preprocessor part This includes geometry and finite element mesh data such as number of elements element connectivity information node coordinates initial order of approximation This preprocessor may per form an automatic division of the domain into elements following some rules specified by user and it may also provide graphical information on elements and nodes The processor part includes generation of the element matrices SE and LE using numerical integration imposition of the boundary conditions assembly of global matrices and solution of the equa tions for the degrees of freedom In the postprocessor part the output data is processed in a desired format for a printout or plotting and secondary variables that are derivable from PREPROCESSOR Read input data Generate the mesh PROCESSOR Calculte the element mat
16. hich eventually are glued together into the globally defined basis functions e in the Galerkin method It is the construction of the basis functions that distinguishes the FEM from other Galerkin approximations We begin our presentation with a discussion of the fundamental notions of the master element the isoparametric element and the finite element space We shall recall the construction of the Galerkin basis functions through the element shape functions and finally introduce the notion of the hp interpolation 3 2 1 1D master element of an arbitary order Geometrically the 1D master element K coincides with the unit interval 0 1 The element space of shape functions X K is identified as polynomials of order p X K P K 3 23 Obviously one can introduce many particular bases than span polynomials of order p In the present implementation we have selected a simple set of hierarchical shape functions Mathematically speaking the basis functions are unions of contributing element shape functions and zero function elsewhere 0 0 5 1 0 0 5 1 0 0 5 1 0 0 5 1 0 05 0 03 0 02 0 01 0 05 0 0 0 5 1 0 0 5 1 Figure 1 1 D Hierarchical shape function defined as follows Me bof X 3 24 POSES Co e a a a The functions admit a simple recursive formula date pet R E 3 25 Xa X1 E X2 O Xi 1 X2 E X1
17. inating the first shape test function we rewrite the first row of the stiffness matrix and the load vector in such a form that would implicitly enforce condition u 0 yo The original element matrices Siu Sig ut Sin L3 DOR 3 45 But S52 zy Snn L Ly l 3 46 Ln get replaced with the modified matrices of the form 1 0 gt 0 D esc irea Bee a 3 47 p p 3 p 4 p Ki K2 K3 K4 Figure 3 Finite element mesh Yo pac 3 48 La Von Here n p 1 is the number of the element degrees of freedom Cauchy Neumann boundary condition at x l Addition of the boundary terms to the bilinear and linear forms results in the following modified element matrices 4 1 t Sin Io MN 3 49 Es dim vd L H 3 50 Ln 3 2 6 Global stiffness matrix and load vector the assembling procedure As global basis functions are constructed by gluing together element shape functions the additivity of integrals 4 implies that the entries of the global matrices are calculated by accumulating the corresponding contributions from element matrices This assembling pro cedure is pivotal in the Finite Element Method We shall illustrate it with a simple mesh consisting of four elements shown in Fig 3 The mesh consists of 4 elements denumerated from the left to the right K K2 X3 and Ky This ordering of elements induces the so called natural order of nodes We begin by listing all nodes of the first element then continue with those node
18. ions with respect to the physical coordinate axe dk 1 dex d h determine the weight w w h get load f f x get material constants a a x b b x c c xi for each d o f ky accumulate for the load vector entries Ly Ly fxiw for each d o f ko accumulate for the stiffness matrix entries dx d d Dk ko Ska sho T a At AUR Tn m z end of the second loop through the d o f end of the first loop through the d o f end of loop through integration points Xk C Xk Xi W 5 4 2 Modification of the element matrices due to boundary conditions In element number Nel element stiffness matrix Sj and load vector Lk Out Sj 4 Ly after BC modification get BC flags for the element vertex nodes for each vertex node i 1 2 CASE Dirichlet boundary get BC data y for each d o f k if k i then Ly y zero out the the row of stiffness matrix Ske 1 Ly Ly Sik y Sig 0 endif end of loop through d o f CASE Neumann boundary get BC data f y accumulate for the stiffness matrix and load vector Si Su B L Li y end of loop through vertex nodes 5 4 3 Assembling global matrices We first have to establish a global denumeration for all basis functions In principle one could follow the numbering of the nodes and then the numbering of the corresponding nodal shape functions However in general such a denumeration may not be optimal from
19. nding exact and approximate solutions The forth option automatic hp refinements will be discussed in the next section Please disregard the call to a testing routine which has been used for debugging 6 Adaptivity 6 1 p refinement unrefinement meshmods pack modord f Modifying order of approximation for an element is easy as it affects only its middle node The memory allocated for the middle node d o f has to be either expanded or shortened depending upon the new order of approximation If the order is increased then the new d o f are initiated with zeros 6 2 h refinement meshmods_pack break f Breaking an element involves creating two new entries in data structure array ELEMS for the element sons and creating one new vertex node and two new middle nodes in array NODES The sons have the same order of approximation as their father their d o f are initiated in routine meshmods_pack inidofh f in such a way that the new representation of the solution will exactly match the one corresponding to the single father element We do not delete the middle node of the father During the h refinement the information about the family tree is stored This includes storing the information on elements fathers and sons The concept is illustrated in Fig 7 6 3 Natural order of elements datastrs pack nelcon f The numbering of elements in the initial mesh from the left to the right and the family tree structure induce the correspondin
20. ng h or p we wish to know bounds on the error measured in the foregoing norms Ordinarily for a single element K these estimates will be of the form 5 llellk lt C 4 64 Pk where C is a constant depending on the data of the problem regularity of the solution hx stands for the element size length and px denotes its order of approximation The exponents r and s are the measure of the rate of convergence of the method with respect to the particular choice of the norm and either h or p refinements In our hp refinment strategy discussed in section 6 we shall use the fact that the h convergence rate is limited by two factors e order of approximation px and e local regularity of the solution expressed in some appropriate norm expressed in terms of derivatives of s 1 order More precisely r min p s 4 65 A convergence rate r lower than order of approximation p indicates a low regularity of the solution Estimates of this type require no information from the actual finite element solution They are known prior to the construction of the solution and called a priori estimates The a priori estimation of errors in numerical methods has long been an enterprise of numerical analysts Such estimates give information on the convergence and stability of various solvers and can give rough information on the asymptotic behavior of errors in calculations as mesh parameters are appropriately varied 4 2 A posteriori err
21. oblem consists of a second order ordinary differential equation a x u x b x u x e x u x f x x 0 1 2 1 accompanied at each of the endpoints x 0 or x l with one of three possible boundary conditions e Dirichlet boundary condition u 0 7 or ull 2 2 e Neumann boundary condition a 0 u 0 7 or a l w l 2 3 e Cauchy boundary condition a 0 u 0 Gou 0 y or alu l Gull y 2 4 The Neumann boundary condition is just a special case of Cauchy boundary condition with constant 0 For discussion of other possible boundary conditions including periodic boundary conditions see 4 The data for the problem consist thus of e coefficients of the differential operator a x b z c z the material constants e right hand side f x the load e boundary conditions data 3 y 2 2 Weak variational formulation of the problem The weak formulation is a reformulation of the strong form and it is from this weak form that the FE approach is established Whenever a smooth classical solution to a problem exists it is also the solution of the weak problem To establish the weak form of the strong form given by equation 2 1 multiply 2 1 by an arbitary so called test function v x and integrate over interval 0 1 We obtain a a u x b x u x e x u x v x dz is f x v x dz 2 5 Next we integrate the first term by parts l l l l au v dx
22. ombination into the variational boundary value problem 3 16 and setting the test functions to the basis functions v em l 1 2 Np we arrive at the algebraic system of equations Find upk k 1 2 N4 such that Nh 3 18 b uo gt gt unkenk eni Hen L 1 2 Nh us k 1 The approximate solution depends only on the space V and is independent of the basis functions eag In order to simplify the notation we shall drop now the approximate space mesh index h remembering that all quantities related to the approximate problem depend upon the index h Finally using the linearity of the bilinear form b u v in u we are led to the following system of linear equations Find uj k 1 N such that N 3 19 Suse a el k 1 Here S denotes the global stiffness matrix 1 deyid d Sy blek ei i a de be ceper dx 3 20 and Live stands for the modified load vector Ld I e b uo ei 2 fe dz yeild f a ade o 9 ei cuger dx The array l Li l ei f fe dx ye 3 22 is called the original load vector Notice the two different meanings of letter l 3 2 Finite Element Method The Finite Element Method is a special case of the Galerkin method and differs from other methods in the way the basis functions are constructed Domain 0 1 is partitioned into disjoint subdomains called finite elements Next for each element K we introduce the cor responding shape functions Xg w
23. or estimation In some cases a more detailed estimate of accuracy can be based on information obtained from the finite element solution itself It is called an a posteriori error estimate and it can be calculated only after the finite element solution has been obtained Interest in a posteriori estimation for finite element methods for elliptic boundary value problems began with the pioneering work of Babuska and Rheinboldt 3 In this note we will present the element implicit residual method introduced in 5 7 1 and applied to a variety of problems in mechanics and physics 4 2 1 The element implicit residual method We introduce the idea of the error estimation using the standard abstract variational formu lation u E uo V 4 l b u v 2 i v Vv V oe Introducing a finite element space V C V we calculate the corresponding finite element solution u by solving the approximate problem Un ug Va 4 l O E e7 The goal is to estimate the residual defined as b up v Uv lIrllvs sup RAUS CM 4 68 vcV Iv where denotes the energy norm lulli a u u 4 69 Here a u v is identified as a symmetric part of b u v that defines a norm For the particular example discussed in the previous section se may select a u v ni alu dx Byu I 4 70 We have to assume then that a x gt 0 and A gt 0 If additionally coefficient c r gt 0 we can include it in the definition of the energy n
24. orm as well a u v A a u cu dx pult 4 71 REMARK 1 Mathematically speaking the residual is a linear functional acting on space V and it must be measured using the dual norm Its choice depends upon the choice of the norm used in the denominator It can be shown 8 that for a class of self adjoint problems where b u v a u v the residual is equal to the error measured in the energy norm We shall represent the residual in the form r us v b us v a Pay Un v lx v Ax v 4 72 ore tp v K where bx u v and lg v are contributions to the bilinear and linear forms from element K respectively and Ax v is an element flux functional We will postulate the following two main assumptions e the element residuals ry uj v are in equilibrium with respect to the finite element space rK Un v bglun v Ik v Ak v 20 Vo e VA K 4 73 e Consistency condition MAx v 20 VveV 4 74 K Here V K denotes the space of element K shape functions possibly incorporating Dirichlet boundary conditions if element K is adjacent to the Dirichlet boundary Next we introduce the local element Neumann problems l Find x V K 4 75 ak x v Ys y Vi V K Here V K H K except for the element adjacent to Dirichlet boundary at x 0 for which V K v H K v 0 201 4 76 We can express now the mesh residual in terms of the element error indicator function
25. phase can begin Solvel f prepares arrays in and iawork calls the prefront routines and then calls the main elimination routines Thus this routine is seen to be the primary driver of the frontal solver The other interface routines are simply for auxilliary purposes For a given element solinl f returns a listing of the destination vectors of the associated modified element solin2 f returns the modified element stiffness matrix and load vector and solout f takes the solution values returned from the frontal solver and inserts them into the data structure for a given node the values of the corresponding dof must be placed into the NODES nod 96dof entry
26. ponding hp interpolant amp 5 defined on master element pci ttnp Do ujX 3 36 J The final hp interpolant over element K is defined as the composition of the master element interpolant amp 5 with the inverse of the element map x pti Uhpl T 2 uix Gn 3 37 Practically that means only that the coefficients u must be determined by solving the appropriate system of equations on the master element Figure 2 Construction of the vertex nodes basis functions 3 2 3 1D hp finite element space Let now interval 0 7 be covered with a FE mesh consisting of disjoint elements K With each element K we associate a possibly different order of approximation p px and element length h hg The element endpoints with coordinates 0 o lt z 4 lt zy lt yy l will be called the verter nodes We define the 1D hp finite element space Xj as the collection of all functions that are globally continuous and whose restrictions to element K live in the element space of shape functions Xp up x u is continuous and u y X K for every element K 3 38 The global basis functions are classified into two groups e the vertex nodes basis functions and e the bubble basis functions The basis function corresponding to a vertex node x is defined as the union of the two adjacent element shape functions corresponding to the common vertex and zero elsewhere The construction is illustrated in
27. rices Apply the boundary conditions Assemble global matrices Solve the system of equations POSTPROCESSOR Calculate quantities of interest Print out and plot results Figure 6 Computer flow diagram of a classical finite element code the solution are computed and printed The results may also be output in a graphical form e g using contour plots 4 Error estimation 4 1 A priori error estimation The error introduced into the finite element solution up because of the approximation of the dependent variable u in an element is inherent to any problem Nh UND Up 5 UhkC hk 4 57 k 1 Here up is the finite element solution and u is the exact solution over the domain We wish to know how the error e u up measured in a meaningful way behaves as the number of elements and or their order of approximation in the mesh is increased There are several ways in which one can measure the difference between any two functions u and up More generally used measures of the difference of the two functions are energy norm L norm and mazimum norm llelg if Jale ce el 4 58 1 2 llellz2 u zl 4 59 0 lello max ez 4 60 The norms listed above are global i e thay apply to th whole domain 0 Analogous definitions hold for any element K EK V a e ce 2 4 61 llel ju lenem TE a 4 62 les max e z 4 63 As we refine the mesh by increasi
28. s eliminated from the system of equations using standard Gaussian operations Thus in the frontal solver approach the operations of assembling and elimination occur simultaneousely The global stiffness matrix never needs to be fully assembled and this leads to the significant savings in memory that has given the frontal solver its popularity Here we will only describe the interface with the frontal solver not the solver itself The interface is constructed via four routines all located in the solverl pack directory solvel f solinl f solin2 f and solout f We will now give an overview of these routines For coding details we refer to the source codes in solverl pack The frontal solution consists of two steps prefront and elimination The prefront re quires two arrays on input in and iawork For each element in contains the number of nodes associated with the corresponding modified element and iawork contains a listing of nicknames for the nodes of the modified element The nicknames are defined as follows for a given node j nick j x 1000 ndof A 80 where ndof is the number of degrees of freedom associated with the node i e is equal 1 for a vertex node and p 1 for a middle node of order p With this information the prefront produces the destination vectors which for a given element denote at what stage of the frontal solution each of its nodes can be eliminated Once this information is constructed the elimination
29. s of the second element that have not been listed yet and so on The order for al a3 a2 a5 a4 al a6 a9 a8 gt 99 P 99 9 9 Ki K2 K3 K4 Figure 4 Natural order of nodes 1 2 2 5 6 4 8 9 10 fi 11 gt S S 79 9 K K2 K3 K4 Figure 5 Natural order of global d o f nodes is depicted in Fig 4 Notice that as the right vertex node of the element is listed before its middle node the natural order of nodes is not equivalent to denumerating simply nodes from the left to the right Finally the natural order of nodes implies the correspond ing natural order of global degrees of freedom obtained by listing the nodal d o f node after node The ordering is depicted in Fig 5 Notice the variable number of degrees of freedom associated with the middle nodes and in particular no degrees of freedom associated with the last middle node at all We are now ready to discuss the assembling procedure We proceed with one element at a time The first element connectivites indicating the global d o f numbers assigned to the local ones are 1 2 3 Consequently after assembling the first element matrices the global matrices will look as follows Si Siz Sis Sh S5 S5 S5 S5 S5 3 51 3 52 Dotted entries indicate zeros Similarly the second element connectivites are 2 4 5 and 6 After assembling the second element matrices we obtain Su S Ss Sh ShctS5h S S5 Sis Sty Sai S3o ie i9 a es Sh Sh S
30. sat isfies the nonhomogeneous Dirichlet data we can make the substitution u ug w where lup to regularity assumptions on the solution w V satisfies the homogeneous Dirichlet boundary conditions and try to determine the perturbation w The corresponding abstract formulation is then as follows Find w V such that 2 15 b w v Uv b ug v Vv EV For the sake of simplicity of presentation we shall use the particular choice of the bound ary conditions throughout the rest of these notes All other cases can be treated in a completely analogous way 3 Finite Element Method 3 1 Galerkin approximation With these preliminaries behind us we are ready to consider Galerkin s method for construct ing approximate solutions to the variational boundary value problem Galerkin s method consists of seeking an approximate solution to variational boundary value problem in a finite dimensional subspace V of space V This procedure leads to the following approximate variational boundary value problem Find up uo Vp such that 3 16 blun Un Lur Vu Vi We introduce a finite set of basis functions e51 55 enn in V that span a finite dimen sional subspace of test functions V in V We then seek the approximate function up uo Vj in the form Ny up x Ugo 5 UhkEhk 3 17 k 1 The unknown coefficients upg k 1 2 Nj are called global degrees of freedom Sub stituting the linear c
31. tempt to choose between h and p refinements Instead after the problem has been solved on the h refined mesh we try to trade the h refined elements for p refinements Towards this goal we first solve the problem on the h refined mesh Next we loop through all just refined elements and compute for each of them the corresponding local numerical h convergence rate of the residual If the rate is optimal i e it is equal to or it exceeds the corresponding order p of approximation for the element we unrefine cluster the element and trade the h refinement for a p refinement i e we increase p to p 1 A rate of convergence below the order of approximation p indicates a lower regularity compare estimate 4 64 and in such a case we leave the h refined element unchanged The formal algorithm looks as follows Note the tollerance factor 9 for each just h refined element if element order px lt 7 then estimate the new element residual by summing up the residuals for its sons compute the numerical convergence rate rg if rg gt 0 9 x py then cluster back the element increase element order from px to px 1 endif endif end of loop through refined elements Consult the meshmods_pack trade f routine for details 6 7 Interactive refinements In many problems we may use our experience on the problem at hand to begin with a better than uniform initial mesh produced by the mesh generator The graph_1D graph1D f routine that displays a gr
32. the error Consequently we shall assume that the data to the problem load f x and boundary condition data y will always be calculated using the exact solution by calling routine data pack ezact f That way we can minimize the number of changes in the code when we want to study a different solution The material data operator coefficients a x x c x and Cauchy boundary data 8 have to be set independently in routines data pack getmat f and getc f Files All files are placed in directory files Besides the input file containing the data for the initial mesh generation discussed earlier the code opens an output file used e g by routine datstrs pack result and an additional file result to be discussed in the next section Both files output and result are automatically opened by the code and need no user s action Running the code Step 1 Define a boundary value problem and select an exact solution that you want to reproduce with the FE code Modify routines data pack getmat f data pack getc f and data_pack exact f accordingly Step 2 Prepare the input file Step 3 Use the provided makefile to compile and link the code Step 4 Type a out to execute the code If the input file is correct and the initial mesh has been generated successfully the code will display the main menu that includes the possibility of solving the problem printing out the content of the data structure arrays and displaying the mesh with the correspo
33. ties has been established the assembling procedure follows the standard algo rithm for each element nel in the mesh calculate element local matrices Aloc Bloc for each nodal d o f i establish element d o f connectivities loc_con nel i nbij nod i 1 end of loop through nodal d o f for each element d o f i determine the connectivity k loc con nel i accumulate for the global load vector Bglob k Bglob k Bloc i for each element d o f j determine the connectivity loc_con nel i accumulate for the global stiffness matrix Aglob k l Aglob k 1 Aloc i 7 end of the second loop through element d o f end of the first loop through element d o f end of loop through elements 5 4 4 Solver Routines assembling the global matrices and solving the corresponding system of equations are to be provided by the user For the sake of presenting an operational code we include in the code an interface with a more complicated frontal solver discussed shortly in the Appendix 5 4 5 Setting up data Files As the main goal of this 1D code is rather academic and focuses more on studying the finite element method than solving practical problems we shall accept a rather unusual way of inputing data Namely we shall assume that we do know the exact solution together with its first and second derivatives The purpose of our finite element computations will be just to compute the finite element approximation and study
34. ule TAS EE dx _ Xn de RM da d dx After switching to the master element coordinate we obtain 1 dXi dE dX dE dXi n Ls S dT K cU de d f cj ORE aa that exa d 3 42 D f 0 As before symbol indicates the composition with the element map e g a z x 3 43 4The support of a function is defined as the closure of a set over which the function takes on values different from zero The evaluation of the integrals in 3 42 is performed using numerical integration In most finite element calculations Gauss quadrature rules are used We note that the Gauss rule of order N will integrate exactly polynomials of degree 2N 1 The explicit formulas for evaluation 3 42 using Gauss quadrature are K LN d s der d Ta d r EE TE sE ote BORO sa e ea E oaeee Se Oa E 3 44 where x zg is the value of element map x calculated at integration point j Notice that for an affine element hy and amp h are independent of integration point j dz K T 3 2 5 Taking into account the boundary conditions Modified element matrices In the case of the first and the last element element stiffness matrix and load vector must be modified to incorporate changes due to the boundary conditions Dirichlet boundary condition at x 0 We use the first element linear shape function 1 premultiplied by yo to construct the lift of the boundary conditions data Instead of elim
35. und in directory files Printing out content of data structure arrays datastrs pack result f The routine prints out the current content of the data structure arrays inlcuding the complete information on elements and nodes It can be conveniently used for debugging the code Graphical output graph_1D graph1D f Th routine displays a graphical representation of the current mesh and plots the corresponding exact and numerical solutions The scale on the right prescribes the color code for different orders of approximation p 1 8 5 4 Processing algorithms We discuss quickly a number of algorithms pivotal in the implementation of any Finite Element Method the calculation of element matrices modification due to the boundary conditions and the assembling procedure Please consult the corresponding routines for details 5 4 1 Evaluation of element matrices elem_pack elem f In element number Nel Out element stiffness matrix Sp kp and load vector Lk determine the element order of approximation and select the corresponding Gauss quadrature determine the element vertex nodes r r geometry d o h xg zr initiate element stiffness matrix 55 4 and load vector Lk for each integration point j evaluate the physical coordinate x zr jh evaluate master element shape functions Y and their derivatives with respect to the master element coordinate dg evaluate the derivatives of the shape funct
36. w Problems with Moving Boundaries Part 1 Variational Principles and A Posteriori Error Estimates Comput Meth Appl Mech Engrg 46 217 251 1984 L Demkowicz J T Oden W Rachowicz and O Hardy Toward a Universal hp Adaptive Finite Element Strategy Part 1 Constrained Approximation and Data Structure Com put Meth Appl Mech Engrg 77 79 112 1989 J T Oden L Demkowicz T Strouboulis and P Devloo Adaptive Method for Problems in Solid and Fluid Mechanics Accuracy Estimates and Adaptive Refinements in Finite Element Computations John Wiley Sons Ltd 249 280 1986 J T Oden L Demkowicz W Rachowicz and T A Westermann Toward a Universal Ap Adaptive Finite Element Strategy Part 2 A Posteriori Error Estimation Comput Meth Appl Mech Engrg 77 113 180 1989 W Rachowicz J T Oden and L Demkowicz Toward a Universal hp Adaptive Finite Element Strategy Part 3 Design of hp Meshes Comput Meth Appl Mech Engrg 77 181 212 1989 A Interface with a frontal solver The frontal solver is a popular choice among direct solvers for finite element codes due to its natural implementation in an element by element scheme In this method a front sweeps through the mesh one element at a time assembling the element stiffness matrices into a global matrix The distinction from the standard assembling procedure is that as soon as all of the contributions for a given dof have been accumulated that dof i
37. y Neumann boundary condition at x l Notice that there is no need to recover Dirichlet boundary condition as it has been simply rewritten into the weak formulation A precise mathematical setting is obtained by introducing the Sobolev space of the first order H 0 1 consisting of functions that are together with their derivatives square inte grable l l H 0 1 fv z I v dz I v dz lt oo 2 9 0 0 Next we identify the subspace of kinematically admissible functions V satisfying the homo geneous Dirichlet boundary condition V v H 0 l v 0 0 2 10 The set of functions satisfying the nonhomogeneous Dirichlet boundary condition has a more complicated algebraic structure of an affine subspace and can be identified as the collection of all sums uo v where ug H 0 1 is a lift of the boundary data and v is an arbitrary function from H 0 I satisfying the homogeneous Dirichlet boundary conditions u H 0 1 u 0 2 49 uo V ug v v V 2 11 The variational formulation can now be written in the form of the abstract variational boundary value problem Find u ug V such that 2 12 b u v lw Vv eV Here l b u v J au v bu v cuv dx Byu I v L 2 13 0 is a bilinear form of arguments u and v and l l v i f x u dx yw l 2 14 is a linear form of test function v Equivalently speaking once we have found a particular function ug H 0 that
Download Pdf Manuals
Related Search