Hierarchical a posteriori error estimation of Bank-Weiser type in the FEniCS Project
Raphaël Bulle, Jack S. Hale, Alexei Lozinski, Stéphane P. A. Bordas, Franz Chouly
HHierarchical a posteriori error estimation of Bank–Weisertype in the FEniCS Project ∗Raphaël Bulle † Jack S. Hale † Alexei Lozinski ‡ Stéphane P. A. Bordas † Franz Chouly § February 9, 2021
Abstract.
In the seminal paper of Bank and Weiser [
Math. Comp. , 44 (1985), pp. 283–301]a new a posteriori estimator was introduced. This estimator requires the solution of a localNeumann problem on every cell of the finite element mesh. Despite the promise of Bank–Weiser type estimators, namely locality, computational efficiency, and asymptotic sharpness,they have seen little use in practical computational problems. The focus of this contributionis to describe a novel implementation of hierarchical estimators of the Bank–Weiser type ina modern high-level finite element software with automatic code generation capabilities. Weshow how to use the estimator to drive (goal-oriented) adaptive mesh refinement and to mixedapproximations of the nearly-incompressible elasticity problems. We provide comparisonswith various other used estimators. An open source implementation based on the FEniCSProject finite element software is provided as supplementary material.
A posteriori error estimation [3] is the de facto tool for assessing the discretization error of fi-nite element method (FEM) simulations, and iteratively reducing that error using adaptive meshrefinement strategies [56].This paper is concerned with the description and justification of an implementation of an errorestimator introduced in the seminal paper of Bank and Weiser [13, Section 6]. In that paper anerror estimate was derived involving the solution of local Neumann problems on a special finiteelement built on nested or hierarchical spaces. Despite its excellent performance and low compu-tational cost, this estimator has seen relatively sparse use in practical computational problems.The overarching goal of this contribution is to provide access to an efficient, generic and extensi-ble implementation of Bank–Weiser type estimators in a modern and widely used finite elementsoftware, specifically, the FEniCS Project [5].
The literature on a posteriori error estimation and adaptive finite element methods is vast, so wefocus on articles on practical software implementations of adaptive finite element methods and ∗ R.B. would like to acknowledge the support of the ASSIST research project of the University of Luxembourg.This publication has been prepared in the framework of the DRIVEN project funded by the European Union’sHorizon 2020 Research and Innovation programme under Grant Agreement No. 811099. F.C.’s work is partiallysupported by the I-Site BFC project NAANoD and the EIPHI Graduate School (contract ANR-17-EURE-0002). † Institute of Computational Engineering, University of Luxembourg, 6 Avenue de la Fonte, 4362 Esch-sur-Alzette,Luxembourg ( [email protected] , [email protected] , [email protected] ) ‡ Laboratoire de Mathématiques de Besançon, UMR CNRS 6623, Université de Bourgogne Franche-Comté, 16route de Gray, 25030 Besançon Cedex, France ( [email protected] ) § Université de Bourgogne Franche-Comté, Institut de Mathématiques de Bourgogne, 21078 Dijon, France( [email protected] ) a r X i v : . [ m a t h . NA ] F e b omparative performance studies.The T-IFISS [21] software package, based on the existing IFISS [34] package, is a finite elementsoftware written in MATLAB/Octave with a focus on a posteriori error estimation and adaptivefinite element methods. Recently [20], T-IFSS has been extended to solve adaptive stochasticGalerkin finite element methods. The stated emphasis of T-IFISS [21] is on being a laboratoryfor experimentation and exploration, and also to enable the rapid prototyping and testing of newadaptive finite element methods. A number of estimation and marking strategies are implementedin T-IFISS, although not the Bank–Weiser estimator we consider in this paper. T-IFISS only worksfor two-dimensional problems and it was never intended to be a high-performance code suitablefor large-scale computations e.g. high-performance computing systems using the Message PassingInterface (MPI).The PLTMG package [12] is one of the oldest open finite element softwares for solving ellipticproblems that is still under active maintenance, and includes many advanced features such as hp -adaptive refinement, a posteriori error estimation, domain decomposition and multigrid precondi-tioning. The a posteriori error estimation is based on a superconvergent patch recovery estimationtechnique introduced in [14]. PLTMG only works in two dimensions and is naturally limited froma usability perspective due to the programming tools available at its inception (Fortran and ANSIC). In [36] an adaptive first-order polynomial finite element method was implemented in a codecalled p1afem using MATLAB. The primary goal was to show how the basic finite element algorithmcould be implemented efficiently using MATLAB’s vectorization capabilities. A standard residualestimator [9] is used to drive an adaptive mesh refinement algorithm. Again, like T-IFISS, p1afemonly works in two dimensions.In [63] a novel methodology for automatically deriving adaptive finite element methods fromthe high-level specification of the goal functional and (potentially non-linear) residual equation wasimplemented in the FEniCS Project. The emphasis of the paper [63], in contrast with the T-IFISStoolbox [21], is on the automatic construction of goal-oriented adaptive finite element methods,without much knowledge required on the part of the user. The implicit residual problems areautomatically localised using bubble functions living on the interior and facets of the cell, andthe dual problem [38] is derived and solved automatically on the same finite element space as theprimal problem, before being extrapolated to a higher-order finite element space using a patch-wise extrapolation operator. In practice the automatically derived estimators seem to be able toeffectively drive adaptive mesh refinement for a range of different PDEs.Explicit residual estimators are also commonly employed by users of high-level finite elementsoftware packages as they can usually be expressed straightforwardly in a high-level form language,e.g. [5, 59]. For example, [41] used the FEniCS Project to implement an explicit residual errorestimator for the Reissner-Mindlin plate problem from [19]. The authors of [31] used the FEniCSProject to implement an explicit residual estimator for elasticity problems within a dual-weightedresidual framework. The dual problem is solved on a higher-order finite element space in order toensure that the weighting by the dual residual solution does not vanish [63]. In [46] the authors usean explicit dual-weighted residual strategy for adaptive mesh refinement of discontinuous Galerkinfinite element methods. In addition, as the name suggests, they can be explicitly computed asthey involve only functions of the known finite element solution and the problem data.In the present work, aside of the Bank–Weiser estimator we will consider an explicit residualestimator [10] named residual estimator in the following, a flux reconstruction based on averagingtechnique estimator [71], referred to as Zienkiewicz–Zhu estimator , and a variant of the Bank–Weiser estimator introduced in [68] and referred to as the bubble Bank–Weiser estimator . Theresidual estimator was proved to be both reliable and (locally) efficient in [68] for any finite elementorder and in any dimension. The proof of reliability and (local) efficiency of Zienkiewicz–Zhuestimator has been derived in [62], for linear finite elements in dimension two and generalised to2ny averaging technique in any dimension in [26] and any finite element order in [15]. The bubbleBank–Weiser estimator was proved to be reliable and locally efficient in [68] for any dimension andany finite element order.A proof of the equivalence between the Bank–Weiser estimator and the exact error was derivedin the original paper [13]. However, this proof requires a saturation assumption [13, 33, 55] askingfor the best approximation with higher order finite elements to be strictly smaller than that oflower order elements and which is known to be tricky to assert in practice. Some progress has beenmade in [55] removing the saturation assumption from the analysis. However, this progress wasmade at the price of restricting the framework to linear polynomial finite elements and dimensiontwo only. The equivalence proof between Bank–Weiser and residual estimators have been improvedby the authors in [22] where it was extended to dimension three.
We show how robust and cheap hierarchical error estimation strategies can be implemented in ahigh-level finite element framework, e.g. FEniCS Project [5], Firedrake [37, 60], freefem ++ [43],Feel ++ [59], GetFEM [61] or the C ++ library Concha [30]. Specifically, the contribution of ourpaper to the existing literature is:• A generic and efficient implementation of the Bank–Weiser estimator in the open sourceFEniCS Project finite element software that works for Lagrange finite elements of arbitrarypolynomial order and in two and three spatial dimensions. The code is released under anopen source (LGPLv3) license [23]. Because the code utilises the existing automatic codegeneration capabilities of FEniCS along with a custom finite element assembly routine, thepackage is very compact (a few hundred lines of code, plus documentation and demos).Additionally, the estimators are be implemented in near mathematical notation using theUnified Form Language, see the Supplementary Material for code snippets.• A numerical comparison of the Bank–Weiser estimator with various estimators we mentionedearlier. We examine the relative efficiency, and their performance within an adaptive meshrefinement loop on various test problems. Unlike [29], we do not aim at running a competitionof error estimators but at stressing out the potential of the Bank–Weiser estimator since, asthe authors of [29] point out, a single error estimation strategy is not sufficient to cover theparticulars of all possible problems.• Relying on results in [17], we show a goal-oriented adaptive mesh refinement algorithm canbe driven by weighted sum of estimators, computed separately on primal and dual problemsdiscretized on the same finite element space. This avoids the extrapolation operation of [63]or the need to compute the dual solution in a higher-order finite element space [18].• Using the same basic methodology as for the Poisson problem, we extend our approach toestimating errors in mixed approximation of nearly incompressible elasticity problems. Thisidea was originally introduced in [3] and is still an active research topic, see e.g. [47] for aparameter-robust implicit residual estimator for nearly-incompressible elasticity. An outline of this paper is as follows; in section 1.4 we outline the main notation and definitionsused in this paper. In sections 2 and 3 we show the derivation of the primal problem and theBank–Weiser error estimator. In section 4 we derive a new method for computing the Bank–Weiser estimator and discuss its implementation in FEniCS. In section 5 we discuss the use of theapproach for various applications such as goal-oriented adaptive mesh refinement and for mixed3pproximations of PDEs. Then, in section 6 we show some results on two and three dimensionalPoisson test problems as well as on linear elasticity problems, before concluding in section 8.
In this section we outline the main notations used in the rest of the paper. Let Ω be an boundedopen domain of R d ( d = 1 , or ), with polygonal/polyhedral boundary denoted by Γ := ∂ Ω .We consider Γ = Γ D ∪ Γ N a partition of the boundary. We denote by n : Γ → R d the outwardunit normal vector along Γ . Let ω be a subset of Ω . For l in R we denote by H l ( ω ) the Sobolevspace of order l . The space H ( ω ) = L ( ω ) is the Lebesgue space of square integrable functionsover ω . The space H l ( ω ) is endowed with the usual inner product ( · , · ) l,ω and norm k·k l,ω . Weomit the subscript l when l = 0 and subscript ω when ω = Ω . We denote H D (Ω) the subspaceof H (Ω) of functions with zero trace on Γ D . We make use of the notation ∂ n v := ∇ v · n for thenormal derivative of a smooth enough function v . For l in R and for a d -dimensional subset ω of Ω , we also define the following vector fields spaces L ( ω ) := (cid:0) L ( ω ) (cid:1) d and H l ( ω ) := (cid:0) H l ( ω ) (cid:1) d , with respective inner products defined as their scalar counterparts, replacing the scalar productby Euclidean inner product or double dot product. The space H D (Ω) is the subspace of H (Ω) offunctions with zero trace on Γ D . From now on, the bold font notation will be reserved to vectorfields. With these notations at hand we can proceed with the rest of the paper. We consider the Poisson problem with mixed Dirichlet and Neumann boundary conditions. Let
Γ = Γ D ∪ Γ N be a partition of the boundary. We apply a Dirichlet boundary condition on Γ D and a Neumann boundary condition on Γ N . Let f ∈ L (Ω) , u D ∈ H / (Γ D ) and g ∈ L (Γ N ) beknown data. We seek a function u : − ∆ u = f in Ω , u = u D on Γ D , ∂ n u = g on Γ N . (1)Problem eq. (1) can be written in an equivalent weak form: Find u ∈ H (Ω) of trace u D on Γ D such that ( ∇ u, ∇ v ) = ( f, v ) + ( g, v ) Γ N , ∀ v ∈ H D (Ω) . (2)The weak problem eq. (2) can be discretized using the Lagrange finite element method. We take amesh T of the domain Ω , consisting of cells T = { T } , facets E = { E } (we call facets the edges indimension two and the faces in dimension three), and vertices N = { χ } . The mesh T is supposedto be regular, in Ciarlet’s sense: h T /ρ T γ, ∀ T ∈ T , where h T is the diameter of a cell T , ρ T the diameter of its inscribed ball, and γ is a positive constant fixed once and for all. The subsetof facets in the interior of the mesh (i.e. those that are not coincident with the boundary Γ ) isdenoted E I . The subset of facets lying on Γ D is denoted E D . The subset of facets lying on Γ N isdenoted E N . The subset of facets lying on the boundary of the domain Γ is denoted E B = E D ∪ E N .Let n + ∈ R d and n − ∈ R d be the outward unit normals to a given edge as seen by two cells T + and T − incident to a common edge E . If we denote P k ( T ) the space of polynomials of order k ona cell T , the continuous Lagrange finite element space of order k on the mesh T is defined by V k := (cid:8) v k ∈ H (Ω) , v k | T ∈ P k ( T ) ∀ T ∈ T (cid:9) . (3)We denote V kD the finite element space composed of functions of V k vanishing on the boundary Γ D . We consider the finite element problem: Find u k ∈ V k such that u k = w on Γ D and: ( ∇ u k , ∇ v k ) = ( f, v k ) + ( g, v k ) Γ N , ∀ v k ∈ V kD , (4)and where w is a discretization of u D on V k . 4 The Bank–Weiser estimator
In this section we derive the general definition of the Bank–Weiser estimator from the equation ofthe error as it was given in the original paper [13]. We also give a concrete example of Bank–Weiserestimator for linear finite elements.
We are interested in estimating the error we commit by approximating the solution u by u k in V kw .We define this error by the function e := u − u k and we want to estimate its norm k∇ e k . Thefirst step towards this will be to derive a new variational problem for which the exact error e isthe solution. For a cell T of the mesh, we introduce the interior residual as r T := ( f + ∆ u k ) | T , (5)and for an edge E , the edge residual J E = if E ∈ E D , [[ ∂ n u k ]] E if E ∈ E I , ( g − ∂ n u k ) | E if E ∈ E N . (6)where the notation [[ v ]] E := v + − v − denotes the jump in the value of the function across an interiorfacet E ∈ E I . Here, v + and v − denote the values of v on the facet E as seen by the two incidentcells T + and T − , respectively. The error function e satisfies what we call the global error equation ( ∇ e, ∇ v ) = X T ∈T ( r T , v ) T + 2 X E ∈E I ( J E , v ) E + X E ∈E N ( J E , v ) E ∀ v ∈ H D (Ω) . (7) We introduce now local finite element spaces in order to derive the finite element approximationof the error. For a cell T of the mesh we define V kT := (cid:8) v k,T ∈ P k ( T ) , v k,T = 0 in (Ω \ T ) ∪ ( T ∩ Γ D ) (cid:9) . (8)A key idea in the Bank–Weiser estimator derivation is to introduce an appropriate finite elementspace for the discretization of error. This non-standard space has two roles. Firstly, for the localproblems involving the cells with facets only in the interior of the domain or on the Neumannboundary, it should remove the constant functions, giving a unique solution. Secondly, and as wewill notice in section 6, solving the local error equation on the finite element space V kT / R does notnecessary lead to an accurate estimation of the error. However, in some cases, the estimation ofthe error can be surprisingly accurate when the space is judiciously chosen. We refer the reader to[1] for a full discussion.Before introducing this non-standard space, we need some more notations. Let k + and k − betwo non-negative integers such that k + > k − > . Let e T be the reference cell fixed once for all(independent from the mesh T ). We denote L e T : V k + e T −→ V k + e T , Im( L e T ) = V k − e T , (9)the Lagrange interpolation operator between the local spaces V k + e T and V k − e T ⊂ V k + e T . Moreover, forany cell T of the mesh, there exists an affine bijection S : e T −→ T e x S ( e x ) =: x (10)5apping e T onto T . From the mapping S we deduce another mapping given by S : V k + T −→ V k + e T v ( x ) ( v )( e x ) := v ( S ( e x )) . (11)If we denote d + the dimension of V k + e T and d − the dimension of V k − e T , given B + e T := { e ϕ , · · · , e ϕ d + } the basis of shape functions of V k + e T and B + T := { ϕ , · · · , ϕ d + } the basis of V k + T , we can always finda mapping S (and a mapping S ) such that S ( ϕ T,i ) = S ( e ϕ T,i ) ∀ i ∈ { , · · · , d + } , (12)We choose S and S so. For a given cell T of the mesh, we define the Lagrange interpolationoperator on T as follows L T := S − ◦ L e T ◦ S . (13)Note, due to (12), the matrix of S in the couple of basis ( B + T , B + e T ) is the identity matrix of size d + × d + . Consequently, if we denote G the matrix of L T in the basis B + T and e G the matrix of L e T in the basis B − e T , we have G = Id − e G Id = e G. (14)For a cell T of the mesh, the local Bank–Weiser space V bw T is defined as the null space of L T , inother words V bw T := ker( L T ) = n v T ∈ V k + T : L T v T = 0 o . (15)With this new space in hands, we can derive a local discrete counterpart of equation eq. (7) onany cell T : Find e T ∈ V bw T such that: (cid:0) ∇ e bw T , ∇ v bw T (cid:1) = (cid:0) r T , v bw T (cid:1) + X E ∈ ∂T (cid:0) J E , v bw T (cid:1) E ∀ v bw T ∈ V bw T . (16)The edge residual J definition takes into account the error on the Neumann boundary data ap-proximation, but it does not take into account the error on the Dirichlet boundary data. As aconsequence, the Bank–Weiser estimator does not take into account the contribution of the Dirich-let boundary data approximation. For a detailed discussion on a priori and a posteriori errorestimation with inhomogeneous Dirichlet boundary conditions see [16]. Finally, on the cell T thelocal Bank–Weiser estimator η bw ,T is defined by η bw ,T := k∇ e bw T k T , (17)where e T is defined in eq. (16) and the global Bank–Weiser estimator by the sum of local estimates η := X T ∈T η ,T . (18) If we assume k = 1 (i.e. we solve eq. (4) using linear finite elements) one can define the space V bw T from the choice of k + = 2 , k − = k = 1 . This example was the case considered in the numericaltests of the original paper [13]. The space V bw T consists of quadratic polynomial functions (in V T )vanishing at the degrees of freedom of the standard linear finite element functions (in V T ) i.e. thedegrees of freedom associated with the vertices of T .6 Algorithms and implementation details
The linear system corresponding to eq. (16) is not accessible in FEniCS this prevent us fromdirectly solving the Bank–Weiser equation. We propose to bypass the problem by constructingthe linear system corresponding to eq. (16) from another linear system derived from finite elementspaces that are accessible directly in FEniCS.
1. We consider the following singular value decomposition (SVD) of GG = U Σ V T , (19)where Σ is a diagonal matrix composed of the singular values of G . The columns of thematrix V are singular vectors of G , associated with singular values. The columns associatedwith the singular values zero span the null space of G . We take the submatrix N made ofthe columns of V spanning the null space of G . Note that, since G does not depend on anycell T , the same property holds for N .2. We build the matrix A + T and vector b + T of the local linear system corresponding to thefollowing variational formulation in the space V k + T , available in FEniCS: (cid:0) ∇ e + T , ∇ v + T (cid:1) = (cid:0) r T , v + T (cid:1) + X E ∈ ∂T (cid:0) J E , v + T (cid:1) E ∀ v + T ∈ V k + T . (20)3. We construct the matrix A bw T and vector b bw T as follow A bw T = N T A + T N and b bw T = N T b + T , (21)where A bw T and b bw T are the matrix and vector which allow to recover the bilinear and linearforms of eq. (16) in a basis of V bw T .4. We solve the linear system A bw T x bw T = b bw T , (22)5. We bring the solution back to V k + T , considering N x bw T , in order to post-process it and computethe local contribution of the Bank–Weiser estimator eq. (17). The global estimator eq. (18)is the square root of the sum of the squared local contributions. We now give more details specific to our implementation in FEniCS of each one of the above steps.1.
Computation of N . This is the key point of our implementation. The operator L T can bewritten as follows: L T : V k + T −→ V k − T −→ V k + T v + ( v + ) (cid:0) G ( v + ) (cid:1) . (23)Then, the matrix G is obtained via the following product G = G G , (24)where G and G are respectively the matrix in the couple of basis ( B + T , B − T ) of the Lagrangeinterpolation operator from V k + T to V k − T , denoted G and the matrix in the same couple of7 ompute G = G G SVD of G to obtain V Extract N from V Computation of N
Compute A + T and b + T Compute A bw T and A bw T Solve eq. (16)Compute localBW estimator
Computation ofthe local estimators Computation ofthe global estimator L oo p o v e r ce ll s Figure 1: Overall process of the Bank–Weiser estimator algorithm.basis of the canonical injection of V k − T into V k + T , denoted G . The matrices G and G can becalculated either using the Finite Element Automatic Tabulator (FIAT) [48] or, as we chooseto do, using the interpolator construction functions of the DOLFIN finite element library [51].The next step consists in computing the unitary matrix V of right singular vectors of G . Thiscomputation is done using the singular value decomposition (SVD) algorithm available in theSciPy library [69]. We can write the matrix V as follows, V = (cid:0) ξ | · · · | ξ d bw | ξ | · · · | ξ d − (cid:1) , (25)where B bw T := { ξ , · · · , ξ d bw } is the set of singular vectors of G corresponding to a zero singularvalue, spanning V bw T and { ξ , · · · , ξ d − } is spanning the supplementary space. The matrix N isthen chosen as the submatrix of V , keeping only the columns from B bw T : N := (cid:0) ξ | · · · | ξ d bw (cid:1) . (26)The linear algebra operations needed to form the submatrix N from V are performed using theNumPy library [66].2. Computation of A + T and b + T . The equation eq. (20) is expressed directly in the Unified FormLanguage (UFL) [6] and efficient C ++ code for calculating the cell local tensors A + T and b + T fora given cell T is then generated using the FEniCS Form Compiler (FFC) [49, 73]. If the cell T has an edge on a Dirichlet boundary E D , the matrix A + T and vector b + T must be modified inorder to enforce the boundary condition.3. Computation of A bw T and b bw T . The matrix A bw T and vector b bw T are constructed using eq. (21).4. Solution of the linear system (22) . The linear system eq. (22) is solved using a partial-pivot LUdecomposition algorithm from the Eigen dense linear algebra library [40].5.
Computation of the Bank–Weiser estimator.
Finally, the solution x bw T is sent back to V k + T using N and the norm of the corresponding function, giving the local estimator eq. (17) is com-puted using standard high-level functions already available within FEniCS. The global estimatoreq. (18) is computed using the information of all the local contributions.8 .3 Additional remarks The custom assembler composed of steps 2.-5. is performed by looping over every cell of the meshand, by virtue of using the abstractions provided by DOLFIN, works in parallel on distributedmemory computers using the Message Passing Interface (MPI) standard. For performance reasonsthese steps have been written in C ++ and wrapped in Python using the pybind11 library so thatthey are available from the Python interface to DOLFIN. In contrast, the first step must only beperformed once since the matrix N is the same for every cell of the mesh.Because we use the automatic code generation capabilities of FEniCS, our approach can bereadily applied to other definitions for the spaces V k + T and V k − T , and to vectorial problems likelinear elasticity, as we will see in the next section.In the numerical results secton we compare several versions of the Bank–Weier estimator andespecially the one we call bubble Bank–Weiser estimator and denote η bT which can be obtainedwith our method by taking V + T as the space V T + Span { ψ T } (the local space of quadratic functionsenriched with the space spanned by the interior bubble function) and V − T as V T . The resultingspace V bw T is spanned by the interior bubble function and the edges bubbles functions of the cell T . Note that this is not the only method we can use to compute the Bank–Weiser estimator, wecould have used e.g. Lagrange multipliers or a penalty method to define the space V bw T . In this section we show a number of applications, including adaptive mesh refinement, goal-oriented estimation and extensions to more complex mixed finite element formulations for thenearly-incompressible elasticity problems.
As well as simply providing an estimate of the global and local error, the estimator can be usedto drive an adaptive mesh refinement strategy. In the following we compare different refinementstrategy all based on the following loop:... → SOLVE → ESTIMATE → MARK → REFINE → ...The loop can be terminated once a given criterion e.g. maximum number of iterations, or globalerror less than a given tolerance, has been reached. A detailed discussion on adaptive refinementmethods can be found in [56]. In the following we expand on the specific algorithms used in ourcase. The weak form eq. (2) is discretized using a standard finite element method implemented withinFEniCS. The resulting linear systems are solved using the appropriate algorithms available withinPETSc [11], e.g. conjugate gradient method preconditioned with Hypre BoomerAMG [35], or directmethods, e.g. MUMPS [7, 8].
The Bank–Weiser estimator η bw is formulated and implemented as described in section 4. Thelocal contributions of the estimator provide an estimate of the local error for each cell in the meshand are subsequently used to mark the mesh. In addition the global estimator can be used todetermine when to stop iterating. 9 .1.3 Mark We have used two distinct marking strategies throughout the results section: the maximum strategyon the three-dimensional test cases and Dörfler strategy on the two-dimensional ones. We followthe presentation in [57]. In the maximum marking strategy [9], a cell is marked if its indicator isgreater than a fixed fraction of the maximum indicator. More precisely, given a marking fraction θ ∈ (0 , , the marked set M ⊂ T is the subset such that: η bw ,T ≥ θ max T ∈T η bw ,T ∀ T ∈ T . (27)In the Dörfler marking strategy [32] (sometimes referred to as the equilibrated marking strategy)enough elements must be marked such that the sum of their estimators is larger than a fixed fractionof the total error. Given a marking fraction θ ∈ (0 , , the marked set M is the subset with minimalcardinality M such that X T ∈M η ,T ≥ θ X T ∈T η ,T . (28)We implement an O ( N log N ) with N := T complexity algorithm for finding the minimumcardinality set by sorting the indicators in decreasing order and finding the cutoff point such thateq. (28) is satisfied. Because of the ordering operation this set is guaranteed to have minimalcardinality. We note that recent work [44, 57] proposes a O ( N ) complexity algorithm for findingthe set with minimum cardinality. We use two-dimensional and three-dimensional variants of the algorithm proposed in [58], some-times referred to as the Plaza algorithm. This algorithm works by subdividing the facets of eachmarked triangle or tetrahedron cell and then subdividing each triangle or tetrahedral cell so thatit is compatible with the refinement on the facets. The algorithm has O ( M ) complexity in thenumber of added mesh vertices M . This algorithm already exists in DOLFIN [51] and was usedfor the numerical results in [63]. In many practical applications it is desirable to control the error in a specific quantity of interest,rather than the (global, i.e. across the entire domain Ω ) energy norm [18]. In this section we showhow the basic Bank–Weiser estimator can be used to control error in a goal functional, rather thanin the natural norm. To do this, we use a weighted marking strategy proposed in [17].Let J : L (Ω) → R be a given linear functional. Associated with J ( u ) and the primal problemeq. (2) is the dual or adjoint problem: Find the dual solution z ∈ H D (Ω) such that ( ∇ v, ∇ z ) = J ( v ) ∀ v ∈ H D (Ω) . (29)The dual problem, like the primal problem, can also be approximated using the finite elementmethod. Find z k ∈ V k such that ( ∇ v k , ∇ z k ) = J ( v k ) = ( c, v k ) + ( h, v k ) Γ ∀ v k ∈ V k . (30)Using Galerkin orthogonality and Cauchy-Schwarz, it follows that |J ( u ) − J ( u k ) | = | ( ∇ ( u − u k ) , ∇ z ) | (31) = | ( ∇ ( u − u k ) , ∇ ( z − z k )) | (32) ≤ k∇ ( u − u k ) kk∇ ( z − z k ) k , (33)10here the inequality holds due to Galerkin orthogonality.Approximating the primal and dual errors k∇ ( u − u k ) k and k∇ ( z − z k ) k with any estimators η u and η z respectively, gives us an estimator for the error in the goal functional | J ( u ) − J ( u k ) | asthe product of η u and η z , thanks to eq. (33): η w := η u η z (34)In addition, if η u and η z are reliable estimators i.e. if there exist two constants C u and C z onlydepending on the mesh regularity such that k∇ ( u − u k ) k ≤ C u η u , and k∇ ( z − z k ) k ≤ C z η z , (35)then, η w is reliable as well | J ( u ) − J ( u k ) | C u C z η w . (36)Note that because the error in the goal functional is bounded by the product of two estimates,the element marking strategy must incorporate information from local indicators for both approx-imations to reduce the error on refinement. There are multiple strategies for doing this in theliterature, see e.g. [54]. We have chosen to implement the weighted goal-oriented (WGO) markingstrategy from [17]. The local WGO estimator is then defined as η w,T := η z η u + η z η u,T + η u η u + η z η z,T ∀ T ∈ T . (37)The marking and refinement using η w,T then follows in exactly the same manner as in the standardadaptive refinement strategy. Our implementation of the Bank–Weiser estimator can be directly applied to mixed formulations of(nearly-incompressible) linear elasticity problems using the results in [47]. In [2] a new a posteriorierror estimator is introduced for mixed formulations of Stokes problems consisting in solving alocal Poisson problem based on the local residuals on each cell. This estimator has been provedto be reliable and efficient in [2] under a saturation assumption. This assumption has been laterremoved in [50]. The reliability and efficiency of the estimator for mixed formulations of linearelasticity is proved in [47] without the need of a saturation assumption. In addition, they showthat the estimator is robust in the incompressible limit.
We consider the problem of linear deformation of an isotropic elastic solid Ω using the Herrmannmixed formulation. We consider the stress tensor σ : Ω → R d × d , the strain tensor ε : Ω → R d × d ,the load f : Ω → R d which belongs to (cid:0) L (Ω) (cid:1) d , the Dirichlet boundary data u D in (cid:0) H / (Γ D ) (cid:1) d ,the Neumann boundary condition (traction) data g in (cid:0) L (Γ N ) (cid:1) d and displacement field u : Ω → R d . The stress and strain tensors are defined by σ := 2 µ ε ( u ) − p Id , (38a) ε ( u ) := 12 (cid:0) ∇ u + ( ∇ u ) T (cid:1) . (38b)where Id is the d × d identity matrix and µ and λ are the Lamé coefficients. The weak form of thislinear elasticity problem reads: find u in H (Ω) of trace u D on Γ D and p in L (Ω) such that µ ( ε ( u ) , ε ( v )) − ( p, div( v )) = ( f , v ) + ( g , v ) Γ N ∀ v ∈ H D (Ω) , (39a) ( q, div( u )) + 1 λ ( p, q ) = 0 ∀ q ∈ L (Ω) . (39b)11he problem given by eqs. (39a) and (39b) admits a unique solution (see e.g. [47]). We introducethe finite element spaces X D ⊂ H D (Ω) and M ⊂ L (Ω) such that X D := (cid:0) V D (cid:1) d , (40)and M := V . Let w be a discretization of u D in X . Considering the stable Taylor–Hood methodof discretization, the mixed finite element approximation of eqs. (39a) and (39b) reads: find u in X D with u = w on Γ D and p in M such that µ ( ε ( u ) , ε ( v )) − ( p , div( v )) = ( f , v ) + ( g , v ) ∀ v ∈ X D , (41a) ( q , div( u )) + 1 λ ( p , q ) = 0 ∀ q ∈ M. (41b)Similarly to eqs. (39a) and (39b) transposed to the discrete context, eqs. (41a) and (41b) have aunique solution. If we denote e := u − u and ε := p − p the discretization error is measured by µ k∇ e T k + k r T k .For a cell T and an edge E the residuals are defined by R T := ( f + div (2 µ ε ( u )) − ∇ p ) | T , (42a) r T := (cid:0) div( u ) + 1 λ p (cid:1) | T , (42b) R E = [[( p Id − µ ε ( u )) n ]] if E ∈ E I , if E ∈ E D , g − ( p Id − µε ( u )) n if E ∈ E N , (42c)Here, once again we derive the a posteriori error estimator from these residuals and a local Poissonproblem, following [47]. Let T be a cell of the mesh, the local Poisson problem read: find e T in V bw T such that µ ( ∇ e T , ∇ v T ) T = ( R T , v T ) T − X E ∈ ∂T ( R E , v T ) E ∀ v T ∈ V bw T . (43)The Poisson estimator is then defined by η := X T ∈T η ,T , (44a) η ,T := 2 µ k∇ e T k T + k r T k T . (44b)This estimator has been proved to be reliable and locally efficient in [47] as well as robust inthe incompressible limit. We illustrate our implementation first on several two dimensional problems as Poisson problemswith solutions of different regularities and with different boundary conditions. Then, we also lookat examples of linear elasticity, and goal-oriented problems. Finally, we treat a three dimensionalexample: a linear elasticity problem on a mesh inspired by a human femur bone. One can findanother example of three dimensional application in [22].We apply different adaptive refinement methods as presented in section 5.1. For each methodwe perform the estimation step with a different estimator among the following: η res the residualestimator, defined in section A, η zz the Zienkiewicz–Zhu estimator, defined in section B. Notethat the Zienkiewicz–Zhu estimator is not defined for quadratic or cubic finite elements nor forlinear elasticity problems, and consequently will be absent from the comparison in these cases (Itis possible to extend the idea of the Zienkiewicz–Zhu estimator to higher-order polynomials via thedefinition of the Scott-Zhang interpolator, see [27, 65]). In addition we compare several versionsof the Bank–Weiser estimator: the bubble Bank–Weiser estimator η bbw defined from the enriched12ubble functions space and η k + ,k − bw for multiple choices of the fine and coarse spaces orders k + and k − .For each one of the following test cases we will first give a comparison of all the refinementstrategies by giving the efficiency of the a posteriori error estimator on the last mesh of thehierarchy, where the efficiency of an estimator η is defined as follows: eff := ηη err , (45)where η err is a higher order approximation of the exact error computed either from the knowledgeof the analytical solution or from a higher-order finite element method on a fine mesh. We consider a 2D L-shaped domain
Ω = ( − , \ [ − , . We solve eq. (1) with f = 0 , Γ D = Γ , u D given by the analytical solution defined below and Γ N = ∅ . In polar coordinates, the exactsolution is given by u exact ( r, θ ) = r / sin (cid:0) / θ + π/ (cid:1) . The exact solution belongs to H / − ε (Ω) for any ε > and its gradient admits a singularity at the vertex of the reentrant corner [39, Chapter5]. L-shaped domains are widely used to test adaptive mesh refinement procedures [53]. In bothlinear and quadratic finite elements all the estimators reach an expected convergence rate ( ≈ − . in the number of degrees of freedom for linear elements and ≈ − for quadratic elements). Thechoice of a posteriori error estimator is not critical for mesh refinement purposes, every estimatorleading to a hierarchy of meshes on which the corresponding errors η err are similar. For brevity wehave not included the convergence plots of these results. Linear elements.
On fig. 2 we can see the initial mesh (top left) used to start the adaptiverefinement strategies. Then, we can see the different refined meshes we obtain after seven refine-ment iterationsAs we can see on fig. 3 the Zienkiewicz–Zhu estimator η zz seems to perform the best in termsof efficiency while the second best estimator is η , . The bubble Bank–Weiser estimator η bbw isoutperformed by almost all the other Bank–Weiser estimators. The residual estimator η res largelyoverestimates the error while the estimators η k + ,k − bw for k − > largely underestimates it, leading topoor error approximations. Among the poor estimators, η , is surprisingly off for linear elementson this test case. This behavior seems to be specific to the L-shaped test cases with linear finiteelements as we will see below. Quadratic elements.
As shown on fig. 4, the best estimator in terms of efficiency is η , which nearly perfectly matches the error η err . We can also notice the very good efficiencies of η , and η , . Once again the Bank–Weiser estimators with k − > drastically underestimate the error.We can notice that the residual estimator is less efficient as the finite element degree increases. We solve eq. (1) on the same two-dimensional L-shaped boundary domain as in section 7.1 butwith different boundary conditions. We consider f = 0 , Γ N = { ( x, y ) ∈ R , x < , y = 0 } and Γ D = Γ \ Γ N . The boundary data are given by g = 0 and u D = u exact = r / sin (cid:0) / θ + π/ (cid:1) .The exact solution belongs to H / − ε (Ω) for any ε > and its gradient has a singularity locatedat the reentrant corner of Γ (see [39, Chapter 5]). As before, each estimator is leading to a a con-vergence rate close to the expected one ( ≈ − . for linear elements, ≈ − for quadratic elements)13igure 2: L-shaped Poisson problem with linear elements: On top left the initial mesh used tostart all the adaptive strategies. From top middle to bottom right, the adaptive meshes obtainedafter seven iterations of refinement strategies steered respectively by η res , η bbw , η zz , η , and η , k − k + ∅ ∅ ∅ ∅ ∅ ∅ η res η bbw η zz η k + ,k − bw and other estima-tors on the last mesh of an adaptively refined hierarchy. k − k + ∅ ∅ ∅ ∅ ∅ ∅ η res η bbw η zz ∅ Figure 4: L-shaped Poisson problem with quadratic elements: efficiencies of η k + ,k − bw and otherestimators on the last mesh of an adaptively refined hierarchy.14 − k + ∅ ∅ ∅ ∅ ∅ ∅ η res η bbw η zz η k + ,k − bw and other estimators on the last mesh of an adaptively refined hierarchy. k − k + ∅ ∅ ∅ ∅ ∅ ∅ η res η bbw η zz ∅ Figure 6: Mixed boundary conditions L-shaped Poisson problem with quadratic elements: efficien-cies of η k + ,k − bw and other estimators on the last mesh of an adaptively refined hierarchy.and the choice of the estimator does not impact the quality of the mesh hierarchy. Linear elements.
First thing we can notice from fig. 5 is that the estimators efficiencies arequite different from those in fig. 3. Most of the Bank–Weiser estimator efficiencies have improved,except when k − > . The Zienkiewicz–Zhu estimator η zz is no longer the most efficient and hasbeen outperformed by η , , η , and η , . The Bank–Weiser estimator η , still performs poorly asin fig. 3, while the residual estimator η res once again largely overestimates the error. Quadratic elements.
As for linear elements, the efficiencies in fig. 6 are very different fromfig. 4, many Bank–Weiser estimators are now underestimating the error. The most efficient esti-mator is η , closely followed by the bubble Bank–Weiser estimator η bbw . As for the previous testcases, the Bank–Weiser estimators with k − > are largely underestimating the error. We solve eq. (1) on a two-dimensional unit square domain
Ω = (0 , with u = u exact on Γ D = Γ ,( Γ N = ∅ ) and f chosen in order to have u ( x, y ) = u exact ( x, y ) = x α , with α > . . In the followingresults we chose α = 0 . . The gradient of the exact solution u admits a singularity along the leftboundary of Ω (for x = 0 ). The solution u belongs to H / − ε for all ε > [45, 52]. Consequently,the value of α determines the strength of the singularity and the regularity of u .Unlike the previous test cases, all the estimators are achieving a convergence rate close to − . instead of − . for linear elements. Moreover, this rate does not improve for higher-order elements(for brevity, the results for higher-order elements are not shown here). The low convergence rateshows how computationally challenging such a problem can be. Once again the choice of estimatoris not critical for mesh refinement purposes. 15 − k + ∅ ∅ ∅ ∅ ∅ ∅ η res η bbw η zz η k + ,k − bw andother estimators on the last mesh of an adaptively refined hierarchy. k − k + ∅ ∅ ∅ ∅ ∅ ∅ η res η bbw η zz ∅ Figure 8: Boundary singularity Poisson problem with linear elements: efficiencies of η k + ,k − bw andother estimators on the last mesh of an adaptively refined hierarchy. Linear elements.
The best estimator in terms of efficiency is η , which slightly overestimatesthe error, closely followed by η , underestimating the error as we can see on fig. 7. Unlike theprevious test case, here the Zienkiewicz–Zhu estimator η zz grandly underestimates the error. Theworst estimator is the residual estimator η res which gives no precise information about the error.We can notice that the poor performance of the estimator η , on the L-shaped test case does notreproduce here. Quadratic elements.
Again, fig. 6 shows that the best estimator is η , closely followed by η , and the bubble estimator η bbw . The residual estimator is getting worse as the finite elementdegree increases. We solve the L-shaped domain problem as described in section 7.1 but instead of controlling errorin the natural norm, we aim to control error in the goal functional J ( u ) = ( c, u ) with c a smoothbump function c (¯ r ) := (cid:26) ε − exp (cid:0) − r (cid:1) ≤ ¯ r < , r ≥ . (46)where ¯ r = (( x − ¯ x ) /ε ) + (( y − ¯ y ) /ε ) , with ε ∈ R a parameter that controls the size of the bumpfunction, and ¯ x ∈ R and ¯ y ∈ R the position of the bumps function’s center. We set ε = 0 . and ¯ x = ¯ y = 0 . . With these parameters the goal functional is isolated to a region close to there-entrant corner.We use the goal-oriented adaptive mesh refinement methodology outlined in section 5.2. We16igure 9: L-shaped goal-oriented Poisson problem with linear elements: On top left the initial meshused to start all the adaptive strategies. From top middle to bottom right, the adaptive meshesobtained after seven iterations of refinement strategies steered by weighted estimators derivedrespectively from η , , η res , η bbw and η zz for both primal and dual problems.use a first-order polynomial finite element method for the primal and dual problem, and the Bank–Weiser error estimation procedure to calculate both η u and η z .The ‘exact’ value of the functional J ( u ) was calculated on a very fine mesh using a fourth-orderpolynomial finite element space and was used to compute higher-order approximate errors for eachrefinement strategy.The weighted goal-oriented strategy refines both the re-entrant corner and the broader regionof interest defined by the goal functional. Relatively less refinement occurs in the regions far awayfrom either of these important areas.The figure fig. 9 shows refined meshes after seven iterations of the weighted goal orientedmethod. We can see that the meshes are mainly refined in the re-entrant corner as well as in theregion on the right top of it where the goal functional focuses. In fig. 10 we show the convergencecurves of some of these adaptive strategies. For each strategy, η u = η z is the estimator specifiedin the legend and η w = η u η z . All the strategies we have tried led to very similar higher-orderapproximate errors. So for the sake of clarity we have replaced the approximate errors by anindicative line computed using a regression from the least squares method (lstsq error), leadingto the line that fits the best the values of the different approximate errors. As we can see, theseadaptive strategies are reaching an optimal convergence rate. Although it is also the case for allthe other strategies we have tried, we do not show the other results for the sake of concision. Inthe left table of fig. 11 we show the efficiencies of the estimators η w where η u = η z = η k + ,k − bw . Onthe right table of fig. 11 we take η u = η z to be the estimators in the left column. As we can seeon the efficiencies are not as good as in section 7.1. The two best estimators are those derivedfrom η , and η , . With η , , they are the only cases where the goal-oriented estimator η w is17 Number of dof10 bbw (primal) bbw (dual) bbw (wgo) zz (primal) zz (dual) zz (wgo) res (primal) res (dual) res (wgo) (primal) (dual) (wgo) lstsq error Figure 10: L-shaped goal-oriented Poisson problem with linear elements: plot comparing conver-gence of some goal-oriented adaptive strategies driven by four different estimators. Expected ratesfor primal and dual problems ( − . ) and goal functional ( − ) shown by triangle markers. Com-parison with an indicative line representing the higher order approximation of the errors of eachstrategy and obtained using least squares method.performing better than the goal-oriented estimator derived from the Zienkiewicz–Zhu estimator.The estimators η w derived from the bubble Bank–Weiser estimator as well as from the residualestimator are poorly overestimating the error. We consider the linear elasticity problem from [28] on the centered unit square domain Ω withhomogeneous Dirichlet boundary conditions on Γ D = Γ ( u D = 0 ). The first Lamé coefficient isset to µ = 100 and the Poisson ratio to ν = 0 . and ν = 0 . . The problem data f is given by f = ( f , f ) with f ( x, y ) = − µπ cos( πy ) sin( πy ) (cid:0) πx ) − (cid:1) ,f ( x, y ) = 2 µπ cos( πx ) sin( πx ) (cid:0) πy ) − (cid:1) . (47)The corresponding exact solution of the linear elasticity problem reads u = ( u , u ) with u ( x, y ) = π cos( πy ) sin ( πx ) sin( πy ) , u ( x, y ) = − π cos( πx ) sin( πx ) sin ( πy ) , (48)the Herrmann pressure is zero everywhere on Ω . In each case we discretize this problem using theTaylor–Hood element and an initial cartesian mesh and we apply our adaptive procedure drivenby the Poisson estimator described in section 5.3.1. We compare the Poisson estimators derivedfrom different Bank–Weiser estimators and the residual estimator.As before, all the refinement strategies are achieving an optimal convergence rate no matter thevalue of ν . fig. 12 shows the results for ν = 0 . . We notice that almost all the Poisson estimatorsderived from Bank–Weiser estimators have a very good efficiency. The best estimator in this case is η , closely followed by η , , η bbw and η , . Although the residual estimator still performs the worst,18 − k + ∅ ∅ ∅ ∅ ∅ ∅ η res η bbw η zz η k + ,k − bw and other estimators on the last mesh of an adaptivelyrefined hierarchy. k − k + ∅ ∅ ∅ ∅ ∅ ∅ η res η bbw η zz ∅ Figure 12: Nearly-incompressible elasticity ( ν = 0 . ) problem with Taylor–Hood elements: effi-ciencies of the Poisson estimators derived from η k + ,k − bw and other estimators on the last mesh of anadaptively refined hierarchy.it is sharper than in all the previous test cases. As we can notice on fig. 13, all the estimators arerobust with respect to the incompressibility constraint. All the efficiencies have slightly increasedand some estimators ( η , and η , ) that where a lower bound of the error previously are now anupper bound. In this test case we consider a linear elasticity problem on a domain inspired by a human femurbone .The goal of this test case is not to provide an accurate description of the behavior of the femurbone but to demonstrate the applicability of our implementation to 3D dimensional goal-orientedproblem with large number of degrees of freedom: the linear elasticity problem to solve on theinitial mesh, using Taylor-Hood element has , degrees of freedom while our last refinementstep reaches , , degrees of freedom.The 3D mesh for analysis is build from the surface model using the C++ library CGAL [4]via the Python front-end pygalmesh. The parameters of the linear elasticity problem for this testcases are the following: the Young’s modulus is set to GPa and the Poisson’s ratio to . (seee.g. [64]). In addition, the load is given by f = (0 , , , the Dirichlet data by u D = 0 on Γ D ( Γ represented as the left dark gray region of the boundary in fig. 14 and g the traction data is definedas g = (0 , , on the center light gray region of the boundary and is constant on the right darkgray region of the boundary g = ( − − , − − , − ) . The femur-shaped domain Ω as well as The STL model of the femur bone can be found at https://3dprint.nih.gov/discover/3dpx-000168 under aPublic Domain license. − k + ∅ ∅ ∅ ∅ ∅ ∅ η res η bbw η zz ∅ Figure 13: Nearly-incompressible elasticity ( ν = 0 . ) problem with Taylor–Hood elements: effi-ciencies of the Poisson estimators derived from η k + ,k − bw and other estimators on the last mesh of anadaptively refined hierarchy.the initial and last meshes are shown in fig. 14. As we can see, the refinement occurs mainly in thecentral region of the femur, where the goal functional J focus. Some artefacts can be seen as stainsof refinement in the central region due to the fact that we use the initial mesh as our geometryand on the left due to the discontinuity in the boundary conditions.In fig. 15 the primal solution is given by the couple ( u , p ) and the dual solution by ( z , κ ) .As we can notice and as expected, the weighted estimator η w converges twice as fast as the primaland dual estimators. In this paper we have shown how the error estimator of Bank–Weiser, involving the solutionof a local problem on a special finite element space, can be mathematically reformulated andimplemented straightforwardly in a modern finite element software with the aid of automatic codegeneration techniques. Through a series of numerical results we have shown that the Bank–Weiseris highly competitive in accurately predicting the total global error and in driving an adaptivemesh refinement strategy. Furthermore, the basic methodology and implementation for the Poissonproblem can be extended to tackle more complex mixed discretizations of PDEs including nearly-incompressible elasticity or Stokes problems.
Acknowledgements
The three-dimensional problem presented in this paper was carried out using the HPC facilities ofthe University of Luxembourg [67].We would like to thank Nathan Sime and Chris Richardson for helpful discussions on FEniCS’sdiscontinuous Galerkin discretizations and PLAZA refinement strategy, respectively. We wouldalso like to thank Prof. Roland Becker for the motivating and fruitful discussions.
Supplementary material
The supplementary material (SM) document contains two snippets showing the implementation ofthe Poisson estimator and the Poisson estimator for the nearly-incompressible elasticity problem.A simplified version of the code (LGPLv3) used to produce the results in this paper is archived at https://doi.org/10.6084/m9.figshare.10732421 . A Docker image [42] is provided in whichthis code can be executed. A port of our approach from the DOLFIN finite element solver to thenew DOLFINX finite element solver https://github.com/fenics/dolfinx is underway at thesame location. 20igure 14: Femur bone linear elasticity problem with Taylor–Hood elements: on the top, the threedifferent regions of the boundary corresponding to different boundary conditions: the left darkgrey region is the non-zero Neumann boundary, the middle light grey region is the zero Neumannboundary and the right dark grey region is the Dirichlet boundary. In the middle, the initial mesh.On the bottom, the last mesh after several steps of adaptive refinement. Number of dof10 S c a l e d e s t i m a t o r s u /(2 || u || + || p ||) (primal) z /(2 || z || + || ||) (dual) w /| J ( u , p )| (wgo) Figure 15: Femur bone linear elasticity problem with Taylor–Hood elements: convergence curvesof the primal, dual and weighted estimators respectively scaled by the norm of the primal solution,dual solution and magnitude of the goal functional evaluated in the primal solution.21 ppendicesA The residual estimators
A.1 Poisson equation
The class of residual estimators, the explicit residual estimator is part of, have been introduced forthe first time in [10]. Let h T be the diameter (see e.g. [65]) of the cell T and h E be the diameterof the face E in dimension three and the length of E in dimension two. The explicit residualestimator [3] on a cell T for the Poisson problems eqs. (2) and (4) is defined as η ,T := h T k f T + ∆ u k k T + X E ∈E I ∩ ∂T h E k [[ ∂ n u k ]] E k E + X E ∈E N ∩ ∂T h E k g E − ∂ n u k k E , (1)where f T := | T | R T f is the average of f on the cell T and g E := | E | R E g is the average of g on theedge E . The global residual estimator reads η := X T ∈T η ,T . (2) A.2 Linear elasticity equations
The residual estimator for the linear elasticity problem eqs. (39a), (39b), (41a) and (41b) is givenby η ,T := ρ T k R T k T + ρ d k r T k T + X E ∈ ∂T ρ E k R E k E , (3)where the residuals R T , r T and R E are respectively defined in eqs. (42a) to (42c) and the constants ρ T , ρ d and ρ E are given by ρ T := h T (2 µ ) − / , ρ d := (cid:0) λ − + (2 µ ) − (cid:1) − , ρ E := h E (2 µ ) − , (4)with h T the diameter of the cell T and h E the length of the edge E . The global estimator reads η := X T ∈T η ,T . (5) B The Zienkiewicz–Zhu estimator
The Zienkiewicz–Zhu estimator is a gradient recovery estimator based on an averaging techniqueintroduced in [71]. This estimator belongs to a general class of recovery estimators, see [24, 25, 70]for recent surveys and a reformulation of the recovery procedure in an H (div) -conforming spacethat has superior performance for problems with sharp interfaces. Despite the fact that somerecovery estimators are now available for higher order finite elements (see for example [72]) weonly consider the original estimator, defined for a piecewise linear finite element framework.Given the finite element solution u ∈ V the numerical flux ρ := ∇ u is a piecewise constantvector field. For each vertex χ ∈ N in the mesh we denote ω χ the domain covered by the unionof cells T having common vertex χ . The recovered flux G ( ρ ) ∈ [ V ] has values at the degrees offreedom associated with the vertices N given by G ( ρ )( χ ) := 1 | ω χ | Z ω χ ρ d x ∀ χ ∈ N . (6)22he local Zienkiewicz–Zhu estimator is then defined as the discrepancy between the recovered fluxand the numerical flux η zz ,T := k G ( ρ ) − ρ k T ∀ T ∈ T . (7)The global Zienkiewicz–Zhu estimator is given by η := X T ∈T η ,T . (8)The code in the supplementary material contains a prototype implementation of the Zienkiewicz–Zhu estimator in FEniCS. We have implemented the local recovered flux calculation in Pythonrather than C ++ , so the runtime performance is far from optimal. References [1]
M. Ainsworth , The influence and selection of subspaces for a posteriori error estimators , Nu-merische Mathematik, 73 (1996), pp. 399–418, https://doi.org/10.1007/s002110050198 .[2]
M. Ainsworth and J. T. Oden , A Posteriori Error Estimators for the Stokes and OseenEquations , SIAM Journal on Numerical Analysis, 34 (1997), pp. 228–245, https://doi.org/10.1137/S0036142994264092 .[3]
M. Ainsworth and J. T. Oden , A Posteriori Error Estimation in Finite Element Analysis ,2011.[4]
P. Alliez, C. Jamin, L. Rineau, S. Tayeb, J. Tournois, and M. Yvinec ,
3D meshgeneration , in CGAL User and Reference Manual, CGAL Editorial Board, 5.1 ed., 2020, https://doc.cgal.org/5.1/Manual/packages.html .[5]
M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson,J. Ring, M. E. Rognes, and G. N. Wells , The FEniCS Project Version 1.5 , Archive ofNumerical Software, Vol 3 (2015), https://doi.org/10.11588/ANS.2015.100.20553 .[6]
M. S. Alnæs, A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells , Unified formlanguage: A domain-specific language for weak formulations of partial differential equations ,ACM Transactions on Mathematical Software, 40 (2014), pp. 1–37, https://doi.org/10.1145/2566630 .[7]
P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent, and J. Koster , A Fully AsynchronousMultifrontal Solver Using Distributed Dynamic Scheduling , SIAM Journal on Matrix Analysisand Applications, 23 (2001), pp. 15–41, https://doi.org/10.1137/S0895479899358194 .[8]
P. R. Amestoy, A. Guermouche, J.-Y. L’Excellent, and S. Pralet , Hybrid schedul-ing for the parallel solution of linear systems , Parallel Computing, 32 (2006), pp. 136–156, https://doi.org/10.1016/j.parco.2005.07.004 .[9]
I. Babuska and M. Vogelius , Feedback and adaptive finite element solution of one-dimensional boundary value problems , Numerische Mathematik, 44 (1984), pp. 75–102, https://doi.org/10.1007/BF01389757 .[10]
I. Babuška and W. C. Rheinboldt , A-posteriori error estimates for the finite elementmethod , International Journal for Numerical Methods in Engineering, 12 (1978), pp. 1597–1615, https://doi.org/10.1002/nme.1620121010 .2311]
S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman,L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C.McInnes, K. Rupp, B. F. Smith, S. Zampini, H. Zhang, and H. Zhang , PETScUsers Manual , Tech. Report ANL-95/11 - Revision 3.7, Argonne National Laboratory, 2016, .[12]
R. E. Bank , PLTMG: A Software Package for Solving Elliptic Partial Differential Equations:Users’ Guide 8.0 , Society for Industrial and Applied Mathematics, Jan. 1998, https://doi.org/10.1137/1.9780898719635 .[13]
R. E. Bank and A. Weiser , Some a posteriori error estimators for elliptic partial differen-tial equations , Mathematics of Computation, 44 (1985), pp. 283–283, https://doi.org/10.1090/S0025-5718-1985-0777265-X .[14]
R. E. Bank, J. Xu, and B. Zheng , Superconvergent Derivative Recovery for LagrangeTriangular Elements of Degree p on Unstructured Grids , SIAM J. Numer. Anal., 45 (2007),pp. 2032–2046, https://doi.org/10.1137/060675174 .[15]
S. Bartels and C. Carstensen , Each averaging technique yields reliable a posteriori errorcontrol in FEM on unstructured grids. Part II: Higher order FEM , Mathematics of Compu-tation, 71 (2002), pp. 971–994, https://doi.org/10.1090/S0025-5718-02-01412-6 .[16]
S. Bartels, C. Carstensen, and G. Dolzmann , Inhomogeneous Dirichlet conditions ina priori and a posteriori finite element error analysis , Numerische Mathematik, 99 (2004),pp. 1–24, https://doi.org/10.1007/s00211-004-0548-3 .[17]
R. Becker, E. Estecahandy, and D. Trujillo , Weighted Marking for Goal-orientedAdaptive Finite Element Methods , SIAM Journal on Numerical Analysis, 49 (2011), pp. 2451–2469, https://doi.org/10.1137/100794298 .[18]
R. Becker and R. Rannacher , An optimal control approach to a posteriori error estimationin finite element methods , Acta Numer., 10 (2001), pp. 1–102, https://doi.org/10.1017/s0962492901000010 .[19]
L. Beirão da Veiga, C. Chinosi, C. Lovadina, and R. Stenberg , A-priori and a-posteriori error analysis for a family of Reissner–Mindlin plate elements , BIT NumericalMathematics, 48 (2008), pp. 189–213, https://doi.org/10.1007/s10543-008-0175-y .[20]
A. Bespalov, D. Praetorius, L. Rocchi, and M. Ruggeri , Goal-oriented error estima-tion and adaptivity for elliptic PDEs with parametric or uncertain inputs , Computer Methodsin Applied Mechanics and Engineering, 345 (2019), pp. 951–982, https://doi.org/10.1016/j.cma.2018.10.041 .[21]
A. Bespalov, L. Rocchi, and D. Silvester , T-IFISS: a toolbox for adaptive FEM com-putation , Computers & Mathematics with Applications, (2020), https://doi.org/10.1016/j.camwa.2020.03.005 .[22]
R. Bulle, F. Chouly, J. S. Hale, and A. Lozinski , Removing the saturation assumptionin Bank–Weiser error estimator analysis in dimension three , Applied Mathematics Letters,107 (2020), p. 106429, https://doi.org/10.1016/j.aml.2020.106429 .[23]
R. Bulle and J. S. Hale , An implementation of the Bank-Weiser error estimator in theFEniCS Project finite element software , (2020), https://doi.org/10.6084/m9.figshare.10732421.v1 . 2424]
Z. Cai, C. He, and S. Zhang , Improved ZZ a posteriori error estimators for diffusion prob-lems: Conforming linear elements , Computer Methods in Applied Mechanics and Engineering,313 (2017), pp. 433–449, https://doi.org/10.1016/j.cma.2016.10.006 .[25]
Z. Cai and S. Zhang , Recovery-Based Error Estimator for Interface Problems: ConformingLinear Elements , SIAM Journal on Numerical Analysis, 47 (2009), pp. 2132–2156, https://doi.org/10.1137/080717407 .[26]
C. Carstensen and S. Bartels , Each averaging technique yields reliable a posteriori errorcontrol in FEM on unstructured grids. Part I: Low order conforming, nonconforming, andmixed FEM , Mathematics of Computation, 71 (2002), pp. 945–969, https://doi.org/10.1090/S0025-5718-02-01402-3 .[27]
C. Carstensen, M. Feischl, M. Page, and D. Praetorius , Axioms of adaptivity , Com-puters & Mathematics with Applications, 67 (2014), pp. 1195–1253, https://doi.org/10.1016/j.camwa.2013.12.003 .[28]
C. Carstensen and J. Gedicke , Robust residual-based a posteriori Arnold–Winther mixedfinite element analysis in elasticity , Computer Methods in Applied Mechanics and Engineer-ing, 300 (2016), pp. 245–264, https://doi.org/10.1016/j.cma.2015.10.001 .[29]
C. Carstensen and C. Merdon , Estimator Competition for Poisson Problems , Journal ofComputational Mathematics, 28 (2010), https://doi.org/10.4208/jcm.2009.10-m1010 .[30]
Concha INRIA Project-Team , Complex Flow Simulation Codes based on High-order andAdaptive methods , https://raweb.inria.fr/rapportsactivite/RA2010/concha/concha.pdf (accessed 17/12/2020).[31] M. Duprez, S. P. A. Bordas, M. Bucki, H. P. Bui, F. Chouly, V. Lleras, C. Lo-bos, A. Lozinski, P.-Y. Rohan, and S. Tomar , Quantifying discretization errors for softtissue simulation in computer assisted surgery: A preliminary study , Applied MathematicalModelling, 77 (2020), pp. 709–723, https://doi.org/10.1016/j.apm.2019.07.055 .[32]
W. Dörfler , A Convergent Adaptive Algorithm for Poisson’s Equation , SIAM Journal onNumerical Analysis, 33 (1996), pp. 1106–1124, https://doi.org/10.1137/0733054 .[33]
W. Dörfler and R. H. Nochetto , Small data oscillation implies the saturationassumption , Numerische Mathematik, 91 (2002), pp. 1–12, https://doi.org/10.1007/s002110100321 .[34]
H. C. Elman, A. Ramage, and D. J. Silvester , IFISS: A Computational Laboratory forInvestigating Incompressible Flow Problems , SIAM Review, 56 (2014), pp. 261–273, https://doi.org/10.1137/120891393 .[35]
R. D. Falgout and U. M. Yang , hypre: A Library of High Performance Preconditioners ,in Comput. Science — ICCS 2002, P. M. A. Sloot, A. G. Hoekstra, C. J. K. Tan, and J. J.Dongarra, eds., no. 2331 in Lecture Notes in Computer Science, Springer Berlin Heidelberg,apr 2002, pp. 632–641, https://doi.org/10.1007/3-540-47789-6_66 .[36] S. Funken, D. Praetorius, and P. Wissgott , Efficient implementation of adaptive P1-FEM in Matlab , Computational Methods in Applied Mathematics, 11 (2011), pp. 460–490, https://doi.org/10.2478/cmam-2011-0026 .[37]
T. H. Gibson, L. Mitchell, D. A. Ham, and C. J. Cotter , Slate: extending Firedrake’sdomain-specific abstraction to hybridized solvers for geoscience and beyond , Geosci. Model Dev.Discuss., (2019), pp. 1–40, https://doi.org/https://doi.org/10.5194/gmd-2019-86 .2538]
M. B. Giles and E. Süli , Adjoint methods for PDEs: a posteriori error analysis andpostprocessing by duality , Acta Numerica, 11 (2002), pp. 145–236, https://doi.org/10.1017/S096249290200003X .[39]
P. Grisvard , Elliptic problems in nonsmooth domains , vol. 24 of Monographs and Studiesin Mathematics, Pitman (Advanced Publishing Program), Boston, MA, 1985.[40]
G. Guennebaud, B. Jacob, and Others , Eigen v3 , 2010, http://eigen.tuxfamily.org (accessed 2020-11-06).[41]
J. S. Hale, M. Brunetti, S. P. Bordas, and C. Maurini , Simple and extensible plate andshell finite element models through automatic code generation tools , Computers & Structures,209 (2018), pp. 163–181, https://doi.org/10.1016/j.compstruc.2018.08.001 .[42]
J. S. Hale, L. Li, C. N. Richardson, and G. N. Wells , Containers for Portable,Productive, and Performant Scientific Computing , Computing in Science & Engineering, 19(2017), pp. 40–50, https://doi.org/10.1109/MCSE.2017.2421459 .[43]
F. Hecht , New development in freefem++ , J. Numer. Math., 20 (2012), pp. 1–14, https://doi.org/10.1515/jnum-2012-0013 .[44]
C. A. R. Hoare , Algorithm 65: find , Communications of the ACM, 4 (1961), pp. 321–322, https://doi.org/10.1145/366622.366647 .[45]
P. Houston, B. Senior, and E. Süli , Sobolev regularity estimation for hp-adaptive finiteelement methods , in Numer. Math. Adv. Appl., Springer Milan, Milano, 2003, pp. 631–656, https://doi.org/10.1007/978-88-470-2089-4_58 .[46]
P. Houston and N. Sime , Automatic Symbolic Computation for Discontinuous GalerkinFinite Element Methods , SIAM Journal on Scientific Computing, 40 (2018), pp. C327–C357, https://doi.org/10.1137/17M1129751 .[47]
A. Khan, C. E. Powell, and D. J. Silvester , Robust a posteriori error estimators formixed approximation of nearly incompressible elasticity , International Journal for NumericalMethods in Engineering, 119 (2019), pp. 18–37, https://doi.org/10.1002/nme.6040 .[48]
R. C. Kirby , Algorithm 839: FIAT, a new paradigm for computing finite element basisfunctions , ACM Transactions on Mathematical Software, 30 (2004), pp. 502–516, https://doi.org/10.1145/1039813.1039820 .[49]
R. C. Kirby and A. Logg , A compiler for variational forms , ACM Transactions on Math-ematical Software, 32 (2006), pp. 417–444, https://doi.org/10.1145/1163641.1163644 .[50]
Q. Liao and D. Silvester , A simple yet effective a posteriori estimator for classical mixedapproximation of Stokes equations , Applied Numerical Mathematics, 62 (2012), pp. 1242–1256, https://doi.org/10.1016/j.apnum.2010.05.003 .[51]
A. Logg and G. N. Wells , DOLFIN: Automated finite element computing , ACM Transac-tions on Mathematical Software, 37 (2010), pp. 1–28, https://doi.org/10.1145/1731022.1731030 .[52]
W. F. Mitchell , A collection of 2D elliptic problems for testing adaptive grid refinementalgorithms , Appl. Math. Comput., 220 (2013), pp. 350–364, https://doi.org/10.1016/j.amc.2013.05.068 . 2653]
W. F. Mitchell , A collection of 2D elliptic problems for testing adaptive grid refinementalgorithms , Applied Mathematics and Computation, 220 (2013), pp. 350–364, https://doi.org/10.1016/j.amc.2013.05.068 .[54]
M. S. Mommer and R. Stevenson , A Goal-Oriented Adaptive Finite Element Method withConvergence Rates , SIAM Journal on Numerical Analysis, 47 (2009), pp. 861–886, https://doi.org/10.1137/060675666 .[55]
R. H. Nochetto , Removing the saturation assumption in a posteriori error analysis , Istit.Lomb. Accad. Sci. Lett. Rend. A, 127 (1993), pp. 67–82 (1994).[56]
R. H. Nochetto, K. G. Siebert, and A. Veeser , Theory of adaptive finite ele-ment methods: An introduction , in Multiscale, Nonlinear Adapt. Approx., R. DeVore andA. Kunoth, eds., Berlin, Heidelberg, 2009, Springer Berlin Heidelberg, pp. 409–542, https://doi.org/10.1007/978-3-642-03413-8_12 .[57]
C.-M. Pfeiler and D. Praetorius , Dörfler marking with minimal cardinality is a lin-ear complexity problem , arXiv1907.13078 [cs, math], (2019), http://arxiv.org/abs/1907.13078 .[58]
A. Plaza and G. Carey , Local refinement of simplicial grids based on the skele-ton , Applied Numerical Mathematics, 32 (2000), pp. 195–218, https://doi.org/10.1016/S0168-9274(99)00022-7 .[59]
C. Prud’homme, V. Chabannes, V. Doyeux, M. Ismail, A. Samake, and G. Pena , Feel++ : A computational framework for Galerkin Methods and Advanced Numerical Methods ,ESAIM: Proceedings, 38 (2012), pp. 429–455, https://doi.org/10.1051/proc/201238024 .[60]
F. Rathgeber, D. A. Ham, L. Mitchell, M. Lange, F. Luporini, A. T. T. Mcrae,G.-T. Bercea, G. R. Markall, and P. H. J. Kelly , Firedrake: Automating the FiniteElement Method by Composing Abstractions , ACM Transactions on Mathematical Software,43 (2017), pp. 1–27, https://doi.org/10.1145/2998441 .[61]
Y. Renard and K. Poulios , GetFEM: Automated FE modeling of multiphysicsproblems based on a generic weak form language , (Submitted), (2020), https://hal.archives-ouvertes.fr/hal-02532422/ .[62]
R. Rodríguez , Some remarks on Zienkiewicz-Zhu estimator , Numerical Methods for PartialDifferential Equations, 10 (1994), pp. 625–635, https://doi.org/10.1002/num.1690100509 .[63]
M. E. Rognes and A. Logg , Automated Goal-Oriented Error Control I: StationaryVariational Problems , SIAM Journal on Scientific Computing, 35 (2013), pp. C173–C193, https://doi.org/10.1137/10081962X .[64]
F. Rupin, A. Saied, D. Dalmas, F. Peyrin, S. Haupert, E. Barthel, G. Boivin, andP. Laugier , Experimental determination of Young modulus and Poisson ratio in cortical bonetissue using high resolution scanning acoustic microscopy and nanoindentation , J. Acoust. Soc.Am., 123 (2008), pp. 3785–3785, https://doi.org/10.1121/1.2935440 .[65]
L. R. Scott and S. Zhang , Finite element interpolation of nonsmooth functions satisfyingboundary conditions , Mathematics of Computation, 54 (1990), pp. 483–483, https://doi.org/10.1090/S0025-5718-1990-1011446-7 .[66]
S. van der Walt, S. C. Colbert, and G. Varoquaux , The NumPy Array: A Structurefor Efficient Numerical Computation , Computing in Science & Engineering, 13 (2011), pp. 22–30, https://doi.org/10.1109/MCSE.2011.37 .2767]
S. Varrette, P. Bouvry, H. Cartiaux, and F. Georgatos , Management of an academicHPC cluster: The UL experience , in 2014 International Conference on High PerformanceComputing & Simulation (HPCS), Bologna, Italy, July 2014, IEEE, pp. 959–967, https://doi.org/10.1109/HPCSim.2014.6903792 .[68]
R. Verfürth , A posteriori error estimation and adaptive mesh-refinement techniques , Jour-nal of Computational and Applied Mathematics, 50 (1994), pp. 67–83, https://doi.org/10.1016/0377-0427(94)90290-9 .[69]
P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Courna-peau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt,M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones,R. Kern, E. Larson, C. J. Carey, I. Polat, Y. Feng, E. W. Moore, J. Vander-Plas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R.Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, andContributors , SciPy 1.0–Fundamental Algorithms for Scientific Computing in Python , Na-ture Methods, 17 (2020), pp. 261–272, https://doi.org/10.1038/s41592-019-0686-2 .[70]
Z. Zhang and N. Yan , Recovery type a posteriori error estimates in finite element methods ,J. Appl. Math. Comput., 8 (2001), pp. 235–251, https://doi.org/10.1007/BF02941963 .[71]
O. C. Zienkiewicz and J. Z. Zhu , A simple error estimator and adaptive procedure forpractical engineering analysis , International Journal for Numerical Methods in Engineering,24 (1987), pp. 337–357, https://doi.org/10.1002/nme.1620240206 .[72]
O. C. Zienkiewicz and J. Z. Zhu , The superconvergent patch recovery (SPR) and adaptivefinite element refinement , Computer Methods in Applied Mechanics and Engineering, 101(1992), pp. 207–224, https://doi.org/10.1016/0045-7825(92)90023-D .[73]
K. B. Ølgaard, A. Logg, and G. N. Wells , Automated Code Generation for Discon-tinuous Galerkin Methods , SIAM Journal on Scientific Computing, 31 (2009), pp. 849–864, https://doi.org/10.1137/070710032 . 28 upplementary materials: Hierarchical a posteriori errorestimation of Bank–Weiser type in the FEniCS Project ∗Raphaël Bulle † Jack S. Hale † Alexei Lozinski ‡ Stéphane P. A. Bordas † Franz Chouly § February 9, 2021
We present here a snippet of DOLFIN Python code showing function to compute the error of aPoisson problem using the Bank–Weiser estimator. from dolfin import *import fenics_error_estimationdef estimate(u_h):"""Bank-Weiser error estimation procedure for the Poisson problem.Parameters-----------u_h: dolfin.FunctionSolution of Poisson problem.Returns-------The error estimate on each cell of the mesh."""mesh = u_h.function_space().mesh() ∗ Raphaël Bulle would like to acknowledge the support of the ASSIST research project of the University ofLuxembourg. This publication has been prepared in the framework of the DRIVEN project funded by the EuropeanUnion’s Horizon 2020 Research and Innovation programme under Grant Agreement No. 811099. † Institute of Computational Engineering, University of Luxembourg, 6 Avenue de la Fonte, 4362 Esch-sur-Alzette,Luxembourg ( [email protected] , [email protected] , [email protected] ) ‡ Laboratoire de Mathématiques de Besançon, UMR CNRS 6623, Université de Bourgogne Franche-Comté, 16route de Gray, 25030 Besançon Cedex, France ( [email protected] ) § Université de Bourgogne Franche-Comté, Institut de Mathématiques de Bourgogne, 21078 Dijon, France( [email protected] ) a r X i v : . [ m a t h . NA ] F e b Construct the Bank-Weiser interpolation operator according to the Indicative snippet of error estimation for linear elasticityequations using Poisson estimator
We give here a snippet of DOLFIN Python code showing function to compute the error of a two-dimensional linear elasticity problem (discretized with Taylor–Hood element) using the Poissonestimator, based on our implementation of the Bank–Weiser estimator. import scipy.linalg as sp.linalgfrom dolfin import *import fenics_error_estimationdef estimate(w_h, mu, lmbda):"""Parameters-----------w_h: dolfin.FunctionSolution of the linear elasticity problem.mu: floatFirst Lamé coefficient.lmbda: floatSecond Lamé coefficient.Returns-------The error estimate on each cell of the mesh."""mesh = w_h.function_space().mesh()u_h = w_h.sub(0)p_h = w_h.sub(1) cs = DirichletBC(X_f, Constant((0., 0.)), 'on_boundary', 'geometric')cs = DirichletBC(X_f, Constant((0., 0.)), 'on_boundary', 'geometric')