Design of Adaptive Finite Element Software

Design of Adaptive Finite Element Software




About this book

During the last years, scientific computing has become an important research branch located between applied mathematics and applied sciences and engineering. Highly efficient numerical methods are based on adaptive methods, higher order discretizations, fast linear and non-linear iterative solvers, multi-level algorithms, etc. Such methods are integrated in the adaptive finite element software ALBERTA. It is a toolbox for the fast and flexible implementation of efficient software for real life applications, based on modern algorithms. ALBERTA also serves as an Environmental-Plant sciences for improving existent, or developing new numerical methods in an interplay with mathematical analysis and it allows the direct integration of such new or improved methods in existing simulation software. The book is accompanied by a full distribution of ALBERTA (Version 1.2) on a CD including an implementation of several model problems. System requirements for ALBERTA are a Unix/Linux Environmental-Plant sciences with C and FORTRAN Compilers, OpenGL graphics and GNU make. These model implementations serve as a basis for students and researchers for the implementation of their own research projects within ALBERTA.



Published by
Published 01 January 2005
Reads 9
EAN13 3540271562
Language English

Legal information: rental price per page €. This information is given for information only in accordance with current legislation.

Report a problem54022842X
In this chapter we describe the implementation of two stationary model prob-lems (the linear Poisson equation and a nonlinear reaction–diffusion equation) and of one time dependent model problem (the heat equation). Here we give an overview how to set up anALBERTAprogram for various applications. We do not go into detail when refering toALBERTAdata structures and functions. A detailed description can be found in Chapter 3. We start with the easy and straight forward implementation of the Poisson problem to learn about the basics ofALBERTA. The examples with the implementation of the nonlinear reaction-diffusion problem and the time dependent heat equation are more involved and show the tools ofALBERTAfor attacking more complex prob-A lems. Removing all L T X descriptions of functions and variables results in the E source code for the adaptive solvers. During the installation ofALBERTA(described in Section 2.4) a subdi-rectoryDEMOwith sources and makefiles for these model problems is cre-ated. The corresponding ready-to-run programs can be found in the files ellipt.c,heat.c, andnonlin.c,nlprob.c,nlsolve.cin the subdirec-toryDEMO/src/Common/. Executable programs for different space dimensions can be generated in the subdirectoriesDEMO/src/1d/,DEMO/src/2d/, and DEMO/src/3d/by callingmake ellipt,make nonlin, andmake heat. Graph-ics output for all problems is generated via a routine void graphics(MESH *mesh, DOF_REAL_VEC *u_h, REAL (*get_est)(EL *)); which shows the geometry given bymesh, as well as finite element func-tion values given byu h, or local estimator values when parameterget est is given, all in separate graphics windows. The source of this routine is DEMO/src/Common/graphics.c, which is self explaining and not described here in detail.
2 Implementation of model problems
2.1 Poisson equation
In this section we describe a model implementation for the Poisson equation d Δu=finΩR, u=gon∂Ω.
This is the most simple elliptic problem (yet very important in many appli-cations), but the program presents all major ingredients for general scalar stationary problems. Modifications needed for a nonlinear problem are pre-sented in Section 2.2.
Fig. 2.1.Solution of the linear Poisson problem and corresponding mesh. The pictures were produced by GRAPE.
Data and parameters described below lead in 2d to the solution and mesh shown in Fig. 2.1. The implementation of the Poisson problem is split into several major steps which are now described in detail.
2.1.1 Include file and global variables
AllALBERTAsource files must include the header filealberta.hwith all ALBERTAtype definitions, function prototypes and macro definitions: #include <alberta.h> For the linear scalar elliptic problem we use four global pointers to data struc-tures holding the finite element space and components of the linear system of equations. These are used in different subroutines where such information cannot be passed via parameters. static const FE_SPACE *fe_space; static DOF_REAL_VEC *u_h = nil; static DOF_REAL_VEC *f_h = nil; static DOF_MATRIX *matrix = nil;
2.1 Poisson equation
fe spacepointer to the actually used finite element space; it is initialized: a by the functioninit dof admin()which is called byGET MESH(), see Section 2.1.4; u h: a pointer to a DOF vector storing the coefficients of the discrete so-lution; it is initialized on the first call ofbuild()which is called by adapt method stat(), see Section 2.1.7; f h: a pointer to a DOF vector storing the load vector; it is initialized on the first call ofbuild(); matrix: a pointer to a DOF matrix storing the system matrix; it is initialized on the first call ofbuild(); The data structureFE SPACEis explained in Section 3.2.14,DOF REAL VECin Section 3.3.2, andDOF MATRIXin Section 3.3.4. Details about DOF adminis-trationDOF ADMINcan be found in Section 3.3.1 and about the data structure MESHfor a finite element mesh in Section 3.6.1.
2.1.2 The main program for the Poisson equation
The main program is very simple, it just includes the main steps needed to implement any stationary problem. Special problem-dependent aspects are hidden in other subroutines described below. We first read a parameter file (indicating which data, algorithms, and solvers should be used; the file is described below in Section 2.1.3). Then we initialize the mesh (the basic geometric data structure), and read the macro triangulation (including an initial global refinement, if necessary). The subdirectoriesMACROin theDEMO/src/*ddirectories contain data for several sample macro triangulations. How to read and write macro triangula-tion files is explained in Section 3.2.16. The macro file name and the number of global refinements are given in the parameter file. Now, the domain’s ge-ometry is defined, and a finite element space is automatically generated via theinit dof admin()routine described in Section 2.1.4 below. A call to graphics()displays the initial mesh. The basic algorithmic data structureADAPT STATintroduced in Section 3.13.1 specifies the behaviour of the adaptive finite element method for sta-tionary problems. A pre–initialized data structure is accessed by the func-tionget adapt stat(); the most important members (adapt>tolerance, adapt>strategy, etc.) are automatically initialized with values from the pa-rameter file; other members can be also initialized by adding similar lines for these members to the parameter file (compare Section 3.13.4). Eventu-ally, function pointers for the problem dependent routines have to be set (estimate,get el est,build,solve). Since the assemblage is done in one step after all mesh modifications, onlyadapt>build after coarsenis used, no assemblage is done before refinement or before coarsening. These addi-tional assemblage steps are possible and may be needed in a more general application, for details see Section 3.13.1.
2 Implementation of model problems
The adaptive procedure is started by a call ofadapt method stat(). This automatically solves the discrete problem, computes the error estimate, and refines the mesh until the given tolerance is met, or the maximal number of iterations is reached, compare Section 3.13.1. Finally,WAIT REALLYallows an inspection of the final solution by preventing a direct program exit with closure of the graphics windows.
int main(int argc, char **argv) { FUNCNAME("main"); MESH *mesh; int n_refine = 0; static ADAPT_STAT *adapt; char filename[100];
/**/ /* first of all, init parameters of the init file */ /**/
/**/ /* get a mesh, and read the macro triangulation from file */ /**/
mesh = GET_MESH("ALBERTA mesh", init_dof_admin, init_leaf_data); GET_PARAMETER(1, "macro file name", "%s", filename); read_macro(mesh, filename, nil); GET_PARAMETER(1, "global refinements", "%d", &n_refine); global_refine(mesh, n_refine*DIM); graphics(mesh, nil, nil);
/**/ /* init adapt structure and start adaptive method */ /**/
adapt = get_adapt_stat("ellipt", "adapt", 2, nil); adapt>estimate = estimate; adapt>get_el_est = get_el_est; adapt>build_after_coarsen = build; adapt>solve = solve;
adapt_method_stat(mesh, adapt); WAIT_REALLY; return(0); }
2.1 Poisson equation
2.1.3 The parameter file for the Poisson equation
The following parameter fileINIT/ellipt.datis used for theellipt.cpro-gram: macro file name: Macro/macro.amc global refinements: 0 polynomial degree: 3
% graphic windows: solution, estimate, and mesh if size > 0 graphic windows: 300 300 300 % for graphics you can specify the range for the values of % discrete solution for displaying: min max % automatical scaling by display routine if min >= max graphic range: 0.0 1.0
solver: 2 % 1: BICGSTAB 2: CG 3: GMRES 4: ODIR 5: ORES solver max iteration: 1000 solver restart: 10 % only used for GMRES solver tolerance: 1.e8 solver info: 2 solver precon: 2 % 0: no precon 1: diag precon % 2: HB precon 3: BPX precon
error norm: estimator C0: estimator C1: estimator C2:
adapt>strategy: adapt>tolerance: adapt>MS_gamma: adapt>max_iteration: adapt>info:
1 % 1: H1_NORM, 2: L2_NORM 0.1 % constant of element residual 0.1 % constant of jump residual 0.0 % constant of coarsening estimate
2 % 0: no adaption 1: GR 2: MS 3: ES 4:GERS 1.e4 0.5 20 8
WAIT: 1 The fileMacro/macro.amcstoring data about the macro triangulation forΩ= d (0,1) can be found in Section 3.2.16 for 2d and 3d. Thepolynomial degree parameter selects the third order finite elements. Bygraphic windows, the number and sizes of graphics output windows are selected. This line is used by thegraphics()routine. For 1d and 2d graphics, the range of function values might be specified (used for graph coloring and height). The solver for the linear system of equations is selected (here: the conjugate gradient solver), and corresponding parameters like preconditioner and tolerance. Parameters for the error estimator include values of different constants and selection of 1 2 the error norm to be estimated (H- orL-norm, selection leads to multipli-cation with different powers of the local mesh size in the error indicators), see Section 3.14.1. An error tolerance and selection of a marking strategy with
2 Implementation of model problems
corresponding parameters are main data given to the adaptive method. Fi-nally, theWAITparameter specifies whether the program should wait for user interaction at additional breakpoints, whenever aWAITstatement is executed as in the routinegraphics()for instance. The solution and corresponding mesh in 2d for the above parameters are shown in Fig. 2.1. As optimal parameter sets might differ for different space dimensions, separate parameter files exist in1d/INIT/,2d/INIT/, and 3d/INIT/.
2.1.4 Initialization of the finite element space
During the initialization of the mesh byGET MESH()in the main program, we have to specify all DOFs which we want to use during the simulation on the mesh. The initialization of the DOFs is implemented in the function init dof admin()which is called byGET MESH(). For details we refer to Sec-tions 3.2.15 and 3.6.2. For the scalar elliptic problem we need one finite element space for the discretization. In this example, we use Lagrange elements and we initialize the degree of the elements via a parameter. The correspondingfe spaceis accessed byget fe space()which automatically stores at the mesh informa-tion about the DOFs used by this finite element space. It is possible to access several finite element spaces inside this function, for instance in a mixed finite element method, compare Section 3.6.2.
void init_dof_admin(MESH *mesh) { FUNCNAME("init_dof_admin"); int degree = 1; const BAS_FCTS *lagrange;
GET_PARAMETER(1, "polynomial degree", "%d", &degree); lagrange = get_lagrange(degree); TEST_EXIT(lagrange)("no lagrange BAS_FCTS\n"); fe_space = get_fe_space(mesh, lagrange>name, nil, lagrange); return; }
2.1.5 Functions for leaf data
As explained in Section 3.2.12, we can “hide” information which is only needed on a leaf element at the pointer to the second child. Such information, which we use here, is the local error indicator on an element. For this elliptic problem we need oneREALfor storing this element indicator. During mesh initialization byGET MESH()in the main program, we have to give information about the size of leaf data to be stored and how to transform
2.1 Poisson equation
leaf data from parent to children during refinement and vice versa during coarsening. The functioninit leaf data()initializes the leaf data used for this problem and is called byGET MESH(). Here, leaf data is one structure struct ellipt leaf dataand no transformation during mesh modifications is needed. The details of theLEAF DATA INFOdata structure are stated in Section 3.2.12. The error estimation is done by the library functionellipt est(), see Section 3.14.1. Forellipt est(), we need a function which gives read and write access to the local element error, and for the marking function of the adaptive procedure, we need a function which returns the local error indicator, see Section 3.13.1. The indicator is stored as theREALmemberestimate ofleaf datastruct ellipt and the functionrw el est()returns for each element a pointer to this member. The functionget el est()returns the value stored at that member for each element.
struct ellipt_leaf_data { REAL estimate; };
one real for the estimate
void init_leaf_data(LEAF_DATA_INFO *leaf_data_info) { leaf_data_info>leaf_data_size = sizeof(struct ellipt_leaf_data); leaf_data_info>coarsen_leaf_data = nil; /* no transformation */ leaf_data_info>refine_leaf_data = nil; /* no transformation */ return; }
static REAL *rw_el_est(EL *el) { if (IS_LEAF_EL(el)) return(&((struct ellipt_leaf_data else return(nil); }
static REAL get_el_est(EL *el) { if (IS_LEAF_EL(el)) return(((struct ellipt_leaf_data else return(0.0); }
2 Implementation of model problems
2.1.6 Data of the differential equation
Data for the Poisson problem are the right hand sidefand boundary values g. For test purposes it is convenient to have access to an exact solution of the problem. In this example we use the function
as exact solution, resulting in
2 10|x| u(x) = e
2 10|x| u(x) =20xe
and 2 210|x| f(x) =Δu(x) =(400|x| −20d) e. d Here,ddenotes the space dimension,ΩR. The functionsuandgrd uare the implementation ofuanduand are optional (and usually not known for a general problem). The functionsgandfare implementations of the boundary values and the right hand side and are not optional.
static REAL u(const REAL_D x) { return(exp(10.0*SCP_DOW(x,x))); }
static const REAL *grd_u(const REAL_D x) { static REAL_D grd; REAL ux = exp(10.0*SCP_DOW(x,x)); int n;
for (n = 0; n < DIM_OF_WORLD; n++) grd[n] = 20.0*x[n]*ux;
return(grd); }
static REAL g(const REAL_D x) { return(u(x)); }
/* boundary values, not optional */
static REAL f(const REAL_D x) /* Delta u, not optional { REAL r2 = SCP_DOW(x,x), ux = exp(10.0*r2); return((400.0*r2  20.0*DIM)*ux); }
2.1 Poisson equation
2.1.7 The assemblage of the discrete system
For the assemblage of the discrete system we use the tools described in Sections 3.12.2, 3.12.4, and 3.12.5. For the matrix assemblage we have to provide an element-wise description of the differential operator. Following the description in Section 1.4.7 we provide the functioninit element()for an initialization of the operator on an element and the functionLALt()for the computation t of det|DFS|ΛAΛon the actual element, whereΛis the Jacobian of the barycentric coordinates,DFSthe the Jacobian of the element parameteriza-tion, andAthe matrix of the second order term. ForΔ, we haveA=idand t det|DFS|ΛΛis the description of differential operator since no lower order terms are involved. For passing information about the JacobianΛof the barycentric coordi-nates and det|DFS|from the functioninit element()to the functionLALt() we use the data structureinfostruct op which stores the Jacobian and the determinant. The functioninit element()calculates the Jacobian and the determinant by the library functionel grd lambda()and the functionLALt() t uses these values in order to compute det|DFS|ΛΛ. Pointers to these functions and to one structureinfostruct op are mem-bers of a structureOPERATOR INFOwhich is used for the initialization of a function for the automatic assemblage of the global system matrix. For more general equations with lower order terms, additional functionsLb0,Lb1, and/orchave to be defined at that point. The full description of the func-tionfill matrix info()for general differential operators is given in Sec-tion 3.12.2. Currently, the functionsinit element()andLALt()only work properly for== DIMDIM OF WORLD . The initialization of theEL MATRIX INFO structure is only done on the first call of the functionbuild()which is called byadapt method stat()during the adaptive cycle (compare Section 3.13.1). By callingdof compress(), unused DOF indices are removed such that the valid DOF indices are consecutive in their range. This guarantees op-timal performance of the BLAS1 routines used in the iterative solvers and admin>size usedis the dimension of the current finite element space. This dimension is printed for information. On the first call,build()also allocates the DOF vectorsu handf h, and the DOF matrixmatrix. The vectoru hadditionally is initialized with zeros and the function pointers for an automatic interpolation dur-ing refinement and coarsening are adjusted to the predefined functions in fe space>bas fcts. The load vectorf hand the system matrixmatrixare newly assembled on each call ofbuild(). Thus, there is no need for interpo-lation during mesh modifications or initialization. On each call ofbuild()the matrix is assembled by first clearing the matrix using the functionclear dof matrix()and then adding element con-tributions byupdate matrix(). This function will callinit element()and LALt()on each element.