*230*
*76*
*40MB*

*English*
*Pages [177]*
*Year 2020*

- Author / Uploaded
- Ali Gerami Matin
- Arzhang Angoshtari

*Table of contents : CoverHalf TitleTitle PageCopyright PageDedicationContentsPrefaceChapter 1: OverviewChapter 2: Mathematical Preliminaries 2.1. Real Numbers 2.2. Functions 2.3. Linear Spaces, Linear Mappings, and Bilinear Forms 2.4. Linear Independence, Hamel Bases, and Dimension 2.5. The Matrix Representation of Linear Mappings and Bilinear Forms 2.6. Normed Linear Spaces 2.7. Functionals and Dual Spaces 2.8. Green’s Formulas Exercises Comments and ReferencesChapter 3: Finite Element Interpolation 3.1. 1D Finite Element Interpolation 3.1.1. The Global Level 3.1.2. The Local Level 3.2 Finite Elements 3.2.1. Simplicial Lagrange Finite Elements of Type (k) 3.2.2. Simplicial Hermite Finite Elements of Type (3) 3.2.3. The Raviart-Thomas Finite Element 3.2.4. The Nedelec Finite Element 3.3. Meshes 3.4. Finite Element Spaces and Interpolations 3.4.1. H1-Conformal Finite Element Spaces 3.4.1.1. Lagrange Elements 3.4.1.2. Hermite Elements 3.4.2. H(div)-Conformal Finite Element Spaces 3.4.3. H(curl)-Conformal Finite Element Spaces 3.4.4. Affine Families of Finite Elements 3.5. Convergence of Interpolations Exercises Computer Exercises Comments and ReferencesChapter 4: Conforming Finite Element Methods for PDEs 4.1. Second-Order Elliptic PDEs 4.2. Weak Formulations of Elliptic PDEs 4.2.1. Dirichlet Boundary Condition 4.2.2. Neumann Boundary Condition 4.2.3. Robin Boundary Condition 4.3. Well-posedness of Weak Formulations 4.4. Variational Structure 4.5. The Galerkin Method and Finite Element Methods 4.5.1. The Stiffness Matrix 4.5.2. Well-posedness of Coercive Discrete Problems 4.5.3. Convergence of Finite Element Solutions 4.6. Implementation: The Poisson Equation 4.6.1. Dirichlet Boundary Condition 4.6.2. Mixed Dirichlet-Neumann Boundary Condition 4.6.3. Robin Boundary Condition 4.7. Time-Dependent Problems: Parabolic Problems 4.7.1. Finite Element Approximations using the Method of Lines 4.7.2. Temporal Discretization 4.7.3. Implementation: A Diffusion Problem 4.8. Mixed Finite Element Methods 4.8.1. Mixed Formulations 4.8.2. Mixed Methods and inf-sup Conditions 4.8.3. Implementation Exercises Computer Exercises Comments and ReferencesChapter 5: Applications 5.1. Elastic Bars 5.2. Euler-Bernoulli Beams 5.3. Elastic Membranes 5.4. The Wave Equation 5.5. Heat Transfer in a Turbine Blade 5.6. Seepage in Embankment 5.7. Soil Consolidation 5.8. The Stokes Equation for Incompressible Fluids 5.9. Linearized Elasticity 5.10. Linearized Elastodynamics: The Hamburg Wheel-Track Test 5.11. Nonlinear Elasticity Exercises Computer ExercisesAppendix A: FEniCS InstallationAppendix B: Introduction to Python B.1. Running Python Programs B.2. Lists B.3. Branching and Loops B.4. Functions B.5. Classes and Objects B.6. Reading and Writing Files B.7. Numerical Python Arrays B.8. Plotting with MatplotlibReferencesIndex*

Finite Element Methods in Civil and Mechanical Engineering

Finite Element Methods in Civil and Mechanical Engineering

A Mathematical Introduction

Dr Arzhang Angoshtari Ali Gerami Matin

First edition published 2021 by CRC Press 6000 Broken Sound Parkway NW, Suite 300, Boca Raton, FL 33487-2742 and by CRC Press 2 Park Square, Milton Park, Abingdon, Oxon, OX14 4RN © 2021 Taylor & Francis Group, LLC CRC Press is an imprint of Taylor & Francis Group, LLC Reasonable efforts have been made to publish reliable data and information, but the author and publisher cannot assume responsibility for the validity of all materials or the consequences of their use. The authors and publishers have attempted to trace the copyright holders of all material reproduced in this publication and apologize to copyright holders if permission to publish in this form has not been obtained. If any copyright material has not been acknowledged please write and let us know so we may rectify in any future reprint. Except as permitted under U.S. Copyright Law, no part of this book may be reprinted, reproduced, transmitted, or utilized in any form by any electronic, mechanical, or other means, now known or hereafter invented, including photocopying, microfilming, and recording, or in any information storage or retrieval system, without written permission from the publishers. For permission to photocopy or use material electronically from this work, access www.copyright.com or contact the Copyright Clearance Center, Inc. (CCC), 222 Rosewood Drive, Danvers, MA 01923, 978750-8400. For works that are not available on CCC please contact [email protected] Library of Congress Cataloging-in-Publication Data Names: Angoshtari, Arzhang, author. | Matin, Ali Gerami, author. Title: Finite element methods in civil and mechanical engineering : a mathematical introduction / Arzhang Angoshtari, Ali Gerami Matin. Description: Boca Raton : CRC Press, 2020. | Includes bibliographical references and index. Identifiers: LCCN 2020030081 (print) | LCCN 2020030082 (ebook) | ISBN 9781138335172 (paperback) | ISBN 9781138335165 (hardback) | ISBN 9780429442506 (ebook) Subjects: LCSH: Finite element method. | Civil engineering--Mathematics. | Mechanical engineering--Mathematics. Classification: LCC TA347.F5 A54 2020 (print) | LCC TA347.F5 (ebook) | DDC 518/.25--dc23 LC record available at https://lccn.loc.gov/2020030081 LC ebook record available at https://lccn.loc.gov/2020030082 Trademark notice: Product or corporate names may be trademarks or registered trademarks, and are used only for identification and explanation without intent to infringe. ISBN: 9781138335165 (hbk) ISBN: 9781138335172 (pbk) ISBN: 9780429442506 (ebk) Typeset in Computer Modern font by KnowledgeWorks Global Ltd. by Arzhang Angoshtari Visit the eResource: www.routledge.com/9781138335165

Dedication To Our Parents

Contents Preface.......................................................................................................................xi Chapter 1

Overview ......................................................................................... 1

Chapter 2

Mathematical Preliminaries............................................................. 5 2.1 2.2 2.3 2.4 2.5

Real Numbers ......................................................................... 5 Functions................................................................................. 5 Linear Spaces, Linear Mappings, and Bilinear Forms ........... 6 Linear Independence, Hamel Bases, and Dimension ............. 9 The Matrix Representation of Linear Mappings and Bilinear Forms ...................................................................... 10 2.6 Normed Linear Spaces.......................................................... 12 2.7 Functionals and Dual Spaces ................................................ 13 2.8 Green’s Formulas.................................................................. 14 Exercises ........................................................................................ 15 Comments and References............................................................. 17 Chapter 3

Finite Element Interpolation.......................................................... 19 3.1

1D Finite Element Interpolation ........................................... 19 3.1.1 The Global Level ..................................................... 19 3.1.2 The Local Level ....................................................... 24 3.2 Finite Elements ..................................................................... 26 3.2.1 Simplicial Lagrange Finite Elements of Type (k) .... 26 3.2.2 Simplicial Hermite Finite Elements of Type (3)...... 30 3.2.3 The Raviart-Thomas Finite Element........................ 31 3.2.4 The Nedelec Finite Element..................................... 33 3.3 Meshes .................................................................................. 34 3.4 Finite Element Spaces and Interpolations............................. 37 3.4.1 H 1 -Conformal Finite Element Spaces ..................... 38 3.4.1.1 Lagrange Elements .................................. 39 3.4.1.2 Hermite Elements .................................... 43 3.4.2 H(div)-Conformal Finite Element Spaces............... 43 3.4.3 H(curl)-Conformal Finite Element Spaces.............. 45 3.4.4 Affne Families of Finite Elements .......................... 45 3.5 Convergence of Interpolations .............................................. 46 Exercises ........................................................................................ 50 Computer Exercises ....................................................................... 51 vii

viii

Contents

Comments and References............................................................. 52 Chapter 4

Conforming Finite Element Methods for PDEs ............................ 55 4.1 4.2

Second-Order Elliptic PDEs ................................................. 55 Weak Formulations of Elliptic PDEs.................................... 56 4.2.1 Dirichlet Boundary Condition.................................. 57 4.2.2 Neumann Boundary Condition ................................ 58 4.2.3 Robin Boundary Condition ...................................... 59 4.3 Well-posedness of Weak Formulations................................. 60 4.4 Variational Structure ............................................................. 61 4.5 The Galerkin Method and Finite Element Methods ............. 62 4.5.1 The Stiffness Matrix ................................................ 63 4.5.2 Well-posedness of Coercive Discrete Problems ...... 64 4.5.3 Convergence of Finite Element Solutions................ 64 4.6 Implementation: The Poisson Equation................................ 66 4.6.1 Dirichlet Boundary Condition.................................. 66 4.6.2 Mixed Dirichlet-Neumann Boundary Condition ..... 70 4.6.3 Robin Boundary Condition ...................................... 74 4.7 Time-Dependent Problems: Parabolic Problems .................. 76 4.7.1 Finite Element Approximations using the Method of Lines .................................................................... 78 4.7.2 Temporal Discretization........................................... 79 4.7.3 Implementation: A Diffusion Problem .................... 79 4.8 Mixed Finite Element Methods ............................................ 82 4.8.1 Mixed Formulations................................................. 83 4.8.2 Mixed Methods and inf-sup Conditions .................. 84 4.8.3 Implementation ........................................................ 85 Exercises ........................................................................................ 90 Computer Exercises ....................................................................... 92 Comments and References............................................................. 93 Chapter 5

Applications................................................................................... 97 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 5.10

Elastic Bars ........................................................................... 97 Euler-Bernoulli Beams ....................................................... 100 Elastic Membranes.............................................................. 104 The Wave Equation............................................................. 105 Heat Transfer in a Turbine Blade........................................ 107 Seepage in Embankment..................................................... 112 Soil Consolidation .............................................................. 117 The Stokes Equation for Incompressible Fluids ................. 122 Linearized Elasticity ........................................................... 126 Linearized Elastodynamics: The Hamburg Wheel-Track Test...................................................................................... 129 5.11 Nonlinear Elasticity ............................................................ 134

ix

Contents

Exercises ...................................................................................... 143 Computer Exercises ..................................................................... 145 Appendix A

FEniCS Installation ..................................................................... 147

Appendix B

Introduction to Python................................................................. 149 B.1 B.2 B.3 B.4 B.5 B.6 B.7 B.8

Running Python Programs.................................................. 149 Lists..................................................................................... 150 Branching and Loops.......................................................... 151 Functions............................................................................. 152 Classes and Objects ............................................................ 152 Reading and Writing Files .................................................. 154 Numerical Python Arrays ................................................... 155 Plotting with Matplotlib...................................................... 157

References ............................................................................................................. 159 Index...................................................................................................................... 161

Preface In the past century, the fnite element method has grown signifcantly and nowadays it is considered as a standard tool for dealing with many engineering and scientifc problems. Numerous computer programs are available today for the implementation of the fnite element method, which relieve one from implementation of the main aspects of the fnite element method from scratch. Signifcant contributions to this feld have been made by engineers and mathematicians. Methods developed by engineers are sometimes suitable for specifc applications or circumstances and cannot be generalized easily. They may lead to confusions when applied to new applications, some of which can be addressed by using rich results already available in the mathematical literature. However, mathematical results are usually not easily accessible to many engineers due to their rigorous nature. The main goal of this book is to provide a modern and concise introduction to the theory and the implementation of the fnite element method for engineers and physicists. The emphasis is on the fundamental aspects of the fnite element method for solving linear partial differential equations. We use a mathematical language to introduce main concepts. Although this approach may look cumbersome at the beginning, it is very useful and powerful not only because it provides a suitable framework to develop, implement, and study fnite element methods but also because it can provide suffcient background and self-confdence to the readers to refer to the mathematical literature. We use a Python based package called FEniCS to implement fnite element methods. FEniCS syntax is close to the mathematical language and it yields concise, readable, and effcient programs. This book is based on the graduate-level course that the frst author taught at the George Washington University in 2016–2019. It is suitable for a one-semester graduate course in engineering and physics. It may be used in a course for advanced undergraduate students in mathematics as well. The reader is assumed to be familiar with elementary calculus and basic programming skills. The book is organized as follows: Chapter 1 contains an overview of the topics that will be discussed in the book. In Chapter 2, mathematical preliminaries are briefy introduced. We only discuss main ideas and results. Although this shallow preliminary chapter is far from being complete, based on our own experience with graduate engineering students, we believe it enables the readers who do not have solid background in mathematics to understand the main ideas and refer to suitable references when necessary. Finite elements, fnite element spaces, and fnite element interpolations are discussed in Chapter 3. In addition to Lagrange and Hermite fnite elements, we also mention the so-called edge and face elements for the curl and div operators. The convergence of fnite element interpolations is briefy discussed as well. xi

xii

Preface

In Chapter 4, we study the application of fnite elements for solving specifc classes of time-independent and time-dependent partial differential equations that arise in many applications. We mention basic results on the well-posedness and the convergence of fnite element methods. A brief discussion of mixed fnite element methods and inf-sup conditions are also given in this chapter. Moreover, FEniCS implementations of fnite element methods for various boundary conditions are provided. In Chapter 5, a selection of applications in civil and mechanical engineering is studied and Python programs for their implementation are provided. Numerous examples, exercises, and computer exercises are given in each chapter. Additional comments and suggested references for further reading are mentioned at the end of the chapters. Several Python programs for the implementation of examples are provided throughout the book. Previous familiarity with Python, although very helpful, is not necessary. Two short appendices on FEniCS installation and Python are provided at the end of the book. The latter is a brief introduction intended for the readers who are not familiar with Python. For the readers convenience, Jupyter notebooks containing all Python programs of each chapter are available online on the companion website of the book. We would like to thank Marmar Mehrabadi for her valuable feedback. We are also grateful to all people in CRC Press who helped us during the preparation of this book. A special thanks goes to our editor Tony Moore for encouraging us to write this book. Arzhang Angoshtari, Ali Gerami Matin Washington, DC May 2020

1 Overview In this introductory chapter, we briefy discuss the basic aspects of the fnite element method for solving differential equations. The goal is to give an overview of the book and to motivate the topics that will be discussed in the following chapters. Only main ideas are discussed and the readers are expected to focus only on the big picture and not on details. Technical details of this chapter are the subject of the following ones. This chapter is not a prerequisite for the rest of the book and the readers may skip this chapter if they wish. Roughly speaking, one may consider the following four stages for the approximation of responses of a system using the fnite element method: (i) Identifying the governing equations, which are assumed to be differential equations in this book; (ii) Deriving weak forms of the governing equations; (iii) Discretizing weak forms by using the Galerkin method to obtain a discrete problem; and (iv) Employing piecewise polynomial fnite element spaces to solve the discrete problem. To be more specifc, let us consider the approximation of the equilibrium position of an elastic membrane Ω shown in Figure 1.1, which is fxed at its boundary ∂ Ω and is under a vertical load f .

Figure 1.1 An elastic membrane under a vertical load: (Left) Undeformed confguration Ω and the load f (dashed arrows); (Right) The vertical displacement u(x, y) at the point (x, y).

Governing Equations The laws of physics imply that the vertical displacement u of the membrane Ω satisfes the boundary value problem −Δ u = f , in Ω, (1.1) u = 0, on ∂ Ω, where Δ is the Laplacian, i.e. Δ u = ∂xx u + ∂yy u. This governing equation is also called Poisson’s equation. Weak Forms The equation (1.1) is called the strong form of the governing equation. To solve the strong form, one should fnd a function which is at least twice differentiable and 1

2

Finite Element Methods in Civil and Mechanical Engineering

satisfes (1.1). Usually, it is hard to directly fnd such a function. Alternatively, it may be easier to fnd a suitable solution candidate frst and then showing that the solution candidate is a real solution by establishing its differentiability. We can obtain a solution candidate for the strong form (1.1) by using one of its weak forms. For example, let w(x, y) be an arbitrary function that vanishes on the boundary. Multiplying the strong form by w and then taking the integral over Ω, one concludes that Z

Z ∂x u ∂x w + ∂y u ∂y w dx dy = f w dx dy,

Ω

for all w.

(1.2)

Ω

A solution u of (1.1) is also a solution of (1.2), however, the converse may not be true in general. Notice that unlike a solution of (1.1) that should be twice differentiable, a solution of (1.2) may be merely once differentiable since only the frst derivatives of u are present in (1.2). The equation (1.2) is called a weak form of (1.1) due to this weaker differentiability requirement. A solution of (1.2) can be considered as a solution candidate for (1.1). For computational purposes, weak forms are usually more suitable than strong forms. The Galerkin Method Since fnding an analytical solution for the weak form is not easy in general, we can try to approximate its solution. The Galerkin method provides a systematic framework for this purpose. More specifcally, given a set of functions {ψ1 , . . . , ψN } called basis functions or global shape functions, one assumes that a solution of the weak form can be approximated by using a function uh of the form N

uh (x, y) = ∑ Ui ψi (x, y), i=1

where Ui , i = 1, . . . , N, are unknown constants. To determine the constants Ui , we assume uh is a solution of the following problem, which is a special case of (1.2) with w = ψ1 , . . . , ψN : Z

Z

∂x uh ∂x ψi + ∂y uh ∂y ψi dx dy =

Ω

f ψi dx dy,

i = 1, . . . , N.

(1.3)

Ω

The above framework for approximating a solution of the weak form (1.2) is called the Galerkin method. The problem (1.3) is called a discrete problem or a discretization of (1.2) and uh is called a discrete solution in the sense that the approximate solution uh is determined by a fnite number of unknowns Ui , i = 1, . . . , N. The discrete problem (1.3) is very suitable for numerical analysis. In fact, it is equivalent to the linear system of equations AN×N · UN×1 = FN×1 ,

(1.4)

3

Overview

where ⎡

⎤

⎢ A(ψ1 , ψ1 ) A(ψ2 , ψ1 ) · · · A(ψN , ψ1 ) ⎥ ⎡ ⎢ ⎥ ⎢ ⎥ ⎢ A(ψ , ψ ) A(ψ , ψ ) · · · A(ψ , ψ ) ⎥ ⎢ ⎢ N 1 2 2 2 2 ⎥ ⎥, U = ⎢ A=⎢ ⎢ ⎢ ⎥ .. .. .. ⎣ ⎢ ⎥ . . . ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ A(ψ1 , ψN ) A(ψ2 , ψN ) · · · A(ψN , ψN ) with A(ψ j , ψi ) =

Z

U1 U2 .. .

⎤

⎡

⎥ ⎢ ⎥ ⎢ ⎥, F = ⎢ ⎦ ⎣

UN

F(ψ1 ) F(ψ2 ) .. .

⎥ ⎥ ⎥, ⎦

F(ψN )

∂x ψ j ∂x ψ j + ∂y ψ j ∂y ψi dx dy,

ZΩ

F(ψi ) =

⎤

(1.5) f ψi dx dy.

Ω

The matrix A is usually called the stiffness matrix and the vector F is called the load vector. In summary, the Galerkin method for approximating a solution of (1.2) results in solving the linear system (1.4) for obtaining the unknown coeffcients Ui . The Finite Element Method We saw that once a set of basis functions {ψ1 , . . . , ψN } is chosen, one can apply the Galerkin method to (1.2) to approximate a solution of (1.1). But how can we determine a suitable choice of basis functions? There are many different approaches for obtaining basis functions which yield different approximation methods such as fnite element methods, fnite-difference methods, spectral methods, etc. We will study fnite element methods in this book. The main aspects of the fnite element method for obtaining basis functions can be summarized as follows. (i) A mesh (also called a triangulation) Th is established over Ω. Each mesh consists of several cells or elements that can be, in principle, of arbitrary shape, see Figure 1.2.

Figure 1.2 A triangular and a rectangular mesh of the rectangular membrane Ω.

(ii) Given a mesh Th , basis functions are chosen such that their restriction to any element of Th are polynomials or “close to” polynomials. The convergence of fnite element methods is closely related to this aspect. Moreover, this leads to simple and effcient computations of stiffness matrices and load vectors.

4

Finite Element Methods in Civil and Mechanical Engineering

(iii) Basis functions are selected to be locally supported, i.e., they are nonzero only on a few elements of Th . For example, given the triangular mesh of Figure 1.2, one may select the basis functions ψ1 and ψ2 such that they are zero everywhere except on a few elements shown in Figure 1.3. Then, the defnition (1.5) of the entries of the stiffness matrix A implies that A(ψ1 , ψ2 ) = A(ψ2 , ψ1 ) = 0. Due to this aspect of the fnite element method, the stiffness matrix A is sparse, i.e., most of the entries of A are zero. The sparsity of A allows us to effciently compute the solution of the linear system (1.4).

Figure 1.3 The shaded elements depict regions where the shape functions ψ1 (left) and ψ2 (right) are nonzero.

Implementation. Based on the Galerkin method, the fnite element method provides a computationally effcient framework for approximating solutions of partial differential equations by using locally supported, piecewise polynomial basis functions. In this book, we will employ FEniCS, a Python package, to implement fnite element methods. An interesting feature of FEniCS is that its syntax is close to the mathematical notion. Figure 1.4 shows a fnite element approximation of a deformation of the elastic membrane of Figure 1.1 under a constant vertical load. In the remainder of this book, we will study the ideas and notions of this chapter in more details.

Figure 1.4 A fnite element approximation of a deformation of the elastic membrane Ω of Figure 1.1 under a constant vertical load: (Left) The underlying mesh; (Right) the fnal confguration of Ω.

2 Mathematical Preliminaries Normed linear spaces and linear operators play a central role in the fnite element method. In this chapter, we study these notions and also fx our notations. We will only mention main ideas in a level suffcient for understanding the discussions of the following chapters. For more examples and discussions on technical details, we refer the reader to the references mentioned at the end of this chapter.

2.1

REAL NUMBERS

The set of real numbers is denoted by R. For any real numbers a, b with a < b, (a, b) represents the open interval {x ∈ R : a < x < b} and [a, b] represents the closed interval {x ∈ R : a ≤ x ≤ b}. By abusing the notation, we also use (a, b) to denote a point in the plane R2 but the meaning is always clear from the context. The empty set is denoted by ∅. Suppose A 6= ∅ is a non-empty subset of R. We say that A is bounded from above if there is a number u such that a ≤ u, for all a ∈ A, and we call u an upper bound of A. We say A is bounded from below if there is a number l such that a ≥ l, for all a ∈ A, and we call l a lower bound of A. The set A is called bounded if it is bounded from below and above. The maximum of A ⊂ R, denoted by max A, is a member of A which is also an upper bound of A. Similarly, the minimum of A ⊂ R, denoted by min A, is a member of A which is also a lower bound of A. For example, max[0, 1] = 1 and min[0, 1] = 0. Notice that the maximum and the minimum of a bounded set may not exist, for example, consider the open interval (0, 1). A basic fact about R allows one to extend the notions of maximum and minimum to all bounded subsets in the following sense: If B ⊂ R is non-empty and is bounded from above, then there is the least upper bound of B, which is called the supremum of B and is denoted by sup B. Similarly, if a nonempty set B is bounded from below, then there is the greatest lower bound of B, which is called the infmum of B and is denoted by inf B. For example, inf(0, 1) = 0, and sup(0, 1) = 1. If max B exits, then max B = sup B, and if min B exists, then min B = inf B.

2.2

FUNCTIONS

Let X and Y be two sets. Recall that a function (also called a mapping or an operator) f : X → Y assigns a unique value f (x) ∈ Y to each x ∈ X. The function f is called a Y -valued function. The set X is called the domain of f and is also denoted by D( f ). A function g is called an extension of f if D( f ) ⊂ D(g), and f (x) = g(x), for all x ∈ D( f ). In this case, the function f is called a restriction of g and we write f = g|D( f ) . Example 2.1. The function g : R → R, g(x) = |x|, is an extension of f : (0, 1) → R, f (x) = x, or equivalently, we have f = g|(0,1) . 5

6

Finite Element Methods in Civil and Mechanical Engineering

The range of f : X → Y is a subset of Y defned as R( f ) = {y ∈ Y : there is x ∈ X with y = f (x)}. If R( f ) = Y , f is called onto or surjective. A function f : X → Y is called one-toone or injective if different points do not take the same values, that is, f (x1 ) = f (x2 ), implies that x1 = x2 . A mapping which is injective and surjective is called bijective. A bijective function f : X → Y is invertible, that is, f admits the inverse function f −1 : Y → X, where for any x ∈ X and y ∈ Y , we have f −1 ( f (x)) = x, and f ( f −1 (y)) = y. Example 2.2. The function f : (0, 1) → R, f (x) = x, is injective but not surjective, while g : (0, 1) → (0, 1), g(x) = x, is bijective. Suppose A ⊂ X and B ⊂ Y . For any arbitrary function f : X → Y , the sets f (A) ⊂ Y and f −1 (B) ⊂ X are defned as f (A) = {y ∈ Y : y = f (x) for some x ∈ A}, and f −1 (B) = {x ∈ X : f (x) ∈ B}. Notice that for defning the set f −1 (B), f does not need to be invertible. Example 2.3. Let f : R → R, f (x) = x2 . Then, f ([0, 2]) = [0, 4], and f −1 ([−2, 1]) = [−1, 1]. Roughly speaking, we say Ω is an open subset of the Euclidean space Rn if Ω is n-dimensional and it does not contain its boundary. The closure Ω of the open subset Ω is the union of Ω and the set of its boundary points ∂ Ω. For example, let the unit ball B1 (0) be the set of all points of R3 such that their distance from the origin 0 is smaller than 1. Then, B1 (0) is an open subset of R3 and B1 (0) is the union of the unit ball and its boundary, which is the unit sphere with radius 1 centered at the origin. Let Ω ⊂ Rn be open. The set of all functions u : Ω → R with continuous derivatives up to order m is denoted by Cm (Ω). For example, suppose Ω ⊂ R2 . Then, C0 (Ω) is the set of continuous functions on Ω and C2 (Ω) is the set of functions u such that u, ∂x u, ∂y u, ∂xx u, ∂xy u, ∂yy u are continuous on Ω, where ∂x u is the partial derivative with respect to x and so on. The set Cm (Ω) is the set of functions u ∈ Cm (Ω) such that u and all of its partial derivatives up to order m can be continuously extended to Ω. Notice that Cm (Ω) ⊂ Cm (Ω) but Cm (Ω) 6⊂ Cm (Ω). Example 2.4. Let Ω = (0, 1), with Ω = [0, 1], and consider f (x) = 1/x, and g(x) = x2 . Then, f ∈ C1 (Ω) but f ∈ / C1 (Ω), while g belongs to both C1 (Ω) and C1 (Ω).

2.3

LINEAR SPACES, LINEAR MAPPINGS, AND BILINEAR FORMS

A linear space, also called a vector space, is a set X together with an addition and a scalar multiplication. Let x, y, z ∈ X and a, b ∈ R. The addition and the scalar multiplication, respectively denoted by x + y and ax, are assumed to satisfy the following conditions: − The addition of any two members of X is another member of X such that x + y = y + x, and x + (y + z) = (x + y) + z;

Mathematical Preliminaries

7

− There is a unique zero member 0 ∈ X with 0 + x = x. Also for any x, there is a unique member −x ∈ X such that x + (−x) = 0; − The scalar multiplication satisfes 1x = x, 0x = 0, and a(bx) = (ab)x; − The addition and the scalar multiplication satisfy a(x + y) = ax + ay, and (a + b)x = ax + bx. For brevity, we only mention a linear space X when its addition and scalar multiplication are known. A simple example of a linear space is the Euclidean space Rn with its standard addition and scalar multiplication. Less trivial linear spaces which are frequently encountered in fnite element methods are discussed in the following examples. Example 2.5. Given an open subset Ω ⊂ Rn , the space Cm (Ω) is a linear space: For any f , g ∈ Cm (Ω), and c ∈ R, the addition and the scalar multiplication are defned pointwise as ( f + g)(x) = f (x) + g(x), and (c f )(x) = c f (x), for all x ∈ Ω. Similarly, one can show that Cm (Ω) is a linear space. Example 2.6. Let Pk (Rn ) denote the space of all polynomials of degree ≤ k in Rn . For example, P2 (R) is the space of polynomials of the form c0 + c1 x + c2 x2 , where ci , i = 0, 1, 2, are arbitrary constants, and P2 (R2 ) is the space of polynomials of the form c00 + c10 x + c01 y + c20 x2 + c11 xy + c02 y2 , where ci j , i, j = 0, 1, 2, are arbitrary constants. By using the pointwise operations similar to those of the previous example, one can show that Pk (Rn ) is a linear space. For any A ⊂ Rn , the set of polynomials on A defned as Pk (A) = {p|A : p ∈ Pk (Rn )}, is also a linear space. L2 Lebesgue Example 2.7 (Lebesgue spaces). For anyR open subset Ω ⊂ Rn , the R 2 2 space is defned as L (Ω) = { f : Ω → R : Ω | f | exists}, where by Ω | f |2 we mean R 2 Ω | f (x1 , . . . , xn )| dx1 · · · dxn . For brevity, we follow this notation for integrals from now on. For example, let Ω = (0, 1), and f (x) = xα . Then, f ∈ L2 (Ω) if α > − 12 . By using the pointwise operations, one can show that L2 (Ω) is a linear space. A vector feld v : Ω → Rn is said to be of L2 -class if its components belong to L2 (Ω). The space of L2 vector felds is denoted by [L2 (Ω)]n and is a linear space. Example 2.8 (Sobolev spaces). The Sobolev space H 1 (Ω) is the space of functions in L2 (Ω) which also have L2 frst-order derivatives, that is, H 1 (Ω) = { f ∈ L2 (Ω) : ∂xi f ∈ L2 (Ω), i = 1, . . . , n}. For example, assuming β > 12 , Ω = (0, 1), and f (x) = xβ , then f ∈ H 1 (Ω). The Sobolev space of vector felds Ω → Rn with components belonging to H 1 (Ω) is denoted by [H 1 (Ω)]n . One can also defne the Sobolev space H m (Ω) for any positive integer m. For example, H 2 (Ω) is the space of functions in L2 (Ω) with frst-order and second-order derivatives also belonging to L2 (Ω). These Sobolev spaces become linear spaces when equipped with the pointwise operations. Example 2.9 (Partly Sobolev classes). Recall that the divergence of a vector feld v : Ω → Rn with components v = (v1 , . . . , vn ) is defned as div v = ∑ni=1 ∂xi vi . The partly Sobolev class associated to the divergence operator is defned as H(div; Ω) = {vv ∈ [L2 (Ω)]n : div v ∈ L2 (Ω)}. This space becomes a linear space with pointwise

8

Finite Element Methods in Civil and Mechanical Engineering

operations. One can also defne the partly Sobolev space H(curl; Ω): Recall that in R3 , the curl operator is given by curl v = (∂x2 v3 − ∂x3 v2 , ∂x3 v1 − ∂x1 v3 , ∂x1 v2 − ∂x2 v1 ). Then, H(curl; Ω) = {vv ∈ [L2 (Ω)]3 : curl v ∈ [L2 (Ω)]3 }, is a linear space. By using the 2-dimensional curl operator curl v = ∂x1 v2 − ∂x2 v1 , the partly Sobolev class H(curl; Ω) can be defned for 2-dimensional vector felds as well. One can show that [H 1 (Ω)]n ⊂ H(div; Ω) ⊂ [L2 (Ω)]n , and [H 1 (Ω)]n ⊂ H(curl; Ω) ⊂ [L2 (Ω)]n . A subset S of a linear space X is called a linear subspace of X if with respect to the addition and the scalar multiplication of X, S is also a linear space. Example 2.10. Any line passing through the origin is a linear subspace of R2 . Example 2.11. Let ∂ Ω be the boundary of an open subset Ω ⊂ Rn . Then, H01 (Ω) = { f ∈ H 1 (Ω) : u|∂ Ω = 0}, is a linear subspace of H 1 (Ω). The polynomial space Pk (Ω) is another subspace of H 1 (Ω). Let X and Y be linear spaces. A mapping L : X → Y is called a linear mapping if for any fnite subset {x1 , . . . , xn } ⊂ X and {c1 , . . . , cn } ⊂ R, we have L(∑ni=1 ci xi ) = ∑ni=1 ci L(xi ). Example 2.12. An m × n matrix M can be considered as a mapping M : Rn → Rm , defned as M(x) = M · x. This mapping is a linear mapping. Example 2.13. For any f ∈ C0 (Ω), consider the mapping L : C0 (Ω) → R, given by R L(u) = Ω f u. Since for any ui ∈ C0 (Ω), and ci ∈ R, i = 1, . . . , n, we have n

L

�

Z n n Z � n c u = f · c u = c f u = i i i i i i ∑ ∑ ∑ ∑ ci L(ui ),

i=1

Ω

i=1

i=1

Ω

i=1

the operator L is a linear operator. Example 2.14. Suppose a ∈ Rm and M is an m × n matrix. A mapping T : Rn → Rm of the form T (x) = a + M · x is said to be an affne mapping. Affne mappings are not linear. The null space of a linear mapping L : X → Y is defned as N(L) = L−1 ({0}), that is, N(L) = {x ∈ X : L(x) = 0}. The null space of L is also called the kernel of L and is denoted by ker L. The null space of L : X → Y is a linear subspace of X and the range R(L) of L is a linear subspace of Y . Moreover, L is injective if and only if ker L = 0.

9

Mathematical Preliminaries

Example 2.15. For the linear mapping L : R3 → R3 given by ⎡ ⎤ 100 L = ⎣0 0 0⎦, 000 the spaces N(L) and R(L) are respectively the y-z plane and the x-axis. Given linear spaces X and Y , a mapping B : X × Y → R is called a bilinear form if for any fnite subsets {x, x1 , . . . , xn } ⊂ X, {y, y1 , . . . , ym } ⊂ Y , and {a1 , . . . , an }, {b1 , . . . , bm } ⊂ R, we have n

n

m

B( ∑ ai xi , y) = ∑ ai B(xi , y), and B(x, ∑ b j y j ) = i=1

i=1

j=1

m

∑ b j B(x, y j ).

j=1

Thus, B(x, y) is a bilinear form if it is linear with respect to each of its entries. A bilinear form B : X × X → R is symmetric if B(x1 , x2 ) = B(x2 , x1 ), for all x1 , x2 ∈ X. Example 2.16. An m × n matrix M can be considered as a bilinear form B : Rn × Rm → R, given by B(x, y) = yT · M · x, where yT is the transpose of the column vector y. Since yT · M · x = xT · M T · y, one concludes that Mn×n is a symmetric bilinear form if and only if the matrix M is symmetric, that is, M = M T . Example 2.17. The mapping B : C0 (Ω) × C0 (Ω) → R, B( f , g) = metric bilinear form.

2.4

R Ω

f g, is a sym-

LINEAR INDEPENDENCE, HAMEL BASES, AND DIMENSION

Suppose a set S, perhaps having infnite members, is a subset of a linear space X. A point s ∈ X is a linear combination of points in S if there are fnite points s1 , . . . , sn in S such that s = ∑ni=1 ci si , ci ∈ R, i = 1, . . . , n. The set of all linear combinations of points of S, also called the linear subspace spanned by S, is denoted by span S. Example 2.18. Let S = {1, x, x2 } ⊂ P2 (R). Then span S = P2 (R). A subset S of a linear space X is called linearly independent if for any fnite collection {s1 , . . . , sn } of elements of S, the only solution of the equation ∑ni=1 ci si = 0, is c1 = · · · = cn = 0. If a subset of X is not linearly independent, it is called linearly dependent. If S is linearly independent, one can show that any nonzero x ∈ span S admits a unique representation x = ∑ni=1 ci si , where si ∈ S, and ci ’s are nonzero numbers, i = 1, . . . , n. A linearly independent subset S of a linear space X is called a (Hamel) basis for X if X = span S. One can show that different bases of X have the same number of members. This number is called the dimension of X and is denoted by dim X. If dim X is fnite, X is called a fnite-dimensional linear space, otherwise, X is called an infnite-dimensional linear space.

10

Finite Element Methods in Civil and Mechanical Engineering

Example 2.19. Consider the set of polynomials S = {1, x, x2 } ⊂ P2 (R). Then S is linearly independent, since let c0 + c1 x + c2 x2 = 0, for any x ∈ R. If at least one of the constant ci is not zero, we have a polynomial on the left side of the equality, which has at most 2 zeros. But this is a contradiction as the left side must be zero for all x ∈ R, and therefore, c0 = c1 = c2 = 0, and S is linearly independent. As span S = P2 (R), we conclude that dim P2 (R) = 3. More generally, one can show that � (n+k)! dim Pk (Rn ) = n+k = k!n! . k Example 2.20. Let P(R) be the space of polynomials of all degrees k ≥ 0 in R and let S = {1, x, x2 , x3 , . . . }. By considering an arbitrary linear combination of fnite members of S and following the approach of Example 2.19, it is straightforward to show that S is linearly independent. Since span S = P(R), we conclude that P(R) is an infnite-dimensional linear space. Since P(R) is a linear subspace of Cm (Ω), Cm (Ω), L2 (Ω), and H 1 (Ω), we conclude that these spaces are also infnitedimensional. Moreover, the partly Sobolev classes H(div; Ω) and H(curl; Ω) are also infnite-dimensional as the space of vector felds with P(R) components is a subspace of these partly Sobolev classes. Let L : X → Y be a linear mapping. The rank-nullity theorem states that dim N(L) + dim R(L) = dim X.

(2.1)

Example 2.21. For the linear mapping L : R3 → R3 of Example 2.15, we have dim N(L) = 2, and dim R(L) = 1 and therefore, dim N(L) + dim R(L) = 3. Example 2.22. Let M : Rn → Rn be a linear mapping. If M is surjective, we have dim R(M) = n, and then (2.1) implies that M is injective as dim N(L) = 0. Conversely, (2.1) implies that if M is injective, then it is also surjective. Thus, a linear mapping Rn → Rn is injective if and only if it is surjective. Example 2.23. As discussed in Example 2.12, a matrix Mm×n can be considered as a linear mapping Rn → Rm . The maximum rank of this linear mapping is m and occurs if it is surjective. In this case, the matrix M is called full rank.

2.5

THE MATRIX REPRESENTATION OF LINEAR MAPPINGS AND BILINEAR FORMS

Suppose X and Y are fnite-dimensional linear spaces with the bases BX = {x1 , . . . , xn }, and BY = {y1 , . . . , ym }, respectively. Any linear mapping L : X → Y can be represented by an m × n matrix as follows: Let L(xi ) = ∑mj=1 l ji y j , where l ji is the j-th component of L(xi ) in the basis BY . Since any x ∈ X can be uniquely written as x = ∑ni=1 ai xi , we can write n

L(x) = L

�

∑ ai xi

i=1

n

n

= ∑ ai L(xi ) = ∑ ai · i=1

i=1

m

�

∑ l ji y j

j=1

m

=

n

�

∑ ∑ l ji ai

j=1 i=1

y j.

11

Mathematical Preliminaries

In the basis BY , L(x) ∈ Y admits a unique representation L(x) = ∑mj=1 b j y j , and therefore, the above relation implies that b j = ∑ni=1 l ji ai , j = 1, . . . , m. This equation can be stated in the matrix form b = [L] · a, where ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ b1 l11 l12 · · · l1n a1 ⎢ b2 ⎥ ⎢ l21 l22 · · · l2n ⎥ ⎢ a2 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ bm×1 = ⎢ . ⎥ , [L]m×n = ⎢ . . . (2.2) . ⎥ , an×1 = ⎢ .. ⎥ . ⎣ .. ⎦ ⎣ .. .. . . .. ⎦ ⎣ . ⎦ lm1 lm2 · · · lmn an bm Notice that a and b respectively contain the components of x in the basis BX and the components of L(x) in the basis BY . The matrix [L] is a matrix representation of L. This matrix representation is not unique and depends on the chosen bases for X and Y . Example 2.24. Consider the linear mapping D : P3 (R) → P2 (R), D(p) = 2 d p(x) dx . We choose the basis B1 = {1, x, x2 , x3 } for P3 (R) and the basis B2 = {1, x, x2 } for P2 (R). Since D(1) = 0, D(x) = 2, D(x2 ) = 4x, and D(x3 ) = 6x2 , the matrix representation of L in the bases B1 and B2 reads ⎡ ⎤ 0200 [D] = ⎣ 0 0 4 0 ⎦ . 0006 Different bases for P3 (R) and P2 (R) lead to different representations for D, see Exercise 2.15. Similarly, one can also write matrix representation of bilinear forms. More specifically, let B : X ×Y → R be a bilinear form. Using the bases BX and BY , one can write n

B(x, y) = B

�

m

∑ ai xi , ∑ b j y j

i=1

j=1

m

=

n

∑ ∑ b j B(xi , y j )ai ,

j=1 i=1

which is equivalent to B(x, y) = bT · [B] · a, where the vector a and b are introduced in (2.2) and the matrix [B]m×n is the matrix representation of the bilinear form in the bases BX and BY given by ⎡ ⎤ B(x1 , y1 ) B(x2 , y1 ) · · · B(xn , y1 ) ⎢ B(x1 , y2 ) B(x2 , y2 ) · · · B(xn , y2 ) ⎥ ⎢ ⎥ [B]m×n = ⎢ ⎥. .. .. .. .. ⎣ ⎦ . . . . B(x1 , ym ) B(x2 , ym ) · · · B(xn , ym ) Notice that the i j-th component of [B] is B(x j , yi ). Example 2.25. Let Ω = (−1, 1) ⊂ R and consider the bilinear form B : P2 (Ω) × R P2 (Ω) → R, B( f , g) = Ω f g. To write the matrix representation of B in the basis

12

Finite Element Methods in Civil and Mechanical Engineering

{1, x, x2 }, we note that Z 1

B(1, 1) =

Z 1

dx = 2, B(1, x) = −1

B(x, x2 ) =

Z 1

xdx = 0, B(1, x2 ) = B(x, x) =

−1

x3 dx = 0, B(x2 , x2 ) =

−1

Z 1

2 x2 dx = , 3 −1

Z 1

2 x4 dx = , 5 −1

and since B is symmetric, we conclude that ⎡ ⎤ 2 0 23 ⎢ ⎥ [B] = ⎣ 0 23 0 ⎦ . 2 2 3 0 5

2.6

NORMED LINEAR SPACES

The notion of norm generalizes the notion of length to linear spaces. More specifcally, a norm on a linear space X is a real-valued function k · k : X → R such that 1. 2. 3. 4.

kxk ≥ 0, for all x ∈ X (positivity); kx + yk ≤ kxk + kyk, for any x, y ∈ X (the triangle inequality); kαxk = |α|kxk, for any x ∈ X and α ∈ R; kxk = 0, if and only if x = 0 (positive defniteness).

A normed linear space (X, k · k) is a linear space X together with a norm k · k on X. For any x ∈ X, kxk is called the norm of x. Example 2.26. A simple example of normed linear spaces is (Rn , k · k), where k · k is the standard norm of Rn , that is, for any x = (x1 , . . . , xn ) ∈ Rn , we have kxk = (x12 + · · · + xn2 )1/2 . Example 2.27. In this example and the next one, we introduce some normed linear spaces which are frequently encountered in the approximation of solutions of differential equations by using fnite element methods. The L2 -norm k · k2 on the Lebesgue space L2 (Ω) of Example 2.7 is given by k f k2 =

Z

| f |2

1 2

, f ∈ L2 (Ω).

Ω

The L2 -norm on [L2 (Ω)]n is componentwise defned, that is, for any v = (v1 , . . . , vn ) ∈ [L2 (Ω)]n , we have � 1 kvk v 2 = kv1 k22 + · · · + kvn k22 2 . For any f ∈ H 1 (Ω), the H 1 -norm k f k1,2 on the Sobolev space H 1 (Ω) of Example 2.8 is given by 1 � k f k1,2 = k f k22 + k∂x1 f k22 + · · · + k∂xn f k22 2 .

13

Mathematical Preliminaries

Finally, for any v = (v1 , . . . , vn ) ∈ [H 1 (Ω)]n , the H 1 -norm on [H 1 (Ω)]n is componentwise defned as � 1 v 1,2 = kv1 k21,2 + · · · + kvn k21,2 2 . kvk Example 2.28. By using the Partly Sobolev classes of Example 2.9, one can defne the normed linear spaces (H(div; Ω), k · kd ) and (H(curl; Ω), k · kc ), where the associated norm of a vector feld v belonging to these spaces is respectively defned as � 2 1 � 2 1 v d = kvk v 2 + kdiv vvk22 2 , and kvk v c = kvk v 2 + kcurl vvk22 2 . kvk The standard norm of the Euclidean space Rn enables us to measure the distance between different points x, y ∈ Rn by using kx − yk. Similarly, we may measure the “distance” between two arbitrary points x, y ∈ X of a normed linear space (X, k · k) by using kx − yk. Two functions f1 , f2 ∈ L2 (Ω) are “close” to each other in L2 (Ω) if k f1 − f2 k2 is small. Similarly, f1 , f2 ∈ H 1 (Ω) are close to each other in H 1 (Ω) if k f1 − f2 k1,2 is small. Suppose f1 and f2 belong to both L2 (Ω) and H 1 (Ω). Then, f1 and f2 may be close in L2 (Ω) but not in H 1 (Ω) since the smallness of k f1 − f2 k2 suggests that values of f1 and f2 are close, however, the smallness of k f1 − f2 k1,2 implies that not only values but also derivatives of f1 and f2 are close, see Figure 2.1. In numerical analysis, norms are employed to measure the difference between approximate solutions and the exact solutions of problems.

Figure 2.1 The two functions of the left panel are “close” in the L2 -norm but not in the H 1 -norm while the functions on the right panel are close in both norms.

2.7 FUNCTIONALS AND DUAL SPACES Let X be a linear space. Any real-valued mapping J : X → R is called a functional. If J is linear, it is called a linear functional. Some physical quantities such as energy can be expressed as functionals. = Ω u2 is a functional. For any Example 2.29. The mapping J : C0 (Ω) → R, J(u) R 0 0 v ∈ C (Ω), the mapping f : C (Ω) → R, f (u) = Ω uv, is a linear functional. R

Let X be a normed linear space. The space of all continuous linear functionals over X is called the dual space of X and is denoted by X 0 or L(X; R). By using the

14

Finite Element Methods in Civil and Mechanical Engineering

pointwise operations, it is easy to show that X 0 is a linear space. It is customary to denote the value of f ∈ X 0 at x ∈ X using the duality brackets as f (x) = h f , xi. Example 2.30 (Dual bases). Let X be a fnite-dimensional normed linear space. Then, X 0 is equal to the space of all linear functionals. Suppose BX = {x1 , . . . , xn } is a basis for X and consider the set of linear functionals BX 0 = { f1 , . . . , fn }, with fi (x j ) = δi j , where the Kronecker delta δi j is 1, if i = j, and 0 otherwise. It is not hard to show that BX 0 is a basis for X 0 : Let g ∈ X 0 be an arbitrary functional. Since BX is a basis for X, any x ∈ X admits a unique representation x = ∑ni=1 ai xi . We can write n

�

g(x) = g

∑ ai xi

i=1

i=1

n

= ∑ g(xi ) fi i=1

n n n � n � n = ∑ g(xi )ai = ∑ g(xi ) ∑ a j δi j = ∑ g(xi ) ∑ a j fi (x j ) i=1

n

�

j=1

i=1

j=1

n

∑ a jx j

j=1

= ∑ g(xi ) fi (x), i=1

which implies that all members of X 0 can be written as a linear combination of the elements of BX 0 , that is, X 0 = span BX 0 . On the other hand, let ∑ni=1 ci fi (x) = 0, for some constants c1 , . . . , cn , and any x ∈ X. In particular, at x = x j , we have 0 = ∑ni=1 ci fi (x j ) = ∑ni=1 ci δi j = c j , j = 1, . . . , n. Thus, BX 0 is also linearly independent. Hence, we conclude that BX 0 is a basis for X 0 and dim X = dim X 0 . The basis BX is said to be the dual basis of BX 0 . In the next chapter, we will see that fnite elements can be defned by using dual bases. Example 2.31. Let X = P1 (K), with K = [0, 1], and consider the set BX 0 = { f1 , f2 } ⊂ X 0 , where f1 (p) = p(0), and f2 (p) = p(1). We show that BX 0 is a basis for X 0 . Let c1 f1 (p)+c2 f2 (p) = 0, for all p ∈ P1 (K). In particular, for p1 (x) = 1−x, and p2 (x) = x, we obtain 0 = c1 f1 (p1 ) + c2 f2 (p1 ) = c1 p1 (0) + c2 p1 (1) = c1 , and 0 = c1 f1 (p2 ) + c2 f2 (p2 ) = c1 p2 (0) + c2 p2 (1) = c2 . Therefore, BX 0 is linearly independent and since dim X 0 = dim X = 2, we conclude that BX 0 is a basis for X 0 . Moreover, since f j (pi ) = δi j , BX = {p1 , p2 } is the dual basis of BX 0 .

2.8

GREEN’S FORMULAS

Let Ω ⊂ Rm and let n = (n1 , . . . , nm ) be the outward unit normal vector feld at the boundary ∂ Ω. The fundamental Green’s formula states that if v ∈ C1 (Ω), then Z

Z

∂xi v = Ω

v ni , i = 1, . . . , m.

(2.3)

∂Ω

By using this formula, one can write other Green’s formulas. For example, given a function v and a vector feld u = (u1 , . . . , um ), one can show that Z Ω

uu·∇v = −

Z

Z

v divu u+ Ω

(vuu)·n. n

(2.4)

∂Ω

where “·” is the standard inner product of Rm given by uu·vv = ∑m i=1 ui vi , and ∇v = (∂x1 v, . . . , ∂xm v) is the gradient of v.

15

Mathematical Preliminaries

Another Green’s formula for the Laplacian Δu = ∑m i=1 ∂xi xi u = div(∇u), can be written as follows: Let m

∂n v = nn·∇v = ∑ ni ∂xi v, i=1

be the normal derivative of v at the boundary ∂ Ω. Then, for any functions u and v, we have Z Z Z ∇u·∇v = − v Δu + v ∂n u. (2.5) Ω

Ω

∂Ω

The above Green’s formulas can be extended to Sobolev spaces based on the following results: 1. If v, u ∈ L2 (Ω) then Ω uv exists; R 2. If u ∈ H 1 (Ω), then u|∂ Ω ∈ L2 (∂ Ω) and ∂ Ω u exists, however, this conclusion generally does not hold if u only belongs to L2 (Ω); R 3. If v ∈ H 1 (Ω) and u ∈ H(div; Ω), then u · n belongs to L2 (∂ Ω) and ∂ Ω (vuu) · n exists. R

By using these results, one can show that Green’s formula (2.4) holds either if u ∈ [H 1 (Ω)]n and v ∈ H 1 (Ω), or if u ∈ H(div; Ω) and v ∈ H 1 (Ω). Similar results hold for (2.5) as well. As will be discussed in the following chapters, Green’s formulas in terms of Sobolev spaces are extensively used in the fnite element method to write weak forms of governing equations.

EXERCISES Exercise 2.1. Suppose B ⊂ R is a non-empty set. Show that: (i) If max B exists, then max B = sup B, and (ii) if min B exists, then min B = inf B. Exercise 2.2. Let Ω ⊂ Rn be an open subset. Is Pk (Ω) a linear subspace of L2 (Ω), Cm (Ω), and Cm (Ω)? Exercise 2.3. Let Ω = (0, 1) × (0, 1) be the unit square and let f (x, y) = (x + y)b . Find values of b such that: (i) f ∈ L2 (Ω); (ii) f ∈ H 1 (Ω); and (iii) f ∈ H 2 (Ω). Exercise 2.4. Show that affne mappings are not linear. Exercise 2.5. Let ∇ f denote the gradient of f ∈ C1 (Ω) and let “·” denoteR the standard inner product of Rn . Show that B : C1 (Ω) × C1 (Ω) → R, B( f , g) = Ω ∇ f ·∇g, is a symmetric bilinear form. Exercise 2.6. Let X and Y be linear spaces. Then, Z = X × Y is also a linear space with the addition (x1 , y1 ) + (x2 , y2 ) = (x1 + x2 , y1 + y2 ), for all (xi , yi ) ∈ Z, i = 1, 2, and the scalar multiplication c(x, y) = (cx, cy), for all c ∈ R. Suppose B : X ×Y → R is a bilinear form. Show that B as the mapping B : Z → R is not a linear mapping. Exercise 2.7. Let L : X → Y be a linear mapping. Show that: (i) ker L is a linear subspace of X; (ii) R(L) is a linear subspace of Y ; and (iii) the linear mapping L is injective if and only if ker L = {0}.

16

Finite Element Methods in Civil and Mechanical Engineering

Exercise 2.8. Show that if a linear mapping L : X → Y is invertible, then its inverse L−1 : Y → X is also a linear mapping. Exercise 2.9. Let S be a subset of a linear space X. Show that span S is a linear subspace of X. Exercise 2.10. By writing a basis, show that: (i) dim P2 (R2 ) = 6; (ii) dim P2 (R3 ) = 10; (iii) dim P3 (R2 ) = 10; and (iv) dim P3 (R3 ) = 20. Exercise 2.11. Let Q2 (R2 ) be the space of all polynomials of degree less than or equal to 2 with respect to each one of the two variables x and y. By deriving a basis for this linear space show that dim Q2 (R2 ) = 9. Also show that Q2 (R2 ) is a linear subspace of Pk (R2 ), k ≥ 4. Exercise 2.12. By generalizing the polynomial space of Exercise 2.11, defne Qk (Rn ) and show that dim Qk (Rn ) = (k + 1)n . This space is useful for defning fnite elements over rectangular elements. Exercise 2.13. Let L : Rn → Rm be a linear mapping. What can we say about the injectivity, the surjectivity, and the bijectivity of L if (i) n > m; (ii) n = m; and (iii) n < m? Justify your answers by using the rank-nullity theorem and give an example for each case. Exercise 2.14. Suppose a non-singular matrix Mk×k (that is, M as mapping Rk → Rk is bijective) is of the form " # An×n Bn×l Mk×k = , Cl×n 0 with k = n + l. Show that the inequality l ≤ n must hold. In Chapters 4 and 5, we will use similar results to study the stability of mixed fnite element methods. (Hint. Use the fact that the bijectivity of M implies the surjectivity of C and then employ Exercise 2.13.) Exercise 2.15. Show that B˜ 1 = {1 + 2x, x, x2 + x3 , x3 } and B˜ 2 = {2, x + 1, x2 } are bases for P3 (R) and P2 (R), respectively. Obtain the matrix representation of the linear mapping of Example 2.24 in these bases. Exercise 2.16. Let Ω be the open Rinterval (−1, 1) and consider the bilinear form B : P3 (Ω) × P3 (Ω) → R, B( f , g) = Ω ∇ f ·∇g. Derive the matrix representation of B in the basis {1, x, x2 , x3 } of P3 (Ω). Exercise 2.17. Let x = (x1 , . . . , xn ) ∈ Rn . Show that kxk = |x1 | + · · · + |xn | is a norm on Rn . (Hint. To prove the triangle inequality, you can use the triangle inequality of real numbers, that is, |a + b| ≤ |a| + |b| for any a, b ∈ R.)

Mathematical Preliminaries

17

Exercise 2.18. Let X = P1 (K), with K = [a1 , a2 ]. Show that BX 0 = { f1 , f2 }, where fi (p) = p(ai ), i = 1, 2, is a basis for X 0 . Obtain the dual basis of BX 0 . In Chapter 3, we will see that BX 0 contains the degrees of freedom of a Lagrange fnite element of degree 1. Exercise 2.19. Let X = P2 (K), with K = [a1 , a2 ]. Show that BX 0 = { f1 , f2 , f3 }, where fi (p) = p(ai ), i = 1, 2, 3, and a3 = (a1 + a2 )/2, is a basis for X 0 and obtain the dual basis of BX 0 . In Chapter 3, we will see that BX 0 contains the degrees of freedom of a Lagrange fnite element of degree 2. Exercise 2.20. Derive Green’s formula (2.4) by using the fundamental Green’s formula (2.3). Exercise 2.21. Derive Green’s formula (2.5) by using Green’s formula (2.4) and the relation Δu = div(∇u).

COMMENTS AND REFERENCES Discussions on real numbers and subsets of the Euclidean space and their properties are available in standard real analysis textbooks such as [3]. A good introduction to linear spaces and linear operators in engineering and science can be found in [11]. A classical text for standard Sobolev spaces is [1]. An introductory discussion on partly Sobolev classes H(div; Ω) and H(curl; Ω) is mentioned in [4, Chapter 2]. A brief discussion of Green’s formulas in the context of solving partial differential equations by using the fnite element method is given in [5, Section 1.2].

Element 3 Finite Interpolation The convergence of fnite element solutions to the true solution of physical models is closely related to the general capability of fnite elements to approximate arbitrary functions. This process of approximating functions with fnite elements is called fnite element interpolation. In this chapter, we introduce several fnite elements and discuss how they can be employed to interpolate functions and vector felds.

3.1 1D FINITE ELEMENT INTERPOLATION To motivate the concept of fnite elements and the main ideas of function interpolation using fnite elements, we begin this chapter by studying piecewise linear interpolations of functions defned on an open interval of R. 3.1.1

THE GLOBAL LEVEL

Consider the unit interval Ω = (0, 1). Given a function f : [0, 1] → R, a simple approach for interpolating (also called approximating) f is to use piecewise linear functions such as the one shown in Figure 3.1. The frst step for obtaining such interpolating functions (also called approximating functions) using fnite elements is to establish a mesh Th on Ω. To this end, we consider a set of points {x0 , x1 , . . . , xN } such that 0 = x0 < x1 < · · · < xN = 1. The mesh Th induced by these points is the collection of subintervals {K0 , . . . , KN−1 } with Ki = [xi , xi+1 ]. The index h of Th refers to the maximum length h of subintervals, that is, h = max{h0 , . . . , hN−1 }, where hi is the length of Ki . The points xi and the subintervals Ki are then called the vertices and the elements (or cells) of Th , respectively. Figure 3.1 suggests that by using the values { f (x0 ), . . . , f (xN )} of f at the vertices of Th , we can uniquely specify a piecewise linear function Ih f with the following properties: (i) Ih f (xi ) = f (xi ), for all vertices xi , and (ii) f |Ki ∈ P1 (Ki ) for all elements Ki . The function Ih f is called the Lagrange interpolant of f of type (1). The values { f (x0 ), . . . , f (xN )} that uniquely determine Ih f are called global degrees of freedom. Here the adjective global means these degrees of freedom are assigned to the whole mesh Th . Notice that Ih f belongs to the set P1 (Th ) that includes all continuous, piecewise linear functions over Th , that is P1 (Th ) = {v ∈ C0 (Ω) : v|Ki ∈ P1 (Ki ), i = 0, . . . , N − 1}.

(3.1)

This set is then called an approximation space as it includes the approximating function Ih f . One can use the fact that P1 (Th ) is a linear space to write Ih f in terms of a 19

20

Finite Element Methods in Civil and Mechanical Engineering

Figure 3.1 A function (the solid line) and its piecewise linear interpolant (the dashed line) associated to the vertices {x0 , . . . , x4 }.

basis of P1 (Th ). More specifcally, consider the hat functions ψi ∈ P1 (Th ), i = 0, . . . , N, given by ⎧ ⎪ h 1 (x − xi−1 ), if x ∈ Ki−1 , ⎨ i−1 1 (3.2) ψi (x) = if x ∈ Ki , (x − x), h ⎪ i i+1 ⎩ 0, otherwise. For the cases i = 0, N, straightforward modifcations of this defnition are needed as shown in Figure 3.2. We have ψ j (xi ) = δi j , i, j = 0, . . . , N.

(3.3)

Also the hat functions are non-zero only on a few elements.

Figure 3.2 The basis functions ψ0 , ψ2 , and ψ4 , which are respectively associated to the vertices x0 , x2 , and x4 .

One can show that {ψ0 , ψ1 , . . . , ψN } is a basis for P1 (Th ), and therefore, dim P1 (Th ) = N + 1. Thus, any v ∈ P1 (Th ) can be uniquely written as v = ∑Nj=0 c j ψ j . Using (3.3), we can write v(xi ) = ∑Nj=0 c j ψ j (xi ) = ci , and therefore, N

v(x) = ∑ v(xi )ψi (x). i=0

For any continuous function f ∈ C0 (Ω), since Ih f ∈ P1 (Th ) and Ih f (xi ) = f (xi ), we conclude that N

Ih f = ∑ f (xi )ψi . i=0

Finite Element Interpolation

21

Thus, the mesh Th induces the interpolation operator Ih : C0 (Ω) → P1 (Th ) that assigns the piecewise linear interpolant Ih f to any continuous function f . In summary, the basis functions ψi on Th allow us to approximate functions. The basis functions ψi are called the global shape functions associated to the global degree of freedom at the vertex xi . The relation (3.3) implies that the value of ψi is 1 at the corresponding vertex and zero at other vertices. A crucial question regarding the Lagrange interpolant Ih f is its ability to approximate f in the following sense: Suppose {Th1 , Th2 , . . . } is a sequence of meshes of Ω such that h1 > h2 > h3 > · · · , and hi → 0 as i → ∞, that is, Thi+1 is more refned than Thi . Let Ih f be the associated Lagrange interpolants of f with h = h1 , h2 , . . . , see Figure 3.3. We want to know if Ih f “converges” to f as h → 0.

Figure 3.3 By increasing the refnement level of meshes, Lagrange interpolants (the dashed lines) provide more accurate approximations of a function (the solid lines).

� � The normed linear spaces L2 (Ω), k · k2 ) and H 1 (Ω), k · k1,2 ) introduced in Example 2.27 provide a suitable framework to answer the above question. Recall that norms can measure the distance between members of linear spaces. Since P1 (Th ) ⊂ H 1 (Ω) ⊂ L2 (Ω), we can interpret the convergence of Ih f to f as h → 0 as follows: Does k f − Ih f k → 0 as h → 0, where k · k can be either k · k2 or k · k1,2 ? We can use the Python based package FEniCS to examine the above question. A brief discussion on FEniCS installation and a short introduction to Python are provided in the appendices. Further references for FEniCS and Python are mentioned at the end of this chapter. The following code calculates the errors k f − Ih f k2 and k f −Ih f k1,2 for f (x) = sin(πx) by using 4 meshes containing 3, 6, 9, and 12 elements of the same size. from dolfin import * # defining the function f = Expression(’sin(pi*x[0])’, degree = 5) # number of divisions of meshes Divisions = [3, 6, 9, 12] # computing interpolants and errors h = [] # max element sizes error_L2, error_H1 = [], [] # initializing errors

22

Finite Element Methods in Civil and Mechanical Engineering

for n in Divisions: # defining the mesh mesh = UnitIntervalMesh(n) # defining the approximation space CGE = FiniteElement("Lagrange", mesh.ufl_cell(), 1) Z = FunctionSpace(mesh, CGE) # interpolating f Interpolant = interpolate(f, Z) # max element size of the mesh h.append(mesh.hmax()) # calculating errors error_L2.append(errornorm(f, Interpolant, norm_type="L2")) error_H1.append(errornorm(f, Interpolant, norm_type="H1")) Let us study the above program. The frst line from dolfin import * imports FEniCS classes such as Expression and UnitIntervalMesh from its DOLFIN library. Alternatively, one may also use from fenics import * FEniCS programs usually start with one of the above lines. The statement f = Expression(’sin(pi*x[0])’, degree = 5) defnes the function f (x) = sin(πx). The formula is expressed as a string with C++ syntax. In 1D, x[0] refers to the x coordinate. In 2D and 3D, the coordinates (x, y) and (x, y, z) are denoted by (x[0],x[1]) and (x[0],x[1],x[2]), respectively. Expressions are interpolated in FEniCS and the degree of polynomials that will be used for this purpose are specifed by the argument degree = 5. The list Divisions = [3, 6, 9, 12] specifes the number of divisions of meshes. In the for loop, interpolants and interpolation errors are calculated for each mesh. A mesh of the unit interval [0, 1] containing n elements of the same size is defned by mesh = UnitIntervalMesh(n) It is easy to obtain mesh attributes. For example, mesh.ufl_cell() and mesh.hmax() respectively denote the shape of the elements and the maximum element size of mesh. The approximation space P1 (Th ), denoted by Z, is then specifed as

23

Finite Element Interpolation

CGE = FiniteElement("Lagrange", mesh.ufl_cell(), 1) Z = FunctionSpace(mesh, CGE) As we will discuss in the next section, this approximation space can be constructed by using Lagrange fnite elements. The frst line assigns a degree 1 Lagrange element to CGE, where as we mentioned earlier, mesh.ufl_cell() denotes the shape of the elements of mesh. The approximation space Z is then defned by using the mesh and the fnite element CGE. The Lagrange interpolant is simply obtained by Interpolant = interpolate(f, Z) Having the function f and Interpolant, the L2 and H 1 errors can be computed using the statement errornorm(f, Interpolant, norm_type) where norm_type is a string specifying the type of error. Figure 3.4 shows the output of the above program: The left panel depicts the interpolants for different meshes and the right panel shows the L2 - and the H 1 -errors versus h.

Figure 3.4 Lagrange Interpolants of type (1) for f (x) = sin(πx) on [0, 1]: (Left) Interpolants for meshes with 3, 6, 9, and 12 uniform elements; (Right) L2 - and H 1 -errors versus h.

It is possible to show that when f is twice differentiable, we have k f − Ih f k2 ≤ C h2 , k f − Ih f k1,2 ≤ C˜ h,

(3.4) (3.5)

where C and C˜ do not dependent on h. These inequalities show that k f − Ih f k → 0 as h → 0, in both the L2 - and the H 1 -norms. The power of h in the inequalities (3.4) and (3.5) is called the convergence rate of the interpolations. Clearly, higher convergence rates imply faster convergence as h → 0. Thus, if f is twice differentiable, then the convergence rate is 2 in the L2 -norm and 1 in the H 1 -norm. To estimate the convergence rates by using the above FEniCS program, we can proceed as follows: Suppose error = C hr , where C is a constant, and let errori = C hri

24

Finite Element Methods in Civil and Mechanical Engineering

and errori−1 = C hri−1 be respectively the errors associated to the mesh hi and the mesh hi−1 . By dividing the expressions for the errors and solving for r, one obtains r=

ln(errori /errori−1 ) . ln(hi /hi−1 )

(3.6)

Notice that the above relation simply implies that the convergence rate is the slope of the log-log diagram of error versus h. By using the outputs h, error_L2, and error_H1 of the above program, we can implement the formula (3.6) as follows: from math import log as ln # log is a dolfin name too rate_L2 = ln(error_L2[-1]/error_L2[-2])/ln(h[-1]/h[-2]) rate_H1 = ln(error_H1[-1]/error_H1[-2])/ln(h[-1]/h[-2]) print(’ L2-convergence rate = %.2f\n H1-convergence rate = %.2f’ % (rate_L2,rate_H1)) For calculating the convergence rates, we used the last two data points that correspond to the meshes with 9 and 12 elements. This program yields the output L2-convergence rate = 2.00 H1-convergence rate = 1.00 This output is consistent with the inequalities (3.4) and (3.5). Also notice that the above convergence rates are the slope of the lines shown in the right panel of Figure 3.4. If f is only one time differentiable, then the convergence rate of Ih f to f is 1 in the L2 -norm and less than 1 in the H 1 -norm. Therefore, the convergence rate depends on the smoothness of f . The convergence rate 2 in the L2 -norm for the Lagrange interpolant of type (1) is called the optimal convergence rate for this type of interpolants in the sense that it is the highest convergence rate due to (3.4). 3.1.2

THE LOCAL LEVEL

Let Th = {K0 , . . . , KN−1 } be a mesh of the unit interval Ω = (0, 1) as discussed in the previous section and suppose Ih f is the Lagrange interpolant of f induced by Th . An interesting feature of Ih f is that it can be constructed locally by using linear interpolants induced by each element. More specifcally, suppose Ii ( f |Ki ) is the Lagrange interpolant of the restriction f |Ki of f on the element Ki = [xi , xi+1 ], see Figure 3.5. The interpolant Ii ( f |Ki ) is called a local interpolant of f in the sense that it is defned only on the element Ki of the mesh Th . We have � Ih f |Ki = Ii ( f |Ki ). Thus, as shown in Figure 3.5, the interpolant Ih f over Th can be constructed by assembling all local interpolants Ii ( f |Ki ) associated to the elements Ki of Th . The local interpolant Ii ( f |Ki ) belongs to the local approximation space P1 (Ki ) that consists of all polynomials of degree ≤ 1 on Ki . Notice that Ii ( f |Ki ) can be

25

Finite Element Interpolation

Figure 3.5 Left: A function f (the solid line) and its Lagrange interpolant Ih f (the dashed line) induced by {x0 , . . . , xN }; Right: The restriction of f to Ki = [xi , xi+1 ] (the solid line) and the local interpolant Ii ( f |Ki ) on Ki (the dashed line).

uniquely determined by the values { f (xi ), f (xi+1 )}. These values at the end points of elements that uniquely specify a member of P1 (Ki ) are called local degrees of freedom in P1 (Ki ). Consider the degree 1 polynomials θi,0 and θi,1 given by θi,0 (xi ) = 1, θi,0 (xi+1 ) = 0, θi,1 (xi ) = 0, θi,1 (xi+1 ) = 1. The polynomials {θi,0 , θi,1 } are called the local shape functions associated to the above local degrees of freedom for P1 (Ki ). Example 2.31 implies that these local shape functions form a basis for P1 (Ki ). One can derive global shape functions (that is, the hat functions of Figure 3.2) by using local shape functions through the fnite element assembly process: For each Ki = [xi , xi+1 ] ∈ Th , we assign the local shape functions θi,0 and θi,1 to xi and xi+1 , respectively. The global function ψi , which is nonzero at the vertex xi of Th , is obtained by noticing that ψi is nonzero only on those elements that contain xi , that is, Ki−1 and Ki . Then, the restriction of ψi to any of these elements is equal to the local shape function assigned to xi . This means that ψi |Ki−1 = θi−1,1 , and ψi |Ki = θi,0 . For example, the global shape function ψ2 in Figure 3.2 is nonzero on K1 and K2 with ψ2 |K1 = θ1,1 , and ψ2 |K2 = θ2,0 . The above discussions imply that the minimum information that we need to locally construct global interpolations over a mesh Th is the triplet (Ki , P1 (Ki ), Σi ) for any Ki ∈ Th , where Σi contains local degrees of freedom. This triplet is called a 1D Lagrange fnite element of type (1). The fact that global results can be obtained by using local information is implemented in the FEniCS code of the previous section by the following segment: CGE = FiniteElement("Lagrange", mesh.ufl_cell(), 1) Z = FunctionSpace(mesh, CGE) Lagrange elements are denoted by Lagrange in FEniCS. The frst line of the above segment defnes CGE as a Lagrange fnite element of type (1), where

26

Finite Element Methods in Civil and Mechanical Engineering

mesh.ufl_cell() refers to Ki , and "Lagrange" and 1 specify the degrees of freedom and the local approximation space. The second line defnes the global approximation space by using the mesh and the fnite element CGE.

3.2

FINITE ELEMENTS

Motivated by the discussions of the previous section, a fnite element is defned to be a triplet (K, P, Σ), where: − The set K, also called element, is a bounded subset of Rd . There is no limitation on the shape of K. For 1D problems, K is a closed interval. In 2D, K is usually a triangle or a rectangle. In 3D, K is usually a tetrahedron, a prism, or a (rectangular) cuboid. − The set P is the local approximation space and is a fnite-dimensional space of polynomials on K that can be Pk (K) or one of its subspaces. − The set Σ is a set containing some information that can uniquely specify any member of P. Mathematically speaking, Σ = {σ1 , . . . , σN } is a basis for the dual space P0 . The members σi , i = 1, . . . , N, are called the local degrees of freedom of the fnite element (K, P, Σ). Notice that the number of degrees of freedom N is the same as dim P. The (local) shape functions or the basis functions of the fnite element (K, P, Σ) are the members of the dual basis of Σ. Therefore, local shape functions {θ1 , . . . , θN } are polynomials belonging to P that satisfy (3.7) σi (θ j ) = δi j , i, j = 1, . . . , N. Given a function f defned on K, the local interpolation operator IK of the fnite element (K, P, Σ) is defned as N

IK f = ∑ σi ( f )θi , i=1

where IK f is also called the P-interpolant of the function f . Since IK f is a linear combination of the local shape functions, it belongs to P. Notice that N

σ j (IK f ) = σ j

�

∑ σi ( f )θi

i=1

N

N

= ∑ σi ( f )σ j (θi ) = ∑ σi ( f )δi j = σ j ( f ), i=1

i=1

that is, the values of degrees of freedom at f and at its P-interpolant are the same. It is straightforward to show that IK p = p, if p ∈ P. We say two fnite elements (K, P, Σ) and (K˜ , P˜ , Σ˜ ) are equal if K = K˜ , P = P˜ , and IK = IK˜ . In the remainder of this section, we introduce a few examples of fnite elements. 3.2.1

SIMPLICIAL LAGRANGE FINITE ELEMENTS OF TYPE (k)

Given a fnite element (K, P, Σ) with dim P = N, suppose there are some points b1 , . . . , bN ∈ K such that σi ( f ) = f (bi ), that is, the degrees of freedom are the values

27

Finite Element Interpolation

of functions at the points {b1 , . . . , bN }. In this case, the points {b1 , . . . , bN } are called the nodes of the fnite element, (K, P, Σ) is called a Lagrange fnite element, and the local shape functions {θ1 , . . . , θN } are also called the nodal basis of P. The associated local interpolation operator is given by IK f = ∑Ni=1 f (bi )θi , that is, Lagrange interpolants are constructed by matching point values at the nodes of Lagrange fnite elements. In this section, a class of Lagrange elements on simplices is introduced. An nsimplex is a subset of Rn defned as follows: A 1-simplex is a closed interval [a0 , a1 ] with vertices {a0 , a1 }. A 2-simplex is a triangle with the vertices {a0 , a1 , a2 }. The face (in 2D also called edge) opposite to the vertex ai is denoted by Fi . The outward unit normal vector of Fi is denoted by ni . A 3-simplex is a tetrahedron with the vertices {a0 , . . . , a3 }. The face Fi and the outward unit normal vector ni are defned similar to 2-simplices. Figure 3.6 shows an n-simplex for n = 1, 2, 3.

Figure 3.6 An n-simplex with the vertices {a0 , . . . , an } for n = 1 (left), n = 2 (center), and n = 3 (right). The face opposite to the vertex ai is denoted by Fi and the vector ni is the outward unit normal vector of Fi .

A useful way for specifying a point x of an n-simplex is through its barycentric coordinates {λ0 , . . . , λn } given by λi (x) x = 1−

(xx − ai ) · ni , (a j − ai ) · ni

(3.8)

where x = (x1 , . . . , xn ), and a j is any vertex in Fi . Notice that λi is the normalized distance of x from the face Fi : We have λi = 1 at ai and λi = 0 at all points of Fi , see Figure 3.7. It is not hard to show that 0 ≤ λi (x) x ≤ 1, and ∑ni=0 λi (x) x = 1. The defnition (3.8) implies that λi (a j ) = δi j , i, j = 0, . . . , n. The point x of an n-simplex with λi (x) x = or the center of gravity of that simplex.

1 n+1 ,

(3.9)

i = 0, . . . , n, is called the barycenter

Example 3.1. The unit 2-simplex is a triangle with the vertices a0 = (0, 0), a1 = (1, 0), and a2 = (0, 1), see Figure 3.7. The barycentric coordinates of a point x = (x1 , x2 ) of the unit 2-simplex are given by λ0 (x1 , x2 ) = 1 − x1 − x2 , λ1 (x1 , x2 ) = x1 , and λ2 (x1 , x2 ) = x2 .

28

Finite Element Methods in Civil and Mechanical Engineering

Figure 3.7 The barycentric coordinate λ0 associated to the vertex a0 = (0, 0) of the unit 2-simplex.

Note that, generally speaking, λi is an affne function of x1 , . . . , xn . Example 3.2. The unit 3-simplex is a tetrahedron with the vertices a0 = (0, 0, 0), a1 = (1, 0, 0), and a2 = (0, 1, 0), and a3 = (0, 0, 1). The defnition (3.8) implies that the barycentric coordinates of a point x = (x1 , x2 , x3 ) of the unit 3-simplex are given by λ0 (x1 , x2 , x3 ) = 1 − x1 − x2 − x3 , and λi (x1 , x2 , x3 ) = xi , i = 1, 2, 3. The simplest Lagrange fnite element on a 1-simplex K = [a0 , a1 ] is the 1D Lagrange fnite element of type (1), also called 1-simplex of type (1), that is introduced in Section 3.1.2. To introduce the simplest Lagrange fnite element on a 2-simplex K, we frst need to show that any polynomial p ∈ P1 (K) can be uniquely determined by its values at the vertices {a0 , a1 , a2 } of K. Notice that dim P1 (K) = 3, and that (3.8) implies x = ∑2i=0 ci λi (x) x = 0, for all x ∈ K. By using that λi ∈ P1 (K), i = 0, 1, 2. Suppose f (x) (3.9), one concludes that f (a j ) = c j . Thus, f (x) x = 0, suggests that c0 = c1 = c2 = 0, x = ∑2i=0 ci λi (x). x Again (3.9) and therefore, {λ0 , λ1 , λ2 } is a basis for P1 (K). Let p(x) implies that p(a j ) = c j . Consequently, any p ∈ P1 (K) can be uniquely determined by its values at the vertices of K with p(x) x = ∑2i=0 p(ai )λi (x). x The simplest Lagrange fnite element on a 2-simplex K with the vertices {a0 , a1 , a2 } is 2-simplex (or triangle) of type (1) (K, P1 (K), ΣK ), where the degrees of freedom ΣK = {σ0 , σ1 , σ2 } are values at the vertices of K given by σi (p) = p(ai ). By abusing the notation, we also write ΣK = {p(a0 ), p(a1 ), p(a2 )}. The relations (3.7) and (3.9) imply that the local shape functions associated to ΣK are {λ0 , λ1 , λ2 }. The local interpolation operator for n-simplex of type (1) with n = 2 is given by n

IK1 f = ∑ f (ai )λi .

(3.10)

i=0

Conventional fnite element diagrams are usually employed to depict fnite elements. These diagrams show the shape of an element and its degrees of freedom. Figure 3.8 depicts the conventional fnite element diagram of 2-simplex of type (1).

Finite Element Interpolation

29

Figure 3.8 The conventional fnite element diagrams of n-simplex of type (1) with n = 1 (left), n = 2 (center), and n = 3 (right).

Similarly, we can defne 3-simplex (or tetrahedron) of type (1) (K, P1 (K), ΣK ), where K is a tetrahedron with vertices {a0 , . . . , a3 } and ΣK = {p(a0 ), . . . , p(a3 )} are values at the vertices of K. The local shape functions are {λ0 , . . . , λ3 } and the associated local interpolation operator is given by (3.10) with n = 3. The conventional fnite element diagram of tetrahedron of type (1) is shown in Figure 3.8. More generally, given an n-simplex K with n = 1, 2, 3, and any integer k ≥ 1, one can defne a Lagrange fnite element (K, Pk (K), ΣK ) called n-simplex of type (k), where Pk (K) is the space of polynomials � of degree k on K and ΣK = {p(b0 ), . . . , p(bN−1 )}, with N = dim Pk (Rn ) = n+k k . The nodes {b0 , . . . , bN−1 } have the barycentric coordinates � j0 jn ,..., , (3.11) k k where j0 , . . . , jn are any choice of non-negative integers such that j0 + · · · + jn = k, and ji < k + 1, i = 0, . . . , n. The local interpolation operator associated to n-simplex of type (k) is denoted by IKk .

Figure 3.9 The conventional fnite element diagrams of n-simplex of type (2) with n = 1 (left), n = 2 (center), and n = 3 (right). Only visible degrees of freedom are shown.

Example 3.3. Let n = 1 and k = 2. Then, the barycentric coordinates of the nodes of 1-simplex of type (2) are given by b0 = (1, 0), b1 = (0, 1), and b2 = ( 12 , 12 ). The conventional fnite element diagram of this element is shown in Figure 3.9. The shape x = λ0 (x) x (1 − 2λ1 (x)), x θ1 (x) x = λ1 (x) x (1 − 2λ0 (x)), x functions can be written as θ0 (x) and θ2 (x) x = 4λ0 (x)λ x 1 (x). x These shape functions are shown in Figure 3.10. Example 3.4. The barycentric coordinates of the nodes of 2-simplex of type (2) are b0 = (1, 0, 0),

b1 = (1/2, 1/2, 0), b2 = (0, 1, 0), b3 = (1/2, 0, 1/2), b4 = (0, 1/2, 1/2), b5 = (0, 0, 1).

30

Finite Element Methods in Civil and Mechanical Engineering

Figure 3.10 Local shape functions of 1-simplex of type (2) on K = [0, 1].

Figure 3.9 shows the conventional fnite element diagram of this fnite element. It is not hard to show that the associated local shape functions are given by θ0 = λ0 (2λ0 − 1), θ1 = 4λ0 λ1 , θ2 = λ1 (2λ1 − 1), θ4 = 4λ1 λ2 , θ5 = λ2 (2λ2 − 1). θ3 = 4λ0 λ2 , Simplicial Lagrange fnite elements are available in FEniCS. For examples, the following code defnes n-simplex of type (2) for n = 1, 2, 3. FE_1simplex = FiniteElement("Lagrange", "interval", 2) FE_2simplex = FiniteElement("Lagrange", "triangle", 2) FE_3simplex = FiniteElement("Lagrange", "tetrahedron", 2) To defne simplicial Lagrange elements in FEniCS, one can also use "CG" instead of "Lagrange" in the above code. 3.2.2 SIMPLICIAL HERMITE FINITE ELEMENTS OF TYPE (3) Suppose degrees of freedom of a fnite element (K, P, Σ) are values of functions at some points {b1 , . . . , bv } ⊂ K together with values of partial derivatives (or directional derivatives) of functions at some points {b¯ 1 , . . . , b¯ d } ⊂ K. The fnite element (K, P, Σ) is then called a Hermit fnite element and the points {b1 , . . . , bv , b¯ 1 , . . . , b¯ d } are called the nodes of the fnite element. In the following, we defne simplicial Hermite fnite elements that involve polynomials of degree 3. Since we will not use 2D and 3D Hermite elements in the remainder of this book, we only discuss the 1D case in some details. Let K = [a0 , a1 ] be a 1-simplex. The Hermite 1-simplex of type (3) is given by (K, P3 (K), ΣK ) with ΣK = {p(a0 ), p0 (a0 ), p(a1 ), p0 (a1 )}, where p0 (x) is the derivative of p(x) with respect to x. The conventional fnite element diagram of this element is shown in Figure 3.11, where circles centered at vertices indicate degrees of freedom associated to derivatives at those vertices. It is not hard to show that the local

31

Finite Element Interpolation

Figure 3.11 The conventional fnite element diagrams of the Hermite n-simplex of type (3) with n = 1 (left), n = 2 (center), and n = 3 (right). Only visible degrees of freedom are shown.

shape functions can be written as 1 (x − a1 )2 (h + 2(x − a0 )) , θ1 (x) = h3 1 θ2 (x) = 3 (x − a0 )2 (h − 2(x − a1 )) , θ3 (x) = h θ0 (x) =

1 (x − a0 )(x − a1 )2 , h2 1 (x − a1 )(x − a0 )2 , h2

(3.12)

where h = a1 − a0 . Figure 3.12 shows these local shape functions on the unit interval.

Figure 3.12 Local shape functions of the Hermite 1-simplex of type (3) on K = [0, 1].

Suppose K is a 2-simplex with the vertices {a0 , a1 , a2 }. The Hermite 2-simplex (or triangle) of type (3) is (K, P3 (K), ΣK ), where ΣK includes values of functions and all their frst-order partial derivatives at the vertices together with the value of functions at the node 13 (a0 + a1 + a2 ). The Hermite 3-simplex (or tetrahedron) of type (3) can be defned in a similar way. The conventional fnite element diagram of these elements are shown in Figure 3.11. Hermite elements are still not available in FEniCS. 3.2.3 THE RAVIART-THOMAS FINITE ELEMENT Let K be a 2-simplex and suppose RT is the space of 2D polynomial vector felds of the form pp(x1 , x2 ) = (c0 + c2 x1 , c1 + c2 x2 ), where c0 , c1 , and c2 are constant. It is not hard to show that RT is a linear space with dim RT = 3. Referring to Figure 3.6, for any p ∈ RT, let σi ( pp) be the “fux” of p across the face Fi , that is, Z

σi ( pp) =

Fi

p · ni ds.

(3.13)

32

Finite Element Methods in Civil and Mechanical Engineering

One can show that the values of these 3 fuxes uniquely specify polynomials belonging to RT. Therefore, (K, RT, ΣK ) with ΣK = {σ0 , σ1 , σ2 }, is a fnite element which is called the (2D) Raviart-Thomas fnite element. The conventional fnite element diagram of this element is shown in Figure 3.13. To motivate the defnition of the local shape functions of this element, let us consider the following example.

Figure 3.13 The conventional fnite element diagram of the 2D (left) and the 3D (right) Raviart-Thomas fnite elements. Arrows normal to faces denote the degrees of freedom associated to faces. Only visible degrees of freedom are shown.

Example 3.5. Let x = (x1 , x2 ), and consider the vector feld θ 0 ∈ RT given by θ 0 (x1 , x2 ) =

1 (xx − a0 ), 2 AK

where AK and a0 are respectively the area and a vertex of the 2-simplex K. Using the notation of Figure 3.6, any point x on the face F0 can be denoted as x = a1 + t(a2 − a1 ), 0 ≤ t ≤ 1. Thus, on F0 we can write θ 0 (x) x · n0 =

1 (a1 − a0 ) · n0 , (a1 + t(a2 − a1 ) − a0 ) · n0 = 2 AK 2 AK

where we used the fact that a2 − a1 is normal to n0 . Since (a1 − a0 ) · n0 is equal to x normal the height h0 of K normal to F0 , we conclude that the component of θ 0 (x) to F0 is constant and is given by θ 0 (x) x · n0 = 1/|F0 |, where |F0 | is the length of F0 . Consequently, one concludes that Z

σ0 (θ θ 0) =

F0

θ 0 · n0 ds =

1 |F0 |

Z

ds = 1. F0

Similarly, one can show that σ1 (θ θ 0 ) = σ2 (θ θ 0 ) = 0. The above example suggests that the local shape functions of the 2D RaviartThomas fnite element are given by θ i (x1 , x2 ) =

1 (xx − ai ), 2 AK

i = 0, 1, 2.

Figure 3.14 shows these local shape functions. The local interpolation operator of the 2D Raviart-Thomas fnite element can be written as 2 Z IKRT v = ∑ v · ni ds θ i . i=0

Fi

33

Finite Element Interpolation

Figure 3.14 The local shape functions of the 2D Raviart-Thomas fnite element.

The 3D Raviart-Thomas fnite element is obtained by using local degrees of freedom similar to (3.13). Figure 3.13 shows the conventional fnite element diagram of this 3D fnite element. More generally, it is also possible to defne 2D and 3D Raviart-Thomas fnite elements of degree k, for any integer k ≥ 1. In FEniCS, the 2D and 3D Raviart-Thomas fnite elements (of degree 1) can be defned as follows: RT_2D = FiniteElement("RT", "triangle", 1) RT_3D = FiniteElement("RT", "tetrahedron", 1) 3.2.4 THE NEDELEC FINITE ELEMENT Suppose K is a 2-simplex and let ti be one of the two unit vectors parallel to the edge Ei of K (in 2D, edges and faces are the same). Let us denote the space of all polynomial vector felds of the form pp(x1 , x2 ) = (c0 + c2 x2 , c1 − c2 x1 ), with c0 , c1 , and c2 being constant, by NE. One can show that NE is a linear space with dim NE = 3. For any p ∈ NE, suppose σi ( pp) is the integral of the tangent component of p along Ei , that is, Z σi ( pp) =

Ei

p · ti ds.

(3.14)

The 2D Ned ´ e´ lec fnite element (of the frst kind) (K, NE, ΣK ) is obtained by assuming ΣK = {σ0 , σ1 , σ2 }. The conventional fnite element diagram of this element is shown in Figure 3.15.

Figure 3.15 The conventional fnite element diagram of the 2D (left) and the 3D (right) ´ e´ lec fnite elements. Arrows parallel to edges denote the degrees of freedom associated to Ned edges. Only visible degrees of freedom are shown.

The local shape functions of this element can be defned as follows: Let R(x1 , x2 ) = (x2 , −x1 ), that is, R(x1 , x2 ) is obtained by a 90 degree rotation of x = (x1 , x2 ) in the clockwise direction. Then, by using an approach similar to that of

34

Finite Element Methods in Civil and Mechanical Engineering

Example 3.5, one can show that the local shape functions are given by θ i (x1 , x2 ) =

R(xx − ai ) �1 , |Ei |ti · R 2 (ai1 + ai2 ) − ai

i = 0, 1, 2,

where i1 and i2 are not equal to i. Figure 3.16 depicts these local shape functions. ´ e´ lec fnite element reads The local interpolation operator of the 2D Ned IKN v

2

=∑

i=0

Z Ei

v · ti ds θ i .

Figure 3.16 The local shape functions of the 2D N´ed´elec fnite element.

Local degrees of freedom of the 3D N´ed´elec fnite element is defned similar to (3.14). Figure 3.15 shows the conventional fnite element diagram of this 3D fnite element. More generally, one can defne 2D and 3D N´ed´elec fnite elements of degree k, for any integer k ≥ 1. In FEniCS, the 2D and 3D N´ed´elec fnite elements (of degree 1) can be defned by the following code: ND_2D = FiniteElement("N1curl", "triangle", 1) ND_3D = FiniteElement("N1curl", "tetrahedron", 1)

3.3 MESHES To construct approximation spaces over a domain Ω by using fnite elements, one needs to employ a mesh of Ω. More specifcally, let Ω be an open, bounded subset of Rn with the continuous boundary ∂ Ω and let Ω = Ω ∪ ∂ Ω. In this book, domains are usually assumed to be polygons for 2D problems and polyhedrons for 3D problems. A mesh of Ω is a set {K1 , . . . , Km }, where 1. Each Ki , called an element or a cell of the mesh, is usually assumed to be a polygonal domain in 2D or a polyhedral domain in 3D that contains its boundary; 2. Ω = K1 ∪ K2 ∪ · · · ∪ Km ; 3. Two different elements either do not intersect or their intersection is a common vertex in 1D, a common vertex or edge in 2D, or a common vertex, edge, or face in 3D. Later, we will see this property of a mesh is important for obtaining approximation spaces over Ω.

Finite Element Interpolation

35

Meshes defned as above are also called geometrically conformal meshes. We consider meshes with simplicial elements in this book, that is, elements are triangles in 2D and tetrahedrons in 3D. These meshes are called triangulations as well. For example, the left panel of Figure 3.17 shows a mesh of a square. The confguration on the right panel of Figure 3.17 is not a mesh since the intersection of each of the two smaller triangles with the bigger one is not an edge of the bigger triangle. We denote the number of vertices, edges, faces, and elements of a mesh by nv , ne , nf , and nel , respectively. Notice that in 2D, edges and faces are the same.

Figure 3.17 An example (left) and a counterexample (right) of a mesh.

Given a mesh {K1 , . . . , Km }, let hKi = diam Ki , be the diameter of the element Ki , that is, hKi is the maximum distance of two different points of Ki . The mesh is usually denoted by Th , where h = max{hK1 , . . . , hKm }. The parameter h refers to the refnement level of the mesh. To numerically approximate solutions, one usually considers a sequence Th1 , Th2 , Th3 , . . . of successively refned meshes of a domain with h1 > h2 > h3 > · · · . This sequence is denoted by {Th }h>0 or simply {Th }, see Figure 3.18.

Figure 3.18 A sequence of meshes of a square.

Meshes may be generated by using a reference element. For example, consider the mesh Th = {K1 , K2 , K3 } shown in Figure 3.19 and let Kˆ be an arbitrary triangle. There are affne transformations fi : Kˆ → Ki , i = 1, 2, 3, that map Kˆ onto Ki with x = vi + Mi xx, where vi is a vector in R2 and Mi is a 2 × 2 matrix. Therefore, the fi (x) affne transformations { f1 , f2 , f3 } allow to generate the mesh Th from the reference element Kˆ . Such meshes are also called affne meshes.

Figure 3.19 An affne mesh can be generated from a reference element Kˆ .

36

Finite Element Methods in Civil and Mechanical Engineering

Notice that if Ω has a curved boundary, it is not possible to completely cover it by meshes that only contain elements with fat faces. In this case, one can use meshes with fat-face elements to approximately cover Ω. One can improve this approximation by using fner meshes at curved boundaries. Another approach is to use meshes with curved-face elements. We will not consider curved-face elements in this book. Meshes of some simple geometries can be directly defned in FEniCS. For example, the following code generates structured meshes of the unit square and the unit cube: Mesh_S = UnitSquareMesh(4,4) Mesh_C = UnitCubeMesh(4,4,4) The numbers inside parenthesis show the number of divisions of each edge. Figure 3.20 shows these meshes.

Figure 3.20 Structured meshes of the unit square and the unit cube.

Unstructured meshes and meshes for more complex geometries can be defned by the mesh generation component of FEniCS called mshr. The following segment generates the unstructured meshes of Figure 3.21. from mshr import * # a rectangle with a rectangular hole Rec1 = Rectangle(Point(0, 0), Point(1, 1)) Rec2 = Rectangle(Point(0.25, 0.25), Point(0.75, .75)) mesh_Rec = generate_mesh(Rec1 - Rec2, 8) # a triangle vertices = [Point(0.0, 0.0), Point(1.0, 0.0), Point(0.5, 1.0)] Tri = Polygon(vertices) mesh_Tri = generate_mesh(Tri, 6) # the unit cube Cube = Box(Point(0, 0, 0), Point(1, 1, 1)) mesh_Cube = generate_mesh(Cube, 7.5) # the unit sphere centered at the origin Spr = Sphere(Point(0, 0, 0), 1.0) mesh_Spr = generate_mesh(Spr, 6)

Finite Element Interpolation

37

Notice that the second argument of the function generate_mesh determines the refnement level of unstructured meshes and does not need to be an integer. Another option for defning meshes in FEniCS is to import meshes that are generated by other programs outside FEniCS.

Figure 3.21 Unstructured meshes of some 2D and 3D domains.

It is easy to recover basic properties of meshes in FEniCS. The following code prints some properties of the mesh of the unit square shown in Figure 3.20. n_el = Mesh_S.num_cells() # number of elements n_f = Mesh_S.num_faces() # number of faces n_v = Mesh_S.num_vertices() # number of vertices h = Mesh_S.hmax() # the maximum diameter of elements print("n_el = %.f, n_f = %.f, n_v = %.f, h = %.3f" % (n_el, n_f, n_v, h)) The output of this code is: n_el = 32, n_f = 32, n_v = 25, h = 0.354

3.4 FINITE ELEMENT SPACES AND INTERPOLATIONS In this section, we use fnite elements to construct approximation spaces. Such approximation spaces are called fnite element spaces. We begin this section by discussing this construction in its full generality suitable for any fnite element. We then apply this construction to fnite elements introduced earlier. The reader is encouraged to read this general discussion again while studying the construction for each specifc element. Let Th = {K1 , . . . , Km } be a mesh of a domain Ω and let {(K, PK , ΣK )}K∈Th be a family of fnite elements associated to the elements of Th . The fnite element space induced by this family of fnite elements is constructed through the fnite element assembly process as follows: 1. Each local degree of freedom of each fnite element (K, PK , ΣK ) is assigned to a geometric entity of the mesh Th such as a vertex, a node, an edge, or a face of Th . Notice that several local degrees of freedom may be assigned to a geometric entity.

38

Finite Element Methods in Civil and Mechanical Engineering

2. The fnite element space Vh associated to Th is defned as the linear space of functions f such that f |K ∈ PK , for all K ∈ Th , and σK,i ( f |K ) = σK 0 , j ( f |K 0 ), if the i-th local degree of freedom σK,i of K and the j-th local degree of freedom σK 0 , j of K 0 are assigned to the same geometric entity of Th . The members of any basis for the linear space Vh are called global shape functions. An important property of the fnite element space Vh is that it admits locally supported global shape functions, that is, bases functions that are zero everywhere on Th except for a few elements. These global shape functions can be constructed by using local shape functions of fnite elements. Before discussing this construction, let us show that fnite element spaces are fnite-dimensional. Let nt be the total number of the geometric entities of Th that were assigned to at least one local degree of freedom in the fnite element assembly process. For any v ∈ Vh , all local degrees of freedom assigned to the same geometric entity of Th take the same value at v. Consequently, any v ∈ Vh can be uniquely specifed by using {g1 , . . . , gnt }, where gi (v) assigned to the i-th geometric entity is the common value of all local degrees of freedom assigned to the i-th geometric entity at v. This argument suggests that the linear space Vh is a fnite-dimensional space with dimVh = nt . One can obtain a set of locally supported global shape functions {ψ1 , . . . , ψnt }, which is the dual basis of {g1 , . . . , gnt } in the sense that gi (ψ j ) = δi j , as follows: Suppose that only 2 local degrees of freedom σK, j and σK 0 ,l are associated to the i-th geometric entity of Th . Let the local shape functions θK, j on K and θK 0 ,l on K 0 satisfy σK, j (θK, j ) = σK 0 ,l (θK 0 ,l ) = 1. Then, ψi is equal to: θK, j on K; θK 0 ,l on K 0 ; and zero on other elements of Th . The same approach also holds if the i-th geometric entity is shared by more than 2-elements. It is straightforward to show the duality relation gi (ψ j ) = δi j . This relation implies that {ψ1 , . . . , ψnt } is a basis for Vh . The (global) interpolation Ih f ∈ Vh of a function f on Ω is then defned as nt

Ih f = ∑ gi ( f )ψi .

(3.15)

i=1

It is not hard to show that on any element K ∈ Th , the global interpolant Ih f satisfes (Ih f )|K = IK ( f |K ), where IK is the local interpolation operator of the fnite element (K, PK , ΣK ). If the fnite element space Vh is a subset of a linear space V , we say Vh is a V conformal space. In the remainder of this section, we will employ the fnite elements introduced in Section 3.2 to derive H 1 -, H(div)-, and H(curl)-conformal fnite element spaces. 3.4.1

H 1 -CONFORMAL FINITE ELEMENT SPACES

Recall that the Sobolev space H 1 (Ω) introduced in Example 2.8 is the space of square integrable functions that also admit square integrable frst-order derivatives. Given a mesh Th of Ω, suppose a function f : Ω → R has the property that f |K ∈ C1 (K), for

Finite Element Interpolation

39

all K ∈ Th . One can show that f ∈ H 1 (Ω) if and only if f ∈ C0 (Ω). In the remainder of this section, we obtain fnite element spaces associated to the simplicial Lagrange fnite elements of type (k) and the simplicial Hermite fnite element of type (3). We use the above result to show that these fnite element spaces are H 1 -conformal. 3.4.1.1

Lagrange Elements

Let Th be a mesh of Ω with simplicial elements and suppose {(K, Pk (K), ΣK )}K∈Th is a family of simplicial Lagrange fnite elements. To obtain the associated fnite element space Pk (Th ) through the fnite element assembly process, we frst assign local degrees of freedom of each K to the corresponding nodes of Th . Then, Pk (Th ) is defned to be the space of functions f such that f |K ∈ Pk (K), for all K ∈ Th , and f |K (bi ) = f |K 0 (bi ), whenever bi is a common node of K and K 0 . Let face of an n-simplex K. We have dim Pk (K) = �n+k F be an (n − 1)-dimensional �n+k−1 , and dim P (F) = . Notice that if p ∈ Pk (K), then p|F ∈ Pk (F). It is not k k k � hard to show that n-simplex of type (k) has n+k−1 nodes on each face of K. Now, k suppose F is a common face of K and K 0 , and let f ∈ Pk (Th ). Then, the requirement that f |K and f |K 0 have the same values at all common nodes on F implies that f |K = f |K 0 on F. Hence, we conclude that Pk (Th ) = { f : f |K ∈ Pk (K), for all K ∈ Th } ∩C0 (Ω). The result mentioned at the beginning of this section then suggests that Pk (Th ) is an H 1 -conformal fnite element space. Let us consider a few examples. Example 3.6. Let Th = {K0 , . . . , KN−1 } be a mesh of an open interval Ω with the vertices {x0 , . . . , xN } as discussed in Section 3.1.1. Also let {(K, P1 (K), ΣK )}K∈Th be a family of 1-simplices of type (1). The conventional fnite element diagram shown in Figure 3.8 tells us that local degrees of freedom are associated to the two vertices of each K. Therefore, two local degrees of freedom are assigned to the internal vertices {x1 , . . . , xN−1 } and only one local degrees of freedom is assigned to x0 and xN . Then, the associated fnite element space P1 (Th ) is the space of continuous, piecewise affne functions over Th defned in (3.1), with dim P1 (Th ) = N + 1, as local degrees of freedom are assigned to N + 1 vertices of Th . Locally supported global shape functions {ψ0 , . . . , ψN } of P1 (Th ) are as follows: The basis function ψi associated to the vertex xi is nonzero only at Ki−1 and Ki , where it is equal to the local shape functions assigned to xi , see Figure 3.2. Notice that in Section 3.1, the global shape functions for P1 (Th ) were obtained similarly. Example 3.7. Let Th be a mesh of a 2D domain and suppose {(K, P2 (K), ΣK )}K∈Th is a family of 2-simplices of type (2). To obtain the associated fnite element space P2 (Th ), the conventional fnite element diagram of 2-simplex of type (2) shown in Figure 3.9 tells us that local degrees of freedom should be assigned to the vertices of Th and the nodes at the middle of the edges of Th , see Figure 3.22. The dimension

40

Finite Element Methods in Civil and Mechanical Engineering

Figure 3.22 Nodes of a 2D mesh that are used in the fnite element assembly process associated to 2-simplices of type (k) with k = 2 (left) and k = 3 (right).

of P2 (Th ) is the number of nodes of Th that are assigned a local degree of freedom, that is, dim P2 (Th ) = nv + ne , where nv and ne are respectively the number of vertices and the number of edges of Th . Locally supported global shape functions {ψ1 , . . . , ψnv +ne } are obtained by using local shape functions. Figure 3.23 shows regions where global shape functions associated to some nodes of the mesh of Figure 3.22 are not zero.

Figure 3.23 Supports of global shape functions induced by 2-simplices of type (2): The global shape function associated to the specifed node is not zero on the shaded region.

Example 3.8. Consider a family of 2-simplices of type (3) on a 2D mesh. The dimension of the associated fnite element space is nv + 2ne + nel , where nel is the number of elements of the mesh, see Figure 3.22. Notice that the interelement continuity of functions belonging to fnite element spaces is closely related to the third property of meshes discussed in Section 3.3. For example, the fnite element assembly process does not result in continuous functions if one considers a family of 2-simplices of type (1) on the confguration on the right panel of Figure 3.17, see also Figure 3.24.

Figure 3.24 The global shape function for a family of 2-simplices of type (1) associated to the vertex b of the counterexample on the right panel of Figure 3.17 is not continuous.

By using locally supported global shape functions {ψ1 , . . . , ψnt } and the general relation (3.15), one can write the (global) Lagrange interpolation operator Ihk

41

Finite Element Interpolation

induced by n-simplices of type (k) as nt

Ihk f = ∑ f (bi )ψi , i=1

where b1 , . . . , bnt are the nodes of the underlying mesh that have been used in the assembly process. In Section 3.1.1, we saw how to employ FEniCS to interpolate functions using 1D Lagrange elements. Interpolations in higher dimensions are very similar. For example, the following code computes Lagrange interpolants of type (1) for f (x, y) = 1 + 2x y2 by using structured and unstructured meshes of the unit square. # defining the function f = Expression(’1.0 + pow(2,x[0])*pow(x[1],2)’, degree = 5) # defining meshes mesh_St = UnitSquareMesh(5,5) # structured mesh Rec = Rectangle(Point(0, 0), Point(1, 1)) mesh_Unst = generate_mesh(Rec, 5) # unstructured mesh # defining approximation spaces CGE = FiniteElement("Lagrange", "triangle", 1) Z_St = FunctionSpace(mesh_St, CGE) Z_Unst = FunctionSpace(mesh_Unst, CGE) # interpolating Interpolant_St = interpolate(f, Z_St) # on structured mesh Interpolant_Unst = interpolate(f, Z_Unst) # on unstructured mesh # saving results in VTK format vtkfile = File(’Data/CGInterpolation_2DSt.pvd’) vtkfile 0 . In the following, we assume that Ω is a polygon or a polyhedron and that {Th }h>0 is a shape-regular family of affne meshes, where shape-regular means there is a constant c independent of h such that hK ≤ c, for all K ∈ Th and all Th , dK where hK is the diameter of K and dK is the diameter of the largest ball inscribed in K. One can show that hdKK ≤ sin2αK , where αK is the smallest angle of K. Thus, {Th }h>0 cannot be shape-regular if elements of Th become too “fat” in the sense that αK → 0 as h → 0, see Figure 3.28.

Figure 3.28 A family of meshes cannot be shape-regular if elements become too fat as h → 0.

47

Finite Element Interpolation

To study the convergence of the Lagrange interpolation operator Ihk , we can use the normed linear spaces (L2 (Ω), k · k2 ) and (H 1 (Ω), k · k1,2 ) as the associated fnite element spaces are H 1 -conformal. In particular, suppose a function f is s + 1 times differentiable with 1 ≤ s ≤ k. Then, one can show that k f − Ihk f k2 ≤ C( f ) hs+1 , k f − Ihk f k1,2 ≤ C˜ ( f ) hs ,

(3.16)

where C( f ) and C˜ ( f ) do not depend on h. The above inequalities show that k f − Ihk f k2 and k f − Ihk f k1,2 tend to zero as h → 0. The power of h in the right side of the above inequalities is called the convergence rate. Notice that the convergence rate depends on the smoothness of f and the norm being used. The maximum convergence rate of Ihk is k + 1 in the L2 -norm and k in the H 1 -norm. The maximum convergence rate in the L2 -norm is called the optimal convergence rate. Example 3.9. If f is not suffciently differentiable, increasing the degree of polynomials k will not necessarily increase the convergence rate. For example, if f is only twice differentiable, the convergence rates of Ih1 f and Ih2 f will be 2 in the L2 -norm. Let v be a 2D or a 3D vector feld with twice differentiable components. Then, one can write kvv − IhRT vvk2 ≤ kvv − IhRT vvkd ≤ D(v) v h, (3.17) N N kvv − Ih vvk2 ≤ kvv − Ih vvkc ≤ D˜ (v) v h, where D(v) v and D˜ (v) v do not depend on h. Notice that since IhRT v ∈ H(div; Ω) and N Ih v ∈ H(curl; Ω), the errors kvv − IhRT vvkd and kvv − IhN vvkc are well-defned. Therefore, the convergence rates of interpolations of twice differentialble vector felds using the Raviart-Thomas and the N´ed´elec fnite elements (of degree 1) is 1 in the L2 , H(div; Ω), and H(curl; Ω) norms. Following the program of Section 3.1.1 and using (3.6), it is straightforward to compute errors and convergence rates in FEniCS. As an example, the following code computes the convergence rates of various interpolations of the vector feld � (3.18) vv(x, y) = (x2 + y2 )b , xy , with b = 1.25, on the unit square. # defining the vector field b_v = 1.25 # parameter of the vector field vf = Expression(("pow(x[0]*x[0]+x[1]*x[1],b)", "x[0]*x[1]"), b = b_v, degree = 5) # number of divisions of meshes Divisions = [3, 6, 9, 12] # computing interpolants and errors

48

Finite Element Methods in Civil and Mechanical Engineering

h = []

# max element sizes

# initializing errors error_L2_Lag1, error_H1_Lag1 = [], [] error_L2_Lag2, error_H1_Lag2 = [], [] error_L2_RT, error_Hdiv = [], [] error_L2_Ned, error_Hcurl = [], []

# # # #

Lagrange type (1) Lagrange type (2) RT Nedelec

for n in Divisions: # defining the mesh mesh = UnitSquareMesh(n,n) # approximation spaces Z_Lag1 = VectorFunctionSpace(mesh, "Lagrange", 1) Z_Lag2 = VectorFunctionSpace(mesh, "Lagrange", 2) Z_RT = FunctionSpace(mesh, "RT", 1) Z_Ned = FunctionSpace(mesh, "N1curl", 1) # interpolants I_Lag1 = interpolate(vf, Z_Lag1) I_Lag2 = interpolate(vf, Z_Lag2) I_RT = interpolate(vf, Z_RT) I_Ned = interpolate(vf, Z_Ned) # max element size of the mesh h.append(mesh.hmax()) # calculating errors error_L2_Lag1.append(errornorm(vf, I_Lag1, norm_type="L2")) error_H1_Lag1.append(errornorm(vf, I_Lag1, norm_type="H1")) error_L2_Lag2.append(errornorm(vf, I_Lag2, norm_type="L2")) error_H1_Lag2.append(errornorm(vf, I_Lag2, norm_type="H1")) error_L2_RT.append(errornorm(vf, I_RT, norm_type="L2")) error_Hdiv.append(errornorm(vf, I_RT, norm_type="Hdiv")) error_L2_Ned.append(errornorm(vf, I_Ned, norm_type="L2")) error_Hcurl.append(errornorm(vf, I_Ned, norm_type="Hcurl")) # computing convergence rates from math import log as ln # log is a dolfin name too # Lagrange of type 1 rate_L2_Lag1 = ln(error_L2_Lag1[-1]/error_L2_Lag1[-2])/ln(h[-1]/h[-2]) rate_H1_Lag1 = ln(error_H1_Lag1[-1]/error_H1_Lag1[-2])/ln(h[-1]/h[-2]) # Lagrange of type 2

49

Finite Element Interpolation

rate_L2_Lag2 = ln(error_L2_Lag2[-1]/error_L2_Lag2[-2])/ln(h[-1]/h[-2]) rate_H1_Lag2 = ln(error_H1_Lag2[-1]/error_H1_Lag2[-2])/ln(h[-1]/h[-2]) # RT rate_L2_RT = ln(error_L2_RT[-1]/error_L2_RT[-2])/ln(h[-1]/h[-2]) rate_Hdiv = ln(error_Hdiv[-1]/error_Hdiv[-2])/ln(h[-1]/h[-2]) # Nedelec rate_L2_Ned = ln(error_L2_Ned[-1]/error_L2_Ned[-2])/ln(h[-1]/h[-2]) rate_Hcurl = ln(error_Hcurl[-1]/error_Hcurl[-2])/ln(h[-1]/h[-2]) print(’Convergence Rates:’) print(’Lagrange of type (1) => L2 = %.2f, % (rate_L2_Lag1,rate_H1_Lag1)) print(’Lagrange of type (2) => L2 = %.2f, % (rate_L2_Lag2,rate_H1_Lag2)) print(’Raviart-Thomas => L2 = %.2f, % (rate_L2_RT,rate_Hdiv)) print(’Nedelec => L2 = %.2f, % (rate_L2_Ned,rate_Hcurl))

H1

= %.2f’

H1

= %.2f’

Hdiv

= %.2f’

Hcurl = %.2f’

The above code works very similar to that of Section 3.1.1. A few comments are in order here. First, notice that following (3.18), the vector feld vf is defned with a parameter b. This parameter is initialized as b = b_v, which is passed as an argument to Expression. After vf is defned, one can access the parameter b of vf by vf.b. For example, vf.b = 2.5, assigns the value 2.5 to b. Also notice that errornorm can compute various norms, which are specifed by the variable norm_type. The output of the above code is: Convergence Rates: Lagrange of type (1) Lagrange of type (2) Raviart-Thomas Nedelec

=> => => =>

L2 L2 L2 L2

= = = =

2.00, 2.99, 1.00, 1.00,

H1 H1 Hdiv Hcurl

= = = =

1.00 1.99 1.00 1.00

These results are consistent with the theoretical values of convergence rates mentioned earlier. From a purely mathematical point of view, higher convergence rates are preferable as they imply faster convergence for small h. In engineering applications, however, it may be computationally very expensive to consider very fne meshes and one may consider only a few meshes in the family {Th }. In such cases where h is not necessarily “small” enough, the inequalities (3.16) and (3.17) suggest that errors also depend on the constants C, C˜ , D, and D˜ . For example, consider Figure 3.29 that shows error versus h diagrams for two cases. For h = h1 , h2 , the errors associated to the case with the higher convergence rate (the dashed curve) are larger than the ones associated to the lower rate case (the solid curve). Consequently, convergence rates may not be of direct practical interest. Nonetheless, as we will discuss in the following chapters,

50

Finite Element Methods in Civil and Mechanical Engineering

computing convergence rates can be useful for debugging programs and for studying the performance of numerical methods.

Figure 3.29 For meshes with not suffciently small h, higher convergence rates do not necessarily imply smaller interpolation errors.

The approach that we considered in this section to study convergence by using a sequence of successively refned meshes is called the h-type approach. Alternatively, we may consider a fxed mesh and study the convergence by increasing the degree of polynomials of fnite elements. This is called the p-type approach. We only consider the h-type approach in this book.

EXERCISES Exercise 3.1. Show that {ψ0 , ψ1 , . . . , ψN } is a basis for P1 (Th ), where the functions ψi are defned in (3.2). Exercise 3.2. Given a fnite element (K, P, Σ), show that if p ∈ P, then IK p = p, that is, the local interpolation operator IK does not change members of P. Exercise 3.3. Verify the barycentric coordinates of the unit 2- and 3-simplices given in Examples 3.1 and 3.2. Exercise 3.4. Let K be a tetrahedron. Show that any p ∈ P1 (K) can be uniquely determined by its values at the vertices of K. Exercise 3.5. For n = 1, 2, 3, and k = 1, verify that the defnition of n-simplex of type (k) based on (3.11) results in n-simplex of type (1) introduced in Sections 3.1.2 and 3.2.1. Exercise 3.6. Show that local shape functions of 1- and 2-simplices of type 2 are those given in Examples 3.3 and 3.4. Exercise 3.7. Write the barycentric coordinates of the nodes of 3-simplex of type (2) and derive the associated local shape functions. The conventional fnite element diagram of this element is shown in Figure 3.9. Exercise 3.8. The 2D Crouzeix-Raviart fnite element is a Lagrange fnite element (K, P1 (K), ΣK ), where K is a 2-simplex and the nodes of this element have the

Finite Element Interpolation

51

barycentric coordinates b0 = (0, 1/2, 1/2), b1 = (1/2, 0, 1/2), and b2 = (1/2, 1/2, 0). Plot the conventional fnite element diagram of this element and derive its local shape functions in terms of barycentric coordinates. Is this element equal to 2-simplex of type (1)? Exercise 3.9. Show that the polynomials given in (3.12) are the local shape functions of the Hermite 1-simplex of type (3). Exercise 3.10. Show that RT, the underlying polynomial space of the RaviartThomas fnite element, is a linear space with dim RT = 3. θ 0 ) = σ2 (θ θ 0 ) = 0. Exercise 3.11. In Example 3.5, show that σ1 (θ Exercise 3.12. Show that NE, the underlying polynomial space of the N´ed´elec fnite element, is a linear space with dim NE = 3. Exercise 3.13. For the N´ed´elec fnite element, show that σi (θ θ ) j = δi j , where the degrees of freedom σi and the local shape functions θ i , i = 1, 2, 3, are introduced in Section 3.2.4. Exercise 3.14. Suppose Ih is the global interpolation operator associated to a mesh Th and the fnite element space Vh . Show that if f ∈ Vh , then Ih f = f . Exercise 3.15. Show that the set of the functions {ψ1 , . . . , ψnv +ne } introduced in Example 3.7 is a basis for P2 (Th ). Exercise 3.16. Determine the dimension of the fnite element space associated to 3-simplices of type (k) on a 3D mesh with k = 2 and k = 3. Exercise 3.17. Determine the dimension of the fnite element space on a 2D mesh associated to a family of 2D Crouzeix-Raviart fnite elements introduced in Exercise 3.8. Is this fnite element space equal to the fnite element space associated to a family of 2-simplices of type (1)? Exercise 3.18. Determine the dimension of the fnite element spaces associated to a family of Hermite 2- and 3-simplices of type (3).

COMPUTER EXERCISES Computer Exercise 3.1. By modifying the program of Section 3.1.1, approximate the L2 and H 1 convergence rates of Largange interpolants of type (1) for the functions (a) f (x) = x0.8 , and (b) f (x) = x0.2 . Justify your results in case convergence rates are not optimal. Computer Exercise 3.2. Plot local shape functions of 1-simplex of type (3) on the unit 1-simplex K = [0, 1].

52

Finite Element Methods in Civil and Mechanical Engineering

Computer Exercise 3.3. Plot local shape functions of 2-simplex of type (2) on the unit 2-simplex. Computer Exercise 3.4. Plot the local shape functions of the 2D Raviart-Thomas fnite element on a 2-simplex with the vertices (0, 0), (2, 1), and (−1, 3). ´ e´ lec fnite Computer Exercise 3.5. Plot the local shape functions of the 2D Ned element on a 2-simplex with the vertices (0, 0), (2, 1), and (−1, 3). Computer Exercise 3.6. Let Ω be the unit square with the vertices v1 = (0, 0), v2 = (1, 0), v3 = (1, 1), and v4 = (0, 1), and let Th = {K1 , K2 } be a mesh of Ω, where K1 and K2 are triangles with the vertices {v1 , v2 , v3 } and {v1 , v3 , v4 }, respectively. Plot the vector feld pp(x, y) = (x2 + 2xy, x + y2 ) on Ω. Also separately plot the fnite element interpolants Ih1 pp, Ih2 pp, IhRT pp, and IhN pp. Which of these interpolants are discontinuous? Which one coincides with pp? Computer Exercise 3.7. Determine the dimension of the fnite element space on a 1D mesh associated to 1-simplices of type (2). Plot the corresponding locally supported global shape functions on a simple mesh consisting of 3 elements. Computer Exercise 3.8. Consider a uniform mesh of Ω = (0, 4) with vertices {0, 1, 2, 3, 4}. Plot global shape functions associated to the vertex 2 induced by a family of: (a) 1-simplices of type (3), and (b) Hermite 1-simplices of type (3). Compare these global shape functions. Computer Exercise 3.9. Consider a sequence of successively refned meshes for the unit square of Computer Exercise 3.6 and plot the errors kvv − Ih1 vvk2 , kvv − Ih2 vvk2 , � b kvv − IhRT vvk2 , and kvv − IhN vvk2 versus h, where vv(x, y) = (x2 + y2 ) 2 , sin(xy) , with b = 2.5. Determine the convergence rate for each case. Computer Exercise 3.10. Repeat Computer Exercise 3.9 with b = 1.5, and b = 0.5. Do the convergence rates coincide with those of Computer Exercise 3.9? Why? Computer Exercise 3.11. Let Ω be the unit cube with the vertices v1 = (0, 0, 0), v2 = (1, 0, 0), v3 = (1, 1, 0), v4 = (0, 1, 0), v5 = (0, 0, 1), v6 = (1, 0, 1), v7 = (1, 1, 1), and v8 = (0, 1, 1). Plot the errors kvv − Ih1 vvk2 , kvv − Ih2 vvk2 , kvv − IhRT vvk2 , and kvv − IhN vvk2 � b versus h, where vv(x, y, z) = (x2 + y2 + z2 ) 2 , sin(xyz), cos(xyz) , with b = 2.5. Determine the convergence rate for each case. Computer Exercise 3.12. Repeat Computer Exercise 3.11 with b = 1.5, and b = 0.5. Do the convergence rates coincide with those of Computer Exercise 3.11?

COMMENTS AND REFERENCES Instructions for installing FEniCS together with several tuturial and solved examples are avalible in the FEniCS Project webpage https://fenicsproject.org/, see

Finite Element Interpolation

53

also Appendix A. The FEniCS Tutorial [9] is a good guide for beginners. A more complete tutorial is the FEniCS Book [10], though some parts of this tutorial are outdated. A good introduction to Python programming for scientifc computations is given by [8], also see Appendix B. The defnition of a fnite element given in Section 3.2 was frst introduced in the classic book of Ciarlet [5]. In this book, a complete study of several types of fnite elements that are commonly used in engineering and science is avaiable. It also discusses curved-face elements and additional numerical errors that are induced by approximating curved boundaries. Another good reference for common fnite elements including the Raviart-Thomas and the N´ed´elec elements is [6]. Proofs of the error estimates (3.16) and (3.17) are also available in [6, Chapter 1].

Conforming Finite Element 4 Methods for PDEs In this chapter, we will study fnite element approximations of solutions of timeindependent and time-dependent classes of second-order partial differential equations (PDEs) that arise in many engineering and scientifc applications. After defning the general form of these PDEs, we introduce the notion of weak forms for boundary value problems associated to these PDEs. We then employ fnite element spaces introduced in the previous chapter to obtain fnite element methods for approximating solutions of weak forms. By considering several examples, we also discuss implementation of PDEs subject to different boundary conditions in FEniCS. Finally, a brief introduction to mixed fnite element methods and their implementation in FEniCS is provided. Approximation methods of this chapter are called conforming fnite element methods in the sense that fnite element spaces containing approximate solutions are subsets of the Sobolev spaces that contain weak solutions of PDEs.

4.1

SECOND-ORDER ELLIPTIC PDES

Let Ω be an open domain in Rn , n = 1, 2, 3, and let ∂ Ω be its boundary. We consider the PDE −div(∇u · D) + bb·∇u + cu = f in Ω, (4.1) where u : Ω → R is the unknown function, D(x) is an n × n matrix at any point x = (x1 , . . . , xn ) of Ω with the components di j (x), and bb(x) = (b1 (x), . . . , bn (x)) is a vector feld. The components di j and bi together with c and f are assumed to be known real-valued functions on Ω. The PDE (4.1) can be written more explicitly as n

−

∑

i, j=1

n

∂x j (di j ∂xi u) + ∑ bi ∂xi u + cu = f in Ω, i=1

and therefore, it is a second-order PDE, that is, the highest order derivative in (4.1) is 2. A PDE of the form (4.1) is called elliptic if for all row vectors w = (w1 , . . . , wn ) ∈ Rn , there exists β > 0 such that w · D(x) · wT ≥ β kwk2 , at any x ∈ Ω, (4.2) √ where wT is the transpose of w and kwk = w · wT is the standard norm of w in Rn . Equivalently, (4.2) can be stated as n

∑

i, j=1

n

di j (x) wi w j ≥ β ∑ wi2 , at any x ∈ Ω. i=1

55

56

Finite Element Methods in Civil and Mechanical Engineering

Elliptic PDEs of the form (4.1) arise in many applications. A few examples include: − For n = 1, b = 0, c = 0, and D(x) = E(x)A(x), one obtains the equilibrium equation ∂x (EA∂x u) + f = 0, of a 1D linearly elastic bar in terms of the displacement u(x), where A(x) and E(x) are respectively the cross section area and Young’s modulus of the bar. − With b = 0, c = 0, and D being the identity matrix I, one recovers Poisson’s equation −Δu = f . This equation has several applications such as torsion of prismatic rods, defections of elastic membranes, incompressible potential fows, and electrostatics. − For D = kI, and c = 0, with k > 0, one obtains the heat transfer equation. − For D = kI, with k > 0, (4.1) reduces to the advection-diffusion equation. The PDE (4.1) has infnitely many solutions in general. To solve this PDE for a specifc solution, one also needs to specify suitable boundary conditions to obtain a well-defned boundary value problem. We will discuss some possible choices of boundary conditions later in this chapter.

4.2

WEAK FORMULATIONS OF ELLIPTIC PDES

For obtaining a solution of the PDE (4.1), we should supplement it with a suitable boundary condition as well. Solving the resulting boundary value problem means fnding a function u : Ω → R that satisfes the PDE (4.1) in Ω and also satisfes the boundary condition on the boundary ∂ Ω. Thus, the solution process involves three basic aspects: Finding a unique function (existence and uniqueness) that is twice differentiable (regularity) and satisfes the boundary value problem. It was observed that the simultaneous achievement of the above three aspects is hard. A common approach is to frst study the existence and the uniqueness of solutions and then their regularity. To this end, alternative formulations of boundary value problems called weak formulations are employed. Usually, it is theoretically and computationally more convenient to study solutions of weak formulations, which are also called weak solutions. If a weak solution is suffciently differentiable (regularity), then it is also a solution of the associated boundary value problem. Therefore, weak solutions are good candidates for solutions of boundary value problems. It should be mentioned that in some applications such as continuum mechanics, weak formulations rather than PDEs are considered to represent the real physics of systems. In the remainder of this section, we write the weak formulation of some boundary value problems associated to (4.1). We do not discuss regularity of weak solutions in this book. Under some mild conditions, one can show that weak solutions of elliptic PDEs are also suffciently differentiable. Studying the regularity of weak solutions of more general PDEs is usually a challenging task.

57

Conforming Finite Element Methods for PDEs

4.2.1

DIRICHLET BOUNDARY CONDITION

The simplest boundary condition for (4.1) is the homogeneous Dirichlet boundary condition u = 0, on ∂ Ω. (4.3) Let v be an arbitrary function that vanishes on the boundary ∂ Ω. By multiplying (4.1) by v and integrating over Ω, one obtains −

Z

v div(∇u · D) +

Ω

Z Ω

(bb·∇u)v +

Z

Z

cuv = Ω

f v.

(4.4)

Ω

Green’s formula (2.4) implies that −

Z

v div(∇u · D) =

Z

Ω

∇u · D · (∇v)T −

Ω

Z

(v∇u · D) · n ,

(4.5)

∂Ω

where (∇v)T is the transpose of the row vector field ∇v and n is the outward unit normal vector field at ∂ Ω. Using the above relation and v|∂ Ω = 0, we can write (4.4) as Z Z Z Z ∇u · D · (∇v)T + (bb·∇u)v + cuv = f v. (4.6) Ω

Ω

Ω

Ω

If u is a solution of (4.1) subject to (4.3), then it satisfies (4.6). Both sides of (4.6) are well-defined for u, v ∈ H01 (Ω), see Example 2.11. Therefore, we can define the following weak formulation of (4.1) subject to the Dirichlet boundary condition (4.3): Find u ∈ H01 (Ω) such that Z

B(u, v) =

f v, for all v ∈ H01 (Ω),

(4.7)

Ω

where the bilinear form B : H01 (Ω) × H01 (Ω) → R is given by B(u, v) =

Z

∇u · D · (∇v)T + (bb·∇u)v + cuv

Ω

=

Z n n Ω

∑

i, j=1

o n di j ∂xi u ∂x j v + v ∑ bi ∂xi u + cuv .

(4.8)

i=1

Notice that (4.1) involves second-order derivatives of u, however, the weak formulation (4.7) only involves first-order derivatives of u and v. A solution of the weak formulation is called a weak solution of the PDE (4.1) subject to (4.3). To show that a weak solution is also a solution of the PDE (also called a strong solution), one should verify that the weak solution is twice differentiable (regularity). Hence, a strong solution is also a weak solution but the converse may not hold, in general. The arbitrary function v in the weak formulation (4.7) is called the test function. The space that contains weak solutions is called the solution space or the trial space. The space that contains test functions is called the test space. The solution and the test spaces of (4.7) are the same and are equal to H01 (Ω).

58

Finite Element Methods in Civil and Mechanical Engineering

In (4.7), the homogeneous Dirichlet boundary condition (4.3) is directly imposed in the solution space H01 (Ω). This boundary condition is then called an essential boundary condition for (4.7). 6 0, then the One may wonder why in (4.7) we assumed that v|∂ Ω = 0. If v|∂ Ω = R term ∂ Ω (v∇u ·D)· n on the right side of (4.5) that contains frst-order derivatives of u should be added to the bilinear form B(u, v) as well. One can show that if u ∈ H 1 (Ω), then u is integrable on the boundary. However, the frst-order derivatives of u are not integrable on the boundary, in general. Consequently, the bilinear form B(u, v) is not well-defned as a function H 1 (Ω) × H 1 (Ω) → R, if v|∂ Ω 6= 0; See also the discussion of Section 2.8 regarding Green’s formulas in terms of Sobolev spaces. A non-homogeneous Dirichlet boundary condition for (4.1) reads u = g, on ∂ Ω,

(4.9)

where g : Ω → R is a known function. Notice that we assumed g is defned on the whole domain and not only on its boundary. If u is a solution of (4.1) subject to (4.9), then w = u − g vanishes at the boundary and is a solution of the problem: Find w ∈ H01 (Ω) such that Z

f v − B(g, v), for all v ∈ H01 (Ω),

B(w, v) =

(4.10)

Ω

where the bilinear form B is defned in (4.8). Thus, the solution of problems subject to non-homogeneous Dirichlet boundary conditions can be reduced to the solution of problems subject to the homogeneous Dirichlet boundary condition. Alternatively, one may consider the following weak formulation for non-homogeneous Dirichlet boundary conditions: Find u ∈ H 1 (Ω) satisfying u|∂ Ω = g such that Z

B(u, v) =

f v, for all v ∈ H01 (Ω).

(4.11)

Ω

The trial and the test spaces of the second weak formulation are not the same. This weak formulation is more suitable in case the function g is only defned at the boundary. 4.2.2

NEUMANN BOUNDARY CONDITION

A Neumann boundary condition for (4.1) is given by (∇u · D) · n = g, on ∂ Ω,

(4.12)

where g : ∂ Ω → R is a known function on the boundary. Notice that if D is the identity matrix I, then (4.12) simply states that the normal derivative ∂n u of u at ∂ Ω is g. By following the approach used for deriving (4.7), one obtains the following weak formulation for the PDE (4.1) subject to (4.12): Find u ∈ H 1 (Ω) such that Z

B(u, v) =

Z

fv+ Ω

∂Ω

gv, for all v ∈ H 1 (Ω),

(4.13)

59

Conforming Finite Element Methods for PDEs

where the bilinear form B(u, v) is the same as (4.8). The solution and the test space of this weak formulation are H 1 (Ω). The derivation of (4.13) is left as an exercise. Notice that the Neumann boundary condition (4.12) is implicitly imposed by the weak formulation and not directly by the solution space. This boundary condition is then called a natural boundary condition for (4.13). It is also possible to have a combination of Dirichlet and Neumann boundary conditions as follows: Suppose that ∂ Ω consists of the two parts ΓD and ΓN and consider the mixed Dirichlet-Neumann boundary condition u = gˆ, on ΓD , (4.14) (∇u · D) · n = g, on ΓN . Let HD1 (Ω) = {v ∈ H 1 (Ω) : v = 0 on ΓD }. Then, a weak formulation for the PDE (4.1) subject to (4.14) can be stated as: Find u ∈ H 1 (Ω) such that u|ΓD = gˆ and Z

B(u, v) =

Z

fv+ Ω

gv, for all v ∈ HD1 (Ω),

(4.15)

ΓN

where the bilinear form B(u, v) is given by (4.8). The derivation of (4.15) is left as an exercise. 4.2.3

ROBIN BOUNDARY CONDITION

Given two functions r : ∂ Ω → R and g : ∂ Ω → R, a Robin boundary condition for the PDE (4.1) is given by ru + (∇u · D) · n = g, on ∂ Ω.

(4.16)

A weak formulation for (4.1) subject to the Robin boundary condition (4.16) can be stated as: Find u ∈ H 1 (Ω) such that Z

B(u, v) +

Z

ruv = ∂Ω

Z

fv+ Ω

gv, for all v ∈ H 1 (Ω),

(4.17)

∂Ω

where the bilinear form B(u, v) is defned in (4.8). The derivation of (4.17) is left as an exercise. Notice that similar to Neumann boundary conditions, Robin boundary conditions are imposed as a natural boundary condition in (4.17). Let us emphasize that the selection of a specifc boundary condition depends on the physical problem that is modeled by (4.1). For example, to model a membrane fxed at its boundary similar to the one discussed in Chapter 1, one can employ the homogeneous Dirichlet boundary condition. It is also possible to have more than one possible choices of boundary conditions for some physical problems.

60

Finite Element Methods in Civil and Mechanical Engineering

4.3 WELL-POSEDNESS OF WEAK FORMULATIONS The weak formulations introduced in Section 4.2 can be considered as special cases of the following abstract problem: Let (X, k·k) be a normed linear space and suppose A : X × X → R and F : X → R are respectively a bilinear form and a linear functional. We consider the abstract problem: Find u ∈ X such that A(u, v) = F(v), for all v ∈ X.

(4.18)

The data X, A, and F for the weak formulations of Section 4.2 are the followings: − The formulation (4.7) for the homogeneous Dirichlet boundary condition: X = H01 (Ω),

Z

A(u, v) = B(u, v),

F(v) =

f v; Ω

− The formulation (4.13) for the Neumann boundary condition: X = H 1 (Ω),

Z

A(u, v) = B(u, v),

Z

F(v) =

fv+ Ω

gv; ∂Ω

− The formulation (4.15) for the mixed Dirichlet-Neumann boundary condition with gˆ = 0: X = HD1 (Ω),

Z

A(u, v) = B(u, v),

F(v) =

Z

fv+ Ω

gv; ΓN

− The formulation (4.17) for the Robin boundary condition: X = H 1 (Ω),

Z

A(u, v) = B(u, v) +

Z

ruv, ∂Ω

F(v) =

Z

fv+ Ω

gv. ∂Ω

By inspecting the derivations of the above weak formulations, we observe that linearity of A(u, v) with respect to u and v stems from different origins: The linearity with respect to u follows form the linearity of the left side of the PDE (4.1) with respect to u while the linearity with respect to v simply follows from the derivation of these weak formulations. This means that if we replace (4.1) with a nonlinear PDE and follow the derivation of Section 4.2, we will obtain a weak formulation of the form A(u, v) = F(v), where A(u, v) is only linear with respect to its second argument v. Our main goal in this section is to study solvability of the abstract problem (4.18). For many problems in engineering and science, it is reasonable to expect that (4.18) has only one solution u for each data F. It is also reasonable to assume the solution u continuously depends on the data F in the sense that “small” changes in F will result in “small” changes in the corresponding solution u. This continuity condition can be stated as: There is α > 0, such that for any F we have kuk ≤ αkFk,

(4.19)

where u is the solution associated to F and k · k denotes appropriate norms. We say the abstract problem (4.18) is well-posed (in the Hadamard sense) if it has a unique

Conforming Finite Element Methods for PDEs

61

solution and if the continuity condition (4.19) holds. Notice that the above defnition of well-posedness is not appropriate for some applications including those with nonunique solutions such as buckling of columns. The Lax-Milgram lemma provides a suffcient condition for the well-posedness of the problem (4.18): Suppose A is X-elliptic, also called coercive, in the sense that there exists γ > 0 such that A(v, v) ≥ γkvk2 , for all v ∈ X.

(4.20)

Then, the Lax-Milgram lemma states that the problem (4.18) is well-posed with α = 1/γ in (4.19). Notice that coercivity is only a suffcient condition for well-posedness and it is possible to have well-posed non-coercive problems as well. One can show that, under some suitable conditions, all the weak formulations of (4.1) with the boundary conditions discussed in Section 4.2 are coercive and consequently, are well-posed. Other important examples of coercive problems that do not belong to the class of elliptic PDEs of the form (4.1) include linearized elasticity and the biharmonic problem. Later we will see that a very useful property of coercive problems is that well-posedness of these problems is automatically inherited by their discretizations.

4.4 VARIATIONAL STRUCTURE It is possible that the abstract problem (4.18) admits a variational structure in the sense that it corresponds to the minimization of a functional. More specifcally, suppose A(u, v) is symmetric, that is, A(u, v) = A(v, u), and consider the functional J : X → R given by 1 J(v) = A(v, v) − F(v). 2

(4.21)

For any t ∈ R and u, v ∈ X, we can write h i t2 J(u + tv) − J(u) = t A(u, v) − F(v) + A(v, v). 2

(4.22)

Therefore, if u solves the abstract problem (4.18) and if A is positive, that is, A(v, v) ≥ 0, for all v ∈ X, we have J(u + tv) − J(u) ≥ 0. In summary, we conclude that if A is symmetric and positive, then a solution of the problem (4.18) is a minimizer of J. Conversely, it is not hard to show that a minimizer of J is a solution of the problem (4.18). We leave the proof as an exercise. When the abstract problem (4.18) corresponds to a minimizer of the functional (4.21), it is also called a variational formulation. Notice that if the bilinear form A is coercive, then it is positive. Thus, the LaxMilgram lemma implies that if the bilinear form A is symmetric and coercive, then the unique solution of (4.18) is the unique minimizer of (4.21). In practice, the functional (4.21) usually represents a specifc type of energy. For example, it represents the stored energy in linearized elasticity.

62

Finite Element Methods in Civil and Mechanical Engineering

Inspection of the bilinear form B(u, v) defned in (4.8) suggests that B(u, v) is symmetric if the vector feld b of the PDE (4.1) vanishes. In this case, it is straightforward to show that the weak formulations of Section 4.2 are also variational formulations. The verifcation of this result is left as an exercise.

4.5 THE GALERKIN METHOD AND FINITE ELEMENT METHODS The Galerkin method provides a framework for approximating solutions of the abstract problem (4.18). Assume Xh is a fnite-dimensional linear subspace of possibly infnite-dimensional, linear space X. To approximate a solution of (4.18) using the Galerkin method, we discretize the problem (4.18) as: Find uh ∈ Xh such that A(uh , vh ) = F(vh ), for all vh ∈ Xh .

(4.23)

Here, uh is the approximate solution, vh is the test function, and Xh is both the solution (or trial) space and the test space. The approximation scheme (4.23) is called a conforming method as Xh ⊂ X. If Xh 6⊂ X, it is called a non-conforming method. If the fnite-dimensional linear space Xh is a fnite element space and Xh ⊂ X, the problem (4.23) is called a conforming fnite element method for approximating a solution of (4.18). As discussed in Section 4.4, if the bilinear form A(u, v) of the abstract problem (4.18) is symmetric, then a solution of (4.18) is also a minimizer of the functional J(v) in X, where J is defned in (4.21). In this case, a solution of the discrete problem (4.23) is a minimizer of J in the space Xh . Approximating a minimizer of J in X by using a minimizer of J in Xh is called the Ritz method. In the sense discussed above, the Ritz method is equivalent to the Galerkin method. To obtain conforming fnite element methods for elliptic PDEs subject to the boundary conditions of Section 4.2, let Th be a simplicial mesh of the underlying domain Ω and let Vh be the fnite element space associated to a family of simplicial Lagrange element of type (k) on Th . Recall that in Section 3.4.1.1, we showed that Vh ⊂ H 1 (Ω). By using suitable data mentioned in Section 4.3 for different boundary conditions, (4.23) leads to conforming fnite element methods for boundary value problems associated to the PDE (4.1), where Xh is selected as follows: − The homogeneous Dirichlet boundary condition: Xh = {vh ∈ Vh : vh = 0 on ∂ Ω};

(4.24)

− Mixed Dirichlet-Neumann boundary condition: Xh = {vh ∈ Vh : vh = 0 on ΓD };

(4.25)

− Neumann and Robin boundary conditions: Xh = Vh . In the remainder of this section, we study well-posedness and convergence of conforming fnite element methods associated to (4.23).

63

Conforming Finite Element Methods for PDEs

4.5.1

THE STIFFNESS MATRIX

We begin by showing that solving (4.23) is equivalent to solving a system of linear equations. Let {ψ1 , . . . , ψN } be a basis for the fnite-dimensional space Xh with dim Xh = N, where Xh is not necessarily a fnite element space. Since uh ∈ Xh , we can write N

uh =

∑ Ujψ j,

(4.26)

j=1

where U1 , . . . ,UN are unknown constants to be determined. It is straightforward to show that uh solves (4.23) if and only if A(uh , ψi ) = F(ψi ),

i = 1, . . . , N.

By using (4.26), the above equations can be written as the linear system [A]N×N · UN×1 = FN×1 ,

(4.27)

where [A] is called the stiffness matrix with the components [A]i j = A(ψ j , ψi ), and U and F are vectors given by ⎡ ⎤ ⎡ ⎤ F(ψ1 ) U1 ⎢ ⎥ ⎢ ⎥ . U = ⎣ ... ⎦ , F = ⎣ ⎦. .. UN

F(ψN )

Notice that the stiffness matrix [A] is the matrix representation of the bilinear form A : Xh × Xh → R in the basis {ψ1 , . . . , ψN } in the sense discussed in Section 2.5. It is straightforward to verify that if A is a symmetric bilinear form, then [A] will be a symmetric matrix. If Xh is a fnite element space, then {ψ1 , . . . , ψN } can be chosen to be the locally supported global shape functions of Xh . If, as for boundary value problems associated to elliptic PDEs, A(u, v) is obtained by integrating terms that involve the multiplication of u and v or their derivatives, the selection of global shape functions will lead to sparse stiffness matrices, that is, most entries of [A] will be zero. The sparsity of [A] is very useful for effciently computing the solution of the linear system (4.27). Example 4.1. Consider the mesh of a square with 8 elements shown in the middle of Figure 3.18. Suppose we use the fnite element space associated to a family of 1-simplex of type (1) to obtain fnite element methods for an elliptic PDE. Then, the size of the stiffness matrix [A] will depend on the selected boundary condition. In particular, the size of [A] is 1 × 1 for the homogeneous Dirichlet boundary condition, 6 × 6 for a mixed Dirichlet-Neumann boundary condition with ΓD being the bottom edge of the square, and 9 × 9 for Neumann and Robin boundary conditions. Example 4.2. Suppose the mesh of Figure 3.23 is employed to obtain a fnite element method for an elliptic PDE using 2-simplex of type (2). Then, the entry of the stiffness matrix corresponding to the global shape functions shown on the left and the right of Figure 3.23 is zero.

64

Finite Element Methods in Civil and Mechanical Engineering

4.5.2 WELL-POSEDNESS OF COERCIVE DISCRETE PROBLEMS Generally speaking, the well-posedness of the abstract problem (4.18) does not imply the well-posedness of the discrete problem (4.23). Coercive problems are exceptions in this regard. More specifcally, suppose the coercivity condition (4.20) holds and let Xh ⊂ X. Then, we can write A(vh , vh ) ≥ γkvh k2 , for all vh ∈ Xh ,

(4.28)

and therefore, A : Xh × Xh → R is also coercive. The Lax-Milgram lemma then tells us that the discrete problem (4.23) is well-posed and has a unique solution. Thus, the well-posedness of conformal discretizations of coercive problems automatically follows from the well-posedness of these problems. The coercivity of A : Xh × Xh → R implies that the stiffness matrix [A]N×N is positive defnite in the sense that vT · [A] · v ≥ 0, for all column vectors v ∈ RN .

(4.29)

Symmetric, positive defnite matrices have positive eigenvalues. 4.5.3

CONVERGENCE OF FINITE ELEMENT SOLUTIONS

Let {Th } be a family of meshes of a domain Ω. For each mesh Th of this family, one can write a fnite element method in the form (4.23) with the solution uh ∈ Xh . Consequently, we obtain a family of approximate solutions {uh } associated to {Th }. We say this family of fnite element solutions converges in the norm k · k to a solution u ∈ X of the problem (4.18) if lim ku − uh k = 0,

h→0

where ku − uh k is the error in the norm k · k. Notice that since Xh ⊂ X, we have u − uh ∈ X and thus, ku − uh k is well-defned when k · k is the norm of X. As will be discussed later, it is also possible to study convergence in a norm other than the norm of the solution space X as long as the other norm of u − uh is also well-defned. To study the convergence of fnite element solutions, let us assume that the bilinear form A : X × X → R is coercive and continuous. This means there are two positive constants γ and ξ such that for all v, w ∈ X, we have A(v, v) ≥ γkvk2 ,

and

|A(v, w)| ≤ ξ kvk kwk.

The coercivity of A implies that the problem (4.18) and the discrete problem (4.23) are well-posed and have unique solutions u and uh , respectively. Since for any wh ∈ Xh , we have A(u, wh ) = F(wh ) = A(uh , wh ), one concludes that A(u − uh , wh ) = 0, for all wh ∈ Xh .

(4.30)

Conforming Finite Element Methods for PDEs

65

Thus, we can write γku − uh k2 ≤ A(u − uh , u − uh ) = A(u − uh , u − uh ) + A(u − uh , uh − vh ) = A(u − uh , u − vh ) ≤ ξ ku − uh k ku − vh k, where we used coercivity in the frst line, the relation (4.30) with wh = uh − vh in the second line, and the linearity of A in its second argument and the continuity of A in the last line. Thus, we conclude that ku − uh k ≤

ξ ku − vh k, for all vh ∈ Xh , γ

which implies that ku − uh k ≤

ξ inf ku − vh k. γ vh ∈Xh

(4.31)

The term inf ku − vh k measures the “ability” of the fnite element space Xh to apvh ∈Xh

proximate u. In particular, we say that the above family of fnite element methods based on (4.23) has the approximability property if for any v ∈ X, we have lim inf ku − vh k = 0.

h→0 vh ∈Xh

The inequality (4.31) implies that a family of discrete solutions associated to a family of conforming fnite element methods for a coercive problem is convergent if that family has the approximability property. This result is called C´ea’s lemma. As mentioned earlier, the weak formulations of elliptic PDEs subject to the aforementioned boundary conditions are coercive problems with X being H 1 (Ω) or one of its linear subspaces. Conforming fnite element methods for these problem can then be obtained by using H 1 -conformal fnite element spaces Xh such as those associated to the n-simplex of type (k). Let Ihk u be an interpolation of the solution u, where the interpolation operator Ihk is associated to the n-simplex of type (k). Using the second inequality of (3.16), we conclude that if {Th } is a shape-regular family of meshes and if u is s + 1 times differentialble with 1 ≤ s ≤ k, we have ku − uh k1,2 ≤

ξ ξ inf ku − vh k1,2 ≤ ku − Ihk uk1,2 ≤ C˜ (u)hs , γ vh ∈Xh γ

(4.32)

where C˜ (u) does not depend on h. Consequently, we have ku − uh k1,2 → 0 as h → 0. Notice that the convergence of fnite element methods follows from the convergence of fnite element interpolations. This is the reason for studying fnite element interpolations before fnite element methods for PDEs. Similar to fnite element interpolations, the power of h in (4.32) is called the convergence rate. The maximum convergence rate of the n-simplex of type (k) in

66

Finite Element Methods in Civil and Mechanical Engineering

the H 1 -norm is k, which occurs if u is k + 1 times differentiable. It is also possible to show that ku − uh k2 ≤ C(u)hs+1 ,

(4.33)

where C(u) does not depend on h. Thus, the above fnite element methods also converge in the L2 -norm. If u is k + 1 times differentiable, one obtains the maximum convergence rate k + 1 in the L2 -norm, which is called the optimal convergence rate.

4.6

IMPLEMENTATION: THE POISSON EQUATION

Now we discuss how to use FEniCS to solve an elliptic PDE subject to different boundary conditions. Consider the 2D Poisson’s equation −Δu = f in Ω,

(4.34)

where Ω is the unit square (0, 1) × (0, 1) and f (x, y) = −4b2 (x2 + y2 )b−1 , with b being a constant. A solution of this PDE is given by ue (x, y) = (x2 + y2 )b .

(4.35)

Throughout this section, we solve (4.34) subject to different boundary conditions such that (4.35) is the solution of the associated boundary value problems. To obtain a known exact solution such as (4.35) for the PDE (4.34), we begin by assuming any arbitrary form for the solution and then obtain the associated function f on the right side of (4.34) by simply plugging the assumed exact solution in the left side of (4.34). This way, one can always construct a test problem with a known solution for a given PDE. Test problems are very useful for debugging codes and for studying the performance of numerical methods since one can calculate the error of approximate solutions. 4.6.1

DIRICHLET BOUNDARY CONDITION

We begin by considering the 2D Poisson equation (4.34) subject to a Dirichlet boundary condition: Find u such that −∂xx u − ∂yy u = f (x, y), in Ω, on ∂ Ω. u = ue , By using (4.11), the weak formulation of the above boundary value problem can be stated as: Find u ∈ H 1 (Ω) satisfying u = ue , on ∂ Ω, such that Z Ω

∇u · ∇v =

Z

f v, for all v ∈ H01 (Ω),

(4.36)

Ω

where “·” on the left side denotes the inner product of vectors. Let Vh be the fnite element space induced by a family of 2-simplices of type (k) on a mesh Th of Ω. By

67

Conforming Finite Element Methods for PDEs

using Vh and its subspace Xh defned in (4.24), we can defne the following conformal fnite element method based on (4.36): Find uh ∈ Vh with uh = Ihk ue , on ∂ Ω, such that Z Ω

∇uh · ∇vh =

Z

f vh , for all vh ∈ Xh .

(4.37)

Ω

Consider a structured simplicial mesh for the unit square with N denoting the number of uniform divisions of each edge. Given N, the degree of fnite element k, and the parameter b of the function f , the function Compute_Dirichlet computes the fnite element solution of (4.37) and its L2 - and H 1 -errors. def Compute_Dirichlet(N, k, b): """ Given mesh division size N, the degree k of Lagrange elements, and the parameter b of the exact solution, this function returns the mesh size h and errors of the finite element solution associated to a Dirichlet boundary condition """ # create mesh and define function space mesh = UnitSquareMesh(N, N) V_h = FunctionSpace(mesh, "Lagrange", k) # the exact solution u_e = Expression(’pow(x[0]*x[0] + x[1]*x[1], b)’, b = b, degree = 3 + k) # define Dirichlet boundary condition def ue_boundary(x, on_boundary): return on_boundary BC = DirichletBC(V_h, u_e, ue_boundary) # define the weak formulation u_h = TrialFunction(V_h) v_h = TestFunction(V_h) f = Expression(’-4*b*b*pow(x[0]*x[0] + x[1]*x[1], b-1)’, b = b, degree = 3 + k) A = inner(grad(u_h), grad(v_h))*dx F = f*v_h*dx # compute finite element solution u_h = Function(V_h) # u_h is redefined solve(A == F, u_h, BC) L2_Error = errornorm(u_e, u_h, norm_type="L2") H1_Error = errornorm(u_e, u_h, norm_type="H1")

68

Finite Element Methods in Civil and Mechanical Engineering

return L2_Error, H1_Error, mesh.hmax() This function works as follows: After defning the mesh mesh, the fnite element space V_h, and the exact solution u_e, the Dirichlet boundary condition is defned by the following segment: def ue_boundary(x, on_boundary): return on_boundary BC = DirichletBC(V_h, u_e, ue_boundary) The function ue_boundary is used to specify the region that the Dirichlet boundary condition will be applied. It returns the boolean value True if x is on the boundary and False otherwise. The argument on_boundary is provided by FEniCS and is True if x is on the boundary of a mesh. Since the Dirichlet boundary condition is applied on the whole boundary in the present problem, ue_boundary can simply return the value of on_boundary. Later, we will see less trivial examples of functions that only specify a specifc portion of the boundary. The function DirichletBC defnes the Dirichlet boundary condition using the value u_e and the function ue_boundary that specifes where the boundary condition should be applied. From the mathematical point of view, the trial and the test spaces of the problem (4.37) are different and are obtained by subjecting Vh to 2 different boundary conditions. In FEniCS, boundary conditions are specifed in the fnal stage of solving discrete equations and not in the defnition of function spaces. Consequently, the defnitions of the trial and the test functions are similar for the problem (4.37): u_h = TrialFunction(V_h) v_h = TestFunction(V_h) After defning the function f , we defne the left and the right sides of the problem (4.37): A = inner(grad(u_h), grad(v_h))*dx F = f*v_h*dx The language employed in FEniCS to express weak forms is called UFL (Unifed Form Language). UFL provides a very interesting feature of FEniCS: The above segment is very similar to the mathematical formula of the problem (4.37). The integration over domains is denoted by dx in UFL. As will be seen later, ds is used for the integration over domains’ boundaries. Having defned the problem (4.37), its solution is computed as: u_h = Function(V_h) # u_h is redefined solve(A == F, u_h, BC)

Conforming Finite Element Methods for PDEs

69

In FEniCS, unknowns of weak formulations are TrialFunction objects and solutions are Function objects. In the mathematical statement (4.37), these concepts are the same and are denoted by uh . After defning the bilinear form A, we do not need u_h as a TrialFunction object anymore. To be consistent with (4.37), we reused the variable name u_h and redefne it as a Function object to store the solution. Given the weak formulation A == F, the function object u_h, and the boundary condition BC, the function solve computes the fnite element solution and stores it in u_h. Finally, the L2 - and H 1 -errors are calculated by using errornorm and the results together with the mesh size mesh.hmax() are returned. By default, solve uses sparse LU decomposition (Gaussian elimination) to solve linear systems. It is possible to change this option and use iterative methods such as preconditioned Krylov solvers, which are more appropriate for large problems. One approach for debugging and checking fnite element programs is to calculate errors and convergence rates associated to a family of meshes. For example, the following program computes the errors and convergence rates by using triangle of type (2), the parameter b = 1.25, and the meshes N = 3, 6, 9, 12. from dolfin import * # number of divisions of meshes Divisions = [3, 6, 9, 12] # computing errors b = 1.25 k = 2 h = [] error_L2, error_H1 = [], []

# # # #

parameter of the exact solution degree of Lagrange element mesh sizes initializing errors

for N in Divisions: Er_L2, Er_H1, hmesh = Compute_Dirichlet(N, k, b) error_L2.append(Er_L2); error_H1.append(Er_H1) h.append(hmesh) # convergence rates from math import log as ln # log is a dolfin name too rate_L2 = ln(error_L2[-1]/error_L2[-2])/ln(h[-1]/h[-2]) rate_H1 = ln(error_H1[-1]/error_H1[-2])/ln(h[-1]/h[-2]) # printing results for i in range(len(Divisions)): print(’N = %2.f, L2_error = %5.2E, H1_error = %5.2E’ % (Divisions[i], error_L2[i], error_H1[i])) print(’L2-convergence rate = %.2f, H1-convergence rate = %.2f’ % (rate_L2,rate_H1))

70

Finite Element Methods in Civil and Mechanical Engineering

This program generates the following output: N = 3, L2_error = 9.63E-04, H1_error = 1.94E-02 N = 6, L2_error = 1.21E-04, H1_error = 4.92E-03 N = 9, L2_error = 3.60E-05, H1_error = 2.20E-03 N = 12, L2_error = 1.52E-05, H1_error = 1.24E-03 L2-convergence rate = 3.00, H1-convergence rate = 1.99 The errors are decreasing and the computed convergence rates are consistent with the theoretical values 3 and 2 suggested by the estimates (4.33) and (4.32), respectively. If the solution u belongs to the polynomial space of the underlying fnite element, then the fnite element solution uh will be equal to u regardless of the refnement level of the underlying mesh. In this case, computed errors should be very small for any mesh, in the order of machine error. This is another approach for designing a test problem to verify fnite element programs. For example, if we assume b = 1 in (4.35), then the exact solution will be a polynomial of degree 2. Therefore, the errors of fnite element solutions induced by triangle of type (k) with k ≥ 2 are expected to be very small. By setting b=1 in the above program, we get the following output, which is consistent with our expectation. N N N N

= 3, = 6, = 9, = 12,

L2_error L2_error L2_error L2_error

= = = =

1.38E-15, 4.54E-15, 9.11E-15, 1.82E-14,

H1_error H1_error H1_error H1_error

= = = =

1.56E-14 3.59E-14 6.16E-14 1.05E-13

4.6.2 MIXED DIRICHLET-NEUMANN BOUNDARY CONDITION Suppose ΓD is the bottom edge of the unit square Ω and let ΓN = Γ1N ∪ Γ2N ∪ Γ3N , where Γ1N , Γ2N , and Γ3N are respectively the right, the top, and the left edges of Ω. We solve (4.34) with the following mixed Dirichlet-Neumann boundary condition: Find u such that ⎧ ⎪ ⎪ −∂xx u − ∂yy u = f (x, y), in Ω, ⎨ on ΓD , u = ue , (4.38) ∂ u = ∂ u , on Γ1N and Γ3N , ⎪ x x e ⎪ ⎩ ∂y u = ∂y ue , on Γ2N . By using (4.15), a weak formulation of the above problem can be stated as: Find u ∈ H 1 (Ω) satisfying u = ue , on ΓD , such that Z

∇u · ∇v =

Z

Ω

Z

fv+ Ω

Z

1 ΓN

v∂x ue +

2 ΓN

v∂y ue −

Z 3 ΓN

v∂x ue , for all v ∈ HD1 (Ω).

Similar to the previous section, let Vh be a fnite element space induced by a family of 2-simplices of type (k) and consider the subspace Xh of Vh defned in (4.25). Then, we can write the following fnite element method for the boundary value problem (4.38): Find uh ∈ Vh with uh = Ihk ue , on ΓD , such that Z Ω

∇uh · ∇vh =

Z

Z

f vh + Ω

1 ΓN

Z

v h ∂x u e +

2 ΓN

vh ∂y ue −

Z Γ3N

vh ∂x ue , for all vh ∈ Xh .

Conforming Finite Element Methods for PDEs

71

The following Python function for solving this fnite element method is the analog of Compute_Dirichlet defned in Section 4.6.1 for Dirichlet boundary conditions. def Compute_DN(N, k, b): """ Given mesh division size N, the degree k of Lagrange elements, and the parameter b of the exact solution, this function returns the mesh size h and errors of the finite element solution associated to a mixed Dirichlet-Nuemann boundary condtion """ # create mesh and define function space mesh = UnitSquareMesh(N, N) V_h = FunctionSpace(mesh, "Lagrange", k) # the exact solution and its first derivatives u_e = Expression(’pow(x[0]*x[0] + x[1]*x[1], b)’, b = b, degree = 3 + k) Dx_u_e = Expression("""2*b*x[0]*pow(x[0]*x[0] + x[1]*x[1], b-1)""", b = b, degree = 3 + k) Dy_u_e = Expression("""2*b*x[1]*pow(x[0]*x[0] + x[1]*x[1], b-1)""", b = b, degree = 3 + k) # marking the boundary using a mesh function boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1) # mark right boundary edges as subdomain 1 class RightEdges(SubDomain): def inside(self, x, on_boundary): tol = 1E-12 # tolerance for comparisons return on_boundary and abs(x[0] - 1) < tol Gamma_R = RightEdges() Gamma_R.mark(boundary_parts, 1) # mark top boundary edges as subdomain 2 class TopEdges(SubDomain): def inside(self, x, on_boundary): tol = 1E-12 # tolerance for comparisons return on_boundary and abs(x[1] - 1) < tol Gamma_T = TopEdges() Gamma_T.mark(boundary_parts, 2) # mark left boundary edges as subdomain 3 class LeftEdges(SubDomain): def inside(self, x, on_boundary):

72

Finite Element Methods in Civil and Mechanical Engineering

tol = 1E-12 # tolerance for comparisons return on_boundary and abs(x[0]) < tol Gamma_L = LeftEdges() Gamma_L.mark(boundary_parts, 3) # define Dirichlet boundary condition def bottom_boundary(x, on_boundary): tol = 1E-12 # tolerance for comparisons return on_boundary and abs(x[1]) < tol BC = DirichletBC(V_h, u_e, bottom_boundary) # define the weak formulation u_h = TrialFunction(V_h) v_h = TestFunction(V_h) f = Expression(’-4*b*b*pow(x[0]*x[0] + x[1]*x[1], b-1)’, b = b, degree = 3 + k) ds = Measure("ds", domain=mesh, subdomain_data=boundary_parts) # boundary integral A = inner(grad(u_h), grad(v_h))*dx F = f*v_h*dx + Dx_u_e*v_h*ds(1) + Dy_u_e*v_h*ds(2) \ - Dx_u_e*v_h*ds(3) # compute finite element solution u_h = Function(V_h) # u_h is redefined solve(A == F, u_h, BC) L2_Error = errornorm(u_e, u_h, norm_type="L2") H1_Error = errornorm(u_e, u_h, norm_type="H1") return L2_Error, H1_Error, mesh.hmax() This function works as follows: After defning the mesh and the function spaces, the exact solution and its frst derivatives are defned. The next step is to mark different parts of boundary for imposing the Nuemann boundary condition. This can be done by defning a function on boundary edges that takes the value i on ΓiN , i = 1, 2, 3. This function is defned by using the MeshFunction object: boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1) The argument mesh.topology().dim()-1 specifes that boundary_parts is defned on parts of the mesh of dimension 1 lower than the dimension of mesh. To assign values to boundary_parts, we use subclasses of the object SubDomain. For example, the following code assigns the value 1 to the right boundary Γ1N , where x = 1:

Conforming Finite Element Methods for PDEs

73

class RightEdges(SubDomain): def inside(self, x, on_boundary): tol = 1E-12 # tolerance for comparisons return on_boundary and abs(x[0] - 1) < tol Gamma_R = RightEdges() Gamma_R.mark(boundary_parts, 1) The condition on_boundary and abs(x[0] - 1) < tol is True if x is on the boundary, that is on_boundary is True, and if x is “very close” to 1. Notice that using == to check the value of foating-point variables is not a good programming practice due to small round-off errors. Instead, we employed abs(x[0] - 1) < tol to check the equality x = 1 with a tolerance tol. After marking Γ1N , Γ2N , and Γ3N , the Dirichlet boundary condition is defned on the bottom edge of the square where y = 0: def bottom_boundary(x, on_boundary): tol = 1E-12 # tolerance for comparisons return on_boundary and abs(x[1]) < tol BC = DirichletBC(V_h, u_e, bottom_boundary) The next step is to defne the weak formulation which involves integrals over the mesh as well as integrals over different parts of the mesh boundary. First, we use boundary_parts to defne the symbol ds for the integration over the boundary: ds = Measure("ds", domain=mesh, subdomain_data=boundary_parts) # boundary integral Then, the symbols ds(1), ds(2), and ds(3) simply denote the integration over Γ1N , Γ2N , and Γ3N , respectively. These symbols are employed similar to dx to defne the weak formulation: A = inner(grad(u_h), grad(v_h))*dx F = f*v_h*dx + Dx_u_e*v_h*ds(1) + Dy_u_e*v_h*ds(2) \ - Dx_u_e*v_h*ds(3) The rest of Compute_DN for solving the fnite element method is similar to the function Compute_Dirichlet. By replacing Compute_Dirichlet with Compute_DN in the program of Section 4.6.1, we can compute errors and convergence rates for the mixed Dirichlet-Neumann boundary condition. For example, the output for k = 2, b = 1.25, and the meshes N = 3, 6, 9, 12, is: N = 3, L2_error = 9.69E-04, H1_error = 1.84E-02 N = 6, L2_error = 1.21E-04, H1_error = 4.80E-03 N = 9, L2_error = 3.58E-05, H1_error = 2.16E-03 N = 12, L2_error = 1.51E-05, H1_error = 1.22E-03 L2-convergence rate = 3.00, H1-convergence rate = 1.98

74

Finite Element Methods in Civil and Mechanical Engineering

We encourage the reader to compare this output with the corresponding output in Section 4.6.1. 4.6.3

ROBIN BOUNDARY CONDITION

Suppose ∂ Ω = Γ0 ∪ Γ1 ∪ Γ2 ∪ Γ3 , where Γ0 , Γ1 , Γ2 , and Γ3 are respectively the bottom, the right, the top, and the left edges of the unit square Ω. Consider the following boundary value problem for (4.34) with a Robin boundary condition: Find u such that ⎧ −∂xx u − ∂yy u = f (x, y), in Ω, ⎪ ⎪ ⎪ ⎪ u on Γ0 , ⎨ − ∂y u = ue − ∂y ue , u + ∂x u = ue + ∂x ue , on Γ1 , (4.39) ⎪ 2, ⎪ u + ∂ u = u + ∂ u , on Γ ⎪ y e y e ⎪ ⎩ u − ∂x u = ue − ∂x ue , on Γ3 . By using (4.17), we can write the following weak formulation for this boundary value problem: Find u ∈ H 1 (Ω) such that Z

∇u · ∇v +

Z

Ω

Z

Z

uv = ∂Ω

fv+ ZΩ

+

Γ2

Γ0

(ue − ∂y ue )v + Z

(ue + ∂y ue )v +

Γ3

Z Γ1

(ue + ∂x ue )v

(ue − ∂x ue )v, for all v ∈ H 1 (Ω).

Let Vh be a fnite element space induced by a family of 2-simplices of type (k). By using Vh , we can write the following fnite element method for the boundary value problem (4.39): Find uh ∈ Vh such that Z

∇uh · ∇vh +

Ω

Z

Z

uh vh = ∂Ω

f vh + ZΩ

+

Z

Γ2

Γ0

(ue − ∂y ue )vh + Z

(ue + ∂y ue )vh +

Γ3

Z Γ1

(ue + ∂x ue )vh

(ue − ∂x ue )vh , for all vh ∈ Vh .

The Python function Compute_Robin for solving the above fnite element method is obtained by some straightforward modifcations of the Python function Compute_DN given in Section 4.6.2: def Compute_Robin(N, k, b): """ Given mesh division size N, the degree k of Lagrange elements, and the parameter b of the exact solution, this function returns the mesh size h and errors of the FE solution associated to a Robin boundary condtion """ # create mesh and define function space mesh = UnitSquareMesh(N, N) V_h = FunctionSpace(mesh, "Lagrange", k) # the exact solution and its first derivatives u_e = Expression(’pow(x[0]*x[0] + x[1]*x[1], b)’, b = b,

Conforming Finite Element Methods for PDEs

degree = 3 + k) Dx_u_e = Expression("""2*b*x[0]*pow(x[0]*x[0] + x[1]*x[1], b-1)""", b = b, degree = 3 + k) Dy_u_e = Expression("""2*b*x[1]*pow(x[0]*x[0] + x[1]*x[1], b-1)""", b = b, degree = 3 + k) # marking the boundary using a mesh function boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1) # mark bottom boundary edges as subdomain 0 class BottomEdges(SubDomain): def inside(self, x, on_boundary): tol = 1E-12 # tolerance for comparisons return on_boundary and abs(x[1]) < tol Gamma_B = BottomEdges() Gamma_B.mark(boundary_parts, 0) # mark right boundary edges as subdomain 1 class RightEdges(SubDomain): def inside(self, x, on_boundary): tol = 1E-12 # tolerance for comparisons return on_boundary and abs(x[0] - 1) < tol Gamma_R = RightEdges() Gamma_R.mark(boundary_parts, 1) # mark top boundary edges as subdomain 2 class TopEdges(SubDomain): def inside(self, x, on_boundary): tol = 1E-12 # tolerance for comparisons return on_boundary and abs(x[1] - 1) < tol Gamma_T = TopEdges() Gamma_T.mark(boundary_parts, 2) # mark left boundary edges as subdomain 3 class LeftEdges(SubDomain): def inside(self, x, on_boundary): tol = 1E-12 # tolerance for comparisons return on_boundary and abs(x[0]) < tol Gamma_L = LeftEdges() Gamma_L.mark(boundary_parts, 3) # define the weak formulation u_h = TrialFunction(V_h)

75

76

Finite Element Methods in Civil and Mechanical Engineering

v_h = TestFunction(V_h) f = Expression(’-4*b*b*pow(x[0]*x[0] + x[1]*x[1], b-1)’, b = b, degree = 3 + k) ds = Measure("ds", domain=mesh, subdomain_data=boundary_parts) # boundary integral A = inner(grad(u_h), grad(v_h))*dx + u_h*v_h*ds F = f*v_h*dx + (u_e - Dy_u_e)*v_h*ds(0) \ + (u_e + Dx_u_e)*v_h*ds(1) + (u_e + Dy_u_e)*v_h*ds(2) \ + (u_e - Dx_u_e)*v_h*ds(3) # compute finite element solution u_h = Function(V_h) # u_h is redefined solve(A == F, u_h) L2_Error = errornorm(u_e, u_h, norm_type="L2") H1_Error = errornorm(u_e, u_h, norm_type="H1") return L2_Error, H1_Error, mesh.hmax() Notice that in the defnition of the weak form, ds is used for the integration over the whole boundary while ds(0), ds(1), ds(2), and ds(3) denote the integration over specifc portion of the boundary. Also notice that the solve command is simply solve(A == F, u_h), as there is no Dirichlet boundary condition. By replacing Compute_Dirichlet with Compute_Robin in the program of Section 4.6.1, we can compute errors and convergence rates for the Robin boundary condition. For example, the output for k = 2, b = 1.25, and the meshes N = 3, 6, 9, 12, is: N = 3, L2_error = 8.64E-04, H1_error = 1.78E-02 N = 6, L2_error = 1.16E-04, H1_error = 4.71E-03 N = 9, L2_error = 3.56E-05, H1_error = 2.13E-03 N = 12, L2_error = 1.54E-05, H1_error = 1.21E-03 L2-convergence rate = 2.92, H1-convergence rate = 1.97 The reader is encouraged to compare this output with the corresponding outputs in Sections 4.6.1 and 4.6.2.

4.7 TIME-DEPENDENT PROBLEMS: PARABOLIC PROBLEMS In this section, we study a class of time-dependent PDEs. Suppose Ω is an open domain in Rn , n = 1, 2, 3, and let t indicate time in a time interval [0, T ]. By using the notation for the elliptic PDE (4.1), we consider the time-dependent PDE ∂t u − div(∇u · D) + bb·∇u + cu = f ,

x ∈ Ω, t ∈ [0, T ],

(4.40)

where u : Ω × [0, T ] → R is the time-dependent unknown function and the given data f (x,t), c(x,t), bb(x,t), and D(x,t) are also time-dependent. We assume that for all

77

Conforming Finite Element Methods for PDEs

w ∈ Rn , there exists β > 0 such that w · D(x,t) · wT ≥ β kwk2 , at any x ∈ Ω, and any t ∈ [0, T ].

(4.41)

The PDE (4.40) is then called parabolic. Example 4.3. The heat equation ∂t u − div(k(x)∇u) = f ,

x ∈ Ω, t ∈ [0, T ],

(4.42)

is an important example of parabolic PDEs, where u(x,t) is the temperature of the point x at time t, k(x) is the thermal conductivity, and f (x,t) is a source term. To solve (4.42), one may specify the temperature g(x,t) at the boundary ∂ Ω for all t and also the initial distribution of temperature at time zero. These can be stated as u(x,t) = g(x,t), x ∈ ∂ Ω, t ∈ [0, T ], (4.43) u(x, 0) = u0 (x), x ∈ Ω. More generally, to solve (4.40), we can impose any of the boundary conditions introduced in Section 4.2 for t ∈ [0, T ] together with the initial condition u(x, 0) = u0 (x), x ∈ Ω. As the following example demonstrates, weak formulations for timedependent parabolic PDEs can be obtained by the approach of Section 4.2. Example 4.4. Consider the PDE (4.40) subject to the boundary and the initial conditions (4.43) with g(x,t) = 0. By multiplying the PDE by an arbitrary v ∈ H01 (Ω) and using Green’s formula, one obtains the following weak formulation: Find u(x,t) such that u(·,t) ∈ H01 (Ω) for any t ∈ [0, T ] and R

R

Ω v∂t u + B(t, u, v) = Ω u(x, 0) = u0 (x),

f (x,t)v, for all v ∈ H01 (Ω),

where B(t, u, v) =

Z

b ∇u · D · (∇v)T + (b·∇u)v + cuv .

Ω

Generally speaking, weak formulations of (4.40) subject to the boundary conditions of Section 4.2 and an initial condition at t = 0 can be stated by the abstract problem: Find u(x,t) such that u(·,t) ∈ X for any t ∈ [0, T ] and R

Ω v∂t u + A(t, u, v) = F(t, v), for all v ∈ X, u(x, 0) = u0 (x),

(4.44)

where A(t, u, v) is linear with respect to u and v and F(t, v) is linear with respect to v. The data X, A, and F are defned similar to those of Section 4.3 by using the timedependent data D(x,t), bb(x,t), c(x,t), and f (x,t). We leave the proof as an exercise.

78

Finite Element Methods in Civil and Mechanical Engineering

4.7.1 FINITE ELEMENT APPROXIMATIONS USING THE METHOD OF LINES Suppose Th is a mesh of Ω and Xh is an H 1 -conformal fnite element space induced by Th with Xh ⊂ X, where X is the trial and the test space of the abstract problem (4.44). We consider the following spatial discretization of (4.44): Find uh (x,t) such that uh (·,t) ∈ Xh for any t ∈ [0, T ] and R Ω vh ∂t uh + A(t, uh , vh ) = F(t, vh ), for all vh ∈ Xh , (4.45) uh (x, 0) = Ih u0 (x), where Ih u0 ∈ Xh is the fnite element interpolation of u0 . The discrete problem (4.45) is equivalent to a linear system of ordinary differential equations (ODEs). To see this, notice that at any t ∈ [0, T ], we have uh (·,t) ∈ Xh , and therefore, we can write N

uh (x,t) =

∑ U j (t)ψ j (x),

(4.46)

j=1

where N = dim Xh , and {ψ1 , . . . , ψN } are the global shape functions of Xh . If uh solves (4.45), then Z

ψi ∂t uh + A(t, uh , ψi ) = F(t, ψi ),

i = 1, . . . , N.

Ω

By using the expansion (4.46), the above relations can be stated as N

∑ ∂t U j

j=1

Z Ω

N

ψ j ψi + ∑ U j A(t, ψ j , ψi ) = F(t, ψi ),

i = 1, . . . , N.

j=1

These equations form a system of coupled ODEs that admits the matrix form MN×N · ∂t UN×1 (t) = −[A]N×N (t) · UN×1 (t) + FN×1 (t), t ∈ [0, T ], (4.47) U(0) = U0 , where the symmetric matrix MN×N is the mass matrix with the components Mi j = Ω ψi ψ j , [A] is the time-dependent stiffness matrix with the components [A]i j = A(t, ψ j , ψi ), and the time-dependent vectors U and F are given by ⎡ ⎤ ⎡ ⎤ U1 (t) F(t, ψ1 ) ⎢ ⎥ ⎢ ⎥ . ... U(t) = ⎣ ⎦ , F(t) = ⎣ ⎦. ..

R

UN (t)

F(t, ψN )

The vector U0 contains the coeffcients of the expansion of the initial value Ih u0 in the basis {ψ1 , . . . , ψN }. By discretizing the system of ODEs (4.47) with respect to t, we can approximate its solution. In summary, to solve the problem (4.44), we frst employed a spatial discretization to obtain (4.47) and then we can employ a temporal discretization to solve the discrete time-dependent problem (4.47). This approach is called the method of lines.

Conforming Finite Element Methods for PDEs

79

4.7.2 TEMPORAL DISCRETIZATION Various fnite-difference approaches are available for the temporal dicretization of the problem (4.45) or equivalently, the time-dependent ODEs (4.47). Here, we discuss two simple approaches, namely, the explicit and the implicit Euler methods. Suppose {t0 ,t1 , . . . ,tm } is a set of points with t0 = 0 < t1 < t2 < · · · < tm = T , and let Δit = ti+1 − ti . Also let uih (x) = uh (x,ti ), that is, uhi is the value of uh at the time step ti . The explicit Euler method seeks the solution of (4.45) iteratively as follows: Let u0h = Ih u0 . Given uih ∈ Xh , fnd ui+1 h ∈ Xh such that 1 Δi t

Z Ω

i i (ui+1 h − uh )vh + A(ti , uh , vh ) = F(ti , vh ), for all vh ∈ Xh .

(4.48)

For the implicit Euler method, A and F are evaluated at ti+1 , that is: Let u0h = Ih u0 . Given uhi ∈ Xh , fnd uhi+1 ∈ Xh such that 1 Δit

Z Ω

i i+1 (ui+1 h − uh )vh + A(ti+1 , uh , vh ) = F(ti+1 , vh ), for all vh ∈ Xh . (4.49)

It is possible to show that the explicit Euler method (4.48) and the implicit Euler method (4.49) converge to the solution of the abstract problem (4.44). Implicit Euler iterations involve solving a linear system and are computationally more expensive than explicit Euler iterations. However, one can use much larger time intervals Δit using the implicit method which implies that fewer time steps are needed for the implicit Euler method. In particular, implicit methods are more suitable for stiff problems, where explicit methods require very fne time steps for accurate approximations. 4.7.3

IMPLEMENTATION: A DIFFUSION PROBLEM

We employ FEniCS to solve the following diffusion problem on the unit square Ω = (0, 1) × (0, 1): Find u(x, y,t) such that ⎧ ⎨ ∂t u − Δu = f , in Ω, t ∈ [0, T ], u = ue , on ∂ Ω, t ∈ [0, T ], (4.50) ⎩ u = ue , at t = 0, with f (x, y,t) = −4b2 (x2 + y2 )b−1 + ct c−1 , where c and b are constants, and ue (x, y,t) = (x2 + y2 )b + t c .

(4.51)

It is straightforward to show that ue (x, y,t) is the solution of (4.50). A weak formulation for this problem reads: Find u(x, y,t) satisfying u(·,t) ∈ H 1 (Ω) for any t ∈ [0, T ], with u = ue , on ∂ Ω, such that R R R 1 Ω v∂t u + Ω ∇u · ∇v = Ω f v, for all v ∈ H0 (Ω), (4.52) at t = 0. u = ue ,

80

Finite Element Methods in Civil and Mechanical Engineering

Let Vh be the fnite element space induced by 2-simplices of type (k) on a mesh Th of Ω and let Xh be the subspace of Vh defned in (4.24). By using the implicit Euler method with a uniform time increment S = ti+1 − ti , we can discretize (4.52) as: Let u0h = Ihk ue at t = 0. Given uhi ∈ Vh , fnd uhi+1 ∈ Vh satisfying uhi+1 |∂ Ω = Ihk ue , such that Z Ω

ui+1 h vh + S

Z Ω

∇ui+1 h · ∇vh = S

Z

Z

uih vh , for all vh ∈ Xh .

f vh + Ω

(4.53)

Ω

The above framework is implemented in the following program. """ This program solves a time-dependent diffusion problem. Given mesh division size N, the degree k of the Lagrange element, the parameters b and c of the exact solution, the final time T, and the time increment S, it computes the L2 error in each time step """ from dolfin import * import numpy as np # parameters b = 1.0; c = 1.0 k = 2 N = 4 T = 20 S = 0.5 error = []

# # # # # #

for the exact solution degree of Lagrange element mesh refinement final time time increment initializing L2-errors

# create mesh and define function space mesh = UnitSquareMesh(N, N) V_h = FunctionSpace(mesh, "Lagrange", k) # the exact solution u_e = Expression(’pow(x[0]*x[0] + x[1]*x[1], b) + pow(t, c)’, b = b, c = c, t = 0, degree = 3 + k) # define Dirichlet boundary condition def ue_boundary(x, on_boundary): return on_boundary BC = DirichletBC(V_h, u_e, ue_boundary) # initial condition u0 = interpolate(u_e, V_h) # define the weak formulation u_h = TrialFunction(V_h) v_h = TestFunction(V_h)

Conforming Finite Element Methods for PDEs

81

f = Expression("""-4*b*b*pow(x[0]*x[0] + x[1]*x[1], b-1) + c*pow(t, c-1)""", b = b, c = c, t = S, degree = 3 + k) A = u_h*v_h*dx + S*inner(grad(u_h), grad(v_h))*dx F = S*f*v_h*dx + u0*v_h*dx # compute solution u_h = Function(V_h) t = S while t 0 such that R qh div vh inf sup Ω ≥ αh . (5.23) qh ∈Ph v ∈V kqh k2 kv vh k1,2 h h In the literature, the above inf-sup condition is called the Babuˇska-Brezzi condition. Example 5.4. The choice of simplex of type (1) for both velocity and pressure leads to unstable mixed methods as it violates (5.23). For any integer k ≥ 2, the choice of simplex of type (k) for velocity and simplex of type (k − 1) for pressure leads to stable mixed method. This stable choice is called the Taylor-Hood element in the literature. A detailed analysis of mixed finite element methods for the Stokes equation is available in [6, Chapter 4] and [4, Chapter 8].

Figure 5.11 The geometry of the problem of flows over a rectangular step.

To implement the boundary value problem (5.20) in FEniCS, we consider a 2D flow over a rectangular step with the geometry shown in Figure 5.11. We assume that f = 0, and that the boundary velocity is zero on the top and the bottom boundaries of the domain. On the left and the right boundaries of the domain, we impose the parabolic velocity profile v (y) = 4U y(H − y), 0 , where H is the height of the H2 domain and U is the maximum velocity. This problem can be implemented using FEniCS as follows. # parameters n = 330 U = 15.0 mu = 1.0 k = 2

# # # #

number of divisions maximum velocity kinematic viscosity degree of polynomial

# generate mesh step = Rectangle(Point(0.4, 0.0), Point(0.8, 0.15)) Rec = Rectangle(Point(0.0, 0.0), Point(2.2, 0.4)) geometry = Rec - step mesh = generate_mesh(geometry, n) # define function spaces T = VectorElement("CG", mesh.ufl_cell(), k)

Applications

125

Q = FiniteElement("CG", mesh.ufl_cell(), k-1) mixed = T * Q W = FunctionSpace(mesh,mixed) # define subdomains inflow = ’near(x[0], 0)’ outflow = ’near(x[0], 2.2)’ topnbottom = ’near(x[1], 0) || near(x[1], 0.4)’ small_rec = """on_boundary && x[0] > 0.4001 && \ x[0] < 0.8001 && x[1] < 0.15001""" # pinpoint for removing the degree of freedom of pressure class PinPoint(SubDomain): def inside(self, x, on_boundary): return x[0] < DOLFIN_EPS and x[1] < DOLFIN_EPS pinpoint = PinPoint() # no-slip boundary condition for velocity at top and bottom noslip = Constant((0.0, 0.0)) bc0 = DirichletBC(W.sub(0), noslip, topnbottom) # no-slip boundary condition for velocity at step bc1 = DirichletBC(W.sub(0), noslip, small_rec) # inflow/outflow on left and right boundary velocity_profile = \ Expression((’4.0*U*x[1]*(0.4 - x[1]) / pow(0.4, 2)’,’0.0’), U = U, degree=4 + k) bc2 = DirichletBC(W.sub(0), velocity_profile, inflow) bc3 = DirichletBC(W.sub(0), velocity_profile, outflow) # pressure at pinpoint zero = Constant((0.0)) bc4 = DirichletBC(W.sub(1), zero, pinpoint, "pointwise") # all boundary conditions bcs = [bc0, bc1, bc2, bc3, bc4] # variational problem (v, p) = TrialFunctions(W) (w, q) = TestFunctions(W) f = Constant((0.0, 0.0)) a = mu*inner(grad(v), grad(w))*dx - div(w) *p*dx - q*div(v)*dx L = inner(f, w)*dx

126

Finite Element Methods in Civil and Mechanical Engineering

# solver V = Function(W) solve(a == L, V, bcs) # extracting sub-functions v, p = V.split() # saving solution in VTK format ufile_pvd = File("Ch5_Stokes/velocity.pvd") ufile_pvd >. In the followings, >>> at the beginning of a line indicates its execution in an interactive environment such a Python shell or a Jupyter notebook. To exit a Python shell type exit(). A Python shell is suitable for examining simple codes. Assume that the fle Test.py contains a Python program possibly containing several thousands of lines. To run this program, one types: $ python3 Test.py Files containing Python programs can be prepared by plain text editors or more advanced editors such as Spyder and PyCharm6 .

B.2

LISTS

A simple list can be defned as: >>> N = [4, 8, 12, 16, 20] The expression N[i-1] specifes the i-th element of a list. For example, the frst and the last elements of N can be printed as follows: >>> print(’first element =’, N[0]) first element = 4 >>> print(’last element =’, N[4]) last element = 20 Notice that the index of Python lists starts at 0. An alternative way to achieve the above results is >>> print(’first element =’, N[-5]) first element = 4 >>> print(’last element =’, N[-1]) last element = 20 6 https://www.jetbrains.com/pycharm/

Introduction to Python

151

It is also possible to access multiple elements of lists. For example, the following code prints the frst three elements of N. >>> print(’first three elements are’, N[0:3]) first three elements are [4, 8, 12] One can add an element to the end of a list: >>> N.append(24) >>> print(N) [4, 8, 12, 16, 20, 24] The number of elements of a list are determined by len(): >>> print(len(N)) 6 It is possible to assign new values to elements of a list as follows: >>> N[0] = -100 >>> print(N) [-100, 8, 12, 16, 20, 24] >>> N[0:3] = [-1, -2, -3] >>> print(N) [-1, -2, -3, 16, 20, 24]

B.3

BRANCHING AND LOOPS

Python has the standard if-else and if-elif-else blocks to branch the fow of programs based on given conditions. Here the keyword elif means “else if”. A simple example is given below. >>> >>> ... ... ... ... ...

A = 21; B = 11 if A > B: print("A is greater than B") elif A < B: print("A is smaller than B") else: print("A is equal to B")

A is greater than B Notice that unlike many programming languages, there is no keywords to specify the end of the block. In Python, lines with the same indentation level are considered to belong to the same block. Therefore, it is very important to pay attention to the indentation of each line when developing codes. The standard for and while loops are also available. For example, the following lines print elements of a list by using these loops.

152

>>> >>> ... ... 1 0.5 -5 >>> >>> ... ... ... 1 0.5 -5

Finite Element Methods in Civil and Mechanical Engineering

M = [1, 0.5, -5] for i in range(len(M)): print(M[i])

i = 0 while i < len(M): print(M[i]) i += 1 # update i

Here i += 1 is equivalent to i = i + 1. Also notice that comments begin with # in Python.

B.4

FUNCTIONS

A function is a tool for grouping statements such that they can be easily called several times in a program. In Python, the defnition of a function begins with def and results are sent back using return. As an example, consider the step function 1, if x ≥ 0, f (x) = 0, otherwise. The following lines defne this function and compute f (0.5). >>> def f(x): ... if x >= 0: ... y = 1 ... else: ... y = 0 ... return y ... >>> f(0.5) 1 As a general rule, all variables defned in a function are local to that function and they will not exist anymore after the function returns its value.

B.5

CLASSES AND OBJECTS

Python is an object-oriented programming language and extensively employs classes. Here we give a simple example of a class in Python. Suppose that we want

153

Introduction to Python

to defne the function � b g(x, y) = x2 + y2 + t c , where a, b, and t represent some physical parameters such as time that may change. To defne this function, one can use the following class G that consists of the variables b,c,t, and also the function value that returns the value of g at (x, y): >>> class G(object): ... def __init__(self, b, c, t): ... self.b = b ... self.c = c ... self.t = t ... def value(self, x, y): ... v = pow(x**2 + y**2, self.b) \ ... + self.t**self.c ... return v ... Suppose b = c = 1, and t = 0. Now, the function g can be defned as an object, also called an instance, of the class G: >>> g = G(1,1,0) To calculate g(2, 3), we use the statement >>> g.value(2,3) 13 Assume that the value of t changes to 1 and we want to compute g(2, 3) with this new value of t. This can be achieved by using the following lines: >>> g.t = 1 >>> g.value(2,3) 14 Functions and variables of a class are commonly called methods and data attributes, respectively. In the above example, G has the method value and the data attributes b,c,t. FEniCS has many classes. To be able to use them, one needs to import them at the beginning of a program: >>> from dolfin import * This way, all classes of FEniCS become available. For example, we can use the class UnitSquareMesh to defne a 10 × 10 mesh of a unit square: >>> mesh = UnitSquareMesh(10, 10)

154

Finite Element Methods in Civil and Mechanical Engineering

The class UnitSquareMesh has several data attributes and methods. To see them, simply type mesh. and press the tab bottom. This will show a long list that begins with >>> mesh. mesh.bounding_box_tree( mesh.cell_name( mesh.cell_orientations( mesh.cells( mesh.color( mesh.coordinates(

mesh.num_entities( mesh.num_entities_global( mesh.num_faces( mesh.num_facets( mesh.num_vertices( mesh.order(

As an example, we can employ the method hmax to calculate the maximum diameter of elements h: >>> mesh.hmax() 0.14142135623730964

B.6

READING AND WRITING FILES

Assume that we have the text fle Data.txt that contains velocities of some points of a mechanical system as follows: Point 1 2 3 4 5 6

Velocity(m/s) 5.5 6.7 6.3 9.1 1.0 3.0

Suppose we want to calculate the average of velocities and save it in another text fle Result.txt. First, we read Data.txt line by line and calculate the average average: >>> data = open(’Data.txt’, ’r’) # open the file for reading >>> data.readline() # skip the first line ’Point\tVelocity(m/s)\n’ >>> Num = 0; sum = 0 >>> for line in data: ... line_list = line.split() # list of strings for each line ... velocity = float(line_list[1]) ... sum += velocity; Num += 1 ... >>> data.close() # close the file >>> average = sum / Num

Introduction to Python

155

>>> average 5.266666666666667 The line data = open(’Data.txt’, ’r’) opens Data.txt for reading. Since the frst line does not contain any velocity, we skip it by using data.readline(). In the for loop, we read the fle line by line. The code line_list = line.split() splits the string line into words separated by blanks and creates the list of strings line_list of these words. For example, in the frst iteration, we have line_list = [’1’, ’5.5’]. To extract the velocity as a number, we select the second element of line_list and convert it to a foating point number by the line velocity = float(line_list[1]). After the for loop, we frst close the data fle and then we compute the average. To write the result, we create Result.txt and open it for writing: >>> result = open(’Result.txt’, ’w’) Then, we write the result and close the fle. >>> result.write(’average velocity = %.8f’ % average) >>> result.close() This will create the text fle Result.txt that contains average velocity = 5.26666667

B.7

NUMERICAL PYTHON ARRAYS

Numerical Python, abbreviated as NumPy, is a Python package designed for scientifc computations. In NumPy, we employ arrays instead of lists to store vectors and matrices. To defne a NumPy array, we frst need to import NumPy. The common statement for this purpose is >>> import numpy as np Numpy arrays can be defned similar to lists. For example, to defne the matrix 4 8 12 M= , 16 20 24 we employ the line >>> M = np.array([[4, 8, 12], [16, 20, 24]]) >>> M array([[ 4, 8, 12], [16, 20, 24]]) Elements of an array are accessed similar to lists:

156

Finite Element Methods in Civil and Mechanical Engineering

>>> M[0][0] 4 >>> M[0,0] 4 NumPy has several tools to defne more complex arrays. For example, an array for a uniform division of the interval [−2, 2] to 4 sub-intervals is obtained by >>> np.linspace(-2,2,5) array([-2., -1., 0., 1.,

2.])

A 2 × 2 zero matrix and a 3 × 2 matrix of ones can be defned as: >>> np.zeros((2,2)) array([[0., 0.], [0., 0.]]) >>> B = np.ones((3,2)) >>> B array([[1., 1.], [1., 1.], [1., 1.]]) Different operations for vectors and matrices are also available. For example, the matrix product M × B of the above matrices is computed as >>> np.dot(M,B) array([[24., 24.], [60., 60.]]) An important feature of Numpy is vectorization, that is, instead of using loops over arrays, many NumPy operations can be directly preformed on arrays. Vectorization is useful for speeding up heavy numerical computations. As an example, suppose that we want to compute f (x) = ex sin(x2 ), where x belongs to a uniform mesh of [0, 1] with 2, 000, 000 vertices. In the following, we preform a standard loop version and a vectorized version of this computation and compare the time spent on each version. >>> import time >>> import math >>> import numpy as np >>> x = np.linspace(0,1,2000000) >>> Time_Start = time.time();\ ... f_for = [math.exp(x[i])*math.sin(x[i]**2) \ ... for i in range(len(x))]; \ ... print("***** for version: %.8f seconds" ... % (time.time() - Time_Start)) ***** for version: 2.00471950 seconds

Introduction to Python

157

>>> Time_Start = time.time();\ ... f_vec = np.exp(x)*np.sin(x**2);\ ... print("***** vectorized version: %.8f seconds" ... % (time.time() - Time_Start)) ***** vectorized version: 0.14114642 seconds

B.8 PLOTTING WITH MATPLOTLIB There are several Python packages for plotting and visualization such as Matplotlib, Mayavi, and Easyviz. Matplotlib is the standard package for plotting curves. The common statement for importing this package is >>> import matplotlib.pyplot as plt We plot two simple examples in this section. First, we plot the function sin(x) on [0, 2π]. This can be achieved by the following lines: >>> import numpy as np >>> x = np.linspace(0, 2*np.pi, 40) >>> y = np.sin(x) >>> plt.plot(x,y) [] >>> plt.show() These lines produce a fgure similar to Figure B.1.

Figure B.1 Plot of sin(x) on [0, 2π].

Next, we simultaneously plot the functions sin(x) and cos(x) on [0, 2π] and add some information to the fgure: >>> x = np.linspace(0, 2*np.pi, 40)

158

Finite Element Methods in Civil and Mechanical Engineering

>>> S = np.sin(x) >>> C = np.cos(x) >>> plt.plot(x, S, ’-bo’, linewidth = ’4’, markersize = ’8’) [] >>> plt.plot(x, C, ’-rs’, linewidth = ’4’, markersize = ’8’) [] >>> plt.xlabel(’x’, fontsize = ’14’) Text(0.5, 0, ’x’) >>> plt.ylabel(’y’, fontsize = ’14’) Text(0, 0.5, ’y’) >>> plt.legend(["sin(x)","cos(x)"], loc = ’best’)

>>> plt.savefig(’Figure_2.png’) >>> plt.show() These lines are self-explanatory and yield a fgure similar to Figure B.2. Moreover, the line plt.savefig(’Figure_2.png’) saves the fgure as Figure_2.png.

Figure B.2 Plots of sin(x) and cos(x) on [0, 2π].

References 1. R.A. Adams and J.J.F. Fournier. Sobolev Spaces. Academic Press, Amsterdam, 2003. 2. A. Angoshtari and A. Gerami Matin. A conformal three-feld formulation for nonlinear elasticity: From differential complexes to mixed fnite element methods. arXiv preprint arXiv:1910.09025, 2020. 3. R.G. Bartle and D.R. Sherbert. Introduction to Real Analysis. John Wiley & Sons, Inc. New York, USA, 2011. 4. D. Boff, F. Brezzi, and M. Fortin. Mixed Finite Element Methods and Applications. Springer-Verlag, Berlin, 2013. 5. P.G. Ciarlet. The Finite Element Method for Elliptic Problems, volume 40. SIAM, Philadelphia, USA, 2002. 6. A. Ern and J. Guermond. Theory and Practice of Finite Elements. Springer-Verlag, New York, 2004. 7. L.C. Evans. Partial Differential Equations. American Mathematical Society, Providence, RI, 2010. 8. H.P. Langtangen. A Primer on Scientifc Programming with Python, volume 6. Springer, Berlin, 2016. 9. H.P. Langtangen and A. Logg. Solving PDEs in Python: The FEniCS Tutorial. Springer, 2017. 10. A. Logg, K.A. Mardal, and G. Wells. Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, volume 84. Springer Science & Business Media, 2012. 11. A.W. Naylor and G.R. Sell. Linear Operator Theory in Engineering and Science. Springer Science & Business Media, 2000. 12. B. Slatkin. Effective Python: 90 Specifc Ways to Write Better Python. Addison-Wesley Professional, 2019.

159

Index abstract problem, 60 coercivity, 61, 64 ellipticity, 61 variational structure, 61 well-posedness, 60, 64 advection-diffusion equation, 56, 108 affne family of fnite elements, 46 affne mapping, 8 affne-equivalent, 45 approximability property, 65 approximation space, 19 Babuˇska-Brezzi condition, 124 barycenter, 27 barycentric coordinates, 27 basis (Hamel), 9 bilinear form, 9 boundary condition, 56 Dirichlet, 57 essential, 58 mixed Dirichlet-Neumann, 59 natural, 59 Neumann, 58 Robin, 59 boundary value problem, 56 bounded set, 5 from above, 5 from below, 5 Box, 36 C´ea’s lemma, 65 CG, 30 conformal space, 38 conforming method, 62 consolidation problem, 117 conventional fnite element diagram, 28 convergence, 46, 64 convergence rate, 23, 47, 65 optimal, 24, 47, 66 curl, 8

degree of freedom global, 19, 38 local, 25, 26 diameter of element, 35 diffusion problem, 79 dimension, 9 DirichletBC, 68, 73, 87 divergence div, 7 dolfin, 22 ds, 68, 73 dual basis, 14 dual space, 13 duality brackets h · , · i, 14 dx, 68 elastic membrane, 56, 104 elasticity linearized, 56, 126, 130 nonlinear, 134 elliptic PDE, 55 error, 46, 64 errornorm, 23 Euler method, 79 explicit, 79 implicit, 79 existence of solution, 56 Expression, 22 FEniCS, 21, 66 installation, 147 fenics, 22 fnite element, 26 assembly process, 25, 37 Crouzeix-Raviart, 50 Hermit, 30 Hermite n-simplex of type (3), 30 Lagrange, 25, 27 ´ e´ lec, 33 Ned node, 27, 30 Raviart-Thomas, 32 161

162

simplex of type (k), 29 space, 37 fnite element method, 62 conforming, 62 mixed, 84, 123 fnite-dimensional, 9 FiniteElement, 22, 41 full rank, 10, 85 function, 5 continuous Cm (Ω), Cm (Ω), 6 bijective, 6 domain, 5 extension, 5 invertible, 6 one-to-one, injective, 6 onto, surjective, 6 range, 6 restriction, 5 functional, 13 FunctionSpace, 22, 41, 87 Gaussian elimination, 69 Gelerkin method, 62 generate_mesh, 36 grad, 68 gradient ∇, 14 Green’s formulas, 14

INDEX

N´ed´elec, 45 Raviart-Thomas, 44 interpolate, 23, 41 interpolation operator, 21 global Ih , 38 local IK , 26 N´ed´elec IKN , IhN , 34, 45 Raviart-Thomas IKRT , IhRT , 32, 44 simplex of type (k) IKk , Ihk , 29, 40 isoparametric family of fnite elements, 46 iterative method, 69 kernel, 8 Kronecker delta δi j , 14 Lagrange, 22, 25, 30 Laplacian Δ, 15, 122 Lax-Milgram lemma, 61 Lebesgue space L2 (Ω), 7 linear combination, 9 linear mapping, 8 linear space, 6 linear subspace, 8 linearly independent, 9 locally supported, 38 lower bound, 5

mass matrix, 78 maximum, 5 Measure, 73 mesh, 19, 34 affne, 35 inf-sup condition, 85 cell, 34 infmum, 5 element, 19, 34 infnite-dimensional, 9 geometrically conformal, 35 initial condition, 77, 92 number of edges ne , 35 initial-boundary value problem, 77, 79, number of elements nel , 35 92, 105 number of faces nf , 35 inner, 68 number of vertices nv , 35 inner product, 14 number of vertices on boundary n∂v , 85 interpolant, 26 global, 38 vertex, 19 Lagrange, 19, 27, 41 MeshFunction, 72 local, 24 method of lines, 78 h-type approach, 50 hat function, 20 heat transfer equation, 56, 77, 82 hyperbolic PDE, 92

163

INDEX

minimum, 5 mixed formulation, 83, 102, 123 MixedElement, 87 mshr, 36 N1curl, 34 nodal basis, 27 non-conforming method, 62 norm, 12 H 1 , k · k1,2 , 12 L2 , k · k2 , 12 k · kc , 13 k · kd , 13 normal derivative ∂n , 15 normed linear space, 12 null space, 8 Numerical Python, 155 NumPy, 155 on_boundary, 68 p-type approach, 50 parabolic PDE, 77 partly Sobolev class, 7 H(curl; Ω), 8 [H(curl; Ω)]n , 135 H(div; Ω), 7 [H(div; Ω)]n , 135 Point, 36 Poisson’s equation, 56, 66, 85, 104 Polygon, 36 polynomial space NE, 33 Pk (Rn ), 7 Qk (Rn ), 16 RT, 31 positive defnite matrix, 64 preconditioned Krylov solver, 69 project, 90 projection, 83, 89 Python programming language, 149 rank-nullity theorem, 10 real numbers, 5 Rectangle, 36

reference element, 35 reference fnite element, 45 regularity of solution, 56 Ritz method, 62 RT, 33 saddle-point variational structure, 84, 123 seepage problem, 112 set closure, 6 open, 6 shape function global, 21, 38 local, 25, 26 shape-regular, 46 simplex, 27 center of gravity, 27 edge, 27 face, 27 vertex, 27 size_t, 72 Sobolev space, 7 H m (Ω), 7 [H 1 (Ω)]n , 7 H01 (Ω), 8 [H01 (Ω)]n , 122 HD1 (Ω), 59 [HD1 (Ω)]n , 127 solution space, 57, 62 solve, 69 span of a set, 9 sparse LU decomposition, 69 sparse matrix, 63 Sphere, 36 split(), 88 stiffness matrix, 63 time-dependent, 78 Stokes equation for fuids, 144 strong solution, 57 SubDomain, 72 supremum, 5 Taylor-Hood element, 124 test function, 57, 62

164

test problem, 66 test space, 57, 62 TestFunction, 68 TestFunctions, 87 the global degrees of freedom, 38 three-point centered-difference formula, 92 implicit, 106, 130 time step, 79 trial space, 57, 62 TrialFunction, 68 TrialFunctions, 87 triangulation, 35 UFL, 68 uniqueness of solution, 56 unit n-simplex, 27, 51 unit ball, 6 unit cube, 36 unit sphere, 6 unit square, 15, 36 UnitCubeMesh, 36 UnitIntervalMesh, 22 UnitSquareMesh, 36 upper bound, 5 variational crimes, 94 variational structure, 61 vector space, 6 VectorFunctionSpace, 42 vectorization, 156 wave equation, 93, 105 weak formulation, 56, 57 weak solution, 56, 57 well-posedness, 60, 64

INDEX