complicated partial differential equations
Proefschrift
ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,
op gezag van de Rector Magnificus prof. ir. K.C.A.M. Luyben, voorzitter van het College voor Promoties,
in het openbaar te verdedigen op dinsdag 3 december 2013 om 12.30 uur door
Kristian Breum ØLGAARD
Master of Science in Civil Engineering, Aalborg Universitet Esbjerg geboren te Ringkøbing, Denemarken
Copromotor: Dr. G. N. Wells
Samenstelling promotiecommissie: Rector Magnificus Voorzitter
Prof. dr. ir. L. J. Sluys Technische Universiteit Delft, promotor Dr. G. N. Wells University of Cambridge, copromotor Dr. ir. M. B. van Gijzen Technische Universiteit Delft
Prof. dr. P. H. J. Kelly Imperial College London
Prof. dr. R. Larsson Chalmers University of Technology Prof. dr. L. R. Scott University of Chicago
Prof. dr. ir. C. Vuik Technische Universiteit Delft
Prof. dr. A. Scarpas Technische Universiteit Delft, reservelid
Copyright © 2013 by K. B. Ølgaard
Printed by Ipskamp Drukkers B.V., Enschede, The Netherlands ISBN 978-94-6191-990-8
Foreword
This thesis represents the formal end of my long and interesting journey as a PhD student. The sum of many experiences over the past years has increased my knowledge and contributed to my personal development. All these experiences originate from the interaction with many people to whom I would like to express my gratitude.
I am most grateful to Garth Wells for giving me the opportunity to come to Delft and to study under his competent supervision. His constructive criticism and vision combined with our nice discussions greatly improved the quality of my research. As the head of the computational mechanics group, Bert Sluys has played a vital role by creating a very nice and supportive working environment where people enjoy a lot of creative freedom. As creativity is key in this research I consider myself lucky to have been part of Bert’s group.
Ronnie Pedersen did a very good job in persuading me to come to Delft for a PhD, and I am happy that he managed to convince me. I am also grateful for enjoying his friendship throughout the years, the good times on the football pitch, and the even better times in ’t Proeflokaal watching football and discussing work and life in general.
A friendly and inspiring working environment is important in order to produce quality work. Therefore, I would like to thank past and present colleagues Rafid Al-Khoury, Roberta Bellodi, Frank Custers, Frank Everdij, Huan He, Cecilia Iacono, Cor Kasbergen, Oriol Lloberas-Valls, Prithvi Mandapalli, Frans van der Meer, Andrei Metrikine, Peter Moonen, Dung Nguyen, Vinh Phu Nguyen, Mehdi Nikbakth, Marjon van der Perk, Frank Radtke, Zahid Shabir, Xuming Shan, Angelo Simone, Mojtaba Talebian, Andy Terrel, Ilse Vegt, Jaap Weerheijm, Sigurd Blöndal, Lars Damkilde, Niels Dollerup, Jens Hagelskjær, Michael Jepsen, Sven Krabbenhøft and Søren Lambertsen. In particular, I would like to thank Frans for the years that we shared the same office and for translating the propositions into Dutch. A special thanks goes to Mehdi, my ‘brother-in-arms’, the only person remaining in the group who was also involved with the FEniCS Project after Garth left for Cambridge and Xuming left for home.
therefore, I would also like to thank all the people in the FEniCS community, in particular my close collaborators from Simula Anders Logg, Martin Alnæs, Marie Rognes and Johan Hake for all the nice discussions, debugging assistance and good ideas. During my PhD, I also had the pleasure of visiting the University of Michigan and in this regard I want to thank Krishna Garikipati, Jake Ostien and his wife Erin for their hospitality during my stay in Ann Arbor.
Outside the office, I enjoyed many hours in the good company of my friends Linda Grimstrup and Lars Freising which definitely improved the quality of my social life a lot. I also want to thank all my former team mates at Vitesse Delft for the many memorable hours on the football pitch trying to learn the secrets behind ‘totaalvoetbal’. Although The Netherlands and Denmark are quite similar in terms of weather, nature and culture it was always nice to receive visitors from home. For this, I would like to thank my friends Kenneth Guldager, Henrik Hansen, Mads Madsen, Christian Meyer, Nick Nørreby and Thomas Sørensen.
Last, but certainly not least I want to thank my parents and my brother and sisters for their encouragement, support, help and visits during my years in Delft. I also wish to thank both of my sons for putting things in perspective which helped me to focus during the last iterations towards finishing this thesis. Of all people, I am most grateful to my wife. I know her patience has been tested to the limit, yet she remained supportive, loving and caring during all the years. For this, and for our sons, I am forever indebted.
The research presented in this thesis was carried out at the Faculty of Civil Engineering and Geosciences at Delft University of Technology. The research was supported by the Netherlands Technology Foundation STW, the Netherlands Organisation for Scientific Research and the Ministry of Public Works and Water Management.
Kristian Breum Ølgaard
Contents
1 Introduction 1
1.1 Research objectives and approach . . . 2
1.2 Outline . . . 4
1.3 The FEniCS Project . . . 5
1.3.1 Simple model problem . . . 6
1.3.2 Unified Form Language . . . 8
1.3.3 FEniCS Form Compiler . . . 11
1.3.4 Unified Form-assembly Code . . . 15
1.3.5 DOLFIN . . . 18
2 FEniCS applications to solid mechanics 29 2.1 Governing equations . . . 30
2.1.1 Preliminaries . . . 30
2.1.2 Balance of momentum . . . 31
2.1.3 Potential energy minimisation . . . 32
2.2 Constitutive models . . . 33
2.2.1 Linearised elasticity . . . 33
2.2.2 Flow theory of plasticity . . . 34
2.2.3 Hyperelasticity . . . 35
2.3 Linearisation issues for complex constitutive models . . . 36
2.3.1 Consistency of linearisation . . . 36
2.3.2 Quadrature elements . . . 38
2.4 Implementations and examples . . . 41
2.4.1 Linearised elasticity . . . 43
2.4.2 Plasticity . . . 43
2.4.3 Hyperelasticity . . . 49
2.4.4 Elastodynamics . . . 52
3 Representations and optimisations of finite element variational forms 57
3.1 Motivation and approach . . . 58
3.2 Representation of finite element tensors . . . 60
3.2.1 Quadrature representation . . . 61
3.2.2 Tensor contraction representation . . . 62
3.3 Quadrature optimisations . . . 64
3.3.1 Eliminate operations on zeros . . . 67
3.3.2 Simplify expressions . . . 69
3.3.3 Precompute integration point constants . . . 71
3.3.4 Precompute basis constants . . . 72
3.3.5 Further optimisations . . . 73
3.4 Performance comparisons of representations . . . 75
3.4.1 Performance for a selection of forms . . . 75
3.4.2 Performance for common, simple forms . . . 80
3.4.3 Performance for forms of increasing complexity . . . 82
3.5 Performance comparisons of quadrature optimisations . . . 86
3.6 Automatic selection of representation . . . 92
3.7 Future optimisations . . . 93
4 Automation of discontinuous Galerkin methods 97 4.1 Extending the framework to discontinuous Galerkin methods . . . . 98
4.1.1 Extending the Unified Form Language . . . 99
4.1.2 Extending the Unified Form-assembly Code . . . 100
4.1.3 Extending the FEniCS Form Compiler . . . 100
4.1.4 Extending DOLFIN . . . 101
4.2 Examples . . . 103
4.2.1 The Poisson equation . . . 104
4.2.2 Steady state advection–diffusion equation . . . 105
4.2.3 The Stokes equations . . . 109
4.2.4 Biharmonic equation . . . 110
4.2.5 Further applications . . . 114
5 Automation of lifting-type discontinuous Galerkin methods 115 5.1 Lifting-type formulation for the Poisson equation . . . 116
5.2 Semi-automated implementation of lifting-type formulations . . . . 117
5.3 Comparison of IP and lifting-type formulations . . . 122
5.4 Future developments . . . 127
6 Strain gradient plasticity 129 6.1 A strain gradient plasticity model . . . 130
6.2 A discontinuous Galerkin formulation for the plastic multiplier . . . 133
6.4 Implementation . . . 137
6.4.1 The predictor step . . . 138
6.4.2 The corrector step . . . 139
6.4.3 Implementing the variational forms . . . 141
6.5 Numerical examples . . . 141
6.5.1 Unit square loaded in shear with strain softening . . . 143
6.5.2 Plate under compressive loading with strain softening . . . . 146
6.5.3 Plate under compressive loading with strain hardening . . . 156
6.5.4 Micro-indentation . . . 160
6.5.5 Computational notes . . . 162
7 Conclusions and future developments 167
References 171 Summary 183 Samenvatting 185 Propositions 187 Stellingen 189 Curriculum vitae 191
1 Introduction
Since the advent of the modern programmable computer in the 1940s, the cost of computing power relative to manpower has decreased significantly. As a con-sequence, high-level programming languages have emerged allowing the imple-mentation of programs in source code using abstractions that are independent of the specific computer architectures on which the program is intended to run. A compiler is then invoked to translate the source code into machine code targeted for the given computer’s central processing unit (CPU). This development has allowed, among other things, researchers and scientists to write programs for investigating and solving various classes of problems numerically.
In engineering, physical phenomena are often described mathematically by partial differential equations (PDEs), and a commonly used method to solve these equations is the finite element method (FEM). Standard finite element software typ-ically provide a problem solving environment for a set of engineering problems using a predefined selection of finite elements. As part of the application program-ming interface (API) a user can often supply subroutines which implement special methods, for instance, the constitutive model in case of a solid mechanics problem. This offers a degree of customisation and flexibility in terms of implementing certain models, but the approach may fall short as the complexity of a model increases.
Strain gradient plasticity is an example of a class of models which can be difficult to implement in traditional finite element software and researchers often resort to implementing their own unique solver targeting a specific model. An implementation involves translating the abstract mathematical representation of the model into source code which can be handled by a compiler, a process which can be tedious, time consuming and error prone. However, by introducing a higher level of abstraction, the burden of this process can be alleviated when it comes to implementing mathematical representations of the FEM for solving PDEs. A possible abstraction consists of a form language for expressing the mathematical formulation of the given problem, and compilers which automatically generate efficient source code from the given mathematical expressions. This thesis is centered around this type of automated mathematical modelling.
1.1
Research objectives and approach
The research presented in this thesis aims at developing concepts, tools and meth-ods which allow researchers and application developers to create efficient solvers for complicated partial differential equations with relatively little effort. Sev-eral software projects aim at providing a flexible framework for solving partial differential equations using the finite element method. These software projects include, among others, traditional finite element libraries and toolboxes such as deal.II (http://www.dealii.org/, Bangerth et al. (2007)), Diffpack (http://www. diffpack.com/, Langtangen (1999)), DUNE (http://www.dune-project.org,
Bas-tian et al. (2008b,a)), GetFEM++ (http://home.gna.org/getfem/), OpenFOAM
(http://www.openfoam.com/) and Cactus (http://cactuscode.org/, Allen et al.
(2000)). However, a bit of ‘hand coding’ is often needed in order to use the above mentioned software. For instance, a user must typically implement (parts of) the assembly algorithm which is cumbersome as the complexity of the problem is increasing. A number of software projects have, therefore, emerged that try to automate the finite element method. These projects include, among others, FINGER (Wang, 1986), Archimedes (Shewchuk and Ghattas, 1993), Symbolic Me-chanics System (Korelc, 1997), GetDP (http://geuz.org/getdp/, Dular et al. (1998)),
FreeFEM++ (http://www.freefem.org/), Sundance (http://www.math.ttu.edu/ ~kelong/Sundance/html/, Long et al. (2010)), Feel++ (http://www.feelpp.org/,
Prud’homme (2006)) and the FEniCS Project (http://fenicsproject.org, Logg
et al. (2012a)). A common feature of these approaches is that they provide a higher level of abstraction for expressing variational forms and thereby lessen the burden on application developers.
The developments presented in this thesis are implemented in various software components of the FEniCS Project which is chosen for a number of reasons. The software is released under an open source license1which makes it possible to obtain and modify the source code. This provides a high degree of freedom and flexibility in terms of implementing advanced models and applications. Furthermore, if the application source code is published, the implementation becomes completely transparent and reproducible, both properties of importance in research. The software contains a problem solving environment that handles the assembly, the application of boundary conditions and the solution of sparse systems of equations. What distinguishes the software from the more conventional finite element packages is that it provides a high degree of mathematical abstraction by implementing a form language for expressing variational forms and relies on form compilers to automatically generate computer code for the local finite element tensor. This approach offers several advantages of which two are of particular interest. Firstly, the time needed to implement, test and debug the code for the local finite element 1All FEniCS core components are licensed under the GNU LGPL version 3 (or any later version) as published by the Free Software Foundation (http://www.fsf.org).
tensor can be reduced. Secondly, various optimisations can be employed by the form compilers during the code generation stage to make the generated code competitive with hand optimised code. The importance of these two advantages is proportional to the complexity of the variational form. Finally, the software is under active development by a growing community which is helpful for receiving feedback when implementing new features and during debugging sessions.
The potential of the FEniCS framework is evident, however, at the time when this work was commenced, the functionality in the FEniCS software was only available for a limited class of problems. For instance, only integration over element interiors was supported which precluded, among other things, discontinuous Galerkin methods from being handled as these methods involve integration over element boundaries. Furthermore, problems like conventional plasticity were not possible to solve because the software could only handle functions coming from a finite element space. Also, the generated code was only efficient for limited classes of problems. The objectives of this work can thus be condensed into the following: extend the automated mathematical modelling framework of FEniCS such that
• discontinuous Galerkin methods can be handled;
• rapid prototyping of advanced models and applications is possible; and • efficiency is maintained also for complex problems in general.
As will be demonstrated in this work, addressing the above three issues has had a significant impact on the range of problems which can be handled in the FEniCS framework and thereby making life easier for researchers and application developers.
A complex application from solid mechanics in the form of a strain gradient plasticity model is considered, as an example, to demonstrate the extensions to the FEniCS framework developed in this work. Strain gradient models are often used to provide regularisation in softening problems and to account for observed size effects at small length scales. An abundance of strain gradient models have been proposed in literature including the models by Aifantis (1984), Gurtin (2004), Fleck and Hutchinson (1997), Fleck and Hutchinson (2001) and Gao et al. (1999) to name a few. The focus in this work is on the class of models involving gradients of fields such as the equivalent plastic strain. An example of such a model is that proposed by Aifantis (1984) which involves the addition of the Laplacian of the equivalent plastic strain to the classical yield condition. A feature of this particular model is that the classical consistency condition leads to a partial differential equation rather than an algebraic equation, as is the case is classical flow theory of plasticity. The partial differential equation is only active in the region undergoing plastic deformations which introduces the difficulty of imposing non-standard boundary conditions on the secondary field on the evolving boundary.
Motivated by the work of Wells et al. (2004) and Molari et al. (2006) who used a discontinuous Galerkin formulation for a strain gradient-dependent damage model, a discontinuous basis can be used to interpolate the secondary field. This provides a natural framework for handling evolving elastic–plastic boundaries and provides local (cell-wise) satisfaction of the yield condition. To satisfy the regularity requirement of the secondary field, a discontinuous Galerkin formulation is used to enforce weak continuity across cell facets. In order to allow the use of a discontinuous constant basis for the secondary field, a so-called lifting-type discontinuous Galerkin formulation, proposed by Bassi and Rebay (1997, 2002), is adopted. A discontinuous constant basis is the natural choice for the secondary field when a linear continuous basis is used for the displacement field. Considering that the formulation involves an additional field variable it is also computationally more efficient if discontinuous constant elements can be used for this particular field.
1.2
Outline
The rest of this chapter contains an overview of the FEniCS Project including details on the components pertinent to the present work. Chapter 2 continues with a demonstration of how to use the FEniCS toolchain for solid mechanics applications. The purpose of this demonstration is twofold. Firstly, it serves as an introduction to the concepts of automated modelling from a solid mechanics point of view, which will give an understanding of how the automated modelling approach can be utilised to also tackle more complex problems. Secondly, the presented models and applications will be used in subsequent chapters, either by extending the models or by using them as a platform for discussing the development of FEniCS components in connection to the work presented in this thesis.
Local finite element tensors can be evaluated using different representations of the tensors. In Chapter 3 the two representations that FFC adopts, the quadrature representation and the tensor contraction representation are presented and comparisons are made between the two representations. Furthermore, optimisation strategies for the quadrature representation are discussed and the performance of these are investigated.
Chapter 4 introduces the extensions implemented in the FEniCS framework to allow a class of discontinuous Galerkin (DG) formulations to be handled in an automated fashion. Building on these abstractions, a semi-automated approach to implementing lifting-type DG formulations is presented in Chapter 5. This chapter also contains a brief comparison, in terms of complexity regarding the implementation and the numerical implications, between a lifting-type formulation and an interior penalty (IP) DG formulation for the Poisson equation.
framework are brought together in an implementation of a lifting-type discontinu-ous formulation for a simple strain gradient plasticity model proposed by Aifantis (Aifantis, 1984). The purpose is to illustrate how researchers and application de-velopers may create solvers for more complex problems on top of the FEniCS software. Finally, in Chapter 7, conclusions are drawn and recommendations for future development related to this work are presented.
1.3
The FEniCS Project
The FEniCS Project is a suite of open source programs for automating the solution of PDEs. The concepts and components which are most important in relation to this work and which will be elaborated on in subsequent chapters are presented. Thus, only a subset of the components in the FEniCS Project is presented. Further details on the components presented here, and other components associated with the FEniCS Project, can be found in the FEniCS book (Logg et al., 2012a) or online athttp://fenicsproject.org. All FEniCS software components, and the
software developed in this work, can be obtained freely athttps://bitbucket.org/ fenics-project2. The FEniCS Project is under continuous development, however,
this presentation and all example code, and software developed and described in this work, is compliant with version 1.0 of the project and its associated components unless stated otherwise.
The majority of developments in this work is implemented in the core com-ponents of FEniCS. However, some of the developments are implemented in the FEniCS Solid Mechanics library3 (Ølgaard and Wells, 2013). In this thesis, several code snippets are presented along with many results from numerical experiments. All example code, and the code which has been used to obtain all the results, can be downloaded fromhttps://bitbucket.org/k.b.oelgaard/ oelgaard-thesis-supporting-material. Note that in order to run the code,
work-ing installations of FEniCS version 1.0 and FEniCS Solid Mechanics version 1.0 are required4.
The procedure of solving PDEs using the FEM can be broken down into the following four steps:
1. Formulate the variational problem of the PDE
2The FEniCS software components have recently moved from Launchpad (https://launchpad.net/
fenics) to Bitbucket. However, as the FEniCS Project is being actively developed the location might
change again in the future. The FEniCS website (http://fenicsproject.org), which is less likely to
move, might be a better starting point for locating the software.
3The FEniCS Solid Mechanics library was formerly known as FEniCS Plasticity (https://launchpad.
net/fenics-plasticity) which focussed solely on plasticity problems. However, to reflect that the
scope of the library has increased to also include more general solid mechanics problems the name was changed during a recent migration from Launchpad to Bitbucket.
4Version 1.0 of the FEniCS Solid Mechanics library can be downloaded fromhttps://bitbucket.
UFC
UFL
FFC
DOLFIN
Figure 1.1: FEniCS toolchain for solving a PDE using the FEM. 2. Discretise the formulation
3. Finite element assembly
4. Solve the global system of equations
Facilities for each of these steps are implemented in separate software components in FEniCS. The relationship between input and output of each component in the FEniCS toolchain for the finite element procedure is shown in Figure 1.1. In short, the variational form of the PDE is expressed in the Unified Form Language (UFL) (Alnæs et al., 2013; Alnæs, 2012), which is given as input to the FEniCS Form Compiler (FFC)5 (Kirby and Logg, 2006, 2007; Logg et al., 2012c; Ølgaard et al., 2008a) that automatically generates efficient C++ code for evaluating the local element tensors. The output from FFC is compliant with the interface defined in Unified Form-assembly Code (UFC) (Alnæs et al., 2009, 2012) and is used by DOLFIN (Logg and Wells, 2010; Logg et al., 2012d), which is the finite element assembler and solver of FEniCS although, in principle, any assembly library which supports UFC can be used.
The key advantage of this modular construction is that it becomes more trans-parent where and how new features and functionality should be implemented. Furthermore, developers and users can pick individual components to form their own applications. In this work for instance, the UFL is augmented with discontinu-ous Galerkin operators6, compiler optimisations are implemented in FFC, while more complex solvers for lifting-type formulations and solid mechanics problems can be implemented on top of the FEniCS toolbox.
1.3.1 Simple model problem
As a model boundary value problem for presenting the FEniCS framework consider the Poisson equation, which for a bodyΩ⊂Rd, where 1≤d≤3, with boundary 5Any compiler that supports UFL as input, and outputs UFC code, can be used instead of FFC in the described toolchain. The Symbolic Form Compiler (Alnæs and Mardal, 2010, 2012), which is also part of FEniCS, is one such example.
6Historically, the DG operators were implemented in the original form language of FFC which was later merged into the richer UFL.
∂Ω and outward unit normal vector n : ∂Ω→Rdreads: −∆u= f inΩ,
u=g onΓD, ∇u·n=h onΓN.
(1.1) Here, u is an unknown scalar field, f is a source term, g is a prescribed value for u on the Dirichlet boundaryΓD, and h is a prescribed value for the outward normal derivative of u on the Neumann boundaryΓN. The boundariesΓD andΓN divide the boundary such thatΓD⊆∂Ω and ΓN=∂Ω\ΓD. To apply the FEniCS framework the problem must be posed as a variational formulation in the following canonical form: find u∈V such that
a(u, v) =L(v) ∀v∈ ˆV, (1.2)
where V is the trial space and ˆV is the test space, a(u, v)and L(v) denote the bilinear and linear forms, respectively. A typical variational form7of (1.1) defines the bilinear and linear forms as:
a(u, v):= Z Ω∇u· ∇v dx (1.3) L(v):= Z Ω f v dx+ Z ΓNhv ds, (1.4)
with the trial and test spaces defined as:
V :=nv∈H1(Ω): v=g onΓDo , (1.5) ˆV :=nv∈H1(Ω): v=0 onΓDo . (1.6) The variational problem in (1.2) must be discretised to compute a finite element solution to the Poisson problem. This is done by using a pair of discrete function spaces for the test and trial functions: find uh∈Vh⊂V such that
a(uh, v) =L(v) ∀v∈ ˆVh⊂ ˆV. (1.7) Thus, after transforming the strong form of the problem into the variational coun-terpart, the FEniCS toolchain, starting with UFL, can be invoked to compute a solution.
UFL code
element = FiniteElement("Lagrange", triangle, 1) u = TrialFunction(element) v = TestFunction(element) f = Coefficient(element) h = Coefficient(element) a = inner(grad(u), grad(v))*dx L = f*v*dx + h*v*ds
Figure 1.2: UFL code for the Poisson problem using continuous-piecewise linear Lagrange polynomials on triangles.
1.3.2 Unified Form Language
In order to compute a solution to the variational problem using the FEM, it is neces-sary to discretise the formulation. The Unified Form Language (UFL) (Alnæs et al., 2013; Alnæs, 2012) enables a user to express the discretisation compactly using a notation which resembles the mathematical notation closely. UFL is implemented as a domain-specific embedded language (DSEL) in Python which, among other things, allow users to define custom operators using all features of the Python programming language when writing UFL code. This section presents the most basic features used throughout in this work, while some of the more advanced func-tionality is presented in subsequent chapters as needed. For a detailed description of the language, refer to Alnæs et al. (2013).
The Poisson problem in (1.7) can be expressed in UFL by the code shown in Figure 1.2. The first line in the code defines the local finite element basis that spans the discrete function space Vhon an element T∈ ThwhereThdenotes the standard triangulation ofΩ. Generally, finite elements are defined in UFL by theirfamily, cellanddegree:
UFL code
element = FiniteElement(family, cell, degree)
which in the given case, in Figure 1.2, means that the basis is a piecewise continuous linear Lagrange triangle. UFL contains a set of predefined finite elementfamily
names, for instance, "Lagrange" as already shown, "Discontinuous Lagrange"
(short name"DG") and"Brezzi-Douglas-Marini"(short name "BDM"). Thecell
argument denotes the polygonal shape of the finite element while the degree
argument denotes the degree of the polynomial space. Although valid cell shapes in UFL are: interval,triangle,tetrahedron,quadrilateralandhexahedron. FFC
Mathematical notation UFL notation A·B dot(A, B) A : B inner(A, B) AB≡A⊗B outer(A, B) AT transpose(A), A.T sym A sym(A) tr A≡Aii tr(A) det A det(A)
Mathematical notation UFL notation
A,i A.dx(i) ∂A ∂xi Dx(A, i) dA dn Dn(A) ∇A grad(A) ∇ ·A div(A)
Table 1.1: (Left) Table of tensor algebraic operators. (Right) Table of differential operators.
value of celland degree depend on the choice of finite elementfamily. It is
important to realise that UFL is only concerned with the abstract operations related to the finite element function spaces; it is left to the form compiler to support the element families, that is, to generate meaningful code for the representation of elements and forms. For mixed finite element methods, product spaces like:
V= [V2×V2]×V1. (1.8)
can easily be generated by either theMixedElementclass or the*operator:
UFL code
V_2 = FiniteElement("Lagrange", triangle, 2) V_1 = FiniteElement("Lagrange", triangle, 1) V = (V_2*V_2)*V_1
W = MixedElement(MixedElement((V_2, V_2)), V_2)
meaning thatVandWare identical. To create a mixed element in which all the
component spaces are identical, theVectorElementcan be used:
UFL code
V = VectorElement(family, cell, degree, dim=None)
wheredimdefaults to the dimension of the given cell unless explicitly specified.
After defining the local finite element basis, the trial function u ∈ Vh, the test function v ∈ Vh and the coefficient functions f , g ∈ Vh can be defined in a straightforward fashion as seen in the code in Figure 1.2. The bilinear and linear forms from (1.3) and (1.4) can then be implemented simply by using the tensor and differential operators defined in UFL, some of which can be seen in Table 1.1. An important thing to note is that the definition of the gradient operator
Mathematical notation UFL notation a b a / b ab a**b,pow(a,b) p f sqrt(f) exp f exp(f) ln f ln(f) |f| abs(f) sign f sign(f)
Mathematical notation UFL notation
cos f cos(f) sin f sin(f) tan f tan(f) arccos f acos(f) arcsin f asin(f) arctan f atan(f)
Table 1.2: (Left) Table of elementary functions. (Right) Table of trigonometric functions.
and not {grad(u)}ij = ∂uj/∂xi. The latter operator is, however, provided in UFL bynabla_grad(A). A similar convention applies to the divergence operator
where nabla_div(A) is provided as an alternative to div(A). In this work, the
operators∇u and∇ ·u follow the UFL definition for the gradient and divergence operators,grad(u)anddiv(u), respectively and should not be confused with the
UFL operatorsnabla_grad(u)andnabla_div(u).
To complete the implementation of the variational forms, integration on the relevant domains must be expressed. In UFL, the integral over the domainR
ΩkI dx is denoted byI*dx(k)while the integral over the exterior boundary of the domain
R
∂ΩkI ds is denoted byI*ds(k)where k is the subdomain number and I is a valid UFL expression. Thus, having completed the implementation of (1.2) in the near-mathematical notation of UFL, the form compiler can be invoked to generate code from the abstract UFL representation.
The last two classes of expressions to be presented in this short introduction to UFL are nonlinear scalar functions and geometric quantities. UFL provides a set of nonlinear scalar functions, presented in Table 1.2, which can be applied to, for instance, scalar valued coefficient functions such as f and g in the Poisson example. It is illegal to apply these functions to any test or trial function as this would render the variational form nonlinear in those arguments. Geometric quantities are related to the local finite element cell T. For instance, the coordinate of the integration point currently being evaluated on T (including its boundary) can be accessed via
cell.x. Other geometric quantities which are particularly useful in relation to this
work are the outward normal to the facet8currently being evaluatedcell.nand the circumradius, the radius of the circumscribed circle of T,cell.circumradius. Basic
usage of nonlinear functions and the integration point coordinate is demonstrated later in Section 1.3.5, while the facet normal and circumradius are frequently used 8A facet is a topological entity of a computational mesh of dimension D−1 (codimension 1) where D is the topological dimension of the cells of the computational mesh. Thus for a triangular mesh, the facets are the edges and for a tetrahedral mesh, the facets are the faces.
when defining discontinuous Galerkin variational forms in Chapter 4. 1.3.3 FEniCS Form Compiler
As shown in Figure 1.1, the FEniCS Form Compiler (FFC) (Kirby and Logg, 2006, 2007; Logg et al., 2012c; Ølgaard et al., 2008a; Ølgaard and Wells, 2010) takes as input a variational form specified in UFL and generates as output C++ code which conforms to the UFC interface, to be described in Section 1.3.4. Central to the finite element method is the assembly of sparse tensors, described in Section 1.3.5, which relies on the computation of the local element tensor AT as well as the local-to-global mapping ιT. Although it is possible to hand code AT and ιT, the process is both tedious and error-prone especially for complex problems. This issue is eliminated by letting FFC generate the code automatically. Introducing a compiler also provides the possibility of applying various optimisation strategies for efficient computation of AT, which would normally not be feasible when developing code by hand. The automated code generation for the general and efficient solution of finite element variational problems is one of the key features of FEniCS.
There are three different interfaces to FFC: a Python interface, a just-in-time (JIT) compilation interface and a command-line interface. Only the latter is presented here while details on the other two interfaces can be found in Logg et al. (2012c). The command-line interface takes a UFL form file or a list of form files as input:
Bash code
$ ffc Poisson.ufl
The form file contains the UFL specification of elements and/or forms, as for instance the code from Figure 1.2 which in this case is saved in the filePoisson.ufl.
The content of a form file is wrapped in a Python script and then executed for further processing in FFC. There exist a number of optional command-line options to control the code generation. Related to this work, the most important options are:
-l language This parameter controls the output format for the generated code. The
default value is “ufc”, which indicates that the code is generated according
to the UFC specification. Alternatively, the value “dolfin” may be used to
generate code according to the UFC format with a small set of additional DOLFIN-specific wrappers.
-r representation This parameter controls the representation used for the
gener-ated element tensor code. There are three possibilities: “auto” (the default),
“quadrature” and “tensor”. FFC implements two different approaches to
code generation. One is based on traditional quadrature and another on a special tensor representation. This will be discussed in Section 3.2. In the case “auto”, FFC will try to select the better of the two representations; that
is, the representation that is believed to yield the best run-time performance for the problem at hand. This issue is addressed in detail in Section 3.6.
-O If this option is used, the code generated for the element tensor is optimised
for run-time performance. The optimisation strategy used depends on the chosen representation. In general, this will increase the time required for FFC to generate code, but should reduce the run-time for the generated code. Note that for very complicated variational forms, hardware limitations can make compilation with some optimisation options impossible. Optimisation strategies are treated in Chapter 3.
As an illustration of the options presented above, the command:
Bash code
$ ffc -l dolfin -r quadrature -O Poisson.ufl
will cause FFC to generate code for the Poisson problem, including DOLFIN wrappers using the quadrature representation with the default optimisation. A list of all available command-line parameters can be seen in FFC manual page by typing ‘man ffc’ on the command-line.
FFC follows the conventional design of a compiler in that it breaks compilation into several sequential stages. The output generated at each stage serves as input for the following stage, as illustrated in Figure 1.3. Introducing separate stages allows development and improvement of each stage to be implemented without affecting other stages of the compilation. Furthermore, adding new stages and dropping existing stages becomes trivial. Each of the stages involved when compiling a form is described in the following. Compilation of elements follow a similar (but simpler) set of stages, and is not described here.
Compiler stage 0: Language (parsing). In this stage, the user-specified form is interpreted and stored as a UFL abstract syntax tree (AST). The actual pars-ing is handled by Python and the transformation to a UFL form object is implemented by operator overloading in UFL.
Input: Python code or.uflfile
Output: UFL form
Compiler stage 1: Analysis. This stage preprocesses the UFL form and extracts form metadata (FormData), such as which elements were used to define the
form, the number of coefficients and the cell type (interval, triangleor tetrahedron). This stage also involves selecting a suitable quadrature scheme
and representation (as discussed earlier) for the form if these have not been specified by the user.
Input: UFL form
Figure 1.3: Compilation of finite element variational forms broken into six se-quential stages: Language, Analysis, Representation, Optimisation, Code gen-eration and Code Format-ting. Each stage gen-erates output based on input from the previous stage. The input/output data consist of a UFL form file, a UFL object, a UFL object and metadata com-puted from the UFL ob-ject, an intermediate rep-resentation (IR), an opti-mised intermediate repre-sentation (OIR), C++ code and, finally, C++ code files (from Logg et al. (2012c)).
Stage 1 Analysis Stage 2 Representation Stage 3 Optimization Stage 4 Code generation Stage 0 Language Stage 5 Code formatting Foo.h / Foo.cpp UFL + metadata IR OIR C++ code Foo.ufl UFL
Compiler stage 2: Code representation. Most of the complexity of compilation is handled in this stage which examines the input and generates all data needed for the code generation. This includes generation of finite element basis func-tions, extraction of data for mapping of degrees of freedom, and generation of the form representation, see Section 3.2, which may involve precomputation of integrals. Both representations available in FFC use tabulated values of finite element basis functions and their derivatives at a suitable set of inte-gration points on the reference element. FFC itself does not generate these values, but relies on the library FIAT (Kirby, 2004, 2012) for the computation of basis functions and their derivatives.
The intermediate representation is stored as a Python dictionary, mapping names of UFC functions to the data needed for generation of the correspond-ing code. In simple cases, likeufc::form::rank, this data may be a simple
number like2. In other cases, likeufc::cell_tensor::tabulate_tensor, the
data may be a complex data structure that depends on the choice of form representation.
Input: preprocessed UFL form and form metadata Output: intermediate representation (IR)
Compiler stage 3: Optimisation. This stage examines the intermediate representa-tion and performs optimisarepresenta-tions. The optimisarepresenta-tion strategy depends on the chosen form representation, see Section 3.3 for optimisations pertinent to the quadrature representation. Data stored in the intermediate representation dictionary is then replaced by new data that encode an optimised version of the function in question.
Input: intermediate representation (IR)
Output: optimised intermediate representation (OIR)
Compiler stage 4: Code generation. This stage examines the optimised intermedi-ate representation and generintermedi-ates the actual C++ code for the body of each UFC function. The code is stored as a dictionary, mapping names of UFC functions to strings containing the C++ code. As an example, the data generated for
ufc::form::rankmay be the string “return 2;”.
This demonstrates the importance of separating stages 2, 3 and 4 as it allows stages 2 and 3 to focus on algorithmic aspects related to finite elements and variational forms, while stage 4 is concerned only with generating C++ code from a set of instructions prepared in earlier compilation stages.
Input: optimised intermediate representation (OIR) Output: C++ code
Compiler stage 5: Code formatting. This stage examines the generated C++ code and formats it according to the UFC format, generating as output one or more
ufc::mesh ufc::function ufc::cell_integral
ufc::cell ufc::finite_element ufc::exterior_facet_integral ufc::form ufc::dofmap ufc::interior_facet_integral
Table 1.3: C++ classes defined in the UFC interface.
.h/.cppfiles conforming to the UFC specification. This is where the actual
writing of C++ code takes place and the stage relies on templates for UFC code available as part of the UFC moduleufc_utils.
Input: C++ code Output: C++ code files
The interface to the code which is generated by FFC is discussed in the following section.
1.3.4 Unified Form-assembly Code
The purpose of Unified Form-assembly Code (UFC) (Alnæs et al., 2009, 2012) is to provide an interface between the problem-specific code generated by form compilers and general-purpose problem solving environments like DOLFIN (described in Section 1.3.5) which implements, among other things, the finite element assembly algorithm. In contrast to other FEniCS components, few changes are made to UFC in order maintain a stable interface between form compilers and DOLFIN. This section gives a brief introduction to the interface, with emphasis on the functions relevant for this work. Furthermore, the UFC numbering convention for mesh entities is discussed.
The UFC interface provides a small set of abstract C++ classes, shown in Ta-ble 1.3 which are commonly used for assembling finite element tensors. The
meshandcellclasses are simple data structures that provide information such
as the geometric dimension and the topological dimension. In addition, the cell
class provides an array of global indices for the mesh entities belonging to the given cell (cell.mesh_entities) and an array of coordinates of the vertices of
the cell (cell.coordinates). The classes function and finite_element define
interfaces for general tensor-valued functions and finite elements respectively. Theformclass defines an interface for assembly of the global tensor
correspond-ing to the given form. This includes functions to createfinite_element,dofmap
and integral objects (ufc::cell_integral, ufc::exterior_facet_integral and ufc::interior_facet_integral) of the variational form.
Of particular interest in relation to this work are thedofmapand integral classes.
The local-to-global degree of freedom mapping on the finite element cell T, ιT, is computed by thedofmap::tabulate_dofsfunction for which the UFC interface is
C++ code /// Tabulate the local-to-global mapping of dofs on a cell
virtual void tabulate_dofs(unsigned int* dofs,const mesh& m, constcell& c) const= 0;
wheredofsis a pointer to an array for the tabulated values on T. UFC only provides
the interface of this function, it is not concerned with computing ιT. The code to compute ιT must be generated by the form compiler. For example, FFC will generate the following code for linear Lagrange elements on triangles.
C++ code /// Tabulate the local-to-global mapping of dofs on a cell
virtual void tabulate_dofs(unsigned int* dofs,const ufc::mesh& m, const ufc::cell& c) const
{
dofs[0] = c.entity_indices[0][0]; dofs[1] = c.entity_indices[0][1]; dofs[2] = c.entity_indices[0][2]; }
Note that FFC associates each degree of freedom with the global vertex number which can be extracted from thecell::entity_indicesarray. For discontinuous
linear Lagrange elements on triangles the generated code is
C++ code /// Tabulate the local-to-global mapping of dofs on a cell
virtual void tabulate_dofs(unsigned int* dofs,const ufc::mesh& m, const ufc::cell& c) const
{
dofs[0] = 3*c.entity_indices[2][0]; dofs[1] = 3*c.entity_indices[2][0] + 1; dofs[2] = 3*c.entity_indices[2][0] + 2; }
because FFC considers all degrees of freedom local to the given element and therefore compute degree of freedom numbers based on the global cell index.
The local finite element tensor is computed inside thetabulate_tensorfunction
which is implemented by all three integral classes although the interface varies slightly. For thecell_integral, the interface is
C++ code /// Tabulate the tensor for the contribution from a local cell
virtual void tabulate_tensor(double* A, constdouble * const * w, const cell& c) const = 0;
where A is a pointer to an array which will hold the values of the local
v0 v1 v2 Vertex Coordinates v0 x = (0, 0) v1 x = (1, 0) v2 x = (0, 1)
Figure 1.4: The UFC reference triangle and the coordinates of the vertices. in the integral. The code which FFC generates for this function varies depend-ing on, for example, the choice of representation and optimisation, issues which are discussed in Chapter 3. (Figures 3.2 and 3.3, on page 63 and 65 respec-tively, show examples of code generated by FFC for this function.) The inter-face forexterior_facet_integral::tabulate_tensoris similar in nature to the
interface forinterior_facet_integral::tabulate_tensorwhich is discussed in
Section 4.1.2 in connection to automation of discontinuous Galerkin methods. The UFC specification also defines a numbering scheme for mesh entities which allows form compilers to access necessary data consistently when generating code, for example, for computing the local tensors and local-to-global mapping as discussed above. Important aspects of this numbering scheme are summarised in the following for triangular cells. Further details on the UFC numbering convention can be found in Alnæs et al. (2012).
The UFC reference triangle, including the coordinates of the three vertices, is shown in Figure 1.4. Mesh entities are identified by the tuple(d, i)where d is the topological dimension of the mesh entity and i is a unique global index of the mesh entity. For convenience, mesh entities of topological dimension 0 are referred to as vertices, entities of dimension 1 are referred to as edges and entities of dimension 2 are referred to as faces. Mesh entities of topological dimension D−1 (codimension 1), with D denoting the topological dimension of the cells of the computational mesh, are referred to as facets. Thus for a triangular mesh, the facets are the edges and for a tetrahedral mesh, the facets are the faces. Following this convention, the vertices of a triangle are identified as v0 = (0, 0), v1 = (0, 1) and v2 = (0, 2), the edges (facets) are e0= (1, 0), e1= (1, 1)and e2= (1, 2), and the cell itself is c0= (2, 0).
The vertices of simplicial cells (intervals, triangles and tetrahedra) are numbered locally based on the corresponding global vertex numbers such that a tuple of increasing local vertex numbers corresponds to a tuple of increasing global vertex numbers. This is illustrated for a simple mesh in Figure 1.5. The remaining mesh entities are numbered within each topological dimension based on a lexicographical
0 1 3 2 v0 v1 v2 v0 v1 v2
Figure 1.5: Local vertex numbering of simplicial mesh based on global vertex numbers.
Entity Incident vertices Non-incident vertices v0= (0, 0) (v0) (v1, v2) v1= (0, 1) (v1) (v0, v2) v2= (0, 2) (v2) (v0, v1) e0= (1, 0) (v1, v2) (v0) e1= (1, 1) (v0, v2) (v1) e2= (1, 2) (v0, v1) (v2) c0= (2, 0) (v0, v1, v2) ∅
Table 1.4: Local numbering of mesh entities on triangular cells.
ordering of ordered tuples of non-incident vertices. For example, the first edge, e0, of a triangle is located opposite vertex v0as shown in Figure 1.6a. The numbering of mesh entities on triangular cells is shown in Table 1.4.
The relative ordering of mesh entities with respect to other incident mesh entities follows by sorting the entities by their indices. Therefore, the pair of vertices incident to edge e0in Figure 1.6a is(v1, v2), not(v2, v1). Due to the vertex numbering convention, this means that two incident simplicial cells will always agree on the orientation of incident subsimplices (for instance facets). This is demonstrated in Figure 1.6b, which shows two incident triangles which agree on the orientation of the common edge. This feature is advantageous when generating code for discontinuous Galerkin methods, as will be demonstrated in Chapter 4. 1.3.5 DOLFIN
Up until now, only the variational form and finite element discretisation has been defined. To obtain a solution to the boundary value problem in (1.1) the computational domain and boundary conditions must be specified which in the
v0 v1 v2
e0
(a) Edges are numbered based on the non-incident vertex. Therefore, e0is located
op-posite vertex v0. 0 1 3 2 v0 v1 v2 v0 v1 v2 e0 e2
(b) Orientation of facets (edges) are defined by the ordered tuple of incident vertices thus e0= (v0, v2)and e2= (v0, v1).
Figure 1.6: Edge numbering and orientation based on sorted tuples of incident and non-incident vertices. As a consequence two incident triangles will always agree on the orientation of the common facet for simplicial cells.
context of FEniCS is handled via a component called DOLFIN, a C++/Python library, which also provides algorithms for finite element assembly and linear algebra functionality to solve the arising system of equations. DOLFIN provides a problem solving environment and is the main user interface to FEniCS. A detailed presentation of DOLFIN is outside the scope of this work but can be found in Logg and Wells (2010) and Logg et al. (2012d).
The necessary DOLFIN functionality to implement a complete solver for the Poisson problem is presented. The intention is to give an impression of the possibilities that are offered by DOLFIN and an understanding of the basic concepts that are developed and used in subsequent chapters. For the model problem under consideration the domainΩ= [0, 1]× [0, 1], in which the source term f = 8π2sin(2πx)sin(2πy)is present, is subjected to homogeneous Dirichlet boundary conditions, g=0 onΓD=∂Ω.
A complete C++ solver for this problem is shown in Figures 1.7 and 1.8. The first line in Figure 1.7 includes the DOLFIN library, while the second line includes the UFC conforming code generated by FFC based on the UFL input for the Poisson problem shown in Figure 1.2. Then follows the definition of the classSourcewhich
is a subclass of theExpressionclass. AnExpressionrepresents a function that can
be evaluated on a finite element space and to suit this purpose it implements an
evalfunction. This function takes as arguments an array ofvalueswhich holds
the return values and an arrayxwhich contains the coordinates of the point where
theExpressionis currently being evaluated. TheSourceclass overloads theeval
C++ code
#include <dolfin.h>
#include "Poisson.h"
using namespace dolfin;
// Source term
classSource : public Expression {
void eval(Array<double>& values, const Array<double>& x) const {
values[0] = 8*pow(DOLFIN_PI, 2)*sin(2*DOLFIN_PI*x[0])*sin(2*DOLFIN_PI*x[1]); }
};
// Sub domain for Dirichlet boundary condition classDirichletBoundary : public SubDomain {
bool inside(const Array<double>& x, bool on_boundary) const {
return on_boundary; }
};
Figure 1.7: Implementation of source term and Dirichlet boundary for the C++ solver for the boundary value problem in (1.1). Program continues in Figure 1.8.
valuesarray.
Next follows the definition of the class DirichletBoundary, a subclass of SubDomain, for the part of the boundary where Dirichlet boundary conditions
are to be applied. TheSubDomainclass implements the functioninsidewhich
eval-uates to true or false depending on whether or not the point given by coordinatesx
is part of the subdomain. In addition to the argumentx, theinsidefunction also
takes the argumenton_boundary, a boolean value, supplied by DOLFIN, which
is true if the pointxis located on ∂Ω. In the given case, the Dirichlet condition
is indeed applied on ∂Ω which means that the overloadedinsidefunction can
simply be implemented by returning theon_boundaryargument.
The remaining part of the C++ solver, themainfunction, is shown in Figure 1.8.
The first line defines the computational mesh and consists of 2048 triangles as the unit square is divided into 32×32 cells and each cell is divided into two 2 triangles. DOLFIN provides functionality for creating simple meshes through the classes: UnitInterval,UnitSquare,UnitCube,UnitCircle,UnitSphere,Interval, RectangleandBoxwhich are useful for testing. For ‘real’ applications, a user can
read a mesh from file in the following way:
C++ code
Mesh mesh("mesh.xml");
provided that the mesh is saved in the DOLFIN XML format. Meshes can be generated by external libraries, such as Gmsh (http://geuz.org/gmsh/), stored in
the Gmsh data format and converted by thedolfin-convertscript to the DOLFIN
XML data format.
Next, the FunctionSpace is defined for the finite element function space Vh
in (1.7). A function space is represented by aMesh, aDofMapand aFiniteElement.
TheDofMapandFiniteElementclasses are generated by FFC based on theelement
definition in Figure 1.2. However, by including the ‘-l dolfin’ option when
compiling the UFL input with FFC:
Bash code
ffc -l dolfin Poisson.ufl
the DOLFIN wrappers are generated, permitting a user to instantiate a function space simply by providing the mesh as argument to the constructor.
The next three lines define an object for the Dirichlet boundary condition u = g = 0 on the boundary ΓD defined by the DirichletBoundary class from
Figure 1.7. The value g=0 is simply represented as a constant.
Then follows the creation of the bilinear and linear forms of the Poisson prob-lem using the function space V as argument. The Poisson::BilinearFormand Poisson::LinearFormclasses are part of the code inPoisson.hgenerated by FFC
C++ code
int main() {
// Create mesh and function space
UnitSquare mesh(32, 32); Poisson::FunctionSpace V(mesh);
// Define boundary condition
Constant g(0.0);
DirichletBoundary boundary; DirichletBC bc(V, g, boundary);
// Define variational forms
Poisson::BilinearForm a(V, V); Poisson::LinearForm L(V); Source f; L.f = f; Constant h(0.0); L.h = h; // Compute solution Function u(V); Matrix A; Vector b; assemble(A, a); assemble(b, L); bc.apply(A, b); solve(A, *u.vector(), b);
// Save solution in PVD format
File file("poisson.pvd"); file << u; // Plot solution plot(u); return 0; }
Figure 1.8: Continuation from Figure 1.7 of C++ code for the Poisson boundary value problem.
and attached to the linear form. The coefficientfis defined by the classSourceas
shown in Figure 1.7, while the coefficienth, the Neumann boundary condition, is
zero in the given case.
AFunction uis then declared to hold the computed solution. TheFunction
class represents a finite element function in V and therefore takes a function space as argument. The functionualso holds a vector of values of the degrees of freedom
associated with the function. A function is evaluated based on linear combinations of basis functions and the values of this vector. This is in contrast to theExpression
class which is evaluated by overloading theevalfunction as seen in Figure 1.7.
To compute a solution foruwhich satisfies the variational problem, defined by
the bilinear and linear formsaandL, the following three steps are applied. Firstly,
the bilinear and linear formsaandLare assembled into theMatrix Aand theVector bby calling the free functionassemblewhich implements an algorithm to assemble
finite element variational forms. The assembly algorithm will be presented later in this section. Secondly, the Dirichlet boundary condition is applied to the linear system of equations using theapplymember function of theDirichletBCobject bc. Thirdly, after applying the boundary condition, the system of equations can be
solved by calling the free functionsolvewhich solves linear systems on the form
Ax=b using the assembled matrixA, the vector of degree of freedom values from uand the assembled vectorbas arguments.
As an alternative to the three steps outlined above, thesolvefunction provides
functionality to solve variational problems in a straightforward fashion namely by: C++ code
solve(a == L, u, bc);
which automatically assembles the system, applies the boundary conditions and solves the linear system which is stored in the functionu.
Finally, the solution is saved in ParaView Data (PVD) format (http://www. paraview.org/) for external post processing and plotted by the built-inplot
com-mand of DOLFIN which enables a quick visual inspection of the computed solution. The computed solution to the Poisson boundary value problem is shown in Fig-ure 1.9.
Python interface
As already mentioned, DOLFIN also provides a Python interface as an alternative to the C++ interface. Most of the Python interface is generated automatically from the C++ interface using SWIG (http://www.swig.org/). In addition, the
Python interface offers seamless integration with UFL and FFC through just-in-time compilation of variational forms and elements which, in combination with the expressiveness of Python, allows solvers to be implemented very compactly. For
Figure 1.9: Computed solution to the Poisson boundary value problem. The warped scalar field u in the figure on the right has been scaled by a factor of 0.5.
this reason, the Python interface to DOLFIN is preferred, whenever feasible, for the examples presented in this thesis.
As an example, the complete solver for the Poisson boundary value problem using the Python interface is shown in Figure 1.10. The code is very similar to the C++ code in Figures 1.7 and 1.8 and the differences are mainly due to the difference in Python and C++ syntax. The two main differences are the definition of the
FunctionSpaceand the definition of the variational forms which are implemented
directly as part of the solver and not in a separate file. Also note that the UFL coordinates have been used to implement the source term f directly as part of the variational formulation. It could also be implemented by subclassing the
Expressionclass and overloading theevalfunction in a similar way to the approach
in the C++ example:
Python code
classSource(Expression):
def eval(self, values, x):
values[0] = 8*pi**2*sin(2*pi*x[0])*sin(2*pi*x[1]) f = Source()
As an alternative, it could be implemented by:
Python code
f = Expression("8*pow(pi,2)*sin(2*pi*x[0])*sin(2*pi*x[1])")
where the string argument to theExpressionclass is given in C++ syntax which is
automatically just-in-time compiled in order to evaluate theExpression. Compared
Python code
from dolfin import *
# Create mesh and define function space
mesh = UnitSquare(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)
# Define Dirichlet boundary
def boundary(x, on_boundary):
return on_boundary
# Define boundary condition
g = Constant(0.0)
bc = DirichletBC(V, g, boundary)
# Define variational problem
u = TrialFunction(V) v = TestFunction(V) x = V.cell().x f = 8*pi**2*sin(2*pi*x[0])*sin(2*pi*x[1]) a = inner(grad(u), grad(v))*dx L = f*v*dx # Compute solution U = Function(V) solve(a == L, U, bc)
# Save solution in PVD format
file = File("poisson.pvd") file << U
# Plot solution
plot(U, interactive=True)
will take place in C++ rather than Python. Assembly algorithm
To conclude this short introduction to the FEniCS Project, the assembly algorithm, implemented in the DOLFINassemblefunction, is presented. The presentation is
given for the assembly of the rank two tensor corresponding to the bilinear form of the Poisson problem in (1.2). A generalisation of the algorithm for multilinear forms is given in Alnæs et al. (2009) and Logg et al. (2012b).
Setting the function space ˆV=V, the tensor A which arises from assembling the bilinear form a is defined by
AI =a φI2, φI1
, (1.9)
where I = (I1, I2)is a multi-index and φk Nk=1is a basis for V. The tensor A is a sparse rank two tensor, a matrix, of dimensions N×N. The matrix A is computed by iterating over the cells of the mesh and adding the contribution from each local cell to the global matrix A. In this case, from (1.3), the local cell tensor ATis defined as: AT,i=aT φiT2, φiT1= Z T∇u· ∇v dx, (1.10)
where i = (i1, i2)is a multi-index AT,i is the ith entry of the cell tensor AT, aTis the local contribution to the form from a cell T∈ ThandnφkTo3
k=1is the local finite element basis for V on T, which is linear Lagrange elements on triangles in this case.
To formulate the assembly algorithm, a local-to-global mapping of degrees of freedom is needed. Let ιT :IT → I denote the collective local-to-global mapping for each T∈ Th
ιT(i) =ι1T(i1), ι2T(i2) ∀i∈ IT, (1.11) where ιjT : [1, 3] → [1, N] denotes the local-to-global mapping for each discrete function space VjandITis the index set
IT= 2
∏
j=1 [1, 3] = (1, 1),(1, 2),(1, 3),(2, 1),(2, 2),(2, 3),(3, 1),(3, 2),(3, 3) . (1.12) That is, ιTmaps a tuple of local degrees of freedom to a tuple of global degrees of freedom. DOLFIN calls the tabulate_tensor and tabulate_dofs functionspresented in Section 1.3.4, in order to compute the local contribution aT and the local-to-global mapping ιTj for each discrete function space from which DOLFIN constructs the collective local-to-global mapping ιT.
The assembly of the matrix A can now be carried out efficiently by iterating over all cells T∈ Th. On each cell T, the cell tensor ATis computed and then added to the global tensor A as outlined in Algorithm 1. The algorithm can be extended Algorithm 1Assembly algorithm.
A=0 forT∈ Th (1) Compute ιT (2) Compute AT (3) Add AT to A according to ιT: fori∈ IT AιT(i) + = AT,i end for end for
to handle assembly over exterior and interior facets, the latter is demonstrated in Section 4.1.4.
2 FEniCS applications to solid mechanics
One of the goals of this work is to tackle complicated solid mechanics models using automated modelling tools. In the previous chapter it was shown how automated modelling could be employed to solve the finite element variational formulation of a Poisson boundary value problem. The Poisson problem provides a simple platform for introducing the concepts behind automated modelling as it is implemented in the FEniCS framework. However, from the simple presentation it may not be immediately clear how more complex problems, like plasticity, can be solved. A natural step is, therefore, to apply the concept of automated modelling to some standard solid mechanics problems. Solid mechanics problems typically involve the standard momentum balance equation, posed in a Lagrangian setting, with different models distinguished by the choice of nonlinear or linearised kinematics, and the constitutive model for determining the stress. The traditional development approach to solid mechanics problems, and traditional finite element codes, places a strong emphasis on the implementation of constitutive models at the quadrature point level. Automated methods, on the other hand, tend to stress more heavily the governing balance equations. Widely used finite element codes for solid mechanics applications provide application programming interfaces (APIs) for users to implement their own constitutive models. The interface supplies kinematic and history data, and the user code computes the stress tensor, and when required also the linearisation of the stress. Users of such libraries will typically not be exposed to code development other than via the constitutive model API.In addition to demonstrating how solid mechanics problems can be solved using automation tools, this chapter presents some of the models that will be further investigated and extended in subsequent chapters. It is not intended as a comprehensive treatment of solid mechanics problems, but should be viewed as a stepping stone towards implementation of classes of plasticity models. The chapter also focuses on some pertinent issues that arise due to the nature of the constitutive models. These issues, and solid mechanics problems in general, have motivated a number of developments in the FEniCS framework.
The common problems of linearised elasticity, plasticity, hyperelasticity and elastic wave propagation are considered. Topics that are addressed in this chapter
via these problems include ‘off-line’ computation of stress updates, linearisation of problems with off-line stress updates, automatic differentiation and time stepping for problems with second-order time derivatives. The presentation starts with the relevant governing equations and the constitutive models under consideration. The important issue of solving and linearising problems in which the governing equation is expressed in terms of the stress tensor (rather than explicitly in terms of the displacement field, or derivatives of the displacement field), and the stress tensor is computed via a separate algorithm is then addressed. These topics are then followed by a number of examples that demonstrate implementation approaches of the described models.
To conclude the chapter, which is primarily based on the work in Ølgaard and Wells (2012a); Ølgaard et al. (2008b), extensions of the FEniCS framework that are particular interesting with respect to solid mechanics problems, and consequently to this work, are summarised.
2.1
Governing equations
2.1.1 PreliminariesThe considered problems will be posed on a polygonal domainΩ ⊂Rd, where 1≤ d ≤3. The boundary of Ω, denoted by ∂Ω, is decomposed into regions ΓD and ΓN such that ΓD∪ΓN = ∂Ω and ΓD∩ΓN = ∅. The outward unit normal vector on ∂Ω will be denoted by n. For time-dependent problems, a time interval of interest I = (0, T]will be considered, where superimposed dots denote time derivatives. The current configuration of a solid body is denoted byΩ; that is, the domainΩ depends on the displacement field. It is sometimes convenient to also define a reference domainΩ0⊂Rdthat remains fixed. For convenience, cases in whichΩ and Ω0coincide at time t = 0 are considered. To indicate boundaries, outward unit normal vectors, and other quantities relative toΩ0, the subscript ‘0’ will be used. When considering linearised kinematics, the domainsΩ and Ω0are both fixed and coincide at all times t. A triangulation of the domainΩ will be denoted byTh, and a triangulation of the domainΩ0will be denoted byTh,0.
The governing equations for the different models will be formulated in the common framework of: find u∈V such that
F(u; w) =0 ∀w∈V, (2.1)
where F : V×V→R is linear in w and V is a suitable function space. If F is also linear in u, then F can be expressed as
