Gascoigne3D Manual

Gascoigne3D_manual

Gascoigne3D_manual

User Manual:

Open the PDF directly: View PDF .
Page Count: 42

Introduction to

Gascoigne 3D

High Performance Adaptive Finite Element Toolkit

Tutorial and manuscript by Thomas Richter

thomas.richter@iwr.uni-heidelberg.de

March 9, 2014

2

Contents

1 Introduction 5

1.1 Installation ..................................... 5

1.2 A minimal example solving a partial differential equation . . . . . . . . . . . 7

2 The parameter file 13

3 Description of the problem 15

3.1 Therighthandside ................................ 15

3.2 BoundaryData................................... 16

3.2.1 Dirichlet boundary data . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.2.2 Neumann conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.2.3 Robinconditions.............................. 19

3.3 Definition of the partial differential equations . . . . . . . . . . . . . . . . . . 19

3.3.1 Nonlinear equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.3.2 Equations on the boundary & Robin boundary data . . . . . . . . . . 22

3.4 Exact Solution and Evaluation of Errors . . . . . . . . . . . . . . . . . . . . . 23

3.5 Systems of partial differential equations . . . . . . . . . . . . . . . . . . . . . 24

3.5.1 Modification of the right hand side . . . . . . . . . . . . . . . . . . . . 24

3.5.2 Modification of Dirichlet data . . . . . . . . . . . . . . . . . . . . . . . 25

3.5.3 Modification for the equation, form and matrix . . . . . . . . . . . . . 25

4 Time dependent problems 27

4.1 Time discretization of the heat equation . . . . . . . . . . . . . . . . . . . . . 27

4.2 Time dependent problem data . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

4.3 Non standard time-discretizations . . . . . . . . . . . . . . . . . . . . . . . . . 30

5 Mesh handling 33

5.1 Definition of coarse meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

5.2 Curvedboundaries ................................. 35

5.3 Adaptive refinement of meshes . . . . . . . . . . . . . . . . . . . . . . . . . . 37

5.3.1 Estimating the Energy-Error . . . . . . . . . . . . . . . . . . . . . . . 37

5.3.2 Picking Nodes for Refinement . . . . . . . . . . . . . . . . . . . . . . . 38

6 Flow problems and stabilization 39

6.1 The function point ................................. 39

6.2 Stabilization of the Stokes system . . . . . . . . . . . . . . . . . . . . . . . . . 39

6.3 Stabilization of convective terms . . . . . . . . . . . . . . . . . . . . . . . . . 40

6.4 LPSStabilization.................................. 41

3

Contents

4

1 Introduction

When referring to files, the following directories are used:

HOME Your home directory

GAS The base-directory of Gascoigne. For the workshop:

GAS/src The source- and include-files of Gascoigne.

Gascoigne is a software-library written in C++. To solve a partial differential equation with

Gascoigne, a small program has to be written using the routines provided by Gascoigne. All

functionality is gathered in the library-file File: libGascoigneStd . This library is linked

to every application.

1.1 Installation

In this section we explain the usage of Gascoigne as a library for own applications.

Gascoigne requires a configuration file:

File: HOME/.gascoignerc

Here, some setting like “where is Gascoigne”, “where are the libraries” are given. The

minimal file looks like:

1#############################################

2SET(GASCOIGNE STD ” ˜/ Softw are / s r c / Ga sc oi gne ” )

3

4#############################################

5#Gascoigne L i b r a r i e s & B i n a ri es

6SET(GASCOIGNE LIBRARY OUTPUT PATH ” ˜/ So ftw are / x86 64 / l i b ” )

7SET(GASCOIGNE EXECUTABLE OUTPUT PATH ” . ” )

You need a file like this in the home directory. To build an application using Gascoigne at

least two files and directories are required:

File: TEST/bin/

File: TEST/src/main.cc

File: TEST/src/CMakeLists.txt

5

1 Introduction

The first directory will be used to store the executable file for this application. The other

two files are the source code of the application and a file to control the compilation of the

application. This could look like:

File: TEST/src/CMakeLists.txt

1INCLUDE($ENV{HOME}/ . g a s c o i g n e r c )

2

3INCLUDE( ${GASCOIGNE STD}/CMakeGlobal . tx t )

4LINK LIBRARIES( ${GASCOIGNE LIBS})

5

6ADD EXECUTABLE( ” Test ” main . cc )

In the first line, the user’s config file is read. The following two lines set global parame-

ters for Gascoigne. The variable GASCOIGNE STD is provided in the config-file, the variable

GASCOIGNE LIBS is given in the file CMakeGlobal.txt. The last line indicates the name of

the application: Test, and lists all code-files to be compiled into the application.

A minimal example for a code-file (it does not really calculate something) looks like this:

File: TEST/src/main.cc

1#include <iostream>

2#include ” s t d l o o p . h”

3

4int main ( int argc , char ∗∗ argv)

5{

6std : : cout << ” Create a Gascoigne Loop ! ” << std : : endl ;

7Gascoigne : : StdLoop Loop ;

8}

For building the application we use the program cmake to generate Makefiles. In the directory

TEST/bin/ type:

1cmake . . / s r c

2make

The first line creates the File: TEST/bin/Makefile . It reads the File: CMakeLists.txt

to put everything together. The second command starts the compilation and linking of the

application and if everything works out, File: TEST/bin/Test should be the executable

file which then can be started.

The command cmake is only needed to create new Makefiles. This is necessary, if new

code-files are added to the application. If only source code is changed, e.g. the file File:

TEST/src/main.cc , you only need to call make in the directory TEST/bin.

6

1.2 A minimal example solving a partial differential equation

1.2 A minimal example solving a partial differential equation

Assume there exists a directory FIRSTEXAMPLE with sub-dirs FIRSTEXAMPLE/bin and

FIRSTEXAMPLE/src. In the src-directory you need the files File: main.cc and File:

CMakeLists.txt . in the bin-directory you call cmake ../src for configuration and then

make for compilation.

File: FIRSTEXAMPLE/src/main.cc

1#include ” s t d l o o p . h”

2#include ” p r o b l e m d e s c r i p t o r b a s e . h”

3#include ”constantrighthandside .h”

4#include ”zerodirichletdata .h”

5#include ” l a p l a c e 2 d . h”

6

7using namespace Gascoigne ;

8

9class Problem : public ProblemDescriptorBase

10 {

11 public :

12 std : : s t r i n g GetName ( ) const {return ” F i r s t Problem” ; }

13 void BasicInit(const ParamFile∗p a r a m f i l e )

14 {

15 GetEquationPointer () = new Laplace2d ;

16 GetRightHandSidePointer ( ) = new OneRightHandSideData ( 1 ) ;

17 GetDirichletDataPointer () = new Z e r o D i r i c h l e t D a t a ;

18

19 ProblemDescriptorBase : : B a s i c I n i t ( p a ra mf i le ) ;

20 }

21 };

22

23 int main ( int argc , char ∗∗ argv )

24 {

25 ParamFile p a r a m f i l e ( ” f i r s t . param” ) ;

26 i f ( argc >1) p a r a m f i l e . SetName ( ar gv [ 1 ] ) ;

27

28 Problem PD;

29 PD. B a s i c I n i t (& pa ra mf i le ) ;

30

31 ProblemContainer PC;

32 PC. AddProblem ( ” l a p l a c e ” , &PD) ;

33

34 StdLoop l o op ;

35 lo op . B a s i c I n i t (& pa ramf i l e , &PC) ;

36 lo op . run ( ” l a p l a c e ” ) ;

37 }

7

1 Introduction

The Problem

Here, the problem to be solved is defined. The Problem gathers all information on the

problem to be solved. This is: which pde, what kind of boundary data, what is the right

hand side, and so on.

By PD.BasicInit() we initialize the problem (will be described later on). This problem is

then added to another class, the ProblemContainer. Here, we can gather different problems,

each together with a keyword, which is laplace in this case. There are applications where

it is necessary to solve different problems in the same program.

The problem set pointers to a large variety of classes. Here, we define the equation to be of

type Laplace2d, an implementation of the two dimensional Laplace equation. For right hand

side and Dirichlet data we define self written classes (will be explained later). In Gascoigne

a very large set of parameters can be set for a problem, all in form of pointers. However these

three classes, Equation,RightHandSide and DirichletData are necessary for every instance

of Gascoigne. The class Problem is derived from the base-class ProblemDescriptorBase:

File: GAS/src/Problem/problemdescriptorbase.h

The equation Laplace2d is part of the Gascoigne-Library:

File: GAS/src/Problem/laplace2d.h

This is explained in detail in Section 3.

The Loop

Then, a Loop is specified, the loop is the class controlling what is happening at run-time.

and what Gascoigne is doing: this is usually a sequence like: initialize, solve a pde, write

the solution to the disk, evaluate functionals, compute errors, refine the mesh and start over.

In BasicInit() we initialize the most basic structure and read data from parameter files.

Here we also pass the ProblemContainer to the different instances.

Finally, the equation is solved by calling loop.run("laplace"), where the keyword laplace

indicates what problem to solve. Let us now have a closer look at this function which is

defined in

File: FIRSTEXAMPLE/src/loop.cc

First, in run(), we define two vectors u,f for storing the solution and the right hand side.

Then, in the for-loop we initialize the problem, solve the problem, write out the solution and

refine the mesh. This is very typical for every program run, since for a reliable simulation it

is necessary to observe the convergence of the solution on a sequence of meshes. Next, we

have a closer look at every step of the inner loop:

8

1.2 A minimal example solving a partial differential equation

1G e tM u l ti L ev e l So l ve r ()−>ReInit(problemlabel );

2G e tM u l ti L ev e l So l ve r ()−>R eI n i tV e ct o r ( u ) ;

3G e tM u l ti L ev e l So l ve r ()−>Re I n i t V ector ( f ) ;

4

5GetSolverInfos()−>GetNLInfo ( ) . c o n t r o l ( ) . mat rixm ustb ebui ld ( ) = 1 ;

The first line initializes the basic ingredients of Gascoigne. In particular, we pass the equa-

tion to be solved and all the necessary data to different parts like “Solver”, “Discretization”.

Then, by GetMultiLevelSolver()->ReInitVector(u) we initialize the vector on the cur-

rent mesh. This basically means that we resize the vector to match the number of degrees

of freedom (which changes with the mesh-size).

Finally, in the last line we tell Gascoigne to assemble the system matrix. This is necessary

whenever the system matrix changes, which is the case if the mesh changes, or if the matrix

depends on the solution uitself.

File: Template1/src/loop.cc

1S o lv e ( u , f , ” R e s u l t s /u” ) ;

2GetMeshAgent()−>global refine ();

By calling Solve(u,f) we solve the PDE. The subsequent line refines the mesh.

File: Template1/src/loop.cc

1s t r i n g Loop : : S ol ve ( V e c t o r I n t e r f a c e& u , V e c t o r I n t e r f a c e& f , s t r i n g

2name)

3{

4G e tM u l ti L ev e l So l ve r ()−>GetS olver ()−>Zero ( f ) ;

5G e tM u l ti L ev e l So l ve r ()−>GetS olver ()−>Rhs ( f ) ;

6

7G e tM u l ti L ev e l So l ve r ()−>GetS olver ()−>SetBoundaryVector( f );

8G e tM u l ti L ev e l So l ve r ()−>GetS olver ()−>SetBoundaryVector(u);

9

10 s t r i n g s t a t us = G et MultiLeve lS olver ()−>

11 S ol v e ( u , f , G e t S o l v e r I n f o s ()−>GetNLInfo ( ) ) ;

12

13 Output ( u , name ) ;

14 }

Here, we first assemble the right hand side fof our problem. Then, we need to initialize

Dirichlet boundary values. Finally, we tell the multigrid solver to solve the equation.

In the last line, we write the solution into a file for visualization.

9

1 Introduction

Output of Gascoigne

After starting the code, Gascoigne first prints out some information on the current mesh,

solver, Discretization and problem data.

Then, for every iteration of the inner loop the output looks like:

...

================== 2 ================ [l,n,c] 5 1089 1024

0: 9.77e-04

M 1: 7.73e-09 [0.00 0.00] - 6.26e-08 [0.012] {3}

[u.00002.vtk]

...

In the first line the current iteration count (2) , the number of mesh-refinemenet levels (5),

the number of mesh nodes (1089) and the number of mesh elements (quads) (1024) is given.

Then, Gascoigne prints control and stastic data for the solution of the algebraic systems.

Gascoigne solves every equation with a Newton-method, even linear ones. This is for reasons

of simplicity as well as having a defect correction method for a better treatment of roundoff

errors. On the left side of the output the convergence history of the newton method is

printed. Here we only need one Newton step. Before this step, the residual of the equation

was 9.77e−4, after this one step the residual was reduced to 7.73e−9. The letter Mindicates

that Gascoigne assembled a new system matrix in this step. The folling two numbers [0.00

0.00] indicate the convergance rates of the newton iteration by printing the factor, by which

the residual was reduced. We only print the first 2 digits, hence Gascoigne gives a 0. The

first number is the reduction rate in the current step, the second number is the average

reduction rate over all Newton steps.

On the right side, the convergence history of the linear multigrid solver is printed. Here,

the values indicate that in newton step 1 we needed {3}iteration of the multigrid solver. In

every of these three steps the error was reduced by a factor of [0.012] and after these three

steps the residual had the value 6.26e-08.

Finally Gascoigne prints the name of the output file. You can visualize it by calling

visusimple u.00002.vtk in the terminal window.

The parameter file

In this directory you can also specify the parameter-file

File: FIRSTEXAMPLE/src/first.param

1// Block Loop

2niter 5

3

4// Block Mesh

5dimension 2

10

1.2 A minimal example solving a partial differential equation

6gridname sq uare . inp

7prerefine 3

8

9// Block BoundaryManager

10 d i r i c h l e t 4 1 2 3 4

11

12 d i r i c h l e t c o m p 1 1 0

13 d i r i c h l e t c o m p 2 1 0

14 d i r i c h l e t c o m p 3 1 0

15 d i r i c h l e t c o m p 4 1 0

16

17 // Block Nix

This minimal example tells Gascoigne the dimension of the problem, which domain to use,

where to apply boundary data.

When starting the program, Gascoigne reads in a parameter file. Here, different parameters

controlling the behaviour of the code can be supplied. For instance in //Block Loop of this

file you can change the values of niter, which tells, how many times Gascoigne solves the

equation and refine the mesh.

In //Block Mesh we supply the mesh file as well as the number of refinement steps realized

before solving the first problem. You can change this parameter in prerefine. Thoe //Block

BoundaryManager will be explained in detail in Section 3.2

11

1 Introduction

12

2 The parameter file

Gascoigne can process a parameter file during runtime. Here, we can set some parameters

and options that can be changed at run-time without need for recompilation of the program.

For the file to be processed, go to the directory containing the parameter file and call

1. . / bin / Test t e s t . param

The data in the param-file is organized in blocks which are read at different points of the

code. A sample parameter file could look like: File: TEST/src/test.param

1// Block S e t t i n g s

2month 10

3year 2009

4name FirstTestOfGascoigne

5

6// Bl ock S o l v e r

7. . .

8

9// Block Nix

Each block has to be initialized by the line

1// Block <blockname>

Note that the slashes in front of the blocks are not treated as a comment here. The parameter

file has to end with a dummy-block, e.g. //Block Nix.

In Gascoigne, the parameter file is stored in the class ParamFile. Reading values from

the parameter-file is done using the classes FileScanner and DataFormatHandler. These

classes are declared in the header File: filescanner.h . Here is a modified File: main.cc

reading values from the parameter file:

File: TEST/src/main.cc

1#include <iostream>

2#include ” s t d l o o p . h”

3#include ” f i l e s c a n n e r . h”

4#include ” p a r a m f i l e . h”

5

6using namespace Gascoigne ;

7using namespace st d ;

8

13

2 The parameter file

9int main ( int argc , char ∗∗ argv)

10 {

11 i f ( argc <2)

12 {c e r r << ” c a l l : T est p a r a m f i l e ! ! ” << end l ;

13 ab ort ( ) ; }

14 ParamFile p a r a m f i l e ;

15 p a r a m f i l e . SetName ( a rgv [ 1 ] ) ;

16

17 int month , year ;

18 s t r i n g name ;

19 DataFormatHandler DFH;

20 DFH. i n s e r t ( ”month” ,&month ) ;

21 DFH. i n s e r t ( ” yea r ” , &year , 2009 ) ;

22 DFH. i n s e r t ( ”name” ,&name , ” t e s t ” ) ;

23 F i l e S c a n n e r FS(DFH, &p a r a m f i le , ” S e t t i n g s ” ) ;

24 cout << month << ”/” << year << ”\t ” << name << end l ;

25 }

To read in one block, we first initialize a DataFormatHandler (line 19). Next, we tell the

DataFormatHandler which keywords to look for using the function insert. The first ar-

gument of the insert function is the keyword we are searching, the second argument the

variable that its value will be written to. Optionally, we can specify a default value as a

third argument which is used if the keyword is not present in the parameter file. Finally, we

initialize an object FileScanner passing the DataFormatHandler, a parameter file and the

block to be read. Note that only this block is read! To read from a different block within

one file, we have to define a second DataFormatHandler.

14

3 Description of the problem

File: GAS/src/Problem/problemdescriptorbase.h

The class ProblemDescriptor provides all information necessary to set up an application.

This includes all problem data like: equation, boundary data, right hand side, initial con-

dition, etc. The most important function is ProblemDescriptor::BasicInit(), where the

problem is defined by setting pointers to the describing classes. A minimal example might

be:

1void Lo calPr oblem Desc ripto r : : B a s i c I n i t ( const ParamFile∗pf )

2{

3GetEquationPointer () = new LocalEquation ;

4GetRightHandSidePointer ( ) = new LocalRightHandSide ;

5GetDirichletDataPointer () = new LocalDirichletData ;

6

7ProblemDescriptorBase : : B a s i c I n i t ( pf ) ;

8}

For an problem to be well defined, at least the equation and the Dirichlet data have to be

specified. All possible pointers to be set are listed in

File: GAS/src/Problem/problemdescriptorbase.h

For the most basic problems Gascoigne provides predefined classes, e.g.

1GetEquationPointer () = new Laplace2d ;

2GetRightHandSidePointer ( ) = new OneRightHanSideData ( 1 ) ;

3GetDirichletDataPointer () = new Z e r o D i r i c h l e t D a t a ;

which defines the Laplace problem with homogeneous Dirichlet data and constant right-hand

side f= 1.

In the following we describe the classes needed to specify the most common ingredients of a

problem.

3.1 The right hand side

File: GAS/src/Interface/domainrighthandside.h

The interface class to specify a function f: Ω →Ras right hand side for the problem is

DomainRightHandSide. Let f(x, y) = sin(x) cos(y) be the right hand side:

15

3 Description of the problem

1class LocalRightHandSide : public DomainRightHandSide

2{

3public :

4int GetNcomp ( ) const {return 1 ; }

5

6std : : s t r i n g GetName ( ) const {return ” Lo cal Right Hand Side ” ; }

7

8double operator ( ) ( in t c , const Vertex2d& v ) const

9{

10 return s i n ( v . x ( ) ) ∗c os ( v . y ( ) ) ;

11 }

12 }

The most important function of this class is double operator()(int c, const Vertex2d&

v) const, here, the value of the right hand side function fat point v∈Ω is returned. The

class Vertex2d is explained in File: GAS/src/Common/vertex.h . For three dimensional

problems the operator double operator()(int c, const Vertex3d& v) const has to be

written.

The second parameter int c as well as the function int GetNcomp() const is only needed

if we solve a system of equations. This is explained in detail in Section 3.5.

A function called GetName() is necessary in every class used in the ProblemDescriptor to

provide a label.

3.2 Boundary Data

In the geometry file, a certain color is assigned to every boundary part of the mesh, see

Section 5. You can visualize meshes (inp-files) with the program visusimple. Here you can

also look at the boundary data, use Select->Scalar. In the parameter file we specify all

boundary colors where Dirichlet boundary conditions are used.

To assign Dirichlet boundary condition to a certain boundary color, the parameter file has

to contain:

1// Block BoundaryManager

2d i r i c h l e t 2 4 8

3d i r i c h l e t c o m p 4 1 0

4d i r i c h l e t c o m p 8 1 0

This definition lets the boundaries with colors 4and 8have Dirichlet condition. The param-

eter dirichlet is given as a vector: the first value is the number of Dirichlet colors, then

the colors follow. If only the color 4is picked for Dirichlet it would look like:

1// Block BoundaryManager

2d i r i c h l e t 1 4

3d i r i c h l e t c o m p 4 1 0

16

3.2 Boundary Data

We have to specify a second parameter dirichletcomp for every boundary color used with

Dirichlet values. Here, we choose the components of the solution where Dirichlet conditions

shall be used. These values are again given as vectors: First the color is specified, then the

number of solution components that have Dirichlet values on this color. Finally the solution

components are listed. For scalar equations with only one component, the last two values

are always 1 0. See Section 3.5 for using Dirichlet colors for systems of partial differential

equations with more than one solution component.

Every boundary color, that is not listed as Dirichlet boundary uses natural boundary condi-

tions given due to integration by parts. For the Laplace Equation

(−∆u, φ)Ω+h∂nu, φiΓ= (∇u, ∇φ)Ω,

this is the homogeneous Neumann condition ∂nu= 0. How to use non-homogeneous Neu-

mann is explained in Section 3.2.2 and details on Robin-boundary conditions are given in

Section 3.2.3.

3.2.1 Dirichlet boundary data

File: src/Problem/diricheltdata.h

Assume, that the colors 4and 8are picked as Dirichlet boundary colors in the parameter

file.

The Dirichlet values to be used for the boundaries with colors specified in the parameter file

are described in the class DirichletData, see File: GAS/src/Interface/dirichletdata.h

Here we show an example how to define a class LocalDirichletData that sets Dirichlet

values

1

2#include ”dirichletdata .h”

3

4class LocalDirichletData : public DirichletData

5{

6public :

7std : : s t r i n g GetName ( ) const {return ” Lo cal D i r i c h l e t Data ; ” }

8void operator ()( DoubleVector& b, const Vertex2d& v , int c o l ) const

9{

10 b . zer o ( ) ;

11 i f ( c o l ==4) b [ 0 ] = 0 . 0 ;

12 i f ( c o l ==8) b [ 0 ] = v . x ( ) + v . y ( ) ;

13 }

14 }

This example sets the boundary value to zero on the boundary with color 4and to the

function g(x, y) = x+yon color 8. The values are not returned as in the case of the right hand

side, but written in the vector b. The value b[i] defines the Dirichlet value of the i-th solution

17

3 Description of the problem

component. For scalar pde’s we only set b[0]. The parameter Vertex2d& v gives the

coordinate and the parameter col gives the color of the node. In three dimensions, the same

operator exists using a Vertex3d& v to indicate the coordinate. In the ProblemDescriptor

we need to set a pointer to this new class: File: TEST/src/problem.h

1class ProblemDescriptor : public ProblemDescriptorBase

2{

3// . . .

4void BasicInit( const ParamFile∗p a r a m f i l e )

5{

6// . . .

7GetDirichletDataPointer () =

8new LocalDirichletData ;

9}

10 // . . .

11 };

3.2.2 Neumann conditions

For non-homogenous Neumann conditions of the type

h∂nu, φiΓN=hgN, φiΓN,

we have to add the term (gN, φ) to the right hand side. In Gascoigne, this Neumann

right hand side is derived from the class BoundaryRightHandSide specified in the File:

GAS/src/Interface/boundaryrighthandside.h . An example for a 2d scalar problem is

given by

1#include ”boundaryrighthandside .h”

2

3class LocalBoundaryRightHandSide : public BoundaryRightHandSide

4{

5public :

6int GetNcomp ( ) const {return 1 ; }

7s t r i n g GetName ( ) const {return ” Lo cal B−RHS” ; }

8

9double operator ( ) ( in t c , const Vertex2d& v , const Vertex2d& n ,

10 int c o l o r ) const

11 {

12 i f ( c o l o r ==0) return 1.0;

13 i f ( c o l o r ==1) return v . x ( ) ∗n . x ( ) + v . y ( ) ∗n . y ( ) ;

14 }

15

16 };

Here, we define Neumann data on two different parts of the boundary with colors 0and 1:

g0(x) = 1, g1(x) = n(x)·x,

18

3.3 Definition of the partial differential equations

where n(x) is the outward unit-normal vector in the point x.

Boundary right hand sides need to be specified in the parameter file in //Block BoundaryManager.

Otherwise these boundary terms are not taken into account: File: TEST/src/test.param

1// Block BoundaryManager

2d i r i c h l e t 1 4

3d i r i c h l e t c o m p 4 1 0

4righthandside 2 0 1

The new parameter righthandside specifies a vector. Here it tells Gascoigne that 2bound-

ary colors have Neumann data. These are the colors 0and 1. In the ProblemDescriptor

we need to set a pointer to this new class: File: TEST/src/problem.h

1class ProblemDescriptor : public ProblemDescriptorBase

2{

3// . . .

4void BasicInit( const ParamFile∗p a r a m f i l e )

5{

6// . . .

7GetBoundaryRightHandSidePointer () =

8new LocalBoundaryRightHandSide ;

9}

10 // . . .

11 };

3.2.3 Robin conditions

Robin boundary data includes very general conditions to be fulfilled on the boundary of the

domain. We can have:

hG(u), φiΓR= 0,

where G(·) can be some operator, e.g. G(u) = ∂nu+u2. Hence, Robin boundary data

means that we have an additional equation that is valid on the boundary. We will explain

this concept in detail after explaining how to specify equations. See Section 3.3.2.

3.3 Definition of the partial differential equations

File: GAS/src/Interface/equation.h

Let the pde to be solved be given in the weak formulation

a(u)(φ)=(f, φ)∀φ,

19

3 Description of the problem

where a(·)(·) is a semi-linear form, linear in the second argument. Gascoigne solves every

problem with a Newton method (also linear problems). With an initial guess u0updates wk

are defined by

a0(uk)(wk, φ)=(f, φ)−a(uk)(φ)∀φ, uk+1 := uk+wk.

The Jacobi matrix is the matrix of the directional derivatives of a(·)(·) in the point ukand

defined by

a0(u)(w, φ) := d

dsa(u+sw)(φ)







s=0

.

For linear problems, it holds

a(u)(w, φ) = a(w)(φ).

To solve, Gascoigne needs to know about the form a(uk)(φ) and its derivative, the matrix

a0(uk)(wk, φ). Form and matrix are given in the class Equation:

1class LocalEquation : public Equation

2{

3public :

4

5int GetNcomp ( ) const ;

6s t r i n g GetName ( ) const ;

7

8void po in t ( double h , const Vertex2d& v ) const ;

9void po in t ( double h , const Vertex3d& v ) const ;

10

11 void Form(VectorIterator b, const FemFunction& U,

12 TestFunction& N) const ;

13 void Matrix ( EntryMatrix& A, const FemFunction& U,

14 const TestFunction& M, const TestFunction& N) const ;

15

16 };

The first function GetNcomp() returns the number of solution components for systems of

partial differential equations. For scalar equations, this function returns 1. GetName() is

a label for the equation. The most important functions are Form(b,U,N), which defines

a(u)(φ) and Matrix(A,U,M,N) which gives the derivative a0(u)(w, φ). The parameters band

Aare the return values, Uis the last approximation u,Nis the test-function φand Mstands

for the direction w. The function point() is meant to set parameters depending on the

mesh size hor on the coordinate v.point is called before each call of Form() or Matrix()

and can be used to set local variables.

For the Laplace equation, the implementation is given by

1void Form(VectorIterator b, const FemFunction& U,

2const , TestFunction& N) const

3{

20

3.3 Definition of the partial differential equations

4b [ 0 ] += U [ 0 ] . x ( ) ∗N. x ( ) + U [ 0 ] . y ( ) ∗N. y ( ) ;

5}

6

7void Matrix ( EntryMatrix& A, const FemFunction& U,

8const TestFunction& M, const TestFunction& N) const

9{

10 A( 0 , 0 ) += M. x ( ) ∗N. x ( ) + M. y ( ) ∗N. y ( ) ;

11 }

The class TestFunction describes the values and derivatives of a discrete function in a

certain point. By N.m() the value is accessed, by N.x(),N.y() and N.z() the directional

derivatives. The class FemFunction is a vector of TestFunctions and used for the solution

function u. For systems of partial differential equations, the index gives the number of the

solution component. For scalar equations, it is always U[0].m().

3.3.1 Nonlinear equations

As example we now consider the nonlinear partial differential equation given by

a(u)(φ)=(∇u, ∇φ)+(h∇u, ∇ui, φ),hx, yi:= Xxiyi.

The form is given by

1void Form(VectorIterator b, const FemFunction& U,

2const , TestFunction& N) const

3{

4b [ 0 ] += U [ 0 ] . x ( ) ∗N. x ( ) + U [ 0 ] . y ( ) ∗N. y ( )

5+ (U [ 0 ] . x ( ) ∗U [ 0 ] . x ( ) + U [ 0 ] . y ( ) ∗U [ 0 ] . y ( ) ) ∗N.m( ) ;

6}

To define the matrix we first build the derivative:

a0(u)(w, φ) = d

dsa(u+sw)(φ)|s=0

=d

ds∇(u+sw),∇φ+d

dsh∇(u+sw),∇(u+sw)i, φ|s=0

= (∇w, ∇φ)+(h∇(u+sw),∇wi, φ)+(h∇w, ∇(u+sw)i, φ)|s=0

= (∇w, ∇φ) + 2(h∇u, ∇wi, φ).

The matrix function is then given by

1void Matrix ( EntryMatrix& A, const FemFunction& U,

2const TestFunction& M, const TestFunction& N) const

3{

4A( 0 , 0 ) += M. x ( ) ∗N. x ( ) + M. y ( ) ∗N. y ( ) ;

5A( 0 , 0 ) += 2 . 0 ∗(U [ 0 ] . x ( ) ∗M. x ( ) + U [ 0 ] . y ( ) ∗M. y ) ∗N.m( ) ;

6}

21

3 Description of the problem

3.3.2 Equations on the boundary & Robin boundary data

Robin boundary data is realized by setting an equation on the boundary. First, we need to

tell Gascoigne to include boundary equations. Herefore, section //Block BoundaryManager

of the parameter file needs to know on what boundary colors we set an equation: File:

TEST/src/test.param

1// Block BoundaryManager

2d i r i c h l e t 1 4

3d i r i c h l e t c o m p 4 1 0

4equ a t i o n 1 1

Here we tell Gascoigne to use the equation on 1color, the color 1. The actual equation is

of type BoundaryEquation defined in the File: GAS/src/Interface/boundaryequation.h

The boundary equation is declared as the equation: you need to specify the Form and the

Matrix. Major difference is that these two function get the outward unit normal vector

and the color of the boundary line in addition. Boundary equations are then defined in the

ProblemDescriptor by setting: File: TEST/src/problem.h

1class ProblemDescriptor : public ProblemDescriptorBase

2{

3// . . .

4void BasicInit( const ParamFile∗p a r a m f i l e )

5{

6// . . .

7GetBoundaryEquationPointer () =

8new LocalBoundaryEquation ;

9}

10 // . . .

11 };

We next give an example of a boundary equation which will realize the Robin data ∂nu+u2=

0, see File: TEST/src/problem.h

1#include ” boundaryequation . h”

2

3

4class LocalBoundaryEquation : public BoundaryEquation

5{

6public :

7

8int GetNcomp ( ) const {return 1 ; }

9s t r i n g GetName ( ) const {return ” Lo cal B−EQ” ; }

10

11 void pointboundary ( double h , const Vertex2d& v , const V ertex2d& n )

12 const

13 {}

14

22

3.4 Exact Solution and Evaluation of Errors

15 void Form(VectorIterator b, const FemFunction& U,

16 TestFunction& N, i nt c o l ) const

17 {

18 i f ( c o l ==1) b [ 0 ] += U [ 0 ] .m( ) ∗U [ 0 ] .m( ) ∗N.m( ) ;

19 }

20

21 void Matrix ( EntryMatrix& A, const FemFunction& U,

22 const TestFunction& M, const TestFunction& N, i nt c o l )

23 const

24 {

25 i f ( c o l ==1) A( 0 , 0 ) += 2 . 0 ∗U[ 0 ] .m( ) ∗M.m( ) ∗N.m( ) ;

26 }

27

28 };

3.4 Exact Solution and Evaluation of Errors

If an analytic solution is known it can be added to the ProblemDescriptor and then be

used for evaluating the error in the L2,H1and L∞norms. The operator in the class

ExactSolution is used to define the solution function. The parameter int c specifies the

solution component for systems of pdes.

File: ../src/problem.h

1#i n c l ud e ” e x a c t s o l u t i o n . h”

2

3class MyExactSolution : public ExactSolution

4{

5public :

6std : : s t r i n g GetName ( ) const {return ”My e xa ct s o l u t i o n ” ; }

7double operator ( ) ( in t c , const Vertex2d& v ) const

8{

9return v . x ( ) ∗v . y ( ) ;

10 }

11 };

12 class Problem : public ProblemDescriptorBase

13 {

14 . . .

15 public :

16 void BasicInit(const ParamFile∗p f )

17 {

18 . . .

19 GetExactSolutionPointer () = new MyExactSolution ;

20 }

21 };

23

3 Description of the problem

The function ComputeGlobalErrors which calculates the L2, L∞and H1norm error has to

be called in the Loop after solving the equation with: File: ../src/myloop.cc

1. . .

2Sol v e (u , f ) ;

3ComputeGlobalErrors(u);

4. . .

3.5 Systems of partial differential equations

One of the main features of Gascoigne is the efficient solution of large systems of partial

differential equations. We call the number of equations in a system (and thus the number

of solution variables) the number of components, denoted by ncomp in Gascoigne. Nearly all

classes have the number of components as member. As example we discuss the two-equation

system in Ω:

(∇u0,∇φ0)+(u0u1, φ0) = (f0, φ0),

(∇u1,∇φ1)+(u0sin(u1), φ1) = (f1, φ1),

u0= 0 on Γ4∪Γ8,

u1= 1 on Γ4,

∂nu1= 0 on Γ8.

Γ8

Γ4

Γ4

Γ8

Ω

Several modifications are necessary:

3.5.1 Modification of the right hand side

See 3.1 for comparison. The right hand side has to know the number of equations, in this

case:

1int GetNcomp ( ) const {return 2 ; }

The operator for specifying the functions f0and f1gets the current component as first

argument. As an example, let us set f0(x, y) = 5 and f1(x, y) = xy:

1double operator ( ) ( in t c , const Vertex2d& v ) const

2{

3i f ( c==0) return 5 ;

4i f ( c==1) return v . x ( ) ∗v . y ( ) ;

5return 0.0;

6}

24

3.5 Systems of partial differential equations

3.5.2 Modification of Dirichlet data

For Dirichlet data two separate issues have to be considered: first the declaration in the

parameter file has to be modified. Second, the function specifying the Dirichlet data has to

be altered. In this example, where the component u0has Dirichlet data on both and u1just

on one boundary color, we set:

1// Block BoundaryManager

2d i r i c h l e t 2 4 8

3d i r i c h l e t c o m p 4 2 0 1

4d i r i c h l e t c o m p 8 1 0

This is to be read as: there are 2 dirichlet colors, the colors 4 and 8. Dirichlet color 4 is

applied to 2 components, components 0 and 1. And so on. Gascoigne uses 0 as index for the

first component. The implementation of the boundary data is as follows:

1class LocalDirichletData : public DirichletData

2{

3public :

4std : : s t r i n g GetName ( ) const {return ” L ocal D i r i c h l e t Data ; ” }

5void operator ()( DoubleVector& b, const Vertex2d& v , int c o l ) const

6{

7b . zer o ( ) ;

8i f ( c o l ==4) b [ 0 ] = 0 . 0 ;

9i f ( c o l ==4) b [ 1 ] = 1 . 0 ;

10 i f ( c o l ==8) b [ 0 ] = 0 . 0 ;

11 }

12 };

Here, the values are sorted into a vector whose components correspond to the components

of the solution u.

3.5.3 Modification for the equation, form and matrix

Like the right hand side, the equation needs to know about the number of solution compo-

nents given in the header file as

1int GetNcomp ( ) const {return 2 ; }

In both functions, Form and Matrix the last approximate solution Uis given as a parameter.

This parameter is a vector and U[c] specifies component cof the solution. The modifications

for the Form are straight-forward, the values for i-th equation are written to the component

b[i]. The test function Nis the same for all equations. We give an example:

1void Form(VectorIterator& b, const FemFunction& U,

2const TestFunction& N) const

3{

25

3 Description of the problem

4b [ 0 ] += U [ 0 ] . x ( ) ∗N. x ( ) + U [ 0 ] . y ( ) ∗N. y ( ) ;

5b [ 0 ] += U [ 0 ] .m( ) ∗U [ 1 ] .m( ) ∗N.m( ) ;

6

7b [ 1 ] += U [ 1 ] . x ( ) ∗N. x ( ) + U [ 1 ] . y ( ) ∗N. y ( ) ;

8b [ 1 ] += U [ 0 ] .m( ) ∗s i n (U [ 1 ] .m( ) ) ∗N.m( ) ;

9}

The Matrix works in a simular way. We assemble a local matrix of size ncomp times ncomp,

where A(i,c) specifies the derivative of equation iwith respect to the component U[c]. In

this example:

A00 =d

ds(∇(u0+sw),∇φ0) + ((u0+sw)u1, φ0)

s=0 = (∇w, ∇φ0)+(wu1, φ0).

A01 =d

ds(∇u0,∇φ0)+(u0(u1+sw), φ0)

s=0 = (u0w, φ0)

A10 =d

ds(∇u1,∇φ1) + ((u0+sw) sin(u1), φ1) = (wsin(u1), φ1),

A11 =d

ds(∇u1,∇φ1)+(u0sin(u1+sw), φ1) = (∇w, ∇φ1)+(u0cos(u1)w, φ1),

and the implementation:

1void Matrix ( EntryMatrix& A, const FemFunction& U,

2const TestFunction& M, const TestFunction& N) const

3{

4A( 0 , 0 ) += M. x ( ) ∗N. x ( ) + M. y ( ) ∗N. y ( ) ;

5A( 0 , 0 ) += M.m( ) ∗U [ 1 ] .m( ) ∗N.m( ) ;

6A( 0 , 1 ) += U [ 0 ] .m( ) ∗M.m( ) ∗N.m( ) ;

7

8A( 1 , 0 ) += M.m( ) ∗s i n (U [ 0 ] .m( ) ) ∗N.m( ) ;

9A( 1 , 1 ) += M. x ( ) ∗N. x ( ) + M. y ( ) ∗N. y ( ) ;

10 A( 1 , 1 ) += U [ 0 ] .m( ) ∗cos (U [ 1 ] .m( ) ) ∗M.m( ) ∗N.m( ) ;

11 }

26

4 Time dependent problems

Gascoigne provides a simple technique for the modeling of time-dependent partial differential

equations like

∂tu+A(u) = fin [0, T ]×Ω,

with a spatial differential operator A(u). Time-discretization is accomplished with Rothe’s

method by first discretizing in time with simple time-stepping methods. For time-dependent

problems, we must provide some additional information:

•We must supply an initial value u(x, 0) = u0(x).

•All problem data like the right hand side, boundary values and differential operator

can depend on the time t∈[0, T ].

•We must describe the time-stepping scheme and choose a time-step size.

4.1 Time discretization of the heat equation

Here, we examplarily demonstrate the time discretization of the heat equation

∂tu−∆u=fin [0, T ]×Ω,

with Dirichlet boundary

u=gon [0, T ]×Ω,

and the initial value

u(·,0) = u0(·) on Ω.

In Gascoigne, time discretization is based on the general θ-scheme

1

θk (um, φ) + a(um, φ)=(f(tm), φ) + 1−θ

θ(f(tm−1), φ) + 1

θk (um−1, φ)−1−θ

θa(um−1, φ),

where θ∈(0,1]. For θ=1

2, this scheme corresponds to the Crank-Nicolson method, for

θ= 1 to the implicit Euler method.

The class StdTimeLoop includes a simple loop for time-stepping problems and replaces the

StdLoop:

27

4 Time dependent problems

1void StdTimeLoop : : run ( const std : : s t r i n g& pro b l em l a be l )

2{

3[ . . . ]

4for ( i t e r =1; i t e r <= n i t e r ; i t e r ++)

5{

6[ . . . ]

7t i m e i n f o . SpecifyScheme ( i t e r ) ;

8TimeInfoBroadcast ();

9

10 // RHS

11 G e tM u l ti L ev e l So l ve r ()−>GetS olver ()−>GetGV( f ) . z e r o ( ) ;

12 G e tM u l ti L ev e l So l ve r ()−>GetS olver ()−>TimeRhsOperator ( f ,u ) ;

13 G e tM u l ti L ev e l So l ve r ()−>GetS olver ()−>TimeRhs ( 1 , f ) ;

14

15 t i m e i n f o . i t e r a t i o n ( i t e r ) ;

16 TimeInfoBroadcast ();

17

18 G e tM u l ti L ev e l So l ve r ()−>GetS olver ()−>TimeRhs ( 2 , f ) ;

19

20 // SOLVE

21 SolveTimePrimal(u, f );

22

23 [ . . . ]

24 }

25 }

Every time-step is treated as a stationary problem. The differential operator evaluated at

the old time step is known in advance and thus, Gascoigne automatically includes it when

assembling the right-hand side, if θ < 1. The function TimeInfoBroadcast() announces the

current time tand the time-step k.

The class StdTimeSolver is used instead of StdSolver. It has to be invoked by the op-

tion solver instat in block //Block MultiLevelSolver of the parameter file. Gascoigne

assembles the time-dependent matrix

1

θk M+A,

where Mis the mass matrix and Athe stiffness matrix which has to be specified as usual

in the class Equation.

1void StdTimeSolver : : AssembleMatrix ( const VectorInterface& gu, double d )

2{

3S td So lv er : : AssembleMatrix ( gu , d ) ;

4

5double s c a l e = d /( dt ∗t h e t a ) ;

6GetMatrix()−>AddM assWithD iff erent Ste ncil ( GetMassMatrix ( ) ,

7GetTimePattern ( ) , s c a l e ) ;

28

4.1 Time discretization of the heat equation

8}

First, the main part regarding a(·,·) is assembled. Then, the part regarding the time-

derivative is added with the proper scaling depending on the time-step kand parameter θ.

The parameters θand the time-step kare provided in the parameter file within the //Block

Loop

1// Block Loop

2scheme Theta

3t h et a 0 . 5

4dt 0 . 0 1

5neuler 0

By the option neuler we can tell Gascoigne to start with a number of neuler initial steps

using the implicit Euler scheme. In this way irregular initial data may be smoothened which

can be necessary for stability reasons. By setting neuler 0, no Euler steps are used and the

loop directly starts with the θ-time-stepping scheme.

The function GetTimePattern() tells Gascoigne, on which solution components the mass-

matrix will act. This TimePattern is specified in the Equation by a further function:

1void HeatEquation : : SetTimePattern ( TimePattern& TP) const

2{

3TP( 0 , 0 ) = 1 . ;

4}

TimePattern TP is a matrix of size ncomp times ncomp. Typically, this matrix is the diagonal

unit-matrix. See Section 4.3 for special choices of the TimePattern.

Time dependent problems usually are initial-boundary value problems and require some

initial data, like

u(x, 0) = u0(x),

at time t= 0. The initial data is specified via the function GetInitialConditionPointer()

in the class ProblemDescriptor. The initial condition can be set in the same way as the

right-hand side:

1class M y I n i t i a l : public DomainRightHandSide

2{

3public :

4std : : s t r i n g GetName ( ) const {return ” I n i t i a l Conditi on ” ; }

5int GetNcomp ( ) const {return 1 ; }

6

7double operator ( ) ( in t c , const Vertex2d& v ) const

8{

9return v . x () ∗v . y ( ) ;

10 }

11 };

12

29

4 Time dependent problems

13 [ . . . ]

14

15 class ProblemDescriptor : public ProblemDescriptorBase

16 {

17 public :

18 std : : s t r i n g GetName ( ) const {return ”Time−Dependent”;}

19 void BasicInit(const ParamFile∗p f )

20 {

21 [ . . . ]

22 GetInitialConditionPointer () = new M y I n i t i a l ;

23 ProblemDescriptorBase : : B a s i c I n i t ( pf ) ;

24 }

25 };

4.2 Time dependent problem data

The problem data like boundary data and right hand side can depend on time. Let us for

example consider the right hand side

f(x, t) = sin(πt)(1 −x2)(1 −y2).

All data classes like DomainRightHandSide,DirichletData or Equation are derived from

the class Application. We can access the current time and the time step via the functions

double GetTime() const and double GetTimeStep() const.

1class RHS : public DomainRightHandSide

2{

3public :

4int GetNcomp ( ) const {return 1 ; }

5s t r i n g GetName ( ) const {return ”RHS” ; }

6

7double operator ( ) ( in t c , const Vertex2d& v ) const

8{

9double t = GetTime ( ) ;

10 return s i n ( M PI ∗t ) ∗( 1 . 0 −v . x ( )∗v . x ( ) ) ∗(1.0 −v . y () ∗v . y ( ) ) ;

11 }

12 };

4.3 Non standard time-discretizations

As example for a non-standard time-depending problem we consider the following system of

partial differential equations

∂tu1+∂tu2−∆u1=f1

∂tu2+u1−∆u2=f2

30

4.3 Non standard time-discretizations

For the implementation of this system, Matrix,Form and SetTimePattern must be given

as:

1void EQ : : Form ( . . . )

2{

3b [ 0 ] += U [ 0 ] . x ( ) ∗N. x ( ) + U [ 0 ] . y ( ) ∗N. y ( ) ;

4

5b [ 1 ] += U [ 0 ] .m( ) ∗N.m( ) ;

6b [ 1 ] += U [ 1 ] . x ( ) ∗N. x ( ) + U [ 1 ] . y ( ) ∗N. y ( ) ;

7}

8

9void EQ : : Matrix ( . . . )

10 {

11 A( 0 , 0 ) += M. x ( ) ∗N. x ( ) + M. y ( ) ∗N. y ( ) ;

12

13 A( 1 , 0 ) += M.m( ) ∗N.m( ) ;

14 A( 1 , 1 ) += M. x ( ) ∗N. x ( ) + M. y ( ) ∗N. y ( ) ;

15 }

16

17 void EQ : : S etTimePat tern ( TimePattern& TP) const

18 {

19 TP( 0 , 0 ) = 1 . ;

20 TP( 0 , 1 ) = 1 . ;

21 TP( 1 , 1 ) = 1 . ;

22 }

31

4 Time dependent problems

32

5 Mesh handling

Mesh handling includes definition of the domain Ω, refinement of meshes, managing of curved

boundaries, input and output of meshes.

The mesh to be used in the calculation is indicated in the parameter file. In the Section

//Block Mesh the mesh to be read is specified:

1// Block Mesh

2

3dimension 2

4gridname mesh . inp

5prerefine 3

Gascoigne distinguishes two types of meshes: inp-files are coarse meshes describing the do-

main Ω. They are user-specified and the file format is explained in Section 5.1. These meshes

are usually as coarse as possible. Refinement is applied during run-time. The parameter

prerefine declares global refinements to be executed before the first solution cycle. It is

also possible to write out and read in locally refines meshes. These meshes are given in the

gup-format:

1// Block Mesh

2

3dimension 2

4gridname mesh . gup

5prerefine 0

During run-time all mesh handling like I/O and refinement is done by the class MeshAgent.

The MeshAgent is a member of the class Loop and is created in the Loop::BasicInit. The

MeshAgent takes care of the mesh refinement hierarchy HierarchicalMesh, it creates the

sequence of multigrid meshes GascoigneMultigridMesh and it finally provides the meshes

on every level used for the computation, the GascoigneMeshes. See the files with the same

names in File: GAS/src/Mesh/* .

If domains with curved boundaries are used, the shape definition is also done in the class

MeshAgent.

5.1 Definition of coarse meshes

Coarse meshes are specified in the inp-format. We indicate the nodes of the mesh, the

elements (that are quadrilaterals in two and hexahedrals in three dimensions) and finally

33

5 Mesh handling

the boundary elements. In two dimensions, these are the lines of quadrilaterals that touch

the boundary, in three dimensions we lists the quadrilaterals that are boundaries of the

domain and of elements. We never list interior lines (or quadrilaterals in three dimensions),

as these are not needed for the calculation. The following example is the coarse mesh for

the L-shaped domain:

811000

00-10

11-10

2-100

3000

4100

5-110

6010

7110

00quad0143

10quad2365

20quad3476

0 4 line 0 1

1 4 line 1 4

2 4 line 4 7

3 4 line 0 3

4 2 line 2 3

5 2 line 2 5

6 8 line 5 6

7 8 line 6 7

−1 10

1

0

−1 0 1

32 4

567

n1n0

n4 n5

n6n7

n2

n3

n3 n2

n0 n1

The file format is very simple: in the row, first the number of nodes in the coarse mesh is

given, followed by the sum of mesh elements plus elements on the boundary. The last three

values are always zero. This example has 8 nodes, 3 quads and 8 lines on the boundary.

After the header line all nodes of the mesh are declared in the format:

nxyz

where nis the index of the node (starting with zero) and xyzare the coordinates. We

always have to give all three coordinates (also in two dimensions).

The nodes are followed by the mesh elements. In two dimensions the format is:

n 0 quad n0 n1 n2 n3

where nis a consecutive index starting with 0. The second parameter is not used, thus

always 0. The mesh element, always quad is indicated by the third entry. The last four

values give the nodes describing the quadrilateral in counter-clockwise sense.

In three spatial dimension, the format is:

34

5.2 Curved boundaries

n 0 hex n0 n1 n2 n3 n4 n5 n6 n7

Here, first the nodes n0 to n3 in the front are given in counter-clock wise sense followed by

the nodes in the back.

Finally, all the boundary elements are given. In two dimensional meshes, these are all lines,

where both nodes are on the boundary. In three dimensions the boundary elements are all

quads with four nodes on the boundary. The format is:

n c line n0 n1

in two dimensions and

n c quad n0 n1 n2 n3

in three dimensions. nis again an iteration index. This index starts with 0. The value c

defines a color for the boundary element. These colors are necessary to distinguish between

different types of boundary conditions, see Sections 3.2.1 and 3.2.2. In two dimensions, the

order of the two nodes n0 n1 is arbitrary. In three dimensions however the four nodes have to

be given in a counter clock-wise sense seen from the outside of the domain! This is necessary

to get a unique definition of the normal vectors on the boundary. Thus, for the hex shown

in Picture above, the inp-file would be:

87000

0000

1100

2110

3010

4001

5101

6111

7011

00hex01234567

04quad0123

14quad1562

24quad3267

34quad4037

44quad5476

54quad4510

5.2 Curved boundaries

A part of the boundary can be assigned a curved boundary function. The idea is easy: Let

the boundary with color col be a curved boundary. We specify a function Rd→Rwhich

has the desired curved boundary as zero-iso-contour. For instance if the boundary Γ of the

domain should have the shape of a circle with radius rand midpoint (mx, my), one possible

function fcircle would be

fcircle(x, y) = (x−mx)2+ (y−my)2−r2.

Whenever an element with boundary nodes on a line with color col is refined, the newly gen-

erated nodes on the boundary line are then pulled onto the zero-iso-contour of the boundary

function. This is done by solving for the root of fcircle with Newton’s method.

35

5 Mesh handling

Adjust Curved BoundaryRefine

In Gascoigne boundary functions are described by the class BoundaryFunction<DIM> in

File: GAS/src/Mesh/BoundaryFunction . The circle is given by:

1class Circle : public BoundaryFunction<2>

2{

3double r , mx , my ;

4public :

5Circle(double r , double mx, double my) : r ( r ) , mx (mx) , my (my)

6{ }

7

8double operator ( ) ( const Vertex2d& v ) const

9{

10 return ( v . x()−mx ) ∗( v . x()−mx )+( v . y ()−my ) ∗( v . y()−my)−r∗r ;

11 }

12 }

The boundary function object has to be passed to the MeshAgent. For this, a derived mesh

agent class has to be specified and created in the Loop. The following minimal steps have to

be taken:

1class LocalMeshAgent : public MeshAgent

2{

3C i r c l e c i r c l e ;

4public :

5LocalMeshAgent ( ) : MeshAgent ( )

6{

7AddShape ( 4 , & c i r c l e ) ;

8}

9};

10

11 class LocalLoop : public StdLoop

12 {

13 public :

14

15 void BasicInit(const ParamFile∗p a r a m f i l e ,

16 const ProblemContainer∗PC,

17 const FunctionalContainer∗FC=NULL)

18 {

36

5.3 Adaptive refinement of meshes

19 GetMeshAgentPointer () = new LocalMeshAgent ( ) ;

20 StdLoop : : B a s i c I n i t ( p aramf i l e , PC, FC ) ;

21 }

22 };

In this example, the boundary with color 4is used as curved boundary.

5.3 Adaptive refinement of meshes

For local refinement, the MeshAgent has the function refine nodes(nodes). The nvector<int>

nodes contains indices of mesh nodes to be refined. If a node is to be refined, all elements

having this node as a corner node will be refined.

File: ../src/myloop.cc

1for ( i t e r =1; i t e r <n i t e r ;++ i t e r )

2{

3. . .

4Sol v e (u , f ) ;

5. . .

6// Estimate Error

7nvector<int>r e f ;

8// p i c k nodes f o r r e fi n em e nt i n r e f

9GetMeshAgent()−>r e f i n e n o d e s ( r e f ) ;

10 }

5.3.1 Estimating the Energy-Error

The class DwrQ1Q2 contains a function for estimation of an energy-error. It can be used in

the following way:

File: ../src/myloop.cc

1#include ”dwrq1q2 . h”

2. . .

3for ( i t e r =1; i t e r <n i t e r ;++ i t e r )

4{

5. . .

6Sol v e (u , f ) ;

7. . .

8// Estimate Error

9DoubleVector eta ;

10 DwrQ1Q2 dwr ( ∗Ge t M ul t iL e v el S o lv e r ()−>Get Solve r ( ) ) ;

11 dwr . EstimatorEnergy ( eta , f , u ) ;

12 . . .

13 }

37

5 Mesh handling

Now, for each node of the mesh the vector ηcontains a measure for the local contribution

to the total energy error.

5.3.2 Picking Nodes for Refinement

There are several methods for picking nodes for refinement. Common to all is to choose nodes

nwith large estimator values eta[n]. The class MalteAdaptor uses a global optimization

scheme for picking nodes:

File: ../src/myloop.cc

1#include ”malteadaptor .h ’ ’

2. . .

3for ( i t e r =1; i t e r <n i t e r ;++ i t e r )

4{

5. . .

6Sol v e (u , f ) ;

7. . .

8// Estimate Error

9DoubleVector eta ;

10 DwrQ1Q2 dwr ( G e t M ul t i Le v el S o lv e r ()−>Get Sol ver ( ) ) ;

11 dwr . EstimatorEnergy ( eta , f , u ) ;

12

13 // Pick Elements

14 nvector<int>refnodes ;

15 MalteAdaptor A( p a r a m f i l e , e t a ) ;

16 A. r e f i n e ( r e f n o d e s ) ;

17

18 // Re fin e

19 GetMeshAgent()−>r e f i n e n o d e s ( r e f n o d e s ) ;

20 }

38

6 Flow problems and stabilization

In order to discretize the Stokes or Navier-Stokes systems one can either use inf-sup- stable

finite element pairs for pressure and velocity which fulfill an inf-sup condition or include

stabilization terms in the equations that ensure uniqueness of the solution and stability.

Gascoigne always uses equal-order elements which are not inf-sup stable, in combination

with different pressure stabilization techniques. If in the Navier-Stokes case the convective

term is dominating (high Reynolds number) further stabilization is required.

6.1 The function point

To set element-wise constants in Gascoigne, the function point is needed:

1void Equation : :

2po in t ( double h , const FemFunction& U, const Vertex2d& v ) const

3{

4

5}

This function is called before every call of Matrix and Form and here the local mesh-size

hand the coordinates of the node vare given as a parameter. That way one can specify

certain criteria for the equation based on location or local mesh width.

To stabilize a given equation or system of equations (like Stokes and Navier-Stokes) one

usually modifies it by adding certain stabilizing terms. These terms often depend on the

local mesh width and have therefore to be defined in this function.

Since the function point is defined as const no class members can be altered. Thus, every

variable which shall be modified here has to be defined as mutable in the header file:

1mutable double variable ;

6.2 Stabilization of the Stokes system

The Stokes equations are given by

∇ · v= 0,−ν∆v+∇p=fin Ω.

39

6 Flow problems and stabilization

The simplest stabilization technique is the addition of an artificial diffusion. Stability is

ensured by modifying the divergence equation to

(∇ · v, ξ) + X

K∈Ωh

αK(∇p, ∇ξ)K,

where αKis an element-wise constant that can be chosen as

αK=α0

h2

K

6ν.

To get the local cell size hKof an element K, we use the function point which is called

before each call of Form and Matrix. In point we can calculate the value αKand save it to

a local variable:

1void St o kes : :

2po in t ( double h , const FemFunction& U, const Vertex2d& v ) const

3{

4al ph aK = a l p h a 0 ∗h∗h / 6. 0/ n u ;

5}

The parameter alpha0 can be read from the parameter file. The variable alphaK has to

be defined as mutable in the header file:

1mutable double alphaK ;

6.3 Stabilization of convective terms

The Navier-Stokes system

∇ · v= 0, v · ∇v−ν∆v+∇p=fin Ω,

needs additional stabilization of the convective term if the viscosity is low. One popular

approach consists of adding diffusion in streamline direction (Streamline Upwind Petrov-

Galerkin (SUPG) method). In combination with the artificial diffusion ansatz for pressure

stabilization, the full stabilized system reads

(∇ · v, ξ) + X

K

(αK∇p, ∇ξ)K= 0,

(v· ∇v, φ) + ν(∇v, ∇φ)−(p, ·∇φ) + X

K

(δKv· ∇v, v · ∇φ)K= (f, φ).

The parameters αKand δKare set to

αK=α06ν

h2

K

+kvkK,∞

hK−1

, δK=δ06ν

h2

K

+kvkK,∞

hK−1

,

40

6.4 LPS Stabilization

where kvkK,∞is the norm of the local velocity vector. The parameters might be chosen as

α0=δ0≈1

2.

For programming the Matrix it is usually not necessary to implement all derivatives. It

should be sufficient to have:

1void Equation : : Matrix ( . . . )

2{

3

4. . .

5// d e r i v a t i v e o f t h e c o nv ec ti on s t a b i l i z a t i o n terms

6A( 1 , 1 ) += d e l t a ∗(U [ 1 ] .m( ) ∗M. x ( ) + U [ 2 ] .m( ) ∗M. y ( ) )

7∗(U [ 1 ] .m( ) ∗N. x ( ) + U [ 2 ] .m( ) ∗N. y ( ) ) ;

8A( 2 , 2 ) += d e l t a ∗(U [ 1 ] .m( ) ∗M. x ( ) + U [ 2 ] .m( ) ∗M. y ( ) )

9∗(U [ 1 ] .m( ) ∗N. x ( ) + U [ 2 ] .m( ) ∗N. y ( ) ) ;

10 . . .

11 }

6.4 LPS Stabilization

Another way to stabilize these equations is the so called Local Projection Stabilization (LPS).

The idea is to suppress local fluctuations by adding the difference between the function itself

and its projection onto a coarser grid.

We use the fluctuation operator

π=id −i2h

where i2his a projection to the coarser grid. The LPS stabilized formulation for the Stokes

system reads

(∇ · v, ξ) + X

K∈Ωh

αK(∇πp, ∇πξ)K= 0

and for the Navier-Stokes systems

(∇ · v, ξ) + X

K

(αK∇πp, ∇πξ)K= 0,

(v· ∇v, φ) + ν(∇v, ∇φ)−(p, ·∇φ) + X

K

(δKv· ∇πv, v · ∇πφ)K= (f, φ).

αKand δKcan be taken from above.

To implement the projection in Gascoigne, one has to use the class LpsEquation instead of

Equation. It provides the additional functions

41

6 Flow problems and stabilization

1void lpspoint(double h , const FemFunction& U,

2const Vertex2d& v ) const {. . . }

3void StabForm ( V e c t o r I t e r a t o r b , const FemFunction& U,

4const FemFunction& UP, const TestFunction& Np)

5const {. . . }

6void StabMatrix ( EntryMatrix& A, const FemFunction& U,

7const TestFunction& Np, const TestFunction& Mp)

8const {. . . }

which have to be used for the stabilizing terms.

The variables UP,Np and Mp stand for the fluctuation operator πapplied to Nand M.Gas-

coigne provides these fluctuation operators if the discretization specified by the keyword

discname in the parameterfile is chosen Q1 or Q2Lps.

Similar as in the case of the SUPG method, it should be sufficient to implement the following

terms in the StabMatrix

1void LpsEquation : : StabMatrix ( . . . )

2{

3

4. . .

5// d e r i v a t i v e o f t h e c o nv ec ti on s t a b i l i z a t i o n terms

6A( 1 , 1 ) += d e l t a ∗(U [ 1 ] .m( ) ∗Mp. x ( ) + U [ 2 ] .m( ) ∗Mp. y ( ) )

7∗(U [ 1 ] .m( ) ∗Np. x ( ) + U [ 2 ] .m( ) ∗Np . y ( ) ) ;

8A( 2 , 2 ) += d e l t a ∗(U [ 1 ] .m( ) ∗Mp. x ( ) + U [ 2 ] .m( ) ∗Mp. y ( ) )

9∗(U [ 1 ] .m( ) ∗Np. x ( ) + U [ 2 ] .m( ) ∗Np . y ( ) ) ;

10 . . .

11 }

42