AFEM USER GUIDE Contents 1. Introduction 1 2. Installing Scilab 2
Contents
1. maxlevels 12 and maxnodes 1000 in the list Notice that the pair dimensions 12 1000 3000 8000 would correspond to the same number of levels and the same maxi mal number of nodes but would allow for the bigger maximal number 3000 of elements on a single mesh and the bigger maximal cumulative number 8000 of elements on all meshes 5 3 Predefined sample equations and discretizations The key word sample_pde allows you to use predefined sample equations and discretizations The corresponding value must be list The length of this list depends on the example as explained below Its first argument always is a string which specifies the example If the list s length is 1 you may only enter this string as value Thus the pairs sample_pde list lisq and sample_pde 1lisq in your parameter list both have the same effect Notice that the use of the keyword sample_pde together with a corresponding value over rides the settings of the keywords diffusion convection reaction force traction robin bdrydispl and solution described in sub section 5 4 If the chosen example also specifies a domain and initial mesh the settings of the keywords sample_mesh init_node init_elem and domain described in subsection 5 5 are also overridden sample_pde allows you to choose among the following examples 10u1sq Laplace equation Au 0 with Dirichlet boundary con dition and exact solution 1 on the square 1 1 subdivided int2. 2 Print on each level all informations on the nodes and elements The default setting for prt1 is 0 The value of p1t1 must be a vector of length seven the default setting being 1 1 1 1 20 10 10 The first component decides whether the meshes are plotted or not with 0 switching plot off and 1 switching plot on The second component determines the plot of the discrete solution as follows 0 no plot 1 contour plot see figure 5 4 2 surface plot with Scilab s built in function plot3d1 see left part of figure 5 5 3 surface plot with Scilab s built in function plot3d see right part of figure 5 5 4 contour plot plus surface plot with Scilab s built in function plot3d1 5 surface plots with Scilab s built in functions plot3d1 and plot3d 6 contour plot plus surface plot with Scilab s built in function plot3d 7 contour plot plus surface plots with Scilab s built in functions plot3d1 and plot3d The third and fourth component of the vector have the same effect for the error indicator and the error respectively if existent The fifth sixth and seventh components of the vector determine the number of 14 R VERFURTH level lines for contour plots of the discrete solution error indicator and error respectively solution of level 4 with 64 elements and 41 nodes 1000 r T T T T T T 1 000 0778 0556 0333 0 111 0 111 0333 0556 0 78 1000 FIGURE 5 4 Contour plot of a functi
3. on the circle x 23 lt 1 subdivided into 8 triangles see left part of figure 5 2 and given by second argument lila Laplace equation Au f with Dirichlet boundary condi tion and exact solution tanh x z on the annulus ig lt x 23 lt 1 subdivided into 16 triangles see figure 5 3 and given by second argument lilrsq Laplace equation Au f with Dirichlet boundary con dition and exact solution Dee E mea on the square 1 1 subdivided into 8 isosceles right angled triangles see left part of figure 5 1 and given by second argument lilrc Laplace equation Au f with Dirichlet boundary con dition and exact solution res Da on the circle 27 z lt 1 subdivided into 8 triangles see left part of figure 5 2 and K given by second argument lilra Laplace equation Au f with Dirichlet boundary con dition and exact solution on the annulus lt Tercera aor i6 x 23 lt 1 subdivided into 16 triangles see figure 5 3 and given by second argument 1ssl Laplace equation Au 0 with Dirichlet boundary con dition and exact solution r3 sin 3y on the L shaped domain 1 1 0 1 x 1 0 subdivided into 6 isosceles right angled triangles see right part of figure 5 1 1lsscs Laplace equation Au 0 with Dirichlet boundary condi I I 2 s tion and exact solution r3sin 3y on the circular sector R VERFURTH rcosy rsing 0 lt r lt 1 0 lt p lt
4. 3m subdivided into 6 triangles see right part of figure 5 2 druisq Diffusion reaction equation Au u 1 with Dirichlet boundary condition and exact solution 1 on the square 1 1 subdivided into 8 isosceles right angled triangles see left part of figure 5 1 and given by second argument e g sample_pde list druisq 100 sets k 100 drbruisq Diffusion reaction equation Au 1 27 1 z2 u 1 27 1 z2 with Dirichlet boundary condition and exact solution 1 on the square 1 1 subdivided into 8 isosceles right angled triangles see left part of figure 5 1 and K given by second argument e g sample _pde list drbruisq 100 sets x 100 drubsq Diffusion reaction equation Au u f with Dirichlet boundary condition and exact solution 1 x 1 73 on the square 1 1 subdivided into 8 isosceles right angled triangles see left part of figure 5 1 and given by second argument e g sample_pde list drubsq 100 sets x 100 drbrubsq Diffusion reaction equation Au K 1 z 1 z2 u f with Dirichlet boundary condition and exact solution 1 x 1 z2 on the square 1 1 subdivided into 8 isosceles right angled triangles see left part of figure 5 1 and given by second argument e g sample _pde list drbrubsq 100 sets amp 100 drilsq Diffusion reaction equation Au Ku f with Dirichlet boundary condition
5. built in sparse linear system solver of Scilab The mesh refinement is either uniform or adaptive with marked edge bisection based on either a standard residual error indicator or a residual error indicator based on edge residuals exclusively or the Zinkiewicz Zhu error estimator in its simplest form We refer to 4 5 for a more detailed description of the mathematical background A detailed description of AFEM and its main function afem are given in sections 4 and 5 below In sections 2 and 3 we first explain how to install Scilab and AFEM respectively 2 INSTALLING SCILAB Scilab is a numerical programming environment comparable to Matlab with a similar syntax It was developed by INRIA France and is currently maintained by the Digiteo foundation in collaboration with INRIA Contrary to Matlab Scilab is completely free of charge There are versions for the operating systems Mac OsX Unix and Windows To use AFEM you first have to get the appropriate version of Scilab from http www scilab org and install it on your computer following the instructions of that web site Useful introductions to Scilab are e g 1 2 3 may be of par ticular interest for those familiar with Matlab 3 INSTALLING AFEM To install AFEM click the link AFEM zip on http www rub de num1 softwareE html or http www rub de num1 software html This installs a zip archive AFEM zip on your computer Unpack the zip archive and put the created
6. see left part of figure 5 1 and given by second ar gument e g sample_pde list dciblsq 0 001 sets 0 001 aidsq Anisotropic diffusion equation 2u _ Pu r r with Dirichlet boundary condition 1 ifa 1 up z lt 1 if z 1 0 else on the square 1 1 subdivided into 8 isosceles right angled triangles see left part of figure 5 1 and given by second argument e g sample_pde list aidsq 0 001 sets 0 001 cbdsq Checker board diffusion equation dAu 0 with d f if lt 0 1 else and Dirichlet boundary condition 0 if z z lt 0 up z lt 1 if z4z gt 0 and z gt 0 1 ifxziz gt 0 and z lt 0 on the square 1 1 subdivided into 8 isosceles right angled triangles see left part of figure 5 1 and given by second 10 R VERFURTH argument e g sample _pde list cbdsq 0 001 sets 0 001 mesh of level 1 with 8 elements and 9 nodes 10 T T T T T r T T 10 08 06 04 02 00 02 04 06 08 10 mesh of level 1 with 6 elements and 8 nodes FIGURE 5 1 Initial mesh for a square left and an L shaped domain right mesh of level 1 with 8 elements and 9 nodes 7 lt 08 Fi in 044 a A if oz N 00 ug N ul A 064 084 i we 10 T T T T T r T T 1 10 08 06 04 02 00 02 04 08 08 10 mesh of level 1 with 6 elements and 8 nodes o T T T T T T r T T 1 10 08 06 04
7. 1 in the console Note that entering afem init_param reftype 1 force myforce will produce the same result while entering afem init_param force 1 reftype myforce will result in an error message since the types of the values do not fit with the keywords 1 is not a function myforce is not a number 5 2 Dimensions The behaviour of afem and its storage requirement is controlled by two keywords maxnodes determines the maximal number of nodes i e vertices of elements Its default value is 10000 The maximal number of elements in a single mesh is automatically set to 2maxnodes the maximal cumulative number of elements over all meshes is automatically set to 6maxnodes maxlevels determines the maximal number of refinement levels Its default value is 21 Instead of prescribing the keywords maxnodes and maxlevels with corresponding values you may also control the behaviour of afem and its storage requirement by entering the keyword dimensions followed 6 R VERFURTH by an integer vector of length four as value The four components of this vector have the following meaning 1 maximal number of levels default is 21 2 maximal number of nodes default is 10000 3 maximal number of elements on a single level default is 20000 4 maximal cumulative number of elements on all levels default is 60000 Thus the pair dimensions 12 1000 2000 6000 in your parameter list has the same effect as the two pairs
8. boolean with the value t indicating a do main with curved boundary The second argument must be a function which is negative inside the domain and positive outside The third ar gument finally must be the gradient of the function given in the second argument The default setting for domain corresponds to the square 1 1 The keyword sample_mesh allows a short hand prescription of the do main and initial mesh Its value must be one of the following strings sq8 criss cross grid for the square 1 1 consisting of 8 isos celes right angled triangles see left part of figure 5 1 16 criss cross grid for the L shaped domain 1 1 0 1 x 1 0 consisting of 6 isosceles right angled triangles see right part of figure 5 1 c8 criss cross grid for the unit circle z x2 lt 1 consisting of 8 triangles see left part of figure 5 2 cs6 criss cross grid for the circular sector rcosy rsiny 0 lt r lt 1 0 lt lt T consisting of 6 triangles see right part of figure 5 2 a16 grid for an annulus with radii 1 and 1 consisting of 16 triangles see figure 5 3 5 6 Error estimation and mesh refinement The following key words which all require a number as a value control the error estimation and adaptive mesh refinement of afem relerrtol Sets the relative error tolerance The adaptive pro cess is terminated once the ratio of the estimated error over the energy norm of the discrete solution is less tha
9. f of the differential equation Then you first enter function y myforce x y 1 2xx 1 x 1 5 x 2 x 2 endfunction in the Scilab console and then use myforce as the corresponding ar gument of afem 4 2 Start Before using AFEM you must make it known to Scilab by first entering mylib lib path to the console Here mylib may be any other name different from the name of any AFEM function any built in function and any of your user defined functions path is a string giving the full path to the AFEM library Suppose for example that you are using Mac OsX that Gargantua is your user name and that you have put the AFEM folder into the subfolder Pantagruel of your documents folder Then the above command must read without the full stop at the end mylib lib users Gargantua documents Pantagruel AFEM If everything is correctly installed Scilab returns the names of all functions in the AFEM library This includes auxiliary functions which are not user relevant and which will not be described below 4 3 Help AFEM comes with a help function afem_help which has an empty argument list and which yields a brief information about the arguments of the main function afem and their use 5 FUNCTION afem 5 1 Overview AFEM s main function is afem It has several argu ments which allow to adjust the parameters of the differential equa tion of its discretization and of the adaptive process and which will be described in mo
10. 02 00 02 04 O08 08 10 FIGURE 5 2 Initial mesh for a circle left and a circular sector right 5 4 Coefficients and right hand side of the differential equa tion The following keywords with functions as corresponding values allow to prescribe the coefficients and right hand side of the differential equation diffusion diffusion A default is A 4 9 convection convection a default is a 0 reaction reaction a default is a 0 force right hand side f default is f 1 traction right hand side g default is g 0 robin function p in Robin condition default is p 0 bdrydispl boundary displacement up default is up 0 AFEM USER GUIDE 11 FIGURE 5 3 Initial mesh for an annulus with radii r and 1 Notice that adopting all default values yields the Laplace equation Au 1 with homogeneous Dirichlet boundary condition Also recall that prescribing the keyword sample_pde described in subsection 5 3 overrides the settings of the above keywords The keyword solution allows to prescribe an eventually known ex act solution of the differential equation Its value must be a vector of length two with its first argument either Zt for a known exact solu tion or Zf for an unknown exact solution and its second argument a built in or user defined function giving the exact solution The default setting is f unit function where unit function corresponds to the function u 1 which is among AFEM s function
11. AFEM USER GUIDE R VERFURTH CONTENTS 1 Introduction 2 Installing Scilab 3 Installing AFEM 4 Using AFEM 4 1 Data structures 4 2 Start 4 3 Help 5 Function afem 5 1 Overview 5 2 Dimensions 5 3 Predefined sample equations and discretizations 5 4 Coefficients and right hand side of the differential equation 10 ap RR G GO ND ND aD 5 5 Initial mesh and domain 11 5 6 Error estimation and mesh refinement 12 5 7 Print and plot output 13 References 14 Index 15 1 INTRODUCTION AFEM is a library of Scilab functions implementing an adaptive linear finite element method for general linear elliptic equations of second order in a bounded connected two dimensional domain with mixed Dirichlet Neumann or Robin boundary conditions div AVu a Vutoau f mQ U Up only n AVu g only pu n AVu g9 ontg The discretization consists of standard linear triangular elements with out any additional stabilization and up replaced by its linear interpo late Irup Date January 5 2012 2 R VERFURTH Find ur Irup S T such that i Vur AVur a Vurur auror Q f purer for f gor ive Q Py holds for all vy S3 T Notice that the domain Q may have a curved boundary In this case the integrals in the discrete problem extend over piecewise linear approxi mations of Q and its boundary The assembly is based on quadrature rules of order one and two The discrete problem is solved using the
12. and exact solution tanh z z2 on the square 1 1 subdivided into 8 isosceles right angled triangles see left part of figure 5 1 and amp given by second argument e g sample_pde list drilsq 100 sets k 100 drilc Diffusion reaction equation Au Ku f with Dirichlet boundary condition and exact solution tanh k ai 23 1 on the circle 2 x3 lt 1 subdivided into 8 triangles see left part of figure 5 2 and given by second argument e g sample_pde list drilc 100 sets amp 100 drila Diffusion reaction equation Au Ku f with Dirichlet boundary condition and exact solution tanh z 23 on the annulus lt 2 z lt 1 subdivided into 16 tri angles see figure 5 3 and x given by second argument e g sample _pde list drila 100 sets amp 100 AFEM USER GUIDE 9 drilch Diffusion reaction equation Au u f with Dirichlet boundary condition and exact solution tanh x x3 tanh K3 tanh K tanh 438 on the circle z x2 lt 1 subdivided into 8 triangles see left part of figure 5 2 with given by second argument e g sample_pde list drilch 100 sets x 100 dciblsq Anisotropic diffusion convection equation u u u u 02 Ax TA ri 7 Or i with Dirichlet boundary condition 1 ifz 1 up z lt 1 ifz 1 0 else on the square 1 1 subdivided into 8 isosceles right angled triangles
13. folder AFEM at any place on your computer that suits you Now you can start AFEM as described in section 4 3 AFEM USER GUIDE 3 4 USING AFEM 4 1 Data structures We refer to 1 2 3 for a detailed introduction to Scilab Here we briefly present those data structures of Scilab that are used for input arguments of AFEM s main function afem Boolean Boolean arguments are transferred to afem by either directly entering Zt for the value true or Zf for the value false or by cre ating an own boolean variable myboolean by typing myboolean t or myboolean f in the Scilab console and then calling afem with myboolean as argument String Strings are enclosed by quotation marks They are transferred to afem by either directly entering name for the string value name or by creating an own string variable mystring by typing mystring name in the Scilab console and then calling afem with mystring as argument Number Numbers may be integers like 2 and 3 or fixed point re als like 3 14 and 2 87 or floating point reals like 0 578D4 0 689D5 0 203D 7 and 0 1763D 5 They are transferred to afem by either directly entering the numerical value or by creating an own variable mynumber by typing mynumber 2 or similar in the Scilab console and then calling afem with mynumber as argument Vector Matriz A vector is either an m x 1 matrix for a column vector or a 1 x n matrix for a row vector Matrices are created by typing the app
14. n relerrtol The default value of relerrtol is 0 01 reftype Determines the refinement type as follows 0 uniform refinement 1 adaptive refinement based on the maximum strategy and edge residuals as error indicator 2 adaptive refinement based on the maximum strategy and edge and element residuals as error indicator 3 adaptive refinement based on the maximum strategy and the Zinkiewicz Zhu indicator The default value for reftype is 0 AFEM USER GUIDE 13 marktol Sets the threshold for the maximum strategy An el ement is marked for refinement if its estimated error is larger than marktol times the maximal estimated error The default value for marktol is 0 5 firstadapt Start the adaptive refinement on this level parti tions with a lesser level number are refined uniformly The default value for firstadapt is 1 5 7 Print and plot output afem has no return values it only prints the required time for the computation In addition it allows to print some information on the partitions and to plot the meshes and discrete solutions and if available error indicators and errors These options are controlled by the keywords prt1 and p1t1 The value corresponding to prt1 is a number which has the following meaning 0 Print on each level the number of nodes degrees of freedom and elements 1 Additionally print the co ordinates of the nodes the global numbers of the nodes and the neighbourhood information of the elements
15. o 8 isosceles right angled triangles see left part of figure 5 1 lubsq Laplace equation Au f with Dirichlet boundary con dition and exact solution 1 z2 1 z2 on the square 1 1 subdivided into 8 isosceles right angled triangles see left part of figure 5 1 11 Laplace equation Au 1 with homogeneous Dirichlet boundary condition the domain and initial mesh must be pre scribed separately AFEM USER GUIDE 7 11sq Laplace equation Au 1 with homogeneous Dirichlet boundary condition on the square 1 1 subdivided into 8 isosceles right angled triangles see left part of figure 5 1 lic Laplace equation Au 1 with homogeneous Dirichlet boundary condition on the unit circle zz x lt 1 subdivided into 8 triangles see left part of figure 5 2 lla Laplace equation Au 1 with homogeneous Dirichlet boundary condition on the annulus lt x 23 lt 1 subdi vided into 16 triangles see figure 5 3 lg Laplace equation Au e7 zi z3 3 with K given by sec ond argument the domain and initial mesh must be prescribed separately lilsq Laplace equation Au f with Dirichlet boundary con dition and exact solution tanh z1 x3 on the square 1 1 subdivided into 8 isosceles right angled triangles see left part of figure 5 1 and given by second argument lilc Laplace equation Au f with Dirichlet boundary con dition and exact solution tanh x 73 4
16. on solution of level 4 with 64 elements and 41 nodes solution of level with 64 elements and 41 nodes N fe oe ra I I S lt FIGURE 5 5 Surface plot of a function with plot3d1 left and plot3d right REFERENCES 1 Micha l Baudin Introduction to Scilab http wiki scilab org Tutorials Programming in Scilab http wiki scilab org Tutorials 3 Eike Rietsch An Introduction to Scilab from a Matlab User s Point of View http wiki scilab org Tutorials 4 R diger Verfiirth Numerik II Finite Elemente http www rub de num1 files lectures NumDgl2 pdf Adaptive Finite Element Methods http www rub de num1 files lectures AdaptiveFEM pdf INDEX sample_pde 6 afem 4 afem_help 4 aidsq 9 a16 12 bdrydispl 10 boolean 3 cbdsq 9 c8 12 convection 10 cs6 12 dciblsgq 9 diffusion 10 dimensions 5 domain 11 drbrubsgq 8 drbruisq 8 drila 8 drilc 8 drilch 9 drilsq 8 drubsq 8 druisq 8 firstadapt 13 force 10 function 4 init_elem 11 init node 11 11 6 lia 7 lic 7 lisq 7 lg 7 lila 7 lilc 7 lilra 7 lilrc 7 lilrsq 7 lilsq 7 list 3 16 12 10u1isq 6 lsscs 7 lss1 7 lubsgq 6 15 marktol 13 matrix 3 maxlevels 5 maxnodes 5 number 3 plt1 13 prtl 13 reaction 10 reftype 12 relerrtol 12 robin 10 sample_mesh 12 solution 11 sq8 12 string 3 traction 10 vector 3
17. re detail in the following subsections Each argument has a default value so that afem may be called with an empty argument list When entering afem in the Scilab console the Poisson equation Au 1 inQ u 0 on AFEM USER GUIDE 5 is solved on the square Q 1 1 with uniform mesh refinement starting from a criss cross grid with eight triangles see left part of figure 5 1 Since this example is not very exciting you probably want to adjust some arguments This is achieved by entering afem init_param keyl valuel keyn valuen in the Scilab console Here key1 keyn is any number of any of the keywords described in the following subsections in any order and valuel valuen are the corresponding values The latter can be booleans strings numbers matrices lists or functions depending on the corresponding keyword Notice that keywords are strings and must be enclosed by quotation marks To clarify this consider two simple examples Suppose that you want to replace the constant right hand side 1 in the above example by the function 1 2x 523 and that you have defined the corresponding function myforce as described at the end of section 4 1 Then you must enter afem init_param force myforce in the console Suppose that you now want to try an adaptive refinement based on the maximum strategy and the error estimator based on edge residuals Then you must enter afem init_param force myforce reftype
18. ropriate sequence of numbers in lexicographic order from top left to bottom right in brackets with columns separated by blanks or commas and rows separated by semicolons Thus the matrices 1 2 G _ and 3 4 5 6 are represented by H 23456 and 1 2 3 4 5 6 A is the transpose of the matrix A A x B is the product of the ma trices A and B ones m n and zeros m n create m x n matrices with all elements equal to 1 and 0 respectively eye n n creates the n x n identity matrix Matrices are transferred to afem by either di rectly entering the matrix or by creating an own matrix valued variable mymatrix and then calling afem with mymatrix as argument List A list is an ordered collection of objects of different types It may be passed directly to afem by entering list obj1 objn for the corresponding argument Here obj1 objn is a sequence of arbi trary length of objects of arbitrary types such as numbers matrices functions lists etc Alternatively you my first create an own variable of type list by entering mylist list obj1 objn 4 R VERFURTH in the Scilab console and then call afem with mylist as argument Function Functions user defined or built in ones are transferred to afem by entering their name for the corresponding argument For a user defined function however you first have to create a corresponding variable To clarify this suppose that you want to use the function 1 2x7 523 as right hand side
19. s 5 5 Initial mesh and domain Unless you use sample_pde in con nection with an example having a predefined initial mesh you must prescribe the initial mesh and the domain Q with the help of the key words init_node init_elem and domain The value corresponding to init_node must be a matrix with three columns The entries of its th row have the following meaning 1 first coordinate of the ith node 2 second coordinate of the ith node 3 type of the ith node with 0 interior node 1 homogeneous Dirichlet condition 2 inhomogeneous Dirichlet condition 1 Neumann condition 2 Robin condition The default setting for init_node corresponds to a criss cross grid with 8 isosceles right angled triangles as depicted in the left part of figure Bal The value corresponding to init_elem must be a matrix with three columns too The entries in the jth column of its ith row gives the 12 R VERFURTH global number of the locally jth node of the ith element The default setting for init elem again corresponds to a criss cross grid with 8 isosceles right angled triangles as depicted in the left part of figure 5 1 The initial mesh must have the property that for every element its first edge is either a boundary edge or is adjacent to the first edge of an other element where by convention the jth edge is opposite to the locally jth vertex The value corresponding to domain must be a list of length three Its first argument must be a
Download Pdf Manuals
Related Search