Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems

An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems Anita Hansbo, Peter Hansbo To cite this version: Anita Hansbo, Peter Hansbo. An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems. Computer Methods in Applied Mechanics and Engineering, Elsevier, 2002, 191 (47-48), pp. 5537-5552. 10.1016/S0045-7825(02)00524-8. hal-01352903 HAL Id: hal-01352903 https://hal.archives-ouvertes.fr/hal-01352903 Submitted on 10 Aug 2016 HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non, lished or not. The documents may come from émanant des établissements d’enseignement et de teaching and research institutions in France or recherche français ou étrangers, des laboratoires abroad, or from public or private research centers. publics ou privés. Distributed under a Creative Commons Attribution| 4.0 International License An unfitted finite element method, based on Nitsches method, for elliptic interface problems a b Anita Hansbo , Peter Hansbo Department of Informatics and Mathematics, University of Trollh€aattan-Uddevalla, Box 957, S-461 39 Trollh€aattan, Sweden Department of Applied Mechanics, Chalmers University of Technology, S-412 96G €ooteborg, Sweden In this paper we propose a method for the finite element solution of elliptic interface problem, usingan approach due to Nitsche. The method allows for discontinuities, internal to the elements, in the approximation across the interface. We show that optimal order of convergence holds without restrictions on the location of the interface relative to the mesh. Further, we derive a posteriori error estimates for the purpose of controllingfunctionals of the error and present some numerical examples. 1. Introduction As a model elliptic interface problem, we consider a stationary heat conduction problem in two di- mensions with a conduction coefficient which is discontinuous across a smooth internal interface. When solvingsuch problems numerically usingthe standard finite element method, one usually takes the dis- continuity of the data into account by enforcingmesh lines along the interface. If this is not done, sub- optimal convergence behaviour will occur, cf. [1,8]. As a motivation for this work, we also have in mind more complicated, time dependent or non-linear, problems where the interface moves with time or duringiteration. In that case, it may be advantageous to use the same mesh on the domain for different, nearby, locations of the interface, since repeated remeshing of the domain to obtain fitted meshes is very costly. We are thus led to study unfitted finite element methods where the interface is allowed to cross the elements. In this paper, we propose an unfitted finite element method, based on a variant of Nitsches method [9], allowingfor discontinuities, internal to the elements, in the approximation across the interface. This method is of optimal order; in particular we show second order convergence in L for appropriately modified piecewise linears on a non-degenerate triangulation. We also consider a posteriori error estimates 1 5538 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 for functionals of the solution, in the spirit of Becker and Rannacher [3], and use these estimates as a basis for adaptively refiningthe mesh. Fitted mesh FE methods for elliptic problems with discontinuous coefficients and homogeneous interface conditions are analysed in [1,6,10]. In [2,4,5], problems with inhomogeneous interface conditions are considered. As for unfitted mesh methods for interface problems, Barrett and Elliott [2] show first order of con- vergence in energy-norm and interior second order L error estimates for a piecewise linear method based on boundary penalty and numerical integration over approximate domains. This approach is close in spirit to the method to be presented here, but contains a consistency error that we avoid. An alternative approach is to construct a FEM-basis where the basis functions fulfill homogeneous interface conditions exactly, and to use the continuous bilinear form (without penalty) to define the method. MacKinnon and Carey [8] use linear basis functions with this property in two dimensions and numerical examples of optimal order of convergence are presented. Li et al. [7] analyse the max-norm interpolation error on cartesian grids of a non-conforming piecewise linear method where the basis func- tions fulfill homogeneous interface conditions exactly, as well as a conforming method based on further subdivision of the triangles which intersect the interface. The latter method is of optimal order in energy norm and numerical examples indicate second order convergence in the max-norm for this method. An outline of the paper is as follows. In Section 2 we formulate the continuous problem that we aim to solve, in Section 3 we define the numerical method used for the approximation, and in Section 4 we prove the approximation properties of the correspondingfinite element spaces. In Section 5 we prove optimal a priori error estimates and in Section 6 we give corresponding a posteriori error estimates that serve as a basis for adaptive mesh refinement. Finally, in Section 7, we give some implementation details and nu- merical examples. 2. Problem formulation and preliminaries Let X be a bounded domain in R , with convex polygonal boundary oX and an internal smooth boundary C dividing X into two open sets X and X . For any sufficiently regular function u in X [ X we 1 2 1 2 define the jump of u on C by ½u :¼ u j  u j , where u ¼ uj is the restriction of u to X . Conversely, for u 1 2 i i i C C X defined in X we identify the pair fu ; u g with the function u which equals u on X . We consider the i 1 2 i i followingstationary heat conduction problem with a discontinuity in the conductivity across C and an inhomogeneous conormal derivative condition on the interface: ðaruÞ¼ f in X [ X ; 1 2 u ¼0on oX; ð1Þ ½u¼0on C; ½ar u¼ g on C: Here n is the outward pointingunit normal to X and r v ¼ n rv. 1 n For a bounded open connected domain D we shall use standard Sobolev spaces H ðDÞ with norm k r;D r 0 and spaces H ðDÞ with zero trace on oD. The inner products in H ðDÞ¼ L ðDÞ is denoted ð Þ . For a 0 D 2 k bounded open set G ¼[ D , where D are open mutually disjoint components of G, we let H ðD [ D Þ i i 1 2 i¼1 denote the Sobolev space of functions in G such that uj 2 H ðD Þ with norm 1=2 k ¼ k k : k;D [D k;D 1 2 i i¼1 2 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5539 1=2 We assume that f 2 L ðXÞ, g 2 H ðCÞ and, for simplicity, that a is constant in X with a > 0. The weak 2 i i form of (1) is as follows: find u 2 H ðXÞ such that aðu; vÞ¼ ðf ; vÞ þðg; vÞ ; 8v 2 H ðXÞ: ð2Þ X C 0 Here aðu; vÞ¼ ðaru; rvÞ is the bilinear form correspondingto the elliptic operator. It is known that this problem has a unique solution which is in H on each subdomain. The followinga priori estimate is valid, see [5]: kuk þkuk 6 Cðkf k þkgk Þ: ð3Þ 1;X 2;X [X 0;X 1=2;C 1 2 Here and below, C and c denote generic constants. 3. The approximation In a standard finite element method, the jump in normal derivative resultingfrom the continuity of the flux, when a 6¼ a , can be taken into account by letting C coincide with mesh lines. We will take an al- 1 2 ternative approach and solve (1) approximately usingpiecewise linear finite elements on a family of con- formingtriangulations T of X which are independent of the location of the interface C. Instead, we shall allow the approximation to be discontinuous inside elements which intersect the interface. We will use the followingnotation for mesh related quantities. Let h be the diameter of K and h ¼ max h . For any element K, let K ¼ K \ X denote the part of K in X .By G :¼fK 2 T : K \ C 6¼;g K2T K i i i h h we denote the set of elements that are intersected by the interface. For an element K 2 G , let C :¼ C \ K h K be the part of C in K. We make the followingassumptions regardingthe mesh and the interface. A1: We assume that the triangulation is non-degenerate, i.e., h =q 6 C 8K 2 T ; K h where h is the diameter of K and q is the diameter of the largest ball contained in K. A2: We assume that C intersects each element boundary oK exactly twice, and each (open) edge at most once. A3: Let C be the straight line segment connecting the points of intersection between C and oK. We as- K;h sume that C is a function of length on C ; in local coordinates K K;h C ¼fðn; gÞ : 0 < n < jC j; g ¼ 0g K;h K;h and C ¼fðn; gÞ : 0 < n < jC j; g ¼ dðnÞg: K K;h Since the curvature of C is bounded, the assumptions A2 and A3 are always fulfilled on sufficiently fine meshes. Thus the assumptions are natural and not very restrictive; they ensure that the curvature of the interface is well resolved by the mesh. h h h We shall seek a discrete solution U ¼ðU ; U Þ in the space V ¼ V  V , where 1 2 1 2 h 1 V ¼f/ 2 H ðX Þ : / j is linear; / j ¼ 0g: i i i K i oX 3 5540 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 Note that functions in V may be discontinuous across C. Since C may intersect two edges of a triangle arbitrarily, the size of the parts K are not fully characterized by the meshsize parameters. To define the method, we will therefore use the function j ¼ðj ; j Þ defined on each element by 1 2 jK j j j ¼ ; jKj where jKj :¼ meas K. Clearly, 06 j 6 1 and j þ j ¼ 1 so that i 1 2 f/g :¼ðj / þ j / Þj 1 2 1 2 C is a convex combination of / ¼ð/ ; / Þ along C. 1 2 The method is defined by the variational problem of finding U 2 V such that a ðU; /Þ¼ Lð/Þ; 8/ 2 V ; ð4Þ where a ðU; /Þ :¼ða rU ; r/ Þ ð½U; far /gÞ ðfar Ug; ½/Þ þðk½U; ½/Þ h i i n n i X [X C C C 1 2 with k sufficiently large (see Lemma 5 below), and Lð/Þ :¼ðf ; /Þ þðj g; / Þ þðj g; / Þ : 2 1 X 1 C 2 C In this method, the conditions at C are satisfied weakly by means of a variant of Nitsches method. With these definitions, we have the followingconsistency relation. Lemma 1. The discrete problem (4) is consistent in the sense that, for u solving (1), a ðu; /Þ¼ Lð/Þ; 8/ 2 V : Proof. We first note that, for u solving(1), g far ug¼ðj þ j Þg far ug j ðg ½ar uÞ n 1 2 n 1 n ¼ j g  j a r u  j a r u þ j a r u  j a r u ¼ j g  a r u ; 2 1 1 n 1 2 2 n 2 1 1 n 1 1 2 n 2 2 2 n 2 and, similarly, g far ug¼ðj þ j Þg far ugþ j ðg ½ar uÞ n 1 2 n 2 n ¼ j g  j a r u  j a r u  j a r u þ j a r u ¼ð1 þ j Þg  a r u ; 1 1 1 n 1 2 2 n 2 2 1 n 1 2 2 n 2 2 1 n 1 so that far ug¼ a r u  j g ¼ a r u þ j g: ð5Þ n 1 n 1 2 2 n 2 1 Since ½u¼ 0, we may use (5) and Greens formula to obtain a ðu; /Þ¼ ðaru; r/Þ ðfar ug; /  / Þ h n X [X 1 2 C 1 2 ¼ðaru; r/Þ ða r u  j g; / Þ þða r u þ j g; / Þ 1 n 1 2 1 2 n 2 1 2 X [X C C 1 2 ¼ ðr ðaruÞ; /Þ þðj g; / Þ þðj g; / Þ ¼ðf ; /Þ þðj g; / Þ þðj g; / Þ ¼ Lð/Þ; 2 1 2 1 X [X 1 C 2 C X 1 C 2 C 1 2 which is the statement of the lemma. 4 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5541 An immediate consequence of Lemma 1 is the condition a ðu  U; /Þ¼ 0; 8/ 2 V ; ð6Þ h h which we will refer to as Galerkin orthogonality. A FE basis for V is easily obtained from a standard FE basis on the mesh by the introduction of new basis functions for the elements that intersect C. Thus, we replace each standard basis function livingon an element that intersects the interface by two new basis functions, namely its restrictions to X and X , re- 1 2 spectively. The collection of basis functions with support in X is then clearly a basis for V , and hence we obtain a basis for V by the identification w ¼ðwj ; wj Þ. If the interface coincides exactly with an element X X 1 2 edge, no new basis functions are introduced on these elements but the approximating functions may still be discontinuous over such an edge. As a consequence, there are six non-zero basis functions on each element that properly intersects C. Further implementation details are considered in Section 7. 4. Approximation property of V Recall that G denotes the set of elements that are intersected by the interface. We will use the following mesh dependent norms: 2 2 kvk :¼ h kvk ; 1=2;h;C K 0;C K2G 2 2 kvk :¼ h kvk ; 1=2;h;C 0;C K2G and 2 2 2 2 jjjvjjj :¼krvk þ kfr vgk þk½vk : 0;X [X 1=2;h;C 1=2;h;C 1 2 We note for future reference that ðu; vÞ 6 kvk kvk : ð7Þ C 1=2;h;C 1=2;h;C 1 2 To show that functions in V approximates functions v 2 H ðXÞ\ H ðX [ X Þ to the order h in the norm h 1 2 jjj jjj, we construct an interpolant of v by nodal interpolants of H -extensions of v and v as follows. 1 2 2 2 Choose extension operators E : H ðX Þ! H ðXÞ such that ðE wÞj ¼ w and i i i kE wk 6 Ckwk 8w 2 H ðX Þ; s ¼ 0; 1; 2: ð8Þ i i s;X s;X Let I be the standard nodal interpolation operator and define I v :¼ðI v ; I v Þ where I v :¼ðI E v Þj : ð9Þ 1 2 i h i i h h;1 h;2 h;i X The followingtheorem is valid. Theorem 2. Let I be an interpolation operator defined as in (9). Then 1 2 jjjv  I vjjj6 C hkvk ; 8v 2 H ðXÞ\ H ðX [ X Þ: A 1 2 h 2;X [X 0 1 2 In the proof of this result, we need to estimate the interpolation error at the interface. To that end, we shall use the followingvariant of a well known trace inequality on a reference element. The crucial fact is that the constant in this inequality is independent of the location of the interface relative to the mesh. 5 5542 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 e e Lemma 3. Map a triangle K 2 G onto the unit reference triangle K K by an affine map and denote by C C the h ~ K K corresponding image of C . Under assumptions A1–A3of Section 3 there exist a constant C, depending on C but independent of the mesh, such that kwk 6 Ckwk kwk ; 8w 2 H ðK KÞ: ð10Þ ~ ~ ~ 0;C C 0;K K 1;K K K K Proof. We start by showingthat 2 2 kwk 6 Cðkwk þkwk kwk Þ: ð11Þ ~ ~ ~ ~ 0;C C 0;C C 0;K K 1;K K ~ ~ K K K K;h Recall that C is the straight line connecting the points of intersection between C and the element K and K;h C ¼fðn; gÞ : 0 < n < jC j; g ¼ dðnÞg: K K;h Assume first that dðnÞ > 0. Since the curvature of the interface is bounded, jd ðnÞj6 CjC j. As the mesh K;h is non-degenerate this implies that on the reference element we may write, using again ðn; gÞ as local co- ordinates, e e C C ¼fðn; gÞ : 0 < n < jC C j; g ¼ d dðnÞg; ~ ~ K K K K;h ~ e e where jd d ðnÞj6 C. We now let D denote the domain bounded by C C ~ and C C ~ and note that by the diver- K K K K;h gence theorem, Z Z Z Z ow 2 1=2 2 2 2 0 2 w dndg ¼ divð0; w Þdndg ¼ w dn þ w ð1 þðd d Þ Þ ds: ð12Þ og ~ ~ D D C C C C ~ ~ K K;h K K As d d is bounded, 2 0 2 1=2 kwk 6 C w ð1 þðd d Þ Þ ds; 0;C C K K C C K K whence (11) follows from (12) usingCauchy–Schwarz  inequality. In a general case where d may switch sign, the same argument may be applied for each part between the e e intersections of C C and C C . ~ ~ K K K K;h It remains to show that the first term on the right in (11) is appropriately bounded. To that end we shall e e e map the triangular part K K and the quadrilateral part K K of K K onto new reference domains. We may assume t q e e that C C ~ intersects K K in ða; 0Þ and in ð0; bÞ, and, by symmetry, that 06 a6 b6 1. K K;h For a ¼ b ¼ 1, the desired trace inequality kwk 6 Ckwk kwk ð13Þ ~ ~ ~ 0;C C 0;K K 1;K K K K;h is valid. For 1=2 < a6 b < 1, we may map the triangular part K K onto the unit reference triangle by a linear map. By the bound from below on a and b, this map is bounded, uniformly in a and b, with uniformly bounded inverse, and hence (13) is valid also in this case. For 1=2 < a < b ¼ 1 the same argument holds, choosingthis time K K as the triangular part which contains the origin. Assume now that a6 1=2. Let ðxx^; yy^Þ¼ Mð~xx; yy~Þ¼ ðyy~; ð1  aÞ ð~xx þ yy~ 1ÞÞ: b e b Then the image K K ¼ MðK K Þ has its corners in ð0; 0Þ, ð1; 0Þ, ð0; 1Þ, PP ¼ðb; ð1  bÞ=ð1  aÞÞ, and there q q holds kwk 6 CðPP Þkwk kwk : ð14Þ ^ ^ ^ 0;C C 0;K K 1;K K ^ q q K K;h 6 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5543 b b An additional argument is needed to show uniformity in PP . Since 06 a6 1=2 and a6 b6 1, PP varies in the domain D D :¼f06 xx ^6 1=2; 1  xx ^6 yy^6 1g[f1=26 xx^6 1; 1  xx^6 yy^6 2ð1  xx ^Þg as a and b vary. Let kwk 0;C C K K;h F ðPP ; w w ^Þ¼ : kwk kwk ^ ^ 0;K Kq 1;K Kq b b b b b We will show that F ðPPÞ¼ sup F ðPP ; w w ^Þ is uniformly bounded. For points RR and SS in D D, assuming 1 ^ w2H ðK K Þ b b without restriction that F ðRRÞ P F ðSSÞ, we have for any w that b b b b b b b b F ðRRÞ F ðSSÞ¼ sup F ðRR; vv^Þ sup F ðSS; vv^Þ6 j sup F ðRR; vv^Þ F ðRR; w w ^Þj þ jF ðRR; w w ^Þ F ðSS; w w ^Þj ¼ I þ II: vv^ vv^ vv^ Given > 0 we may choose w w ^ such that I 6 =2. Note that F ðRR; vv^Þ is continuous for fixed vv^ since the only b b b dependence of RR lies in the domains of integration. We may thus take jRR  SSj small enough so that II 6 =2. b b b Hence F ðPP Þ is continuous on the compact set D D, and thus (14) holds uniformly in PP . Finally, since M is bounded, uniformly in a and b, with uniformly bounded inverse, (13) follows and the proof is complete. Proof of Theorem 2. Recall that K ¼ K \ X and let v ¼ E v denote the extension of v to X. By a standard i i i i i interpolation estimate we obtain krðv  I v Þk ¼ krðv  I v Þk 6 krðv  I v Þk 6 Chkv k : i i h h h;i 0;K i i 0;K i i 0;K i 2;K i i Summingover all triangles that intersect X , it follows by (8) that 2 2 2 2  2 krðv  I v Þk 6 Ch kv k 6 Ch kv k : ð15Þ i i i h;i 0;X i 2;K 2;X i i K\X 6¼; Next we consider the jumps on the interface. Since the mesh is non-degenerate, it follows from Lemma 3, scaled by the map from the reference triangle, that 2 2 2 1 2 1 h kwk 6 Cðh kwk þkwk Þ; 8w 2 H ðKÞ: K 0;C K 0;K 1;K Hence it follows, usingagain a standard interpolation estimate, that X X 2 2 2 1  1  1 h k½v  I vk 6Ch kv  I v k ¼ Ch kv  I v k i i h K h 0;C K h;i 0;C K i i 0;C K K K i i X X 2 2 2 2     2 6C ðh kv  I v k þkv  I v k Þ6Ch kv k : h h K i i 0;K i i 1;K K i 0;K i i Summingthe contributions from K 2 G , we get from (8) that k½v  I vk 6 Ch kv k 6 Chkvk : ð16Þ h 1=2;h;C i 2;[K2G 2;X [X h 1 2 i¼1 Finally, Lemma 3 applied to r w and scalinggives 2 2 2 2 2 h kr wk 6 Cðkwk þ h kwk Þ; 8w 2 H ðKÞ; K n 0;C 1;K K 2;K whence similar arguments as above yield kr ðv  I v Þk 6 Chkv k : ð17Þ n i i i h;i 1=2;h;C 2;X Since j < 1, the theorem now follows from (15)–(17). 7 5544 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5. A priori error estimates We will first show coercivity of the discrete form, for which purpose we will need the followinginverse inequality. Lemma 4. For / 2 V , the following inverse inequality holds: 2 2 kfr /gk 6 C kr/k : n I 1=2;h;C 0;X [X 1 2 Proof. Since / 2 V is linear on K , we have h i jC j jC jjK j 2 2 K 2 K i 2 2 2 2 h kj r / k 6 h j jC jjr/ j ¼ h j kr/ k ¼ h kr/ k 6 Ckr/ k : K i n i K K i K i K i i 0;C i i 0;K 0;K 0;K K i 2 i i jK j jKj In the last step above we have used that jC j6 h , jK j6 h , and, since the mesh is non-degenerate, K K i jKj P ch . The result follows by summation over the elements. Lemma 5. The discrete form a ð Þ is coercive on V , i.e., a ðv; vÞ P Cjjjvjjj 8v 2 V ; provided k is chosen sufficiently large. It is also continuous, i.e., a ðu; vÞ6 Cjjjujjj jjjvjjj 8u 2 V ; 8v 2 V : Proof. Continuity of the discrete form follows directly from the definitions. To prove coercivity, we use (7) to find that for any > 0 2 1=2 2 1=2 a ðv; vÞ¼ka rvk  2ð½v; far vgÞ þkk ½vk h n 0;X [X C 0;C 1 2 2 2 1=2 1=2 P ka rvk  2kfar vgk k½vk þkk ½vk 0;X [X 1=2;h;C 1=2;h;C 0;C 1 2 2 2 2 1=2 P ka rvk  kfar vgk þ k  k½vk : 0;X [X 1=2;h;C 0;C 1 2 K K2G It then follows from Lemma 4 that 1 1 2C max a I X 2 2 1=2 1=2 a ðv; vÞ P ka rvk þ  ka rvk 0;X [X 0;X [X 1 2 1 2 2 2 2 2 þ kfar vgk þ k  k½vk : 1=2;h;C 0;C K2G Taking  ¼ 4C max a, coercivity follows if kj ¼ ch where c > 4C max a. I X I X K K Theorem 6. Under assumptions A1–A3of Section 3, and for U solving (4) and u solving (1), the following a priori error estimates hold: jjju  Ujjj6 Chkuk ð18Þ 2;X [X 1 2 and ku  Uk 6 Ch kuk : ð19Þ 0;X 2;X [X 1 2 8 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5545 Proof. For any v 2 V , jjju  Ujjj6 jjju  vjjj þ jjjv  Ujjj. Further, by Lemma 5 and orthogonality, we have that jjjU  vjjj 6 Ca ðU  v; U  vÞ¼ Ca ðu  v; U  vÞ6 Cjjju  vjjj jjjU  vjjj; h h and it follows that jjju  Ujjj6 Cjjju  vjjj 8v 2 V : Taking v ¼ I u and invokingthe interpolation result of Theorem 2, (18) follows. For (19) we use a duality argument. Define z ¼ðz ; z Þ by 1 2 ða rz Þ¼ e in X ; i ¼ 1; 2; i i i i z ¼0on oX \ oX ; i i ð20Þ ½z¼0on C; ½ar z¼0on C; where e ¼ u  U . By Greens formula and (5) with u ¼ z, g ¼ 0, we have that i i i kek ¼ ðr ðarzÞ; eÞ ¼ðarz; reÞ ða r z ; e Þþða r z ; e Þ 1 n 1 1 2 n 2 2 0;X X [X X 1 2 ¼ðarz; reÞ ðfar zg; ½eÞ ¼ a ðz; eÞ n h X C since ½z¼ 0. Thus, usingthe symmetry of a ð Þ and applyingthe orthogonality relation (6) and Theorem 2, we find that kek ¼ a ðz  I z; eÞ6 Cjjjz  I zjjj jjjejjj6 Chkzk jjjejjj: ð21Þ h h h 0;X 2;X [X 1 2 Finally, by the elliptic regularity result (3), we have kzk 6 Ckek , whence the estimate (19) follows 2;X [X 0;X 1 2 from (21) and (18). 6. A posteriori error estimates In this Section, we prove a posteriori error estimates and formulate adaptive algorithms for the finite element method (4), followingBecker and Rannacher [3]. We will consider control of linear functionals jðeÞ of the error, and define the local and global estimators as 2 2 2 2 1 2 1=2 E ðUÞ¼ ðh kf þr ðarUÞk þ h k½a rU k þ h kg ½arUk þ h k½Uk Þ ð22Þ K K i i K K 0;K [K 0;oK 0;C K 0;C K K 1 2 and 1=2 EðUÞ¼ h E ðUÞ : ð23Þ K2T We then have the followinga posteriori error estimate. Theorem 7. For a continuous linear functional jð Þ on L ðXÞ, let J 2 L ðXÞ be defined by Riesz’ representation 2 2 theorem, i.e., jð Þ :¼ðJ; Þ . Then there is a positive constant C such that jðeÞ6 CEðUÞkJk : ð24Þ 0;X 9 5546 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 Proof. Let z be the solution to the problem ða rz Þ¼ J in X ; i ¼ 1; 2; i i i z ¼0on oX \ oX ; i i ð25Þ ½z¼0on C; ½ar z¼0on C: We first note that ðe; JÞ ¼ðare; rzÞ ða r z ; e Þ þða r z ; e Þ : 1 n 1 1 2 n 2 2 X X [X C C 1 2 Now, since ½u¼ 0, we see that, by (5), ða r z ; e Þ þða r z ; e Þ ¼ðU ; a r z Þ ðU ; a r z Þ ¼ð½U; far zgÞ: 1 n 1 1 2 n 2 1 1 1 n 1 2 2 n 2 n C C C C Thus, jðeÞ¼ðe; JÞ ¼ðare; rzÞ þð½U; far zgÞ: ð26Þ X X [X 1 2 Now, take Z 2 V . From Galerkin orthogonality we then have a ðe; ZÞ¼ 0, and since ½z¼ 0, ½u¼ 0 on the h h interface, we get 0 ¼ a ðe; ZÞ¼ ðare; rZÞ ð½e; far ZgÞ ðfar eg; ½ZÞ þðk½e; ½ZÞ h n n X [X C C C 1 2 ¼ðare; rZÞ þð½U; far ZgÞ þðfar eg; ½z  ZÞ þðk½U; ½z  ZÞ : ð27Þ n n X [X C C C 1 2 Denote by n the outward pointingunit normal to K. Subtracting(27) from (26) and integratingby parts we get jðeÞ¼ðare; rðz  ZÞÞ þð½U; far ðz  ZÞgÞ ðfar eg; ½z  ZÞ ðk½U; ½z  ZÞ n n X [X C C C 1 2 ¼ ðf r ðarUÞ; z  ZÞ þð½U; far ðz  ZÞgÞ ðfar eg; ½z  ZÞ ðk½U; ½z  ZÞ n n K [K C C C 1 2 K2T ð½an rU; z  ZÞ þða r e ; z  Z Þ ða r e ; z  Z Þ : K 1 n 1 1 1 2 n 2 2 2 oKnC C C K2T We now note that ða n re ; z  Z Þ ða n re ; z  Z Þ ðfar eg; ½z  ZÞ 1 1 1 1 2 2 2 2 n C C C ¼ða r e ; z  Z Þ ða r e ; z  Z Þ ðj a r e þ j a r e ; ½z  ZÞ 1 n 1 1 1 2 n 2 2 2 1 1 n 1 2 2 n 2 C C C ¼ðj ½ar e; z  Z Þ þðj ½ar e; z  Z Þ 2 n 1 1 1 n 2 2 C C ¼ðj ðg ½ar UÞ; z  Z Þ þðj ðg ½ar UÞ; z  Z Þ 2 n 1 1 1 n 2 2 C C ¼ ðð1  j Þðg ½ar UÞ; z  Z Þ j n j j j¼1 and thus we find that jðeÞ¼ ðf r ðarUÞ; z  ZÞ  ð½an rU; z  ZÞ K [K oKnC 1 2 K2T X X X þ ð½U; far ðz  ZÞg  k½z  ZÞ þ ðð1  j Þðg ½ar UÞ; z  Z Þ : ð28Þ n j n j j C C K K K2G K2G j¼1 h h 10 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5547 Further, by Cauchy–Schwarz inequality, and choosing kj ¼ ch with c > 4C max a , I X i K K 2 2 X X X X jðeÞ6 q x þ q x ; ð29Þ K;j S;j K;j S;j K2T j¼1 K2G j¼1 h h where q ¼ h kf þr ðarUÞk ; x ¼ h kz  Zk ; K K;1 K;1 0;K [K K 0;K [K 1 2 1 2 1=2 1=2 q ¼ h k½an rUk ; x ¼ h kz  Zk ; K K;2 K;2 K 0;oKnC K 0;oKnC 1=2 1=2 q ¼ h k½uk ; x ¼ h kfar ðz  ZÞg  c½z  Z=h k S;1 n K S;1 K 0;C K 0;C K K and 1=2 1=2 q ¼ h kg ½ar Uk ; x ¼ h max kz  Z k : n S;2 j j S;2 K 0;C K 0;C K K Note that the right hand side of (29) is a weighted sum of local residuals q , q . To conclude the proof K;i S;i we bound the weights in (29) by choosing Z ¼ I z, applyingthe interpolation estimates (15)–(17) and finally the regularity estimate kzk 6 CkJk ; 2;X [X 0;X 1 2 whence (24) follows. 7. Implementation and numerical examples Recall that a FE basis for V is obtained from a standard FE basis on the mesh by replacingeach standard basis function livingon an element that intersects the interface by two new basis functions, namely its restrictions to X and X , respectively. Both these new basis functions are, however, represented in the 1 2 implementation using the same nodes from the original triangulation. The points of intersection between the element edges and the interface are not used to represent the new basis function, and the geometry of the interface and the element parts does not come into play until integrating the terms in the bilinear form. One may thus alternatively think of the approximation as beingdefined on two overlappingmeshes formed by triangles covering the subdomains. In order to implement the discontinuous approximation, we first determined the set G of triangles 0 00 intersected by C. For each K 2 G , we assigned two identical copies K and K . We assumed that K, with nodes fi; j; kg was split by C into a triangle and a quadrilateral; the case of a split into two triangles was handled by creating a quadrilateral with one side of zero length. We assigned to K the triangular part and 00 00 0 0 to K the quadrilateral part. Thus, on K we created a new node i on the far side of C, and on K we created 0 0 0 two new nodes, k and j , on the other side of C, see Fig. 1. To ensure continuity across the edges of K and K (away from C), we also checked if the new nodes had already been created by the same process on the neighboring elements. After having completed this process, we thus had two independent meshes, one completely covering X , and the other completely covering X . The elements crossed by C had been dou- 1 2 bled, but coincided geometrically. To numerically evaluate a ð Þ and Lð Þ, we used the following strategy. The triangles that were not crossed by C were handled in the usual way. On the elements crossed by C, we used centroid quadrature to 11 5548 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 j′ k′ K′ K′′ i′ Fig. 1. The split of a triangle used in the implementation of the proposed method. evaluate all terms, both on the quadrilateral side and on the triangular side. On the interface, we used two- point Gaussian quadrature. All contributions were then assembled usingthe old and new nodes defined 0 0 0 by the splittingprocess. We emphasize that the new nodes i , j , and k are to be considered convenient support points for the definition of a continuous, piecewise linear, approximation rather than nodes in the standard finite element sense. The solution at these points is computed but is outside the domain of interest. 7.1. Example 1 We considered solutions to the ordinary differential equation d du du du i 1 2 a ¼ 1; ½uð1=2Þ ¼ 0; a ð1=2Þ¼ a ð1=2Þ: i 1 2 dx dx dx dx The domain is ð0; 1Þ, with an interface at x ¼ 1=2. While this is a one-dimensional problem, we solved it numerically in 2D on the domain ð0; 1Þð0; 1Þ, with zero Neumann boundary conditions at y ¼ 0and y ¼ 1. The equation has a closed-form solution, given by 2 2 ð3a þ a Þx x a  a þð3a þ a Þx x 1 2 2 1 1 2 u ðxÞ¼  ; u ðxÞ¼  : 1 2 2 2 4a þ 4a a 2a 4a þ 4a a 2a 1 2 1 1 2 2 1 2 We chose a ¼ 1=2, a ¼ 3 and performed a numerical convergence test for different approaches: a 1 2 standard unfitted FE method, a fitted standard FE method, and the proposed unfitted method. The con- vergence results are given in Fig. 2, and show the suboptimal behaviour of a standard unfitted method, and the small difference in computational error between the proposed method and a fitted method. We remark that an exact correspondence cannot be obtained, since the meshes will not be exactly the same for the fitted and unfitted methods. In Fig. 3 elevations of the solutions obtained with the different unfitted methods are shown on the final mesh consistingof 2048 elements. 12 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5549 -3.5 Nitsche -4 Unfitted standard FEM Fitted standard FEM -4.5 -5 -5.5 -6 -6.5 -7 -7.5 -8 -8.5 -5 -4 -3 -2 -1 0 1 Log(h) Fig. 2. L -norm convergence of different methods applied to the interface problem. Fig. 3. Elevation of the discontinuous approximation obtained from the proposed method (left) and the continuous standard FE approximation (right) on the final mesh. 7.2. Example 2 Here we considered a less trivial two-dimensional example (from [7]). The exact solution is given by ; if r6 r ; uðx; yÞ¼ 2 2 2 r r r 0 0 þ ; if r > r ; a a a 2 2 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 where r :¼ x þ y , and, on the domain ð0; 1Þð0; 1Þ, we chose r ¼ 3=4, a ¼ 1, and a ¼ 1000. The 0 1 2 boundary conditions were symmetry boundaries at x ¼ 0 and y ¼ 0 and Dirichlet boundary conditions correspondingto the exact solution at x ¼ 1 and y ¼ 1. The right-hand side yielding this exact solution is f ¼4. In Fig. 4, we give the elevation of the approximate solution on the last mesh in a sequence. The corresponding L - and maximum norm convergence is given in Fig. 5. We note in particular that the order of convergence in maximum norm seems to be approximately equal to the (second order) L -norm con- vergence, though we have no theoretical results to back this up. Log(L2 error) 5550 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.8 0.6 0 0.4 1 0.2 0.9 0.8 0.7 0.6 0.5 0.4 0 0.3 0.2 0.1 Fig. 4. Elevation of the approximate solution. L2 error -4.5 Maximum norm error -5 -5.5 -6 -6.5 -7 -7.5 -8 -5 -4 -3 -2 -1 0 Log(h) Fig. 5. L -norm and maximum norm convergence. 7.3. Example 3 The third example was solved on the domain ð0; 1Þð0; 1Þ, with zero Dirichlet boundary conditions. Centred in the domain is an ellipsoidal inclusion with conduction parameter a ¼ 1; outside of the ellipse min we set a ¼ 6. In Fig. 6 we show the first mesh and associated solution. We used an adaptive algorithm max correspondingto control of the L ðXÞ-error, but made no attempt to tune the interpolation constants or solve a dual problem. Thus, the example only gives an indication of how the adapted meshes appear in such a case. For each successive mesh we refined the third of the elements containingthe highest element contribution to the total error as defined by (29). In Fig. 7 we give the last adapted mesh in a sequence, together with an elevation of the corresponding solution. Log(error) A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5551 Fig. 6. First mesh and elevation of the approximate solution. Fig. 7. Last adaptively refined mesh and elevation of the approximate solution. 8. Concluding remarks In this paper, we have introduced and analysed a new method for elliptic interface problems on unfitted meshes, i.e., meshes which are independent of the location of the interface. Unlike the standard unfitted finite element method, the proposed approach leads to optimal convergence rates. This have been shown in the model situation of piecewise linear approximations in two dimensions. In future work, we will address the general situation of higher order polynomial approximations in three dimensions. Other extensions under consideration include Stefan problems and fictitious domain type simulations. References [1] I. Babusska, The finite element method for elliptic equations with discontinuous coefficients, Computing5 (1970) 207–213. [2] J.W. Barrett, C.M. Elliott, Fitted and unfitted finite-element methods for elliptic equations with smooth interfaces, IMA J. Numer. Anal. 7 (1987) 283–300. [3] R. Becker, R. Rannacher, A feed-back approach to error control in finite element methods: basic analysis and examples, East- West J. Numer. Math. 4 (1996) 237–264. [4] J.H. Bramble, J.T. King, A finite element method for interface problems in domains with smooth boundaries and interfaces, Adv. Comput. Math. 6 (1996) 109–138. [5] Z. Chen, J. Zhou, Finite element methods and their convergence for elliptic and parabolic interface problems, Numer. Math. 79 (1998) 175–202. [6] M. Feistauer, V. Sobot ııkov aa, Finite element approximation of nonlinear problems with discontinuous coefficients, M2AN Math. Model. Numer. Anal. 24 (1990) 457–500. 15 5552 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 [7] Z. Li, T. Lin, X. Wu, New Cartesian grid methods for interface problems using the finite element formulation, NCSU Technical Report: CRSC-TR99-12, Center for Research in Scientific Computation, North Carolina State University, 1999. [8] R.J. MacKinnon, G.F. Carey, Treatment of material discontinuities in finite element computations, Int. J. Numer. Meth. Eng. 24 (1987) 393–417. [9] J. Nitsche, U Uber ein Variationsprinzip zur Loosungvon € Dirichlet–Problemen bei Verwendungvon Teilr €aaumen, die keinen Randbedingungen unterworfen sind, Abh. Math. Univ. Hamburg 36 (1971) 9–15. [10] A. Z Zen ııssek, The finite element method for nonlinear elliptic equations with discontinuous coefficients, Numer. Math. 58 (1990) 51–77. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Computer Methods in Applied Mechanics and Engineering Unpaywall

An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems

Hansbo, Anita; Hansbo, Peter
Computer Methods in Applied Mechanics and Engineering – Nov 1, 2002
Download PDF
17 pages

Loading next page...
 
/lp/unpaywall/an-unfitted-finite-element-method-based-on-nitsche-s-method-for-eBr9XLLP70

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
Unpaywall
ISSN
0045-7825
DOI
10.1016/s0045-7825(02)00524-8
Publisher site
See Article on Publisher Site

Abstract

An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems Anita Hansbo, Peter Hansbo To cite this version: Anita Hansbo, Peter Hansbo. An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems. Computer Methods in Applied Mechanics and Engineering, Elsevier, 2002, 191 (47-48), pp. 5537-5552. 10.1016/S0045-7825(02)00524-8. hal-01352903 HAL Id: hal-01352903 https://hal.archives-ouvertes.fr/hal-01352903 Submitted on 10 Aug 2016 HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non, lished or not. The documents may come from émanant des établissements d’enseignement et de teaching and research institutions in France or recherche français ou étrangers, des laboratoires abroad, or from public or private research centers. publics ou privés. Distributed under a Creative Commons Attribution| 4.0 International License An unfitted finite element method, based on Nitsches method, for elliptic interface problems a b Anita Hansbo , Peter Hansbo Department of Informatics and Mathematics, University of Trollh€aattan-Uddevalla, Box 957, S-461 39 Trollh€aattan, Sweden Department of Applied Mechanics, Chalmers University of Technology, S-412 96G €ooteborg, Sweden In this paper we propose a method for the finite element solution of elliptic interface problem, usingan approach due to Nitsche. The method allows for discontinuities, internal to the elements, in the approximation across the interface. We show that optimal order of convergence holds without restrictions on the location of the interface relative to the mesh. Further, we derive a posteriori error estimates for the purpose of controllingfunctionals of the error and present some numerical examples. 1. Introduction As a model elliptic interface problem, we consider a stationary heat conduction problem in two di- mensions with a conduction coefficient which is discontinuous across a smooth internal interface. When solvingsuch problems numerically usingthe standard finite element method, one usually takes the dis- continuity of the data into account by enforcingmesh lines along the interface. If this is not done, sub- optimal convergence behaviour will occur, cf. [1,8]. As a motivation for this work, we also have in mind more complicated, time dependent or non-linear, problems where the interface moves with time or duringiteration. In that case, it may be advantageous to use the same mesh on the domain for different, nearby, locations of the interface, since repeated remeshing of the domain to obtain fitted meshes is very costly. We are thus led to study unfitted finite element methods where the interface is allowed to cross the elements. In this paper, we propose an unfitted finite element method, based on a variant of Nitsches method [9], allowingfor discontinuities, internal to the elements, in the approximation across the interface. This method is of optimal order; in particular we show second order convergence in L for appropriately modified piecewise linears on a non-degenerate triangulation. We also consider a posteriori error estimates 1 5538 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 for functionals of the solution, in the spirit of Becker and Rannacher [3], and use these estimates as a basis for adaptively refiningthe mesh. Fitted mesh FE methods for elliptic problems with discontinuous coefficients and homogeneous interface conditions are analysed in [1,6,10]. In [2,4,5], problems with inhomogeneous interface conditions are considered. As for unfitted mesh methods for interface problems, Barrett and Elliott [2] show first order of con- vergence in energy-norm and interior second order L error estimates for a piecewise linear method based on boundary penalty and numerical integration over approximate domains. This approach is close in spirit to the method to be presented here, but contains a consistency error that we avoid. An alternative approach is to construct a FEM-basis where the basis functions fulfill homogeneous interface conditions exactly, and to use the continuous bilinear form (without penalty) to define the method. MacKinnon and Carey [8] use linear basis functions with this property in two dimensions and numerical examples of optimal order of convergence are presented. Li et al. [7] analyse the max-norm interpolation error on cartesian grids of a non-conforming piecewise linear method where the basis func- tions fulfill homogeneous interface conditions exactly, as well as a conforming method based on further subdivision of the triangles which intersect the interface. The latter method is of optimal order in energy norm and numerical examples indicate second order convergence in the max-norm for this method. An outline of the paper is as follows. In Section 2 we formulate the continuous problem that we aim to solve, in Section 3 we define the numerical method used for the approximation, and in Section 4 we prove the approximation properties of the correspondingfinite element spaces. In Section 5 we prove optimal a priori error estimates and in Section 6 we give corresponding a posteriori error estimates that serve as a basis for adaptive mesh refinement. Finally, in Section 7, we give some implementation details and nu- merical examples. 2. Problem formulation and preliminaries Let X be a bounded domain in R , with convex polygonal boundary oX and an internal smooth boundary C dividing X into two open sets X and X . For any sufficiently regular function u in X [ X we 1 2 1 2 define the jump of u on C by ½u :¼ u j  u j , where u ¼ uj is the restriction of u to X . Conversely, for u 1 2 i i i C C X defined in X we identify the pair fu ; u g with the function u which equals u on X . We consider the i 1 2 i i followingstationary heat conduction problem with a discontinuity in the conductivity across C and an inhomogeneous conormal derivative condition on the interface: ðaruÞ¼ f in X [ X ; 1 2 u ¼0on oX; ð1Þ ½u¼0on C; ½ar u¼ g on C: Here n is the outward pointingunit normal to X and r v ¼ n rv. 1 n For a bounded open connected domain D we shall use standard Sobolev spaces H ðDÞ with norm k r;D r 0 and spaces H ðDÞ with zero trace on oD. The inner products in H ðDÞ¼ L ðDÞ is denoted ð Þ . For a 0 D 2 k bounded open set G ¼[ D , where D are open mutually disjoint components of G, we let H ðD [ D Þ i i 1 2 i¼1 denote the Sobolev space of functions in G such that uj 2 H ðD Þ with norm 1=2 k ¼ k k : k;D [D k;D 1 2 i i¼1 2 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5539 1=2 We assume that f 2 L ðXÞ, g 2 H ðCÞ and, for simplicity, that a is constant in X with a > 0. The weak 2 i i form of (1) is as follows: find u 2 H ðXÞ such that aðu; vÞ¼ ðf ; vÞ þðg; vÞ ; 8v 2 H ðXÞ: ð2Þ X C 0 Here aðu; vÞ¼ ðaru; rvÞ is the bilinear form correspondingto the elliptic operator. It is known that this problem has a unique solution which is in H on each subdomain. The followinga priori estimate is valid, see [5]: kuk þkuk 6 Cðkf k þkgk Þ: ð3Þ 1;X 2;X [X 0;X 1=2;C 1 2 Here and below, C and c denote generic constants. 3. The approximation In a standard finite element method, the jump in normal derivative resultingfrom the continuity of the flux, when a 6¼ a , can be taken into account by letting C coincide with mesh lines. We will take an al- 1 2 ternative approach and solve (1) approximately usingpiecewise linear finite elements on a family of con- formingtriangulations T of X which are independent of the location of the interface C. Instead, we shall allow the approximation to be discontinuous inside elements which intersect the interface. We will use the followingnotation for mesh related quantities. Let h be the diameter of K and h ¼ max h . For any element K, let K ¼ K \ X denote the part of K in X .By G :¼fK 2 T : K \ C 6¼;g K2T K i i i h h we denote the set of elements that are intersected by the interface. For an element K 2 G , let C :¼ C \ K h K be the part of C in K. We make the followingassumptions regardingthe mesh and the interface. A1: We assume that the triangulation is non-degenerate, i.e., h =q 6 C 8K 2 T ; K h where h is the diameter of K and q is the diameter of the largest ball contained in K. A2: We assume that C intersects each element boundary oK exactly twice, and each (open) edge at most once. A3: Let C be the straight line segment connecting the points of intersection between C and oK. We as- K;h sume that C is a function of length on C ; in local coordinates K K;h C ¼fðn; gÞ : 0 < n < jC j; g ¼ 0g K;h K;h and C ¼fðn; gÞ : 0 < n < jC j; g ¼ dðnÞg: K K;h Since the curvature of C is bounded, the assumptions A2 and A3 are always fulfilled on sufficiently fine meshes. Thus the assumptions are natural and not very restrictive; they ensure that the curvature of the interface is well resolved by the mesh. h h h We shall seek a discrete solution U ¼ðU ; U Þ in the space V ¼ V  V , where 1 2 1 2 h 1 V ¼f/ 2 H ðX Þ : / j is linear; / j ¼ 0g: i i i K i oX 3 5540 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 Note that functions in V may be discontinuous across C. Since C may intersect two edges of a triangle arbitrarily, the size of the parts K are not fully characterized by the meshsize parameters. To define the method, we will therefore use the function j ¼ðj ; j Þ defined on each element by 1 2 jK j j j ¼ ; jKj where jKj :¼ meas K. Clearly, 06 j 6 1 and j þ j ¼ 1 so that i 1 2 f/g :¼ðj / þ j / Þj 1 2 1 2 C is a convex combination of / ¼ð/ ; / Þ along C. 1 2 The method is defined by the variational problem of finding U 2 V such that a ðU; /Þ¼ Lð/Þ; 8/ 2 V ; ð4Þ where a ðU; /Þ :¼ða rU ; r/ Þ ð½U; far /gÞ ðfar Ug; ½/Þ þðk½U; ½/Þ h i i n n i X [X C C C 1 2 with k sufficiently large (see Lemma 5 below), and Lð/Þ :¼ðf ; /Þ þðj g; / Þ þðj g; / Þ : 2 1 X 1 C 2 C In this method, the conditions at C are satisfied weakly by means of a variant of Nitsches method. With these definitions, we have the followingconsistency relation. Lemma 1. The discrete problem (4) is consistent in the sense that, for u solving (1), a ðu; /Þ¼ Lð/Þ; 8/ 2 V : Proof. We first note that, for u solving(1), g far ug¼ðj þ j Þg far ug j ðg ½ar uÞ n 1 2 n 1 n ¼ j g  j a r u  j a r u þ j a r u  j a r u ¼ j g  a r u ; 2 1 1 n 1 2 2 n 2 1 1 n 1 1 2 n 2 2 2 n 2 and, similarly, g far ug¼ðj þ j Þg far ugþ j ðg ½ar uÞ n 1 2 n 2 n ¼ j g  j a r u  j a r u  j a r u þ j a r u ¼ð1 þ j Þg  a r u ; 1 1 1 n 1 2 2 n 2 2 1 n 1 2 2 n 2 2 1 n 1 so that far ug¼ a r u  j g ¼ a r u þ j g: ð5Þ n 1 n 1 2 2 n 2 1 Since ½u¼ 0, we may use (5) and Greens formula to obtain a ðu; /Þ¼ ðaru; r/Þ ðfar ug; /  / Þ h n X [X 1 2 C 1 2 ¼ðaru; r/Þ ða r u  j g; / Þ þða r u þ j g; / Þ 1 n 1 2 1 2 n 2 1 2 X [X C C 1 2 ¼ ðr ðaruÞ; /Þ þðj g; / Þ þðj g; / Þ ¼ðf ; /Þ þðj g; / Þ þðj g; / Þ ¼ Lð/Þ; 2 1 2 1 X [X 1 C 2 C X 1 C 2 C 1 2 which is the statement of the lemma. 4 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5541 An immediate consequence of Lemma 1 is the condition a ðu  U; /Þ¼ 0; 8/ 2 V ; ð6Þ h h which we will refer to as Galerkin orthogonality. A FE basis for V is easily obtained from a standard FE basis on the mesh by the introduction of new basis functions for the elements that intersect C. Thus, we replace each standard basis function livingon an element that intersects the interface by two new basis functions, namely its restrictions to X and X , re- 1 2 spectively. The collection of basis functions with support in X is then clearly a basis for V , and hence we obtain a basis for V by the identification w ¼ðwj ; wj Þ. If the interface coincides exactly with an element X X 1 2 edge, no new basis functions are introduced on these elements but the approximating functions may still be discontinuous over such an edge. As a consequence, there are six non-zero basis functions on each element that properly intersects C. Further implementation details are considered in Section 7. 4. Approximation property of V Recall that G denotes the set of elements that are intersected by the interface. We will use the following mesh dependent norms: 2 2 kvk :¼ h kvk ; 1=2;h;C K 0;C K2G 2 2 kvk :¼ h kvk ; 1=2;h;C 0;C K2G and 2 2 2 2 jjjvjjj :¼krvk þ kfr vgk þk½vk : 0;X [X 1=2;h;C 1=2;h;C 1 2 We note for future reference that ðu; vÞ 6 kvk kvk : ð7Þ C 1=2;h;C 1=2;h;C 1 2 To show that functions in V approximates functions v 2 H ðXÞ\ H ðX [ X Þ to the order h in the norm h 1 2 jjj jjj, we construct an interpolant of v by nodal interpolants of H -extensions of v and v as follows. 1 2 2 2 Choose extension operators E : H ðX Þ! H ðXÞ such that ðE wÞj ¼ w and i i i kE wk 6 Ckwk 8w 2 H ðX Þ; s ¼ 0; 1; 2: ð8Þ i i s;X s;X Let I be the standard nodal interpolation operator and define I v :¼ðI v ; I v Þ where I v :¼ðI E v Þj : ð9Þ 1 2 i h i i h h;1 h;2 h;i X The followingtheorem is valid. Theorem 2. Let I be an interpolation operator defined as in (9). Then 1 2 jjjv  I vjjj6 C hkvk ; 8v 2 H ðXÞ\ H ðX [ X Þ: A 1 2 h 2;X [X 0 1 2 In the proof of this result, we need to estimate the interpolation error at the interface. To that end, we shall use the followingvariant of a well known trace inequality on a reference element. The crucial fact is that the constant in this inequality is independent of the location of the interface relative to the mesh. 5 5542 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 e e Lemma 3. Map a triangle K 2 G onto the unit reference triangle K K by an affine map and denote by C C the h ~ K K corresponding image of C . Under assumptions A1–A3of Section 3 there exist a constant C, depending on C but independent of the mesh, such that kwk 6 Ckwk kwk ; 8w 2 H ðK KÞ: ð10Þ ~ ~ ~ 0;C C 0;K K 1;K K K K Proof. We start by showingthat 2 2 kwk 6 Cðkwk þkwk kwk Þ: ð11Þ ~ ~ ~ ~ 0;C C 0;C C 0;K K 1;K K ~ ~ K K K K;h Recall that C is the straight line connecting the points of intersection between C and the element K and K;h C ¼fðn; gÞ : 0 < n < jC j; g ¼ dðnÞg: K K;h Assume first that dðnÞ > 0. Since the curvature of the interface is bounded, jd ðnÞj6 CjC j. As the mesh K;h is non-degenerate this implies that on the reference element we may write, using again ðn; gÞ as local co- ordinates, e e C C ¼fðn; gÞ : 0 < n < jC C j; g ¼ d dðnÞg; ~ ~ K K K K;h ~ e e where jd d ðnÞj6 C. We now let D denote the domain bounded by C C ~ and C C ~ and note that by the diver- K K K K;h gence theorem, Z Z Z Z ow 2 1=2 2 2 2 0 2 w dndg ¼ divð0; w Þdndg ¼ w dn þ w ð1 þðd d Þ Þ ds: ð12Þ og ~ ~ D D C C C C ~ ~ K K;h K K As d d is bounded, 2 0 2 1=2 kwk 6 C w ð1 þðd d Þ Þ ds; 0;C C K K C C K K whence (11) follows from (12) usingCauchy–Schwarz  inequality. In a general case where d may switch sign, the same argument may be applied for each part between the e e intersections of C C and C C . ~ ~ K K K K;h It remains to show that the first term on the right in (11) is appropriately bounded. To that end we shall e e e map the triangular part K K and the quadrilateral part K K of K K onto new reference domains. We may assume t q e e that C C ~ intersects K K in ða; 0Þ and in ð0; bÞ, and, by symmetry, that 06 a6 b6 1. K K;h For a ¼ b ¼ 1, the desired trace inequality kwk 6 Ckwk kwk ð13Þ ~ ~ ~ 0;C C 0;K K 1;K K K K;h is valid. For 1=2 < a6 b < 1, we may map the triangular part K K onto the unit reference triangle by a linear map. By the bound from below on a and b, this map is bounded, uniformly in a and b, with uniformly bounded inverse, and hence (13) is valid also in this case. For 1=2 < a < b ¼ 1 the same argument holds, choosingthis time K K as the triangular part which contains the origin. Assume now that a6 1=2. Let ðxx^; yy^Þ¼ Mð~xx; yy~Þ¼ ðyy~; ð1  aÞ ð~xx þ yy~ 1ÞÞ: b e b Then the image K K ¼ MðK K Þ has its corners in ð0; 0Þ, ð1; 0Þ, ð0; 1Þ, PP ¼ðb; ð1  bÞ=ð1  aÞÞ, and there q q holds kwk 6 CðPP Þkwk kwk : ð14Þ ^ ^ ^ 0;C C 0;K K 1;K K ^ q q K K;h 6 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5543 b b An additional argument is needed to show uniformity in PP . Since 06 a6 1=2 and a6 b6 1, PP varies in the domain D D :¼f06 xx ^6 1=2; 1  xx ^6 yy^6 1g[f1=26 xx^6 1; 1  xx^6 yy^6 2ð1  xx ^Þg as a and b vary. Let kwk 0;C C K K;h F ðPP ; w w ^Þ¼ : kwk kwk ^ ^ 0;K Kq 1;K Kq b b b b b We will show that F ðPPÞ¼ sup F ðPP ; w w ^Þ is uniformly bounded. For points RR and SS in D D, assuming 1 ^ w2H ðK K Þ b b without restriction that F ðRRÞ P F ðSSÞ, we have for any w that b b b b b b b b F ðRRÞ F ðSSÞ¼ sup F ðRR; vv^Þ sup F ðSS; vv^Þ6 j sup F ðRR; vv^Þ F ðRR; w w ^Þj þ jF ðRR; w w ^Þ F ðSS; w w ^Þj ¼ I þ II: vv^ vv^ vv^ Given > 0 we may choose w w ^ such that I 6 =2. Note that F ðRR; vv^Þ is continuous for fixed vv^ since the only b b b dependence of RR lies in the domains of integration. We may thus take jRR  SSj small enough so that II 6 =2. b b b Hence F ðPP Þ is continuous on the compact set D D, and thus (14) holds uniformly in PP . Finally, since M is bounded, uniformly in a and b, with uniformly bounded inverse, (13) follows and the proof is complete. Proof of Theorem 2. Recall that K ¼ K \ X and let v ¼ E v denote the extension of v to X. By a standard i i i i i interpolation estimate we obtain krðv  I v Þk ¼ krðv  I v Þk 6 krðv  I v Þk 6 Chkv k : i i h h h;i 0;K i i 0;K i i 0;K i 2;K i i Summingover all triangles that intersect X , it follows by (8) that 2 2 2 2  2 krðv  I v Þk 6 Ch kv k 6 Ch kv k : ð15Þ i i i h;i 0;X i 2;K 2;X i i K\X 6¼; Next we consider the jumps on the interface. Since the mesh is non-degenerate, it follows from Lemma 3, scaled by the map from the reference triangle, that 2 2 2 1 2 1 h kwk 6 Cðh kwk þkwk Þ; 8w 2 H ðKÞ: K 0;C K 0;K 1;K Hence it follows, usingagain a standard interpolation estimate, that X X 2 2 2 1  1  1 h k½v  I vk 6Ch kv  I v k ¼ Ch kv  I v k i i h K h 0;C K h;i 0;C K i i 0;C K K K i i X X 2 2 2 2     2 6C ðh kv  I v k þkv  I v k Þ6Ch kv k : h h K i i 0;K i i 1;K K i 0;K i i Summingthe contributions from K 2 G , we get from (8) that k½v  I vk 6 Ch kv k 6 Chkvk : ð16Þ h 1=2;h;C i 2;[K2G 2;X [X h 1 2 i¼1 Finally, Lemma 3 applied to r w and scalinggives 2 2 2 2 2 h kr wk 6 Cðkwk þ h kwk Þ; 8w 2 H ðKÞ; K n 0;C 1;K K 2;K whence similar arguments as above yield kr ðv  I v Þk 6 Chkv k : ð17Þ n i i i h;i 1=2;h;C 2;X Since j < 1, the theorem now follows from (15)–(17). 7 5544 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5. A priori error estimates We will first show coercivity of the discrete form, for which purpose we will need the followinginverse inequality. Lemma 4. For / 2 V , the following inverse inequality holds: 2 2 kfr /gk 6 C kr/k : n I 1=2;h;C 0;X [X 1 2 Proof. Since / 2 V is linear on K , we have h i jC j jC jjK j 2 2 K 2 K i 2 2 2 2 h kj r / k 6 h j jC jjr/ j ¼ h j kr/ k ¼ h kr/ k 6 Ckr/ k : K i n i K K i K i K i i 0;C i i 0;K 0;K 0;K K i 2 i i jK j jKj In the last step above we have used that jC j6 h , jK j6 h , and, since the mesh is non-degenerate, K K i jKj P ch . The result follows by summation over the elements. Lemma 5. The discrete form a ð Þ is coercive on V , i.e., a ðv; vÞ P Cjjjvjjj 8v 2 V ; provided k is chosen sufficiently large. It is also continuous, i.e., a ðu; vÞ6 Cjjjujjj jjjvjjj 8u 2 V ; 8v 2 V : Proof. Continuity of the discrete form follows directly from the definitions. To prove coercivity, we use (7) to find that for any > 0 2 1=2 2 1=2 a ðv; vÞ¼ka rvk  2ð½v; far vgÞ þkk ½vk h n 0;X [X C 0;C 1 2 2 2 1=2 1=2 P ka rvk  2kfar vgk k½vk þkk ½vk 0;X [X 1=2;h;C 1=2;h;C 0;C 1 2 2 2 2 1=2 P ka rvk  kfar vgk þ k  k½vk : 0;X [X 1=2;h;C 0;C 1 2 K K2G It then follows from Lemma 4 that 1 1 2C max a I X 2 2 1=2 1=2 a ðv; vÞ P ka rvk þ  ka rvk 0;X [X 0;X [X 1 2 1 2 2 2 2 2 þ kfar vgk þ k  k½vk : 1=2;h;C 0;C K2G Taking  ¼ 4C max a, coercivity follows if kj ¼ ch where c > 4C max a. I X I X K K Theorem 6. Under assumptions A1–A3of Section 3, and for U solving (4) and u solving (1), the following a priori error estimates hold: jjju  Ujjj6 Chkuk ð18Þ 2;X [X 1 2 and ku  Uk 6 Ch kuk : ð19Þ 0;X 2;X [X 1 2 8 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5545 Proof. For any v 2 V , jjju  Ujjj6 jjju  vjjj þ jjjv  Ujjj. Further, by Lemma 5 and orthogonality, we have that jjjU  vjjj 6 Ca ðU  v; U  vÞ¼ Ca ðu  v; U  vÞ6 Cjjju  vjjj jjjU  vjjj; h h and it follows that jjju  Ujjj6 Cjjju  vjjj 8v 2 V : Taking v ¼ I u and invokingthe interpolation result of Theorem 2, (18) follows. For (19) we use a duality argument. Define z ¼ðz ; z Þ by 1 2 ða rz Þ¼ e in X ; i ¼ 1; 2; i i i i z ¼0on oX \ oX ; i i ð20Þ ½z¼0on C; ½ar z¼0on C; where e ¼ u  U . By Greens formula and (5) with u ¼ z, g ¼ 0, we have that i i i kek ¼ ðr ðarzÞ; eÞ ¼ðarz; reÞ ða r z ; e Þþða r z ; e Þ 1 n 1 1 2 n 2 2 0;X X [X X 1 2 ¼ðarz; reÞ ðfar zg; ½eÞ ¼ a ðz; eÞ n h X C since ½z¼ 0. Thus, usingthe symmetry of a ð Þ and applyingthe orthogonality relation (6) and Theorem 2, we find that kek ¼ a ðz  I z; eÞ6 Cjjjz  I zjjj jjjejjj6 Chkzk jjjejjj: ð21Þ h h h 0;X 2;X [X 1 2 Finally, by the elliptic regularity result (3), we have kzk 6 Ckek , whence the estimate (19) follows 2;X [X 0;X 1 2 from (21) and (18). 6. A posteriori error estimates In this Section, we prove a posteriori error estimates and formulate adaptive algorithms for the finite element method (4), followingBecker and Rannacher [3]. We will consider control of linear functionals jðeÞ of the error, and define the local and global estimators as 2 2 2 2 1 2 1=2 E ðUÞ¼ ðh kf þr ðarUÞk þ h k½a rU k þ h kg ½arUk þ h k½Uk Þ ð22Þ K K i i K K 0;K [K 0;oK 0;C K 0;C K K 1 2 and 1=2 EðUÞ¼ h E ðUÞ : ð23Þ K2T We then have the followinga posteriori error estimate. Theorem 7. For a continuous linear functional jð Þ on L ðXÞ, let J 2 L ðXÞ be defined by Riesz’ representation 2 2 theorem, i.e., jð Þ :¼ðJ; Þ . Then there is a positive constant C such that jðeÞ6 CEðUÞkJk : ð24Þ 0;X 9 5546 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 Proof. Let z be the solution to the problem ða rz Þ¼ J in X ; i ¼ 1; 2; i i i z ¼0on oX \ oX ; i i ð25Þ ½z¼0on C; ½ar z¼0on C: We first note that ðe; JÞ ¼ðare; rzÞ ða r z ; e Þ þða r z ; e Þ : 1 n 1 1 2 n 2 2 X X [X C C 1 2 Now, since ½u¼ 0, we see that, by (5), ða r z ; e Þ þða r z ; e Þ ¼ðU ; a r z Þ ðU ; a r z Þ ¼ð½U; far zgÞ: 1 n 1 1 2 n 2 1 1 1 n 1 2 2 n 2 n C C C C Thus, jðeÞ¼ðe; JÞ ¼ðare; rzÞ þð½U; far zgÞ: ð26Þ X X [X 1 2 Now, take Z 2 V . From Galerkin orthogonality we then have a ðe; ZÞ¼ 0, and since ½z¼ 0, ½u¼ 0 on the h h interface, we get 0 ¼ a ðe; ZÞ¼ ðare; rZÞ ð½e; far ZgÞ ðfar eg; ½ZÞ þðk½e; ½ZÞ h n n X [X C C C 1 2 ¼ðare; rZÞ þð½U; far ZgÞ þðfar eg; ½z  ZÞ þðk½U; ½z  ZÞ : ð27Þ n n X [X C C C 1 2 Denote by n the outward pointingunit normal to K. Subtracting(27) from (26) and integratingby parts we get jðeÞ¼ðare; rðz  ZÞÞ þð½U; far ðz  ZÞgÞ ðfar eg; ½z  ZÞ ðk½U; ½z  ZÞ n n X [X C C C 1 2 ¼ ðf r ðarUÞ; z  ZÞ þð½U; far ðz  ZÞgÞ ðfar eg; ½z  ZÞ ðk½U; ½z  ZÞ n n K [K C C C 1 2 K2T ð½an rU; z  ZÞ þða r e ; z  Z Þ ða r e ; z  Z Þ : K 1 n 1 1 1 2 n 2 2 2 oKnC C C K2T We now note that ða n re ; z  Z Þ ða n re ; z  Z Þ ðfar eg; ½z  ZÞ 1 1 1 1 2 2 2 2 n C C C ¼ða r e ; z  Z Þ ða r e ; z  Z Þ ðj a r e þ j a r e ; ½z  ZÞ 1 n 1 1 1 2 n 2 2 2 1 1 n 1 2 2 n 2 C C C ¼ðj ½ar e; z  Z Þ þðj ½ar e; z  Z Þ 2 n 1 1 1 n 2 2 C C ¼ðj ðg ½ar UÞ; z  Z Þ þðj ðg ½ar UÞ; z  Z Þ 2 n 1 1 1 n 2 2 C C ¼ ðð1  j Þðg ½ar UÞ; z  Z Þ j n j j j¼1 and thus we find that jðeÞ¼ ðf r ðarUÞ; z  ZÞ  ð½an rU; z  ZÞ K [K oKnC 1 2 K2T X X X þ ð½U; far ðz  ZÞg  k½z  ZÞ þ ðð1  j Þðg ½ar UÞ; z  Z Þ : ð28Þ n j n j j C C K K K2G K2G j¼1 h h 10 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5547 Further, by Cauchy–Schwarz inequality, and choosing kj ¼ ch with c > 4C max a , I X i K K 2 2 X X X X jðeÞ6 q x þ q x ; ð29Þ K;j S;j K;j S;j K2T j¼1 K2G j¼1 h h where q ¼ h kf þr ðarUÞk ; x ¼ h kz  Zk ; K K;1 K;1 0;K [K K 0;K [K 1 2 1 2 1=2 1=2 q ¼ h k½an rUk ; x ¼ h kz  Zk ; K K;2 K;2 K 0;oKnC K 0;oKnC 1=2 1=2 q ¼ h k½uk ; x ¼ h kfar ðz  ZÞg  c½z  Z=h k S;1 n K S;1 K 0;C K 0;C K K and 1=2 1=2 q ¼ h kg ½ar Uk ; x ¼ h max kz  Z k : n S;2 j j S;2 K 0;C K 0;C K K Note that the right hand side of (29) is a weighted sum of local residuals q , q . To conclude the proof K;i S;i we bound the weights in (29) by choosing Z ¼ I z, applyingthe interpolation estimates (15)–(17) and finally the regularity estimate kzk 6 CkJk ; 2;X [X 0;X 1 2 whence (24) follows. 7. Implementation and numerical examples Recall that a FE basis for V is obtained from a standard FE basis on the mesh by replacingeach standard basis function livingon an element that intersects the interface by two new basis functions, namely its restrictions to X and X , respectively. Both these new basis functions are, however, represented in the 1 2 implementation using the same nodes from the original triangulation. The points of intersection between the element edges and the interface are not used to represent the new basis function, and the geometry of the interface and the element parts does not come into play until integrating the terms in the bilinear form. One may thus alternatively think of the approximation as beingdefined on two overlappingmeshes formed by triangles covering the subdomains. In order to implement the discontinuous approximation, we first determined the set G of triangles 0 00 intersected by C. For each K 2 G , we assigned two identical copies K and K . We assumed that K, with nodes fi; j; kg was split by C into a triangle and a quadrilateral; the case of a split into two triangles was handled by creating a quadrilateral with one side of zero length. We assigned to K the triangular part and 00 00 0 0 to K the quadrilateral part. Thus, on K we created a new node i on the far side of C, and on K we created 0 0 0 two new nodes, k and j , on the other side of C, see Fig. 1. To ensure continuity across the edges of K and K (away from C), we also checked if the new nodes had already been created by the same process on the neighboring elements. After having completed this process, we thus had two independent meshes, one completely covering X , and the other completely covering X . The elements crossed by C had been dou- 1 2 bled, but coincided geometrically. To numerically evaluate a ð Þ and Lð Þ, we used the following strategy. The triangles that were not crossed by C were handled in the usual way. On the elements crossed by C, we used centroid quadrature to 11 5548 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 j′ k′ K′ K′′ i′ Fig. 1. The split of a triangle used in the implementation of the proposed method. evaluate all terms, both on the quadrilateral side and on the triangular side. On the interface, we used two- point Gaussian quadrature. All contributions were then assembled usingthe old and new nodes defined 0 0 0 by the splittingprocess. We emphasize that the new nodes i , j , and k are to be considered convenient support points for the definition of a continuous, piecewise linear, approximation rather than nodes in the standard finite element sense. The solution at these points is computed but is outside the domain of interest. 7.1. Example 1 We considered solutions to the ordinary differential equation d du du du i 1 2 a ¼ 1; ½uð1=2Þ ¼ 0; a ð1=2Þ¼ a ð1=2Þ: i 1 2 dx dx dx dx The domain is ð0; 1Þ, with an interface at x ¼ 1=2. While this is a one-dimensional problem, we solved it numerically in 2D on the domain ð0; 1Þð0; 1Þ, with zero Neumann boundary conditions at y ¼ 0and y ¼ 1. The equation has a closed-form solution, given by 2 2 ð3a þ a Þx x a  a þð3a þ a Þx x 1 2 2 1 1 2 u ðxÞ¼  ; u ðxÞ¼  : 1 2 2 2 4a þ 4a a 2a 4a þ 4a a 2a 1 2 1 1 2 2 1 2 We chose a ¼ 1=2, a ¼ 3 and performed a numerical convergence test for different approaches: a 1 2 standard unfitted FE method, a fitted standard FE method, and the proposed unfitted method. The con- vergence results are given in Fig. 2, and show the suboptimal behaviour of a standard unfitted method, and the small difference in computational error between the proposed method and a fitted method. We remark that an exact correspondence cannot be obtained, since the meshes will not be exactly the same for the fitted and unfitted methods. In Fig. 3 elevations of the solutions obtained with the different unfitted methods are shown on the final mesh consistingof 2048 elements. 12 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5549 -3.5 Nitsche -4 Unfitted standard FEM Fitted standard FEM -4.5 -5 -5.5 -6 -6.5 -7 -7.5 -8 -8.5 -5 -4 -3 -2 -1 0 1 Log(h) Fig. 2. L -norm convergence of different methods applied to the interface problem. Fig. 3. Elevation of the discontinuous approximation obtained from the proposed method (left) and the continuous standard FE approximation (right) on the final mesh. 7.2. Example 2 Here we considered a less trivial two-dimensional example (from [7]). The exact solution is given by ; if r6 r ; uðx; yÞ¼ 2 2 2 r r r 0 0 þ ; if r > r ; a a a 2 2 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 where r :¼ x þ y , and, on the domain ð0; 1Þð0; 1Þ, we chose r ¼ 3=4, a ¼ 1, and a ¼ 1000. The 0 1 2 boundary conditions were symmetry boundaries at x ¼ 0 and y ¼ 0 and Dirichlet boundary conditions correspondingto the exact solution at x ¼ 1 and y ¼ 1. The right-hand side yielding this exact solution is f ¼4. In Fig. 4, we give the elevation of the approximate solution on the last mesh in a sequence. The corresponding L - and maximum norm convergence is given in Fig. 5. We note in particular that the order of convergence in maximum norm seems to be approximately equal to the (second order) L -norm con- vergence, though we have no theoretical results to back this up. Log(L2 error) 5550 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.8 0.6 0 0.4 1 0.2 0.9 0.8 0.7 0.6 0.5 0.4 0 0.3 0.2 0.1 Fig. 4. Elevation of the approximate solution. L2 error -4.5 Maximum norm error -5 -5.5 -6 -6.5 -7 -7.5 -8 -5 -4 -3 -2 -1 0 Log(h) Fig. 5. L -norm and maximum norm convergence. 7.3. Example 3 The third example was solved on the domain ð0; 1Þð0; 1Þ, with zero Dirichlet boundary conditions. Centred in the domain is an ellipsoidal inclusion with conduction parameter a ¼ 1; outside of the ellipse min we set a ¼ 6. In Fig. 6 we show the first mesh and associated solution. We used an adaptive algorithm max correspondingto control of the L ðXÞ-error, but made no attempt to tune the interpolation constants or solve a dual problem. Thus, the example only gives an indication of how the adapted meshes appear in such a case. For each successive mesh we refined the third of the elements containingthe highest element contribution to the total error as defined by (29). In Fig. 7 we give the last adapted mesh in a sequence, together with an elevation of the corresponding solution. Log(error) A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 5551 Fig. 6. First mesh and elevation of the approximate solution. Fig. 7. Last adaptively refined mesh and elevation of the approximate solution. 8. Concluding remarks In this paper, we have introduced and analysed a new method for elliptic interface problems on unfitted meshes, i.e., meshes which are independent of the location of the interface. Unlike the standard unfitted finite element method, the proposed approach leads to optimal convergence rates. This have been shown in the model situation of piecewise linear approximations in two dimensions. In future work, we will address the general situation of higher order polynomial approximations in three dimensions. Other extensions under consideration include Stefan problems and fictitious domain type simulations. References [1] I. Babusska, The finite element method for elliptic equations with discontinuous coefficients, Computing5 (1970) 207–213. [2] J.W. Barrett, C.M. Elliott, Fitted and unfitted finite-element methods for elliptic equations with smooth interfaces, IMA J. Numer. Anal. 7 (1987) 283–300. [3] R. Becker, R. Rannacher, A feed-back approach to error control in finite element methods: basic analysis and examples, East- West J. Numer. Math. 4 (1996) 237–264. [4] J.H. Bramble, J.T. King, A finite element method for interface problems in domains with smooth boundaries and interfaces, Adv. Comput. Math. 6 (1996) 109–138. [5] Z. Chen, J. Zhou, Finite element methods and their convergence for elliptic and parabolic interface problems, Numer. Math. 79 (1998) 175–202. [6] M. Feistauer, V. Sobot ııkov aa, Finite element approximation of nonlinear problems with discontinuous coefficients, M2AN Math. Model. Numer. Anal. 24 (1990) 457–500. 15 5552 A. Hansbo, P. Hansbo / Comput. Methods Appl. Mech. Engrg. 191 (2002) 5537–5552 [7] Z. Li, T. Lin, X. Wu, New Cartesian grid methods for interface problems using the finite element formulation, NCSU Technical Report: CRSC-TR99-12, Center for Research in Scientific Computation, North Carolina State University, 1999. [8] R.J. MacKinnon, G.F. Carey, Treatment of material discontinuities in finite element computations, Int. J. Numer. Meth. Eng. 24 (1987) 393–417. [9] J. Nitsche, U Uber ein Variationsprinzip zur Loosungvon € Dirichlet–Problemen bei Verwendungvon Teilr €aaumen, die keinen Randbedingungen unterworfen sind, Abh. Math. Univ. Hamburg 36 (1971) 9–15. [10] A. Z Zen ııssek, The finite element method for nonlinear elliptic equations with discontinuous coefficients, Numer. Math. 58 (1990) 51–77.

Journal

Computer Methods in Applied Mechanics and Engineering – Unpaywall

Published: Nov 1, 2002

There are no references for this article.