Table Of ContentNumerical integration over implicitly defined domains for
higher order unfitted finite element methods∗
Maxim A. Olshanskii∗ Danil Safin∗
6
1
0
2
Abstract
n
a The paper studies several approaches to numerical integration over a domain defined
J
implicitlybyanindicatorfunctionsuchasthelevelsetfunction. Theintegrationmethods
2 are based on subdivision, moment–fitting, local quasi-parametrization and Monte-Carlo
2 techniques. As an application of these techniques, the paper addresses numerical solution
of elliptic PDEs posed on domains and manifolds defined implicitly. A higher order un-
]
A fitted finite element method (FEM) is assumed for the discretization. In such a method
the underlying mesh is not fitted to the geometry, and hence the errors of numerical inte-
N
grationovercurvilinearelementsaffecttheaccuracyofthefiniteelementsolutiontogether
. with approximation errors. The paper studies the numerical complexity of the integration
h
t procedures and the performance of unfitted FEMs which employ these tools.
a
m
1 Introduction
[
1
Numerical approaches for solving PDEs that integrate the underlying geometric information,
v
such as isogeometric analysis [14], are in the focus of research over the last decade. In the iso-
2
8 geometric analysis, the functions used for geometry representation are also employed to define
1 functionalspacesintheGalerkinmethod. Severalotherapproachestrytomakethegridgenera-
6 tion and geometry description independent. Unfitted finite element methods, such as immersed
0
boundary methods [16], extended FEM [4, 9], cut FEM [5], or trace FEM [19], use a sufficient
.
1 regular background grid, but account for geometric details by modifying the spaces of test and
0
trial functions or the right-hand side functional.
6
The accuracy of unfitted FEM depends on several factors. These are the approximation
1
: property of basic finite element space, the accuracy of underlying geometry recovering, and the
v
error introduced by numerical integration. The present paper first reviews the analysis of two
i
X higher order finite element methods: one is an unfitted FEM for the Neumann problem in a
r boundedcurvilineardomain, anotheroneisanarrow-bandFEMforanellipticPDEposedona
a closedsmoothmanifold. Inbothcases,thedomain(avolumeorasurface)isgivenimplicitlyby
a discrete level-set function. Practical implementation of these methods (as well as many other
unfitted FEM) leads to the following problem: Given a simplex K ∈ RN, a smooth function f
defined on K and a polynomial φ of degree q such that |∇φ | ≥ c > 0 on K, evaluate the
h h 0
integral
(cid:90)
I (φ ,f):= fdx, with Q={x∈K : φ (x)>0}. (1)
K h h
Q
1Department of Mathematics, University of Houston, Houston, Texas 77204-3008
(molshan,dksafin)@math.uh.edu
∗PartiallysupportedbyNSFthroughtheDivisionofMathematicalSciencesgrants1315993and1522252.
1
The paper focuses on several numerical approaches to problem (1) in R2. In the context of
numericalsolutionofPDEs, f istypicallyapolynomial; q+1istheorderofgeometryrecovery.
We note that for q = 1 one has to integrate f over a polygon (or polyhedron). Then for a
polynomialfunctionf anexactnumericalintegrationisstraightforwardthroughsubdividingthe
polygon(polyhedron)intoafinitenumberoftriangles(tetrahedra)andapplyingstandardGauss
quadratures on each triangle. However, for q > 1 the problem appears to be less trivial and
building an exact quadrature rule for I (φ ,f) does not look feasible. This can be realised by
K h
considering the simple 1D example with K =(0,1) and f ≡1. Solving (1) becomes equivalent
to finding the root of φ ∈(0,1). The latter problem is resolved (only) in radicals for 2≤q ≤4
h
and is well known to have no general algebraic solution for q >4 by the Abel theorem.
The problems of building numerical quadratures for the implicitly defined volume integrals
(cid:82)
(1) and the implicitly defined surface integrals fds, with S = {x ∈ K : φ (x) = 0} have
S h
beenalreadyaddressedintheliteratureandseveraltechniqueshavebeenappliedinthecontext
of XFEM and other unfitted FE methods. One straightforward approach consists of employing
thesmearedHeavisidefunctionH ,see,e.g.,[22]. Thenfortheregularizedproblem,oneapplies
ε
a standard Gaussian quadrature rule on the simplex K with weights ω and nodes x :
i i
(cid:90) (cid:90)
(cid:88)
I (φ ,f)= fH(φ )dx≈ fH (φ )dx≈ ω f(x )H (φ (x )).
K h h ε h i i ε h i
K K i
However, for a general superposition of K and the zero level set of φ , the smearing leads to
h
significant integration errors which are hard to control. Another numerical integration tech-
nique is based on an approximation of Q by elementary shapes. Sub-triangulations or quadtree
(octree) Cartesian meshes are commonly used for these purposes. On each elementary shape a
standard quadrature rule is applied. The sub-triangulation is often adaptively refined towards
the zero level of φ . The approach is popular in combination with higher order extended FEM
h
for problems with interfaces, see, e.g., [1, 17, 8], and the level-set method [15, 13]. Although
numerically stable, the numerical integration based on sub-partitioning may significantly in-
crease the computational complexity of a higher order finite element method, since the number
of function evaluations per K scales with h−p for some p > 0 depending on the order of the
FEM. In several recent papers [18, 21, 10] techniques for numerical integration over implicitly
defineddomainsweredevisedthathaveoptimalcomputationalcomplexity. Themoment–fitting
method from [18] uses polynomial divergence free basis of vector function to approximate the
integrandandfurtherreducethevolumeintegralstoboundaryintegralstofindtheweightsofa
quadratureformulabyaleast-squarefittingprocedure. Werecallthemoment–fittingmethodin
section3. ForthecasewhenK isahyper-rectangle,theapproachin[21]convertstheimplicitly
given geometry into the graph of an implicitly defined height function. The approach leads to
a recursive algorithm on the number of spatial dimensions which requires only one-dimensional
root finding and one-dimensional Gaussian quadrature. In [10], a zero level-set of φ is ap-
h
proximated by higher order interface elements. These elements are extended inside K so that
K is covered by regular and curvilinear simplexes. For curvilinear simplexes a mapping to the
reference simplex is constructed. Further standard Gauss quadratures are applied.
Inthispaper,wedevelopanapproachfor(1)basedonthelocalquasi-parametrizationofthe
zero level set of φ . Similar to [21] the zero level set is treated as a graph of an implicitly given
h
function. With the help of a 1D root finding procedure the integration over Q is reduced to the
integration over regular triangles and the recursive application of 1D Gauss quadratures. The
technique is also related to the method of local parametrization for higher order surface finite
element method of [11]. We compare the developed method with several other approaches for
thenumericalintegrationinthecontextofsolvingpartialdifferentialequationsindomainswith
curvilinear boundaries and over surfaces. For the comparison purpose we consider the method
2
of sub-triangulation for Q, the moment–fitting method, and the Monte-Carlo method.
The rest of the paper is organized as follows. We first recall unfitted finite element method
forsolvinganellipticPDEinadomainwithcurvilinearboundaryandanellipticPDEposedon
a surface. For the surface PDE we use the method from [20] of a regular extension to a narrow
band around the surface. The error analysis of these unfitted FE methods is also reviewed.
Furtherinsection3wediscussthemethodsfornumericalintegrationof (1), whichfurtherused
to build the FEM stiffness matrices. Section 4 collects the result of numerical experiments.
Section 5 concludes the paper with a few closing remarks.
2 Unfitted FEM
We assume that Ω is an open bounded subset in RN, N = 2,3, with a boundary Γ, which is
a connected C2 compact hypersurface in RN. In this paper we apply unfitted FE methods to
elliptic equations posed in Ω and on Γ. As model problems, let us consider the Poisson and the
Laplace–Beltrami problems:
−∆u+αu=f in Ω,
∂u (2)
=0 on Γ,
∂n
and
−∆ u+αu=g on Γ, (3)
Γ
with some strictly positive α∈L∞(Ω) or α∈L∞(Γ), respectively.
2.1 Preliminaries
Todefinefiniteelementmethods,weneedaformulationofthesurfacePDE(3)basedonnormal
extension to a narrow band. First, we introduce some preliminaries. Denote by Ω a domain
d
consisting of all points within a distance from Γ less than some d>0:
Ω ={x∈R3 : dist(x,Γ)<d}.
d
Let φ:Ω →R be the signed distance function, |φ(x)|:=dist(x,Γ) for all x∈Ω . The surface
d d
Γ is the zero level set of φ:
Γ={x∈R3 : φ(x)=0}.
We may assume φ<0 on the interior of Γ and φ>0 on the exterior. We define n(x):=∇φ(x)
for all x ∈ Ω . Thus, n is the normal vector on Γ, and |n(x)| = 1 for all x ∈ Ω . The Hessian
d d
of φ is denoted by H:
H(x)=D2φ(x)∈R3×3 for all x∈Ω .
d
Forx∈Γ,thenon-zeroeigenvaluesofH(x)aretheprincipalcurvatures. Hence,onecanchoose
such sufficiently small positive d = O(1) that I−φH is uniformly positive definite on Ω . For
d
x∈Ω denote by p(x) the closest point on Γ. Assume that d is sufficiently small such that the
d
decomposition x = p(x)+φ(x)n(x) is unique for all x ∈ Ω . For a function v on Γ we define
d
its extension to Ω :
d
ve(x):=v(p(x)) for all x∈Ω .
d
Thus, ve is the extension of v along normals on Γ.
We look for u solving the following elliptic problem
−divµ(I−φH)−2∇u+αeµu=geµ in Ω ,
d
∂u (4)
=0 on ∂Ω ,
∂n d
3
with µ=det(I−φH). The Neumann condition in (4) is the natural boundary condition. The
following results about the well-posedness of (4) and its relation to the surface equations (3)
have been proved in [20]:
(i) For g ∈ L2(Γ), the problem (4) has the unique weak solution u ∈ H1(Ω ), which satisfies
d
(cid:107)u(cid:107) ≤C(cid:107)ge(cid:107) , with a constant C depending only on α and Γ;
H1(Ωd) L2(Ωd)
(ii) For the solution u to (4) the trace function u| is an element of H1(Γ) and solves a weak
Γ
formulation of the surface equation (3).
(iii) The solution u to (4) satisfies (∇φ)·(∇u)=0. Using the notion of normal extension, this
can be written as u=(u| )e in Ω ;
Γ d
(iv) Additionally assume Γ ∈ C3, then u ∈ H2(Ω ) and (cid:107)u(cid:107) ≤ C(cid:107)ge(cid:107) , with a
d H2(Ωd) L2(Ωd)
constant C depending only on α, Γ and d.
2.2 FEM formulations
Let Ωbulk ⊂ RN, N = 2,3, be a polygonal (polyhedral) domain such that Ω ⊂ Ωbulk for
problem (2) and Ω ⊂Ωbulk for problem (3). Assume we are given a family {T } of regular
d h h>0
triangulations of Ωbulk such that max diam(T) ≤ h. For a triangle T denote by ρ(T) the
T∈Th
diameter of the inscribed circle. Denote
β = sup diam(T)/ inf ρ(T). (5)
T∈Th T∈Th
For the sake of presentation, we assume that triangulations of Ωbulk are quasi-uniform, i.e., β is
uniformly bounded in h.
It is computationally convenient not to align (not to fit) the mesh to Γ or ∂Ω . Thus, the
d
computational domain Ω approximates Ω or Ω and has a piecewise smooth boundary which
h d
is not fitted to the mesh T .
h
Let φ be a continuous piecewise, with respect to T , polynomial approximation of the
h h
surface distance function φ in the following sense:
(cid:107)φ−φ (cid:107) +h(cid:107)∇(φ−φ )(cid:107) ≤chq+1 (6)
h L∞(Ω) h L∞(Ω)
withsomeq ≥1(forproblem(2)agenericlevel-setfunctionφcanbeconsidered, notnecessary
a signed distance function). Then one defines
Ω ={x∈R3 : φ (x)<0} for problem (2)
h h
(7)
Ω ={x∈R3 : |φ (x)|<d } for problem (3),
h h h
with d = O(h), d ≤ d. Note, that in some applications the surface Γ may not be known
h h
explicitly and only a finite element approximation φ to the distance function φ is known.
h
Otherwise, one may set φ :=J (φ), where J is a suitable piecewise polynomial interpolation
h h h
operator. Estimate (6) is reasonable, if φ is a polynomial of degree q and φ∈Cq+1(Ω ). The
h d
latter is the case for Cq+1-smooth Γ.
The space of all continuous piecewise polynomial functions of a degree r ≥1 with respect to
T is our finite element space:
h
V :={v ∈C(T ) : v| ∈P (T) ∀T ∈T }. (8)
h h T r h
4
The finite element method for problem (2) reads: Find u ∈V satisfying
h h
(cid:90) (cid:90)
[∇u ·∇v +αeu v ] dx= fev dx ∀v ∈V , (9)
h h h h h h h
Ωh Ωh
where αe,fe are suitable extensions of α and f to Ωbulk. For problem (3), the finite element
method is based on the extended formulation (4) and consists of finding u ∈V that satisfies
h h
(cid:90) (cid:90)
(cid:2)(I−φ H )−2∇u ·∇v +αeu v (cid:3) µ dx= gev µ dx ∀v ∈V . (10)
h h h h h h h h h h h
Ωh Ωh
It should be clear that only those basis functions from V contribute to the finite element
h
formulations (9), (10) and are involved in computations that do not vanish everywhere on
Ω . Error estimates for the unfitted finite element methods (9), (10) depend on the geometry
h
approximation and the order of finite elements. Let Γ∈Cr+1 and assume u∈Hr+1(Ω) solves
the Neumann problem (2) and u ∈V solves (9). Then it holds
h h
(cid:107)ue−u (cid:107) +h(cid:107)ue−u (cid:107) ≤C(hr+1+hq+1), (11)
h L2(Ωh) h H1(Ωh)
where a constant C is independent of h, r is the degree of the finite element polynomials, and
q+1 the geometry approximation order defined in (6).
FortheNeumannproblemwithα=0andacompatibilityconditiononthedata,theestimate
in (11) is proved in [3] subject to an additional assumption on Ω . In that paper it is assumed
h
that ∂Ω and Γ match on the edges of K ∈T intersected by Γ. This is not necessary the case
h h
for Ω defined implicitly from the discrete level function φ . For implicitly defined domains,
h h
the convergence result in (11) follows from the more recent analysis in [12], where a suitable
mapping of Ω on Ω was constructed.
h
For the FE formulation (10) of the Laplace-Beltrami problem (3), approximations to φ and
H are required. If Γ is given explicitly, one can compute φ and H and set φ =φ, H =H and
h h
µ =det(I−φ H ) in (9). Otherwise, if the surface Γ is known approximately as, for example,
h h h
the zero level set of a finite element distance function φ , then, in general, φ (cid:54)= φ and one
h h
has to define a discrete Hessian H ≈ H and also set µ = det(I−φ H ). A discrete Hessian
h h h h
H can be obtained from φ by a recovery method, see, e.g., [2, 23]. Assume that some H is
h h h
provided and denote by p≥0 the approximation order for H in the (scaled) L2-norm:
h
|Ωh|−12(cid:107)H−Hh(cid:107)L2(Ωh) ≤chp, (12)
where |Ω | denotes the area (volume) of Ω .
h h
The convergence of the finite element method (9) is summarized in the following result
from [20]. Let Γ ∈ Cr+2, f ∈ L∞(Γ), and assume u ∈ W1,∞(Γ)∩Hr+1(Γ) solves the surface
problems (3) and u ∈V solves (9). Then it holds
h h
(cid:107)u−u (cid:107) ≤C(hr+hp+1+hq), (13)
h H1(Γ)
where a constant C is independent of h, and r ≥1, p≥0, q ≥1 are the finite elements, Hessian
recovery, and distance function approximation orders defined in (8), (6) and (12), respectively.
Numericalexperimentsin[20]showthat(cid:107)u−u (cid:107) typicallydemonstratesaoneorderhigher
h L2(Γ)
convergence rate than the H1(Γ) norm of the error.
Note that both error estimates (11) and (13) assume exact numerical integration. In the
introduction, we discussed that for q > 1 the exact numerical integration is not feasible. The
errorofthenumericalquadratureshouldbeconsistentwiththefiniteelementinterpolationand
geometric error to ensure that the FEM preserves the optimal accuracy.
5
Remark 1 For the purpose of improving algebraic properties of an unfitted finite element
method a stabilization procedure was suggested in [6] for elliptic equations posed in volumetric
domains and further extended to surface PDEs in [7]. The procedure consists of adding a
special term penalizing the jump of the solution gradient over the edges (faces) of the triangles
(tetrahedra) cut by ∂Ω .
h
Let Th := {K ∈ T : meas (K ∩∂Ω ) > 0} (Th is the set of all elements having non-
Γ h N−1 h Γ
empty intersection with the boundary of the numerical domain Ω ). By Fh denote the set of
h Γ
all edges in 2D or faces in 3D shared by any two elements from Th. Define the term
Γ
(cid:90)
(cid:88)
J(u,v)= σ n ·∇u n ·∇v .
F F F
F(cid:74) (cid:75)(cid:74) (cid:75)
F∈Fh
Γ
Here n ·∇u denotesthejumpofthenormalderivativeofuacrossF;σ ≥0arestabilization
F E
param(cid:74)eters(in(cid:75) ourcomputationswesetσ =1). Theedge-stabilizedtracefiniteelementreads:
F
Find u ∈V such that
h h
a (u ,v )+J(u ,v )=f (v ), (14)
h h h h h h h
for all v ∈ V . Here a (u ,v ) and f (v ) are the bilinear forms and the right-hand side
h h h h h h h
functional corresponding to the finite element methods (9) or (10).
ForP continuousbulkfiniteelementmethodsonquasi-uniformregulartetrahedralmeshes,
1
theoptimalordersofconvergencefor(14)wereprovedin[6]. Inournumericalstudieswetested
thestabilizedformulation(14)withhigherorderelements. Weobservedverysimilarconvergence
results for the formulations with and without J(u ,v ) term, including sub-optimal/irregular
h h
behaviour with moment-fitting quadratures and optimal with other integration techniques. For
thesystemsoflinearalgebraicequationsweuseexact‘backslash’solvesineithercase. Therefore,
we shall not report results for (14) in addition to finite element formulations (9) and (10).
3 Numerical integration
For a simplex K ⊂ RN, K ∈ T we are interested in computing integrals over Q = {x ∈
h
K : φ (x) > 0} with a certain accuracy O(hm), i.e. for a sufficiently smooth f we look for a
h
numerical quadrature I (φ ,f) such that
h,K h
|I (φ ,f)−I (φ ,f)|≤chm, m=min{q,r}+N, (15)
K h h,K h
with a constant c uniformly bounded over K ∈ T . We restrict ourselves with the two-
h
dimensional case, N =2.
As a pre-processing step we may compute a simple polygonal approximation Q to the
K
curvilinear integration domain Q. To find the polygonal subdomain Q , we invoke a root
K
finding procedure (several iterations of the secant method are used in our implementation) to
find all intersection points of the zero level set of φ (x) with the edges of K. Q is defined as a
h K
convex hull of these intersection points and the vertices of K lying in Ω . The integration of a
h
polynomial function f over Q can be done exactly using Gaussian quadratures on a (macro)
K
triangulation of QK. The integration over Q(cid:101) =Q(cid:52)QK has to be done approximately. We can
write
(cid:90) (cid:90) (cid:90)
fdx= fdx+ sign(φ )fdx.
h
Q QK Q(cid:101)
For the numerical integration over Q(cid:101) and Q, we consider several approaches. We start with the
most straightforward technique, the Monte-Carlo method.
6
Figure 1: Local subdivisions of cut triangles.
FortheMonte-Carlomethod weseedM (cid:29)1points{x }usingauniformrandomdistribution
i
over a narrow rectangular strip S containing Q(cid:101). Further compute
(cid:90) |S|(cid:88)M
sign(φ )fdx≈ g(x ), g(x)=f(x)sign(φ (x))χ (x). (16)
h M i h Q(cid:101)
Q(cid:101) i=1
(cid:114)
(cid:12) (cid:12)
Further we calculate the variance |S| (cid:12)((cid:80)M g(x ))2/M −(cid:80)M (g(x )2)(cid:12)/M. If the variance
(cid:12) i=1 i i=1 i (cid:12)
exceeds a predefined threshold ε = O(hm), then we increase the set of points used and update
the sum on the right-hand side of (16). Noting |S| = O(h3), a conservative estimate gives
M (cid:39)O(h3−2m).
IntheMonte-Carlomethodthenumberoffunctionevaluationsperagridelementintersected
by ∂Ω is too far from being optimal. A sub-triangulation method below allows to decrease the
h
number of function evaluations per cell and still constitutes a very robust approach.
In the sub-triangulation algorithm we construct a local triangulation of Q. This is done by
finding O(h3−m) points on the curvilinear boundary ∂Ω . In our implementation, these points
h
arefoundasintersectionsofauniformraycorntailoredtoabasispointon∂Q ,seeFigure1for
K
the example of how such local triangulations were constructed for cut triangles K and K (left
1 2
plot) and a cut triangle K (right plot) of a bulk triangulation (FE functions are integrated over
thegreenarea). Further,theintegraloveracutelementiscomputedasthesumofintegralsover
the resulting set of smaller triangles. The approach can be viewed as building a local piecewise
linear approximation of φ with some h(cid:48) =O(hm(cid:48)), with m(cid:48) =min{q,r}, i.e. h(cid:48) =O(h2) for P
h 2
elements and h(cid:48) =O(h3) for P elements.
3
The number of function evaluations per a triangle intersected by ∂Ω is O(h3−m), which is
h
better than with the Monte-Carlo method, but still sub-optimal. The first integration method
delivering optimal complexity we consider is the Moment-Fitting method from [18].
In the Moment-Fitting method, one first defines a set of points {x }, i = 1,...,M, for a
i
given cell K intersected by ∂Ω . The choice of the points can be done for a reference triangle
h
and is independent on how ∂Ω intersects K: {x } can be regularly spaced, come from a
h i
conventionalquadraturescheme,orevenrandomlydistributed. Theboundaryoftheintegration
domain Q consists of straight edges E and the curvilinear part I, cf. Figure 2. Let n ,
k k
n = (∇φ )/|∇φ | to be the unit outward normal vectors for each E and I respectively.
I h h k
7
Themoment–fittingmethodcalculatesquadra-
ture weights for the chosen points {x } by
i
finding the least-square solution to the sys-
tem
M (cid:90)
(cid:88)
ω g (x )= g dx (17)
i j i j
i=1 Q
for a given set of basis functions {g }, j =
j
1,...,K. For example, for P2 FEM we use
the basis
G ={g }={1, x, y, x2, xy, y2} (18)
j
The number M of points {x } is recom-
i
mended in [18] to exceed the number of basis
functions, i.e., K <M holds.
Figure 2: Integration domain Q for the The integrals on the right-hand side of
moment–fitting. (17) are evaluated approximately by resort-
ing to divergence-free basis functions. For
divergence-free basis functions, one applies the divergence theorem to reduce dimensions of
the integrals and relate integration over implicit surface to simple line integrals. In 2D, the
div-free basis complementing (18) is given by
(cid:26) 1 0 0 x y y2 2xy x2 0 (cid:27)
F ={f }=
j 0, 1, x, −y, 0, 0, −y2, −2xy, x2
Since ∂Q is a closed piecewise smooth curve, for f ∈F it holds
j
0=(cid:90) divf (x)dx=(cid:90) f ·n ds=(cid:90) f · ∇φh ds+(cid:88)(cid:90) f ·n ds
j j ∂Q j |∇φ | j k
Q ∂Q I h k Ek
Hence, we obtain
(cid:90) f · ∇φh ds=−(cid:88)(cid:90) f ·n ds
j |∇φ | j k
I h k Ek
Integrals over any interval E on the right side can be computed with a higher order Gauss
k
quadrature rule. This allows to build a quadrature for the numerical integration over implicitly
givencurvilinearedgeI basedontheinteriornodes{x }. Tothisend,onecalculatestheweights
i
{v } for the set of nodes {x } by solving the following system:
i i
(cid:88)M f (x )· ∇φh(xi) v =(cid:90) f · ∇φh ds, j =1,...,K.
j i |∇φ (x )| i j |∇φ |
i=1 h i I h
For area integration, take the second set of functions related to G as divh =2g :
j j
(cid:26) x x2/2 xy x3/3 x2y/2 xy2 (cid:27)
H={h }=
j y xy y2/2 x2y xy2/2 y3/3
The divergence theorem gives
2(cid:90) g dx=(cid:90) h ·n ds=(cid:90) h · ∇φh ds+(cid:88)(cid:90) h ·n ds
j j ∂Q j |∇φ | j k
Q ∂Q I h k Ek
8
pp
11
qˆ
1
n
φ(x)=c
0 q
1
p
1
Q(cid:101)
q2 qˆ2
p
2
p
2
Now the right-hand side values in (17) are (approximately) computed using the above identity
and the previously computed surface quadrature rule for the numerical integration over I:
2(cid:90) g dx≈(cid:88)M h (x )· ∇φh(xi) v +(cid:88)(cid:90) h ·n ds
j j i |∇φ (x )| i j k
Q i=1 h i k Ek
This algorithm can be extended to higher order quadratures by expanding the function sets G,
F and H; or to integration over a higher dimensional curvilinear simplex by adding a further
moment–fitting step.
The complexity of the moment–fitting integration is optimal, e.g. O(1) of function evalua-
tions per triangle. However, the weights computed by the fitting procedure are not necessarily
all non-negative. There is no formal prove of the resulting quadrature accuracy as in (15) or
the consistency order. Furthermore, in experiments we observe that the finite element methods
with stiffness matrices assembled using moment–fitting can be less stable compared to applying
other numerical integration techniques discussed here. Below we consider another integration
method of optimal complexity.
Consider the curvilinear remainder Q(cid:101) and denote by p1, p2 the intersection points of ∂Ωh
with Q . If there are more than 2 such points, then the calculations below should be repeated
K
for each of the simply-connected component of Q(cid:101). Choose points {qi} on (p1,p2) as a nodes
of a Gaussian quadrature with weights {ω }. Consider the normal vector n for the line passing
i
through p and p . For each point q , one finds the point qˆ on the boundary such that
1 2 i i
qˆ =q +αn, α∈R, and φ (qˆ )−c =0, where c is the φ -level value for ∂Ω , i.e. c =0 for
i i h i 0 0 h h 0
the FEM (9) and c =±d for the FEM (10). Secant method finds qˆ up to machine precision
0 h i
within a few steps. Further we employ the same 1D quadrature rule to place points {r } on
ij
each segment (qi, qˆi). The integral over Q(cid:101) is computed through
(cid:90)
(cid:88) (cid:88)
sign(φ )f(x)dx≈I (f)=|p −p | |q −qˆ |sign(φ (q ))ω f(r )ω . (19)
h h,Q(cid:101) 1 2 i i h i i ij j
Q(cid:101) i j
TheweightsforintegratingoverQmaybecombinedinω =|p −p ||q −qˆ |sign(φ (q ))ω ω
ij 1 2 i i h i i j
(cid:80) (cid:80)
to write the final quadrature formula f(r )ω . We remark that for problem (10), the
i j ij ij
9
factor sign(φ (q )) in (19) is replaced by sign(φ (q )+d ) or sign(d −φ (q )) depending on
h i h i h h h i
the level set of ∂Ω .
h
Similartothemoment–fittingmethod,thecomplexityofthenumericalintegrationisoptimal,
i.e. O(1) of function evaluations per triangle. All weights ω are positive, and the accuracy
ij
analysis of numerical integration of a smooth function over Ω, which invokes the constructed
quadrature for handling boundary terms, is straightforward and outlined below.
We estimate the error of integration of a sufficiently smooth f over Ω. The integration
uses a conventional quadrature scheme for interior cells of T and the quadrature (19) for the
h
curvilinear remainders of cut cells. This composed numerical integral is denoted by I (f).
h,K
Assumethe2DquadratureoverregulartriangleshasO(hm)accuracy. Bythetriangleinequality,
we have
(cid:12)(cid:90) (cid:12) (cid:12)(cid:90) (cid:90) (cid:12) (cid:12)(cid:90) (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)
(cid:12) fdx−Ih,Ω(f)(cid:12)≤(cid:12) fdx− fdx(cid:12)+(cid:12) fdx−Ih,Ω(f)(cid:12)
(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)
Ω Ω Ωh Ωh
We apply the co-area formula to estimate the first term
(cid:12)(cid:90) (cid:90) (cid:12) (cid:90) (cid:90)
(cid:12) (cid:12)
(cid:12) fdx− fdx(cid:12)≤ |f|dx≤(cid:107)f(cid:107)L∞ 1dx
(cid:12) (cid:12)
Ω Ωh
(Ω(cid:52)Ωh) (Ω(cid:52)Ωh)
+(cid:107)φh(cid:90)−φ(cid:107)L∞ (cid:90)
≤(cid:107)f(cid:107) |∇φ|dxdt≤Chq(cid:107)f(cid:107) .
L∞ L∞
−(cid:107)φh−φ(cid:107)L∞{φ=t}
For the second term we estimate
(cid:12)(cid:90) (cid:12) (cid:12)(cid:90) (cid:12)
(cid:12) (cid:12) (cid:88) (cid:12) (cid:12)
(cid:12) fdx−Ih,Ω(f)(cid:12)≤ (cid:12) fdx−Ih,K(f)(cid:12),
(cid:12) (cid:12) (cid:12) (cid:12)
Ωh K∈Th K
where I (f) is a quadrature we use to integrate f over K ∩Ω . Decomposing the mesh into
h,K h
interiorcellsTint andcutcellsTΓ (thoseintersectedby∂Ω )andapplyingtriangleinequalities,
h h h
we obtain:
(cid:12)(cid:90) (cid:12) (cid:12)(cid:90) (cid:12) (cid:12)(cid:90) (cid:12)
(cid:12) (cid:12) (cid:88) (cid:12) (cid:12) (cid:88) (cid:12) (cid:12)
(cid:12) fdx−Ih,Ω(f)(cid:12)≤ (cid:12) fdx−Ih,K(f)(cid:12)+ (cid:12) fdx−Ih,K(φh,f)(cid:12).
(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)
Ωh K∈Tint K K∈TΓ K∩Ωh
h h
The interior cells are integrated with the error Chm by a conventional method:
(cid:12)(cid:90) (cid:12)
(cid:88) (cid:12)(cid:12) fdx−Ih,K(f)(cid:12)(cid:12)≤Chm.
(cid:12) (cid:12)
K∈Tint K
h
For Q=K∩Ωh =QK ∪Q(cid:101)K, Q(cid:101)K =Q(cid:52)QK and polygonal QK defined earlier, we have
(cid:12)(cid:90) (cid:12) (cid:90)
(cid:88) (cid:12)(cid:12)(cid:12) fdx−Ih,K(φh,f)(cid:12)(cid:12)(cid:12)≤Chm+ (cid:88) | sign(φh)fdx−IQ(cid:101)K(f)|.
K∈TΓ K∩Ωh K∈TΓ Q(cid:101)K
h h
The estimate below assumes that a Gaussian 1D quadrature with P nodes is used in the
construction of (19). For each interval (p ,p ) we introduce the local orthogonal coordinate
1 2
system (s,t) = x(s)+tn, where s : (0,|p −p |) → (p ,p ) parameterizes the interval. The
2 1 1 2
10
