Efficient quadrature rules for finite element discretizations of nonlocal equations

In this paper we design efficient quadrature rules for finite element discretizations of nonlocal diffusion problems with compactly supported kernel functions. Two of the main challenges in nonlocal modeling and simulations are the prohibitive computational cost and the nontrivial implementation of discretization schemes, especially in three-dimensional settings. In this work we circumvent both challenges by introducing a parametrized mollifying function that improves the regularity of the integrand, utilizing an adaptive integration technique, and exploiting parallelization. We first show that the "mollified" solution converges to the exact one as the mollifying parameter vanishes, then we illustrate the consistency and accuracy of the proposed method on several two- and three-dimensional test cases. Furthermore, we demonstrate the good scaling properties of the parallel implementation of the adaptive algorithm and we compare the proposed method with recently developed techniques for efficient finite element assembly.

Authors

• 5 publications
• 6 publications
• 1 publication
• 20 publications
04/20/2016

An algorithm for the optimization of finite element integration loops

We present an algorithm for the optimization of a class of finite elemen...
12/10/2019

A Condensed Constrained Nonconforming Mortar-based Approach for Preconditioning Finite Element Discretization Problems

This paper presents and studies an approach for constructing auxiliary s...
02/05/2021

Error analysis of a decoupled finite element method for quad-curl problems

Finite element approximation to a decoupled formulation for the quad–cur...
05/04/2020

03/31/2021

Parallel two-scale finite element implementation of a system with varying microstructures

We propose a two-scale finite element method designed for heterogeneous ...
05/21/2020

A cookbook for finite element methods for nonlocal problems, including quadrature rules and approximate Euclidean balls

The implementation of finite element methods (FEMs) for nonlocal models ...
05/09/2018

The main computing tasks of a finite element code(FE) for solving partia...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Nonlocal equations have become the model of choice in applications where the global behavior of the system is affected by long-range forces at small scales. In particular, these equations are preferable to partial differential equations (PDEs) in presence of anomalous behavior, such as superdiffusion and subdiffusion, multiscale behavior, and discontinuities or irregularities in the solution that cannot be captured by classical models. For these reasons, nonlocal models are currently employed in several scientific and engineering applications including surface or subsurface transport

[7, 8, 24, 48, 49], fracture mechanics [31, 34, 51], turbulence [33, 42], image processing [9, 16, 28, 35] and stochastic processes [10, 18, 37, 39, 40].

The most general form [19] of a nonlocal operator for a scalar function is given by

 Lu(x)=2∫Rn(u(y)−u(x))γ(x,y)dy,

where , the kernel, is a compactly supported function over , the ball of radius centered at . We refer to as horizon or interaction radius; this quantity determines the extent of the nonlocal interactions and represents the length scale of the operator. The integral form allows one to catch long-range forces within the length scale and reduces the regularity requirements on the solutions. It also highlights the main difference between nonlocal models and PDEs, i.e., the fact that interactions can occur at distance, without contact. The kernel depends on the application and determines the regularity properties of the solutions; the choice of its parameters or functional form is among the most investigated open questions in the current nonlocal literature [42, 11, 20, 21, 30, 43, 44, 55, 57, 56]. In this work we limit ourselves to smooth integrable kernels since the treatment of more complex functions is not germane to the issues investigated in this paper, as clarified later on.

The integral nature of the operator poses several modeling and numerical challenges including the treatment of nonlocal interfaces [3, 12], the prescription of nonlocal boundary conditions [15, 23] and the design of efficient discretization schemes and numerical solvers [2, 13, 17, 22, 47, 50, 54]. In fact, the numerical solution of nonlocal equations becomes prohibitively expensive when the ratio between the interaction radius and the discretization size increases. Even though the nonlocal literature offers several examples of meshfree, particle-type discretizations [14, 45, 46, 50]

, in this paper we focus on finite element (FE) methods. This allows us to easily deal with nontrivial domains, achieve high-order accuracy, and use mesh adaptivity. Furthermore, the nonlocal vector calculus

[26] provides a means for a rigorous stability and convergence analysis of variational methods as it allows us to analyze nonlocal diffusion problems in a similar way as elliptic PDEs [27].

When cast in a variational form, the nonlocal problem associated with the operator results in a bilinear form characterized by the following double integral

 ∫Ω∪Γ∫(Ω∪Γ)∩Bδ(x)(u(y)−u(x))(φ(y)−φ(x))γ(x,y)dydx

where we explicitly reported the domain of integration in the inner integral. Here, and are the domain of interest and the corresponding “nonlocal boundary” and is an appropriate test function. Thus, the variational setting introduces further computational challenges. Not only do we have to numerically evaluate a double integral, but the integrand function is discontinuous, due to the compact support of and to the fact that, in FE settings, the tests functions are also compactly supported.

The paper by D’Elia et al.[22] thoroughly describes the challenges associated with nonlocal FE discretizations and proposes approximation techniques for efficient and accurate implementations. In particular, the authors introduce “approximate balls” that facilitate the assembly procedure by substituting the Euclidean ball with suitable polygonal approximations. They also suggest a set of quadrature rules for the outer and inner integration and analyze the convergence properties of the resulting scheme. With the same spirit, in this work, we propose an alternative way to efficiently evaluate the integral above by circumventing the issue of integrating a truncated function. The key idea of this paper is the introduction of a mollifier [41] to approximate the discontinuous kernel function; by doing so, the new, approximate, and parameterized kernel is a smooth function for which standard Gaussian quadrature rules can be employed over every element without compromising their accuracy. Additionally, we introduce adaptive quadrature rules for the numerical integration of the outer integral. In fact, contrary to intuition, sophisticated integration techniques for the outer integral are required in order to prevent the quadrature error from exceeding the FE one [22].

The main contributions of this paper are

• The introduction of a parametrized, smooth, approximate kernel, by means of a mollifier, that yields a smooth integrand over every element. This allows us to avoid the tedious and impractical task of determining the intersections of the ball with the elements, and hence represents a major advantage of our method in three-dimensional simulations.

• The design of adaptive quadrature rules for the outer integral and of a parallel algorithm for efficient simulations.

• The theoretical proof and numerical illustration of the convergence of the approximate, mollified solution to the analytic one as the mollified kernel approaches .

• The numerical illustration of the convergence of the mollified solution to the exact one as we refine the mesh and a numerical study of the dependence of the convergence behavior with respect to the parameters.

• The demonstration via two-dimensional and three-dimensional numerical tests of the scalability of our algorithm.

Outline of the paper

In the following section we define the notation that is used throughout the paper and recall important results on nonlocal calculus. In Section 3 we introduce the mollifier function and the associated approximate, parametrized weak form of the nonlocal diffusion problem. We also analyze the convergence of the solution of the latter to the original weak solution. In Section 4 we describe the nonlocal FE discretization, with special focus on the assembly procedure, and briefly recall its challenges. In Section 5 we introduce adaptive quadrature rules for the numerical integration of the outer integral. In Section 6 we illustrate our theoretical results via two- and three-dimensional tests. We also discuss a parallel implementation of the FE assembly procedure and show the corresponding scaling results. Finally, we summarize our contributions in Section 7.

2 Preliminaries

Let be open and bounded, . Given some , we define the interaction domain of as the set of all points not in that are within a distance from points in , i.e.

 Γ={y∈Rn∖Ω:|x−y|≤δfor somex∈Ω}, (1)

see Figure 1 for a graphical example in . Note that depends on even if it is not explicitly indicated. Let be an integrable nonnegative symmetric kernel that is also radialeeeFor a discussion on nonpositive kernels and nonsymmetric kernels, see [38] and [18], respectively., namely and . We also assume that has bounded support over the ball of radius centered at , i.e. . For a scalar function we define the nonlocal Laplacian as

 Lu(x)=2∫Rn(u(y)−u(x))γ(x,y)dy. (2)

The strong form of a nonlocal Poisson problem is then given by: for , and , find such that

 {−Lu(x)=f(x),x∈Ωu(x)=g(x),x∈Γ (3)

where the second condition in (3) is the nonlocal counterpart of a Dirichlet boundary condition for PDEs and it is referred to as Dirichlet volume constraintfffFor definition and analysis of Neumann volume constraints we refer to [25] and for its numerical treatment we refer to, e.g., [23].. Such condition is required [27] to guarantee the well-posedness of (3). The weak form of the Poisson problem is obtained by multiplying the first equation in (3) by a test function in and by applying the nonlocal first Green’s identity [26]. This yields

 0 =∫Ω(−Lu−f)φdx (4) =∬(Ω∪Γ)2(u(y)−u(x))(φ(y)−φ(x))γ(x,y)dydx−∫Ωf(x)φ(x)dx,

Then, the weak form of the nonlocal diffusion problem reads as follows. For and , find such that

 A(u,v)=F(v),∀v∈V0,subject tou=g inΓ, (5)

where

 A(u,φ) =∬(Ω∪Γ)2(u(y)−u(x))(φ(y)−φ(x))γ(x,y)dydx, (6) F(φ) =∫Ωf(x)φ(x)dx,

and where the function spaces are defined as

 V ={φ∈L2(Ω∪Γ):|||φ|||<∞%andv|Rn∖Ω=0} (7) V0 ={φ∈V:φ|Γ=0}, VΓ ={p:Γ→R:∃φ∈Vsuch that% φ|Γ=p}.

Here, the energy semi-norm is defined as

 |||φ|||2=A(φ,φ), (8)

the space is the dual space of and is a nonlocal trace space. Note that since the kernel is integrable and translation invariant, the energy semi-norm is a norm in the constrained space and satisfies a Poincaré inequality [27]. Furthermore, by construction, the bilinear form defines an inner product on and it is continuous and coercive with respect to the energy norm . Finally, the latter is equivalent to the norm; this allows us to establish an equivalence relationship between and . Together with the continuity of , these facts yield the well-posedness of the weak form (5) [27].

As commonly done in the PDE context, we recast the problem in by simply rewriting the solution as , where and is an extension of to zero into , known in the FE framework as a lifting function. Thus, equation (5) can be rewritten in terms of as follows

 ∬(Ω∪Γ)2(w(x)−w(y))(φ(x)−φ(y))γ(x,y)dydx=∫Ωf(x)φ(x)dx +∬(Ω∪Γ)2(˜g(x)−˜g(y))(φ(x)−φ(y))γ(x,y)dydx,∀φ∈V0, (9)

or, equivalently,

 A(w,φ)=˜F(φ)∀φ∈V0. (10)

The latter is useful for implementation purposes as it allows us to automatically take into account the presence of a non-homogeneous Dirichlet volume constraint.

3 Weak form approximation

We introduce a parametrized approximation of the bilinear form defined in (6) with the purpose of obtaining a weak problem that is computationally less challenging. Specifically, the approximated bilinear form is associated with a parametrized kernel function that is still integrable, radial, and compactly supported, but not discontinuous in . This fact makes the numerical integration of the inner integral in, e.g., (5), a much simpler task, compared to the case of discontinuous kernel functions.

For simplicity of exposition, we rewrite the “exact” kernel as

 γ(x,y)=Cδη(x,y)X(y∈Bδ(x)) (11)

where is a scaling constant that guarantees that the nonlocal operator associated with is such that as . Clearly, by definition, . Given , we approximate the bilinear form defined in (6) with the parametrized bilinear form obtained by replacing the kernel with

 γϵ(x,y)=Cδ,ϵη(x,y)μδ,ϵ(x,y), (12)

where is an appropriately scaled mollifier function. Inspired by the mollifier function introduced in [41], for we define as the following radial function

 μδ,ϵ(|x−y|)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩1 for 0≤|x−y|<δ−εξ((δ−ε)−|x−y|ε)%forδ−ε≤|x−y|≤δ+ε0 for |x−y|>δ+ε (13) ξ(r)=(128256+315256r−420256r3+378256r5−180256r7+35256r9),

where, for given , is such that the nonlocal operator , associated with , converges to as . The constant is also such that it converges to as . Furthermore, it follows from the definition of that

 limϵ→0μδ,ϵ(x,y)=X(y∈Bδ(x)),

which, together with the property of , implies that converges pointwise to as . Note that the support of the mollifier is bigger than the one of the original kernel function as it corresponds to . The parametrized bilinear form is therefore defined as follows

 Aϵ(u,φ)=∬(Ω∪Γ)2(u(y)−u(x))(φ(y)−φ(x))γϵ(x,y)dydx, (14)

and consequently the approximate weak formulation of the nonlocal volume-constrained problem now reads: find such that

 Aϵ(uϵ,φ)=F(φ)∀φ∈V0subject to uϵ=g on Γ. (15)

Note that the weak formulation above is defined over the same function spaces of (5). This is allowed because the parametrized kernel belongs to the same class of kernels as . As a consequence, problem (15) is also well-posed.

3.1 Convergence of the approximate weak solution

We find a bound for the energy norm of the difference between solutions of the weak form (5) and (15). First, we note that by subtracting (5) from (15) we obtain

 A(u,φ)=Aϵ(uϵ,φ)∀φ∈V0. (16)

Our ultimate goal is to find a bound for , or, equivalently, for , being and solutions of (5) and (15) respectively. We first consider a generic test function; equality (16) implies

 |A(u−uϵ,φ)| =|A(u,φ)−A(uϵ,φ)|=|Aϵ(uϵ,φ)−A(uϵ,φ)| =∣∣ ∣ ∣∣∬(Ω∪Γ)2(uϵ(x)−uϵ(y))(φ(x)−φ(y))(γϵ(x,y)−γ(x,y))dydx∣∣ ∣ ∣∣.

By expanding the product, we have

 ∣∣ ∣ ∣∣∬(Ω∪Γ)2(uϵ(x)−uϵ(y))(φ(x)−φ(y))(γϵ(x,y)−γ(x,y))dydx∣∣ ∣ ∣∣ ≤ ∬(Ω∪Γ)2|uϵ(x)φ(x)||γϵ(x,y)−γ(x,y)|dydx + ∬(Ω∪Γ)2|uϵ(y)φ(y)||γϵ(x,y)−γ(x,y)|dydx + ∬(Ω∪Γ)2|uϵ(x)φ(y)||γϵ(x,y)−γ(x,y)|dydx + ∬(Ω∪Γ)2|uϵ(y)φ(x)||γϵ(x,y)−γ(x,y)|dydx.

By switching the order of integration and renaming dummy variables in the second and fourth terms above, we obtain

 2∫Ω∪Γ|uϵ(x)φ(x)|∫Ω∪Γ|γϵ(x,y)−γ(x,y)|dydx (17) + 2∫Ω∪Γ|φ(x)|∫Ω∪Γ|uϵ(y)||γϵ(x,y)−γ(x,y)|dydx ≤ 2c1(ϵ)∫Ω∪Γ|uϵ(x)φ(x)|dx + 2∫Ω∪Γ|φ(x)|∥uϵ∥L2(Ω∪Γ)∥γϵ−γ∥L2(Ω∪Γ)dx ≤ 2c1(ϵ)∥uϵ∥L2(Ω∪Γ)∥φ∥L2(Ω∪Γ)+2c2(ϵ)|Ω∪Γ|12∥uϵ∥L2(Ω∪Γ)∥φ∥L2(Ω∪Γ),

where we used the Cauchy-Schwarz inequality for the outer and inner integral for the first and second term, respectively, and where

 c1(ϵ) =maxx∈Ω∪Γ∫Bδ+ϵ(x)∩(Ω∪Γ)|γϵ(x,y)−γ(x,y)|dy≤∫Bδ+ϵ(0)|γϵ(0,y)−γ(0,y)|dy, c22(ϵ) =maxx∈Ω∪Γ∫Bδ+ϵ(x)∩(Ω∪Γ)(γϵ(x,y)−γ(x,y))2dy≤∫Bδ+ϵ(0)(γϵ(0,y)−γ(0,y))2dy.

We recall that, by definition, pointwise as ; thus, as , for . Furthermore, thanks to the properties of and , the integrals above are well-defined.

To obtain the final estimate, we consider

and recall that for the kernels considered in this work the energy norm is equivalent to the norm. In particular there exists a positive constant such that . Thus, we have the following estimate

 ∥u−uϵ∥2L2(Ω∪Γ) ≤Ceq|||u−uϵ|||2 =CeqA(u−uϵ,u−uϵ) ≤Ceqk(ϵ)∥uϵ∥L2(Ω∪Γ)∥u−uϵ∥L2(Ω∪Γ),

where is obtained from the constants in (17). We finally conclude that

 ∥u−uϵ∥L2(Ω∪Γ)≤Ceqk(ϵ)∥uϵ∥L2(Ω∪Γ), (18)

where the constant is such that as .

4 Finite element formulation

In this section we introduce a FE discretization of problem (10), highlight the associated computational challenges, and describe how the formulation introduced in the previous section helps circumventing them. Let be a shape-regular triangulation of into finite elements ; the latter can either be triangles and/or quadrilaterals in two dimensions and tetrahedra and/or hexaheadra in three dimensionsgggSee [22] for a description of appropriate triangulation techniques for nonlocal problems.. The parameter represents the size of the triangulation and corresponds to the larger element diameter. Also, let be a finite dimensional subspace of of dimension , proportional to , and let be a basis for . In this work we consider Lagrange basis functions over the triangulation . Thus, we can write the FE solution of equation (10) as . By using this expression and , equation (10) reduces to the algebraic system

 AW=F, (19)

where

is the vector whose components are the degrees of freedom of the numerical solution

, is such that , and is the stiffness matrix with entries

 Aij=A(φi,φj)=∬(Ω∪Γ)2(φi(x)−φi(y))(φj(x)−φj(y))γ(x,y)dydx. (20)

4.1 Circumventing computational challenges

The computation of the entries of the stiffness matrix raises several diverse challenges. In this work we specifically focus on the challenges related to the presence of the indicator function in the definition of the kernel. Other challenges, such as the presence of singularities in fractional-type kernels or peridynamics kernels are not considered here. We point out that our method can be combined with any technique that takes into account the presence of the singularity. In fact, while the singularity is located at the center of the ball, the issues considered in this paper arise at the boundary. To make our description clear, we rewrite (20) by explicitly indicating the domain of integration, i.e.

 Aij =Cδ∫Ω∪Γ∫(Ω∪Γ)∩Bδ(x)(φi(x)−φi(y))(φj(x)−φj(y))η(x,y)dydx (21) =CδNL∑l=1NL∑k=1∫El∫Ek∩Bδ(x)(φi(x)−φi(y))(φj(x)−φj(y))η(x,y)dydx,

where we split the integrals over the elements with the purpose of using composite quadrature rules. In fact, global quadrature rules used, e.g., over the ball for the inner integration are not convenient due to the basis functions’ bounded support [22].

It is evident that the first challenge that one has to face is the integration over partial elements: when the element is not fully contained in the ball, standard quadrature rules such as Gauss quadrature rules defined over are not suitable due to the presence of the discontinuity induced by the indicator function. Thus, it is necessary to determine the intersection regions and define quadrature rules there. This task, while affordable in two dimensions, becomes extremely complex and impractical in three dimensions. Furthermore, when is a Euclidean ball, such regions are curved so that appropriate approximations or quadrature rules for curved domains must be taken into account [22]. A key observation is that these issues do not arise in case of smooth kernel functions, e.g. functions that do not abruptly jump to zero outside of , but that approach zero smoothly. This would allow the use of quadrature rules defined over the whole element , circumventing the issue of determining intersections or integrating over curved regions.

The parametrized kernel introduced in Section 3 is such that the transition to zero happens smoothly (as an example, for constant kernel functions , the kernel function is a piece-wise polynomial in . Thus, the inner integration can be performed over the whole element , using accurate enough quadrature rules, without worrying about the presence of a discontinuity. We then propose to solve the approximate, parametrized problem

 Aϵ(wh,ϵ,φi)=˜F(φi),∀i=1,…Nh, (22)

for which the entries of the stiffness matrix, that, with an abuse of notation, we still denote by , read

 Aij =Cδ,ϵNL∑l=1NL∑k=1∫El∫Ek(φi(x)−φi(y))(φj(x)−φj(y))η(x,y)μδ,ϵ(x,y)dydx. (23)

By avoiding the problem of determining intersecting elements, this approach makes three-dimensional implementation a much simpler task.

Remark 1

The convergence of the solution to the continuous solution depends on both the discretization parameter and the mollifying parameter . An adaptive quadrature procedure, introduced in the following section, will further contribute to the overall approximation error, as we discuss and illustrate in Section 6.

As already pointed out, the use of the mollifier, in place of the characteristic function, removes the difficulty of integrating discontinuous functions. The transition region of the mollifier has thickness

We recall that we denote by the stiffness matrix corresponding to the parametrized bilinear form . It is convenient to rewrite its entries as , where each term is given by

 A11ij =∫Ω∪Γ∫Ω∪Γγϵ(x,y)φi(x)φj(x)dydx, (24) A12ij =−∫Ω∪Γ∫Ω∪Γγϵ(x,y)φi(x)φj(y)dydx, (25) A21ij =−∫Ω∪Γ∫Ω∪Γγϵ(x,y)φi(y)φj(x)dydx, (26) A22ij =∫Ω∪Γ∫Ω∪Γγϵ(x,y)φi(y)φj(y)dydx. (27)

The following proposition allows us to express only as a sum of two of the terms above, as we show in Corollary 1.

Proposition 1

Let and , and let be a symmetric function, then

 ∫Ω∪Γ∫Ω∪Γg(x,y)f1(x)f2(y)dydx=∫Ω∪Γ∫Ω∪Γg(x,y)f1(y)f2(x)dydx. (28)
Proof 1
 ∫Ω∪Γ∫Ω∪Γg(x,y)f1(x)f2(y)dydx = ∫Ω∪Γ∫Ω∪Γg(x,y)f1(x)f2(y)dxdy reversing the order = ∫Ω∪Γ∫Ω∪Γg(x,y)f1(y)f2(x)dydx renaming variablesand using the symmetry of g
Corollary 1

The entries of the stiffness matrix satisfy the following equality

 Aij=2A11ij+2A12ij=2A21ij+2A22ij

.

Proof 2

The proof follows from the definition of and Proposition 1. Namely,

Let and denote sets of quadrature points and associated weights representing two different composite quadrature rules for the numerical integration over the region . To preserve the equality between and , and between and we use and to numerically evaluate the inner and outer integrals as follows:

 A11ij=∑q1∈Q1∑q2∈Q2γϵ(xq1,xq2)φi(xq1)φj(xq1)wq2wq1, (29) A12ij=∑q1∈Q1∑q2∈Q2γϵ(xq1,xq2)φi(xq1)φj(xq2)wq2wq1, (30) A21ij=∑q2∈Q2∑q1∈Q1γϵ(xq2,xq1)φi(xq1)φj(xq2)wq1wq2, (31) A22ij=∑q2∈Q2∑q1∈Q1γϵ(xq2,xq1)φi(xq1)φj(xq1)wq1wq2. (32)

According to Corollary 1, only two terms among the ones above need to be computed. We choose to compute and in (31) and (32). We adopt an adaptive scheme for the outer quadrature and Gaussian composite quadrature rules for the inner integration over the elements that intersect the ball. These choices are empirical, i.e. they have been guided by numerical experiments that showed that for and using an adaptive quadrature for the outer integral is more efficient than using it for the inner one.

For every we define

 Kl={m∈{1,…,NL}:∥x−y∥ℓ∞≥δ+ε,∀x∈El,∀y∈Em}, (33)

and , i.e. the complement of in . For and let denote the subset of composed of those quadrature points and weights obtained by only considering the quadrature points (and associated weights) that lie within element . Then, integrals (31) and (32) can be rewritten as

 A21ij=NL∑l=1∑q2l∈Ql2∑m∈Jl∑q1m∈Qm1γϵ(xq2l,xq1m)φmi(xq1m)φlj(xq2l)wq1mwq2l, (34) A22ij=NL∑l=1∑q2l∈Ql2∑m∈Jl∑q1m∈Qm1γϵ(xq2l,xq1m)φmi(xq1m)φlj(xq1m)wq1mwq2l. (35)

Note that in the equations above the terms in the sum are nonzero only when the support of the basis functions intersects the elements. An adaptive quadrature with midpoint refinement is adopted for . The pseudo-code that describes the adaptivity algorithm is reported in Algorithm 1: we employ recursive calls with input arguments , , , , and . In each call and are fixed parameters, and represent the minimum and the maximum level of refinement, with . is the current level of refinement. In the initial call and are the ones defined above and , while in the recursive calls these three arguments are subject to changes as described below.

• If , then is split in sub-elements using the midpoint rule and for each of them the adaptive integration function is called again increasing by one and using the same index set .

• If integration is performed on for the outer integral and on each element indexed by for the inner integral, without any further refinement for . The numerical integration is performed using standard Gauss Legendre quadrature rules both for the outer and inner integrals.

• If , from the index set two new index sets and are extracted, for which is either integrated or further refined. For any in , the maximum distance from to is computed. If this distance is less than , then is added to , otherwise the minimum distance from to is computed. If this distance is less than then is added to the index set . If is non-empty, integration is performed on for the outer integral and on each element indexed by for the inner integral. If is non-empty, then is split in sub-elements and for each of them the adaptive integration function is called again increasing by one and using as index set.

Algorithm 1
function Adaptive Integration(, , , , )
if  then
split into new sub-elements
for   do
Adaptive Integration(, , , , )
end for
else if  then
Integration(, )
else

for all  do
if  then
else if  then
end if
end for
if  then
Integration(
end if
if  then
split into new sub-elements
for   do
Adaptive Integration(, , , , )
end for
end if
end if
end function

5.1 Approximate maximum and minimum distances between elements

Evaluating the exact distances between two elements can be computationally expensive, especially for unstructured three-dimensional meshes. Hence, in practice, we use conservative distances that are simple to compute in place of the maximum and the minimum. Namely, we first loop over the nodes of each element to find the minimum and maximum coordinates in each dimension, denote them by , and by , . Here and are the vectors containing the minimum and maximum coordinates of the bounding box containing . Similarly, and are the vectors containing the minimum and maximum coordinates of the bounding box containing . For each dimension evaluate the two quantities