High order unfitted finite element methods on level set domains using isoparametric mappings

https://doi.org/10.1016/j.cma.2015.12.005Get rights and content

Abstract

We introduce a new class of unfitted finite element methods with high order accurate numerical integration over curved surfaces and volumes which are only implicitly defined by level set functions. An unfitted finite element method which is suitable for the case of piecewise planar interfaces is combined with a parametric mapping of the underlying mesh resulting in an isoparametric unfitted finite element method. The parametric mapping is constructed in a way such that the quality of the piecewise planar interface reconstruction is significantly improved allowing for high order accurate computations of (unfitted) domain and surface integrals. We present the method, discuss implementational aspects and present numerical examples which demonstrate the quality and potential of this method.

Introduction

In the recent years unfitted finite element methods have drawn more and more attention. The possibility to handle complex geometries without the need for complex and time consuming mesh generation is very appealing. The methodology of unfitted finite element methods, i.e. methods which are able to cope with interfaces or boundaries which are not aligned with the grid, has been investigated for different problems. For boundary value problems with non-aligned boundaries methods such as penalty methods  [1], [2], the fictitious domain method  [3], [4], the immersed boundary method  [5], and other unfitted finite element methods  [6] have been developed. For unfitted interface problems extended finite element methods (XFEM) have been developed in (among others)  [7], [8], [9], [10], [11]. Also partial differential equations on surfaces which are not meshed have been considered using the trace finite element method (TraceFEM)  [12].

In the community of unfitted finite element methods mostly piecewise linear discretizations are considered. One major issue in the design and realization of high order methods is the problem of numerical integration on domains which are only prescribed implicitly, for instance as the zero level of a scalar function, the so-called level set function. The use of standard integration rules (which ignore the cut position on cut elements) is not a good option as the integrand does not provide the necessary smoothness.

The purpose of this paper is the presentation of a new approach which allows for high order numerical integration on domains prescribed by level set functions. The approach is new, robust and fairly simple to implement. The method is geometry-based and can be applied to unfitted interface or boundary value problems as well as to partial differential equations on surfaces.

We briefly discuss the state of the art in the literature to put the new approach into context. One approach that is often used in order to realize numerical integration on implicit domains consists of two steps: the approximation of the interface with a piecewise planar interface and a tessellation algorithm to divide the subdomains of a cut element into simple geometries, on which standard quadrature rules can be applied. A famous example for quadrilaterals and cubes is the marching cube algorithm  [13]. For simplices a detailed discussion of this approach can be found in (among others)  [14, Chapter 5], [15] for triangles and tetrahedra and in  [16], [17] also for 4-prisms and pentatopes (4-simplices). Many simulation codes which apply unfitted finite element discretizations, e.g.  [18], [19], [20], [21], [22] make use of such strategies. In order to increase the accuracy of this integration one often combines this approach with subdivisions and adaptivity. Especially on octree-based meshes this can be done very efficiently  [23]. However, this tessellation approach is only second order accurate (w.r.t. the finest subtriangulation) which complicates the realization of high order methods.

One approach to solve this problem is the tailoring of quadrature points and weights which provide high order accurate integration rules for implicit domains. Such a construction of points and weights is based on the idea of fitting integral moments, cf.  [24], [25]. Although this results in accurate integration rules it has two shortcomings. First, the computation of quadrature rules is fairly involved. This aspect is typically outwayed by the resulting high accuracy. Secondly, the major problem, quadrature weights are in general not positive. This can lead to stability problems as positivity of mass or stiffness matrices in finite element formulations can no longer be guaranteed.

For special cases satisfactory answers to the question of high order numerical integration strategies have been found which allow for an implementation of high order unfitted finite element methods. We mention a few interesting approaches. For quadrilaterals and hexahedra in  [26] a numerical integration algorithm is presented which can achieve arbitrary high order accuracy and guarantee positivity of integration weights at the same time. The approach is based on the idea of interpreting the interface as a piecewise graph over hyperplanes. In  [27] an unfitted boundary value problem is considered. Instead of aiming at a high order accurate geometrical approximation of the boundary a correction in the imposition of boundary values is applied in order to recover a high order method. For the discretization of partial differential equations on surfaces, in  [28] a parametric mapping of the interface from a piecewise planar representation to a high order representation based on quasi-normal fields is applied. Although the approach cannot be carried over straight-forwardly to the case where also the integration on sub-domains is important, that approach has important similarities to the approach presented in this paper.

In the literature of extended finite element methods (XFEM), there exist other approaches which are based on a tessellation approach and aim at improving the resulting piecewise planar approximations by means of a parametric mapping of the sub-triangulation   [29], [30], [31]. These approaches are typically technically involved, especially in more than two dimensions. Moreover, ensuring robustness of these approaches is difficult.

The approach presented in this paper is similar to the approaches in  [28], [29], [30] in that it is also based on a piecewise planar interface which is significantly improved using a parametric mapping. The important difference is, that we consider a parametric mapping of the underlying mesh rather than the sub-triangulation or the interface. According to the mapping of the mesh the use of isoparametric finite element is natural.

At this point, we would also like to mention the paper  [32] where the construction of a mesh deformation, which is also used in combination with isoparametric finite elements, is specifically designed to align a given mesh to a given interface position. The goal in that paper is to obtain a mesh that is conforming to a given interface without changing the mesh topology while keeping the quality of the mesh. Our approach is in a similar virtue.

In this paper we present a strategy for the numerical approximation of domains implicitly defined by an approximate signed distance function ϕ, the level set function. The approach is based on the assumption that a numerical strategy to deal with interfaces stemming from a piecewise linear level set function exactly is available.

In the context of unfitted finite element methods integrals of the form Sfdx have to be computed for S being an interface or a domain which is only implicitly defined through the level set function ϕ. Consider the case S=Ωi=Ω{ϕ<0} for a polygonal domain Ω which has a suitable triangulation T. As an exact evaluation of the integral Ωifdx is in general practically not feasible we consider the approximation of Ωi with the domain Ωi,hΩ{Ihϕ<0} where Ih is the continuous piecewise linear interpolation of ϕ. To Ωi,h an explicit representation can be found on which quadrature rules can be applied element by element: ΩifdxΩi,hfdxTTiωif(xi). Here ωi and xi are the integration weights and points which depend on the cut configuration on T. The accuracy of this approach is limited by the approximation quality of Ωi,h which is only second order accurate. By a transformation of the mesh which is represented by an explicit mapping Ψh, parametrized by a finite element function, the piecewise planar interface is mapped approximately onto the zero level set of a corresponding high order accurate level set function, cf. Fig. 1.1.

According to this mapping we have that Ψh(Ωi,h) is a high order accurate approximation to Ωi which, by construction, still has an explicit representation. The integral can hence be approximated as ΩifdxΨh(Ωi,h)fdxTTiωi|det(Ψh(xi))|f(Ψh(xi)), where the integration weights and points are the same as in (1.1). In contrast to (1.1) the accuracy of the quadrature in (1.2) is no longer bounded to second order accuracy but essentially depends on Ψh. The finite element spaces (used in a discretization of an (unfitted) PDE problem) have to be adapted correspondingly which renders a resulting method an isoparametric finite element method: Let Wh be the finite element space corresponding to the piecewise planar interface approximation with Γ1{Ihϕ=0}. Then, after transformation with Ψh, the appropriate isoparametric finite element space is Whk{φhΨh1|φhWh}. Note that this implies that shape functions are not necessarily polynomials on the mapped domain.

From a computational point of view only two things change compared to the case of a piecewise planar interface. First, a suitable (approximated) mapping Ψh has to be constructed. A major contribution of this paper is the discussion of this. Secondly, the deformation has to be considered in the definition of the shape functions and the quadrature of deformed elements. The treatment of the latter aspect is well-known and is not discussed here.

Owing to the combination of a tessellation algorithm for a piecewise planar interface and a parametric mapping, we obtain the four most important features of this approach:

  • An explicit high order accurate geometrical approximation of the exact interface.

  • Guaranteed positiveness of quadrature weights on the interface and in the sub-domains.

  • Unchanged cut topology compared to the piecewise planar case.

  • Easy integration into existing codes as the approach builds on a piecewise planar interface.

Together, this renders the method highly accurate, robust, efficient and fairly simple to implement. We again note, that the improved geometry approximation is based on a parametric mapping of the initial mesh, s.t. a combination with isoparametric finite elements is crucial.

The main purpose of this study is the presentation of the new approach which allows for an explicit high order geometry approximation of domains which are implicitly defined as level sets. We focus on the discussion of an efficient construction of a suitable mesh deformation and related implementational aspects. Although this discussion goes beyond simple concepts in that it addresses important details of the method, we do not aim at a thorough error analysis, yet. This is topic of ongoing research and will be published in a different paper.

The paper is structured in the following way. The starting point for the method is the ability to deal with interfaces which are piecewise planar. As we mainly consider simplicial meshes this piecewise planar interface typically corresponds to a (continuous) piecewise linear level set function. We consider a corresponding decomposition strategy as given and do not discuss it here. Instead, we refer to the literature instead, cf.  [14, Chapter 5] or  [17, Chapter 4].

The crucial component of the new approach is the construction of an explicit mapping Ψh which is suitable for implementation. The construction of Ψh is discussed in Section  2 and consists of several steps. We first characterize the desired properties of the transformation Ψh for a general case. Afterwards we restrict to an important special case with simplicial elements and a given level set function which is an approximate signed distance function. For this case we present an explicit construction of a suitable mapping.

Numerical examples demonstrating the quality of the explicit geometrical representation obtained by this new approach are presented in Sections  3 Numerical examples I: Geometry approximation, 4 Numerical examples II: a high order unfitted isoparametric finite element method for an interface problem. While the examples in Section  3 focus on the accuracy of the geometry approximation, in Section  4 a high order unfitted (isoparametric) finite element formulation for an interface problem is considered which proves the practical use of the method.

Section snippets

Construction of the mesh transformation Ψh

We start with the formulation of properties that we demand from a suitable transformation Ψh. Let Γ1 be a piecewise planar interface and [Vhk]d the space of continuous vector-valued piecewise polynomial functions of degree k. We seek for a transformation Ψh[Vhk]d, s.t. Ψh=argminΨhVdist(Ψh(Γ1),Γ)  or at least  dist(Ψh(Γ1),Γ)chk+1 with V[Vhk]d the subset of admissible transformations which fulfill the following constraints:

  • Homeomorphy: Ψh,Ψh1[C0(Ω)]d.

  • Shape regularity: Consider an

Preliminaries

In the following we consider complex geometries described by level set functions and their approximation with the approach proposed before. For the barrier parameter we choose γ=0.1 (independent of the degree k) and consider different polynomial degrees k. We note that the case k=1 corresponds to Ψh(x)=x, the case of a piecewise linear interface. To solve the one-dimensional nonlinear problem in (2.6b) we choose a (relative) tolerance of 10−14 which resulted in only 23 iterations for each

Elliptic interface problem: The disc

Problem description. We consider the problem from  [17, section 2.5.1.4], where an elliptic interface problem of the following form is considered. div(αu)=finΩi,i=1,2,αun=0onΓ,βu=0onΓ,u=gDonΩ. Here αi, βi, i=1,2 are domain-wise constants, (α1,α2)=(2,1), (β1,β2)=(1,3/2). We consider a circular interface Γ with radius R=0.6 as interface Γ in the square domain [1,1]2. The data gD and f is chosen such that the solution to (4.1) is u(x)={α2U(r(x))+β2,xΩ1α1U(r(x))+β1,xΩ2  with  U(r)=cos(πr

Conclusion

A new unfitted finite element method with a high order geometrical approximation of level set domains is presented and discussed in detail. The method is efficient and easy to implement. Numerical examples reveal its potential in handling complex geometries robust and highly accurate.

As an example for the potential of the method, the method has been combined with a rather standard Nitsche-XFEM discretization for an unfitted interface problem and optimal order of convergence has been observed.

Acknowledgments

The author would like to express his appreciation to Arnold Reusken for his valuable and constructive feedback on a first draft of this manuscript. The author also greatly appreciates the fruitful discussion with Joachim Schöberl on the topic of mesh transformations and shape regularity in the context of this study.

References (36)

  • J.W. Barrett et al.

    Finite element approximation of the Dirichlet problem using the boundary penalty method

    Numer. Math.

    (1986)
  • P. Bastian et al.

    An unfitted finite element method using discontinuous Galerkin

    Internat. J. Numer. Methods Engrg.

    (2009)
  • T.-P. Fries et al.

    The extended/generalized finite element method: an overview of the method and its applications

    Internat. J. Numer. Methods Engrg.

    (2010)
  • S. Groß et al.

    A finite element based level set method for two-phase incompressible flows

    Comput. Vis. Sci.

    (2006)
  • R. Massjung

    An unfitted discontinuous Galerkin method applied to elliptic interface problems

    SIAM J. Numer. Anal.

    (2012)
  • M.A. Olshanskii et al.

    A finite element method for elliptic equations on surfaces

    SIAM J. Numer. Anal.

    (2009)
  • W.E. Lorensen et al.

    Marching cubes: A high resolution 3d surface construction algorithm

  • T.A. Nærland

    Geometry decomposition algorithms for the Nitsche method on unfitted geometries

    (2014)
  • Cited by (166)

    • A high order geometry conforming immersed finite element for elliptic interface problems

      2024, Computer Methods in Applied Mechanics and Engineering
    • A CutFEM method for phase change problems with natural convection

      2024, Computer Methods in Applied Mechanics and Engineering
    • Fast immersed boundary method based on weighted quadrature

      2023, Computer Methods in Applied Mechanics and Engineering
    View all citing articles on Scopus
    View full text