Adaptively weighted numerical integration over arbitrary domains

Computers and Mathematics with Applications 67 (2014) 1682–1702

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Adaptively weighted numerical integration over arbitrary domains Vaidyanathan Thiagarajan ∗ , Vadim Shapiro Spatial Automation Laboratory, University of Wisconsin - Madison, United States

article

info

Article history: Received 4 October 2013 Received in revised form 10 January 2014 Accepted 1 March 2014 Available online 13 April 2014 Keywords: Numerical integration Quadrature Cubature Moment fitting equations Shape sensitivity analysis Meshfree analysis

abstract In adaptively weighted numerical integration, for a given set of quadrature nodes, order and domain of integration, the quadrature weights are obtained by solving a system of suitable moment fitting equations in least square sense. The moments in the moment equations are approximated over a simplified domain that is homeomorphic to the original domain, and then are corrected for the deviation from the original domain using shape sensitivity analysis. Using divergence theorem, the moments reduce to integrals over the boundary of the simplified domain. The proposed approach supports accurate and efficient computation of quadrature weights for integration of a priori unknown functions over arbitrary 2D and 3D solid domains. Experimental results (2D) indicate that adaptively weighted integration compares favorably with more traditional approaches. Because the adaptively weighted integration avoids excessive domain subdivision, it is useful in many applications and meshfree analysis in particular. Published by Elsevier Ltd.

1. Introduction 1.1. Motivation Numerical integration is a fundamental computation that is central to many problems in scientific and engineering computing, in large part because they require numerical integration of partial differential equations. For example, in structural analysis, numerical integration is used to assemble stiffness and mass matrices, as well as computing integral properties of physical fields such as average stresses and displacements. Formally, the problem is to numerically evaluate the integral:

 Ω

f (x) dΩ

(1)

where f : Ω → R is an integrable function defined over an arbitrary domain1 Ω ⊂ Rd (d = 2 or 3) that is typically represented as either 2D or 3D solid model [1]. For the purposes of this paper, we will assume that the integrand function f is continuous.2 If function f is known a priori, then the domain integral (1) can usually be reduced to a boundary integral

∗

Corresponding author. Tel.: +1 518 698 8380. E-mail addresses: [email protected] (V. Thiagarajan), [email protected] (V. Shapiro).

1 By ‘arbitrary domain’ we mean any 2D or 3D closed regular set [1] whose boundary is an orientable manifold. 2 For dealing with discontinuous functions refer to [2,3] and references therein. http://dx.doi.org/10.1016/j.camwa.2014.03.001 0898-1221/Published by Elsevier Ltd.

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

1683

Fig. 1. Illustration of excessive boundary cell fragmentation in quadtree decomposition.

via application of the divergence theorem [4–8]. For the rest of the paper, we will assume that f is an arbitrary (continuous) function that is not known a priori. Numerical integration over simple regular domains (interval, simplex, cube, and other ‘‘cells’’) is a classical widely studied topic that is relatively well understood [9,10]. Integration over more complex domains, such as those arising in most engineering applications, usually require that domain Ω is approximated by a union of simple cells, so that the integration can be carried out in a cell by cell fashion. In order to avoid integration errors, the union of cells must closely approximate the domain—either by cells that conform to the boundary ∂ Ω or by non-conforming cells that are recursively decomposed to cover the original domain. The first approach leads to a difficult (theoretically and computationally) problem of meshing. The latter case, illustrated in Fig. 1, leads to excessive fragmentation near domain’s boundary, resulting in significant computational cost and loss of accuracy. The goal of this paper is to devise an accurate and efficient method (called Adaptively Weighted Numerical Integration method) that avoids domain meshing and excessive fragmentation. 1.2. Overview and outline Given an arbitrary domain Ω with a set of appropriately chosen quadrature nodes and order of integration, we compute the quadrature weights by solving a system of linear moment fitting equations [11–13] for an appropriate set of basis functions in the least square sense. Setting up the moment-fitting equations involves computing integrals of the basis functions (known as moments) over an arbitrary geometric domain. This task itself is non-trivial, because it either entails some kind of domain decomposition, or can be reduced to repeated boundary integration as was proposed in [14] for bivariate domains. Adaption of the latter approach avoids excessive domain fragmentation, but leads to a significant computational overhead, particularly in 3D. We overcome this challenge by first computing the moments over a simpler domain Ω0 , usually a polygon/polyhedron, that is homeomorphic to and is a reasonable approximation (in a Hausdorff distance sense) of the original domain Ω . The computed moments are then corrected for the deviation from Ω based on first-order Shape Sensitivity Analysis (SSA) [15]. The moments over approximate domain are further reduced to boundary integrals by the application of divergence theorem [4–8]. The approximate moment fitting equations, thus obtained, can then be easily solved for quadrature weights that adapt to the original geometric domain—hence the name Adaptively Weighted Numerical Integration (AW). The resulting weights are not exact, but are accurate enough to integrate functions over any arbitrary domain when they are well approximated by the chosen set of basis functions. This integration scheme could be advantageously employed in cell decomposition methods, for example using quadtrees [16] or octrees [17], to reduce fragmentation. In Fig. 2(b), we show a simple comparison to demonstrate this. Here, the AW method proposed in this paper and Geometrically Adaptive (GA) integration method proposed by Luft et al. [18], are compared for integrating a 5th-degree polynomial (f = x2 y2 + x2 y3 + x3 + 100x + 10y + 2) over a wavy domain (see Fig. 2(a)) using quadtree with 3-point integration rule in each of the cells. The analytical integral computed using MATLAB’s symbolic toolbox is used as the reference for comparison. The graph in Fig. 2(b) clearly indicates that the proposed AW method is far more efficient than GA in terms of speed and memory, as it requires fewer quadtree subdivisions to achieve a desired accuracy. Additionally, the AW method is straightforward to implement and possesses a number of additional attractive features discussed in Section 5.2. The rest of the paper is organized as follows. Section 2 summarizes the relevant background material on geometrically adaptive methods of integration over piecewise linear and arbitrary domains. Section 3 summarizes the classical moment fitting approach and its relevance to the current formulation. The moment approximation for arbitrary domains is derived in Section 4. In Section 5, we outline the algorithm for adaptively weighted numerical integration (AW) and discuss its properties, including its running time. Section 6 describes experimental results and compares performance of AW method with other methods. The conclusion in Section 7 includes discussion of open issues that warrants further research.

1684

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

Fig. 2. (a) Wavy domain (b) Accuracy comparison between GA [18] and the proposed AW method.

2. Background Ideally, one would like to have a numerical integration method that is provably accurate and efficient. The first requirement implies that the method should avoid ad hoc and heuristic steps. The main challenge with the second requirement lies in the integrator’s necessity to adapt to both the integrand and geometric domain. The issues related to integrand adaptivity is well understood (e.g., see [19]). On the other hand, geometric adaptivity is far from complete and is the main focus of our research. In particular, the desire of efficiency suggests that the geometric adaptivity should require no or minimal subdivision. Hence, in this section, we will survey relevant ideas and literature on numerical integration techniques with emphasis on geometric adaptivity. 2.1. Numerical integration over polygonal and polyhedral domains Numerical integration over a regular n-dimensional interval, i.e. a box, is well understood (see for e.g. [20]). Stroud [21] describes cubature rules for a number of different simple domains (such as hexagon, octahedron, cubical shell, etc.) in his book four decades ago. Since then, many additional specialized cubature formulas have been developed for a variety of shapes and spaces [22–24]. More recently, Dasgupta [8] proposed a method to integrate symbolically integrable functions in the context of polygonal/polyhedral finite elements. Here, application of divergence theorem (once for polygonal domain and twice for polyhedral domain) converts the 2D/3D integral into an 1D integral over the edges of the polygon/polyhedron. The resulting integral is then efficiently evaluated using 1D Gauss Quadrature rules. For a priori unknown functions, Sommariva and Vianello [25] used Green’s integral formula with thin-plate splines basis functions to obtain a meshless cubature method for convex, nonconvex and even multiply connected polygonal domains. Further, Xiao and Gimbutas [13] devised a numerical algorithm for constructing efficient, high-order quadratures in two and higher dimensions by using moment fitting equations along with node elimination scheme. This algorithm was successfully applied to integrate general functions over triangle and square [13]. Mousavi et al. [26,3] used the same algorithm to construct efficient quadrature rules for bivariate polynomials to integrate functions over convex and non-convex polygons. Specifically, Mousavi and Sukumar [2] constructed efficient quadrature for integration over convex polygons and polyhedra based on moment fitting equations and monomial basis functions. They used Lasserre’s method [27,28] for the integration of homogeneous functions arising in the moment computations. The authors also observed that integration of the known basis functions may be transformed to boundary integrals using divergence theorem, at least in principle. In practice, repeated and accurate integration over the boundary of an arbitrary

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

1685

Fig. 3. Hierarchical partitioning (quadtree) (a) Gray: interior cells; Red: boundary cells (b) Boundary cells are randomly excluded from integration. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

domain is computationally expensive. We shall see, however, in Sections 4 and 5 that moment fitting and divergence theorem, when suitably applied, leads to efficient computation of weights in the AW method. 2.2. Adaptation by hierarchical partitioning The boundary of an arbitrary domain is a piecewise smooth curve in 2D or surface in 3D. Numerical integration over such domains must adapt to the domain’s boundary. This can be achieved by meshing—subdividing the domain into a union of cells that conform to the domain’s boundary and then applying well established quadrature rules to each of the cells. Recall that our goal is to avoid mesh generation. However, in principle, the boundary may be polygonized, but none of the proposed methods scale to handle polygons or polyhedra of such complexity. A common approach to geometric adaptation is hierarchical partitioning. In its simplest form, a quadtree [16] or octree [17] subdivision breaks the domain into non-conforming cells recursively. This kind of subdivision primarily leads to two types of integration cells: interior cells and boundary cells (see Fig. 3(a)). Interior cells can be easily handled using ordinary lattice rules (also called the Cartesian product rule). Several approaches to dealing with boundary cells have been proposed. For example: (a) The simplest method is to randomly include/exclude boundary cells (see Fig. 3(b)) and then use ordinary lattice rules on the randomly included cells [29]. (b) Alternatively, one could map the boundary cells to smoothly deformed rectangles [30] or other canonical regions [31] and then use ordinary lattice rules on it. (c) More recently, convenient space parametrization of boundary cells was found to produce good results [32,33,18]. (d) Monte-Carlo based methods [19] can be used to accurately integrate the boundary cells. (e) Finally, level-set based approaches [34,35] to deal with boundary cells. All of the proposed approaches to dealing with boundary cells are problematic. Specifically, referring to the above, (a) lacks accuracy, (b) does not generalize to 3D, (b) and (d) are computationally expensive, (b) requires function evaluation outside of the domain, (e) is specific to a particular representation of the domain and (c) involves heuristic steps. One of the main applications of the AW method proposed in this paper is in accurately integrating arbitrary functions arising in the boundary cells of such non-conforming decomposition, both in 2D and in 3D, without the need for excessive subdivision or unrealistic requirements. 2.3. Adaptation by characteristic function The problem of geometric adaptivity can be reduced in to integrand adaptivity over regular rectangular/box domain. Thus, well known integrand adaptivity techniques (over regular domains) can be used to compute the integral over arbitrary domains. To be precise, let f be the function to be integrated over the arbitrary domain Ω ⊂ D, where D is a regular domain. Let us define a characteristic (Heaviside) function H (x): H (x) =



0 1

: x ̸∈ Ω :x∈Ω

 so that Ω

f (x)dΩ =



f (x)H (x)dD. D

Now we have transformed the original integration domain Ω of the integral (Fig. 4(a)) into a regular box domain D (Fig. 4(b)). Thus, the required integral can be computed by sampling the modified function f (x)H (x) over the extended domain D using the ordinary lattice rules. Although this transformation is mathematically perfectly legitimate and sensible, it poses a serious problem when evaluated numerically owing to the discontinuity of the Heaviside function (H (x)). This method is particularly popular with immersed boundary methods [36]. A recent survey [37] discusses a number of strategies employed in Finite Cell Method [38] in dealing with the discontinuity of the Heaviside function. These include local meshing [39–41], heuristic modification of integration weights [42,37], smoothed step function approach [43,44], and equivalent polynomial approximation combined with moment fitting approach [45,37]. These methods help alleviate the problem but at additional computational cost and without solving the underlying fundamental issue.

1686

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

Fig. 4. Integrand adaptivity (a) Original integration domain (b) Extended integration domain.

Fig. 5. Reference and deformed domains for SSA.

2.4. Adaptation by shape sensitivity analysis We could also try to evaluate integral (1) in two steps: first integrate function f (x) over another domain Ω0 that approximates Ω , and then correct the result by the difference between the two integrals. To the best of our knowledge, this approach has not been considered in the context of integration over arbitrary domains. One popular approximation is based on the application of first-order Shape Sensitivity Analysis (SSA) [15] from the shape optimization [15] literature. As a first order approximation, consider a reference domain Ω0 that is homeomorphic to and in the neighborhood of the original domain Ω . In other words, there exists a continuous bijective mapping T (X, t ) (with continuous inverse) that deforms the reference domain into the original domain. In this context, the original domain Ω could also be referred to as the deformed domain. For simplicity, let us choose a homeomorphic piecewise linear domain as our reference domain (see for e.g. Fig. 5). On applying SSA, and assuming that f can be suitably extended outside of Ω , we get the following approximation to the integral:

 I = Ω

fdΩ ≈

 Ω0

fdΩ0 +

 Γ0

fVN dΓ0 ,

(2)

where VN is the normal component of the design velocity that is defined formally in Section 4.2. Thus the original integral over Ω was approximated using SSA over a simpler domain Ω0 . Furthermore, at least in principle, the integrals over the polygonal/polyhedral domain (Ω0 ) can be computed using the quadrature rules via the moment fitting equations as explained in [13,26] (or see Section 3). Although this method would give a reasonable approximation to the original integral, it suffers from a serious drawback. In order to evaluate Eq. (2), function f (x) must be sampled over the polygonal/polyhedral boundary Γ0 . This either assumes that Γ0 lies completely within the domain, or requires extension of the function to outside the domain (Ω ), usually based on application specific heuristic arguments as in [46]. To circumvent this rather severe restriction, the boundary term in Eq. (2) could be written in domain form using divergence theorem:

 I ≈ Ω0

[f + ∇ · (f V)]dΩ0 ,

(3)

reducing the original integration over Ω to that over approximate simplified domain Ω0 . Although this eliminates the need for function extension, it does require the domain velocity computation which is computationally expensive, ambiguous and/or not easily implemented [47]. These difficulties perhaps explain why SSA has not been used for integration purposes until now. Note, however, that SSA can be used without any problems when the integrand f is defined everywhere in Rd , as is the case with basis functions used in moment fitting. We will rely on this observation in Section 4 in formulating the AW method described in Section 5.

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

1687

3. Moment fitting equations A quadrature in Rd is a formula of the form: n 

wi f (xi ) ≈

i =1

 Ω

W (x)f (x)dΩ

(4)

where Ω ⊂ Rd is the integration region, f is an integrand defined on Ω , and W (x) is a non-negative weight function. A quadrature rule is determined by points xi ∈ Rd that are usually called quadrature nodes, and the quadrature weights wi [13]. A standard technique for constructing quadrature rules is to solve the following system of moment fitting equations [11–13]:



W (x)b1 (x)dΩ



   Ω   b1 (x1 )    W (x)b2 (x)dΩ   b2 (x1 )  Ω  =  ···     ··· bm (x1 )   W (x)bm (x)dΩ

b1 (x2 ) b2 (x2 )

··· bm (x2 )

··· ··· ··· ···

b1 (xn ) w1 b2 (xn )  w2 





. · · ·  · · · wn bm (xn )

(5)

Ω

The same equations may be written in a more compact form:

{M}m×1 = [A]m×n {W}n×1

(6)

with



W (x)b1 (x)dΩ



  Ω      W (x)b2 (x)dΩ  {M} =  Ω      ···   W (x)bm (x)dΩ

(7)

Ω

3 where {bi }m i=1 is the set of basis functions and {M} is the vector of moments defined over the domain of integration Ω . The quadrature nodes and the corresponding weights {xi , wi }ni=1 can be determined numerically by solving the above set of moment fitting equations. The resulting quadrature rule can be used to integrate any function that is in the function space spanned by these basis functions, assuming that the integral of basis functions and solution to moment fitting equations were computed exactly. Each integration point in d-dimensions contribute d + 1 unknowns: d coordinate components and m ⌉ could serve as an estimate for the number of points required to integrate m basis functions in d one weight. Thus, ⌈ d+ 1 dimensions [2]. The moment equations are not easy to solve as they form a system of non-linear equations. However, if we fix the position of the integration points as in [2], the equations become linear and can be solved easily. Knowing {M}, one could solve the above set of equations for weights using linear least squares as:

{W}n×1 = [AĎ ]n×m {M}m×1 (8) Ď Ď where [A ] is the Moore–Penrose pseudoinverse [49]. One popular stable algorithm to numerically compute [A ] is based on the QR factorization [50]. The moment fitting equations become linear only because we choose the quadrature nodes a priori. When the domain of integration is a convex polygon/polyhedron, there are a number of ways to choose the quadrature nodes as suggested in [2]. For multivariate domains, [48] suggests using approximate Fekete points, though it is not clear how to generate such points efficiently in an arbitrary domain. We will take a common approach [19] of adapting the tensor-product rule to arbitrary domains with quadrature nodes restricted to lie within the domain by appropriate scaling as illustrated in Fig. 6. 4. Approximation of moments We now set the weight W (x) = 1 in Eq. (4), and focus on evaluating Ω 1.f (x)dΩ via computation of moment fitting equations. We will construct the moment equations by first computing {M} over a simplified domain Ω0 that is homeomorphic to and is a reasonable approximation to the original domain Ω , and then correcting the value of the computed moments using shape sensitivity analysis. Note that the basis functions bi are defined throughout Rd and can be evaluated at points in or outside the original domain Ω . Below, we explain these steps in detail. We will denote position vectors of points on the original domain Ω by lowercase letters (such as x ∈ Ω ) and simplified domain Ω0 by uppercase letters (such as X ∈ Ω0 ).



3 Popular choices of basis functions include the bivariate polynomials xp yr (p + r ≤ o) and trivariate polynomials xp yr z s (p + r + s ≤ o) in 2D and 3D respectively (for the given integration order o). Other recommended choices of basis functions include Legendre and Chebyshev polynomials [48].

1688

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

Fig. 6. One-dimensional lattice rules are applied along the vertical line segments that are themselves allocated using the lattice rule in the horizontal direction.

4.1. Integration over polygonal/polyhedral domains We begin by replacing the integral of the basis function bi

 Mi =

Ω

bi (x) dΩ

(9)

by an integral of the same function over a piecewise linear (polygonal or polyhedral) domain Ω0 (simplified domain). The latter can be evaluated using divergence theorem and symbolic integration as explained in [8], because the basis functions bi are preselected and known a priori. 2 For the  purpose of illustration, let us assume that the simplified domain is the polygon (Ω0 ⊂ R ). Our objective is to simplify Ω bi (X)dΩ0 to an integral over the polygonal boundary. In order to apply the divergence theorem, we first require 0

a vector function φi = φXi (X , Y )ˆi + φYi (X , Y )ˆj such that the following is satisfied:

 Ω0

bi (X , Y )dXdY =

 Ω0

∇ · (φi )dXdY .

It can be easily proved [8] that for φi (X , Y ) = theorem to the RHS of Eq. (10) yields:

 Ω0

bi (X , Y )dXdY =

 Γ0

(10)



bi (X , Y )dX ˆi + 0ˆj the above is satisfied. Now, application of divergence

φXi NX dΓ0

(11)

where NX = N · ˆi is the X -component of the normal to the boundary. For a polygon, the normal is constant over its edges. ne k Hence, for a polygon with ne -edges (i.e. Γ0 = ∪k= 1 E0 ) one could rewrite the above more explicitly as:

 Ω0

bi (X , Y )dXdY =

ne   k=1

E0k

It is to be noted that φXi (X , Y ) =

φXi (X , Y ) =

X (m+1) Y n

(m+1)

φXi (X , Y )NXk dE0k . 

(12)

bi (X , Y )dX can be obtained symbolically (For example, if bi (X , Y ) = X m Y n , then

).

Analogously, the integral over the polyhedral domain Ω0 ⊂ R3 can be derived to be an integral over its faces as [8]:

 Ω0

bi (X , Y , Z )dXdYdZ =

nf   k =1

F0k

χXi (X , Y , Z )NXk dF0k

(13)

where χXi (X , Y , Z ) = bi (X , Y , Z )dX and F0k being the kth polygonal face of the polyhedron with nf faces. As each of these integrals is an integral over a polygon, Eq. (12) could in turn be applied to Eq. (13) to further simplify the above to a line integral over the edges of the polyhedron. However, it is more advantageous to not simplify Eq. (13) further for reasons that will become apparent in Section 4.3. Thus, in general if Ω0 ⊂ Rd (d = 2 or 3) is the polygonal/polyhedral domain of n∗ k integration with boundary Γ0 ⊂ Rd−1 comprising of n∗ edges/faces (i.e. Γ0 = ∪k= 1 Γ0 ), then:



 Ω0

bi (X)dΩ0 =

n∗   k=1

Γ0k

βXi (X)NXk dΓ0k

where βXi (X) = φXi (X , Y ) for polygons and βXi (X) = χXi (X , Y , Z ) for polyhedrons.

(14)

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

1689

4.2. Shape sensitivity analysis (SSA) We will take the simplified domain Ω0 as our reference domain, and estimate the value of the integral (9) over the original domain Ω as the value of the integral (14) plus a correction term C due to deformation of the reference domain Ω0 into the domain Ω . To estimate the correction term, we will use the well-known first-order Shape Sensitivity Analysis (SSA) [15], under the assumption that the piecewise linear reference domain Ω0 is homeomorphic to and approximates the deformed (original) domain Ω (see Fig. 5). Let T be a homeomorphic mapping that transforms the reference domain in to the deformed domain i.e. x = T(X, t ). Thus, x ∈ Ω is a function of the shape parameter t which denotes the amount of shape change in the design variable direction such that t = 0 represents the reference domain Ω0 [15]. If T(X, t ) is assumed to be regular enough in the neighborhood of t = 0, then it can be expanded using the Taylor series around the initial mapping point T(X, 0) as [15]: x = T(X, t )

= T(X, 0) + t

∂T (X, 0) + · · · ∂t

= X + t Vˆ (X) + · · ·

(15)

  ˆ (X) = ∂ T (X, 0) = dx  where V ∂t dt 

is the design velocity vector that is defined over the reference domain Ω0 . First-order

t =0

shape sensitivity is developed by considering only the first two terms of the above Taylor series. Thus, ignoring the higher order terms, T is assumed to be a homeomorphic linear mapping in t given by:

ˆ (X), x = T(X, t ) = X + t V

X ∈ Ω0 .

(16)

Having defined the mapping, the sensitivity of integral (9) to shape perturbation follows directly from the definition of SSA (see Appendix A for details):

      = d  = ∇ · (bi (X)Vˆ (X)) dΩ0 . b ( x ) d Ω i  dt t =0 dt Ω0 Ω t =0  Thus, Mi ≈ Ω bi (X) dΩ0 + C with the correction term C defined as 0    ˆ C =t ∇ · (bi (X)V(X))dΩ0 = ∇ · (bi (X)V(X))dΩ0 = bi (X)VN (X)dΓ0 , dMi 

Ω0

Ω0

(17)

(18)

Γ0

ˆ is the design velocity and VN (X) = V(X) · N(X) is the normal component of the design velocity at a point X on where V = t V the boundary Γ0 of the reference (piecewise linear) domain Ω0 . 4.3. Expression for {M} Thus, the final approximation of Mi is given as a sum:

 Mi =

Ω

bi (x)dΩ ≈

 Ω0

bi (X)dΩ0 + C .

(19)

Replacing the first term by an equivalent boundary integral (14) and the second term by its equivalent expression from Eq. (18), we obtain an approximation to the moment over arbitrary domain Ω :

 n∗     Mi  ≈ [βXi (X)NXk + bi (X)VN (X)]dΓ0k for i ∈ {1, 2, . . . , m}. k Ω

k=1

(20)

Γ0

The resulting boundary integral can be evaluated accurately and efficiently. For 2D domains, approximate {M} could be obtained from Eq. (20) simply by employing the usual Gaussian quadrature over the edges of the polygon. Similarly, integration over 3D domains bounded by triangles or convex quadrilaterals is straightforward. For more general 3D domains, we proceed in two steps. First, we compute the appropriate weights for a chosen set of quadrature nodes for each of the polygonal faces of the approximating polyhedron by solving the moment fitting equations (Eq. (8)) with exact {M} given by Eq. (14). Then, using the determined weights for the polygonal faces, we obtain approximate {M} for the arbitrary 3D domain by numerically evaluating Eq. (20) over the faces of the polyhedron. Using this approximate {M}, the approximate weights for the arbitrary domain are obtained by solving the moment fitting equations i.e. Eq. (8). The second term in the integrand of Eq. (20) is the shape sensitive term that accounts for the deviation of the original domain Ω from the polygonal/polyhedral domain Ω0 . In other words, this term ensures that the weights determined using the moment fitting equations (Eq. (8)) adapt to the domain of integration. In fact, the exact moments for a polygonal/polyhedral

1690

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

domain can be recovered from Eq. (20) by simply omitting the second term. In general, as the polygonal/polyhedral approximation (Γ0 ) approaches (in Hausdorff sense) the boundary of the original domain (Γ ), the second terms goes to 0, and the first term approaches the exact moments for the original domain. Note that the second term requires the computation of VN . For our purposes, it is sufficient to assume that the design ˆ (X) is piecewise continuous over every edge/face Γ0k of Ω0 . Then, the perturbation of the domain by distance velocity vector V α in the direction normal to the boundary is x − X = α N(X), which implies that

ˆ (X) · N(X) = t VN (X) = V (X) · N(X) = t V

x−X t

· N(X) = α N(X) · N(X) = α.

(21)

Thus, VN (X) is simply the shortest distance (α) from X to the original boundary (Γ ) as measured in the normal direction. 5. Adaptively weighted numerical integration (AW) all pieces together to arrive at adaptively weighted (AW) integration procedure to evaluate  We are now ready to put d f ( x ) d Ω over any Ω ∈ R ( d = 2 or 3). Ω 5.1. Algorithm The procedure follows essentially the same steps for 2D and 3D domains, and produces a set of integration nodes and weights that can be used to integrate any sufficiently smooth function over domain Ω . For the purpose of analysis, we assume that the shape is bounded by a degree k curve/surface and that k is known a priori (or at least can be estimated a posteriori). Conceptually, the procedure consists of the following three subtasks that we describe below in detail. Initialization involves setting the data structures required for solving the system of moment equations (8). Note that the dimensions of the matrix [AĎ ] depend on the choice of basis functions and on the number of integration nodes, both of which depend on the required order of integration. Specifically, initialization consists of the following steps: 1. Determine the order of integration (o) based on the function to be integrated. Choose a set of basis functions {bi }m i=1 of order up to o. A simple choice is the trivariate polynomials xp yr z s (3D) or bivariate polynomials xp yr (2D). 2. For the chosen order of integration (o), determine the number4 and location of quadrature nodes. One simple scheme is to use the tensor-product rule ensuring all points are inside the domain as explained in [18] (for example, see Fig. 6). 3. Choose a piecewise linear domain Ω0 (polygon in 2D or polyhedron in 3D) that is homeomorphic to and is a reasonable approximation of the original domain Ω . There are many ways to accomplish this and we will see some examples in the following section. Compute moments The left hand side of Eq. (6) is a vector of m moments, one moment for each basis function. The adaptively weighted evaluation procedure follows the development in Section 4, namely: 4. Determine the quadrature nodes and weights for each edge (2D) or face (3D) of the simplified domain Ω0 . If the edge/face approximates a surface of degree 2k−1 in the original boundary Γ , choose the integration order to be (o+k) (for example, k = 2 for conic sections and quadric surfaces). 5. Using the basis functions selected in step 1, evaluate the approximation to {M} over the original domain Ω from Eq. (20), by summing contribution over individual edges (2D) or faces (3D) of the approximate domain Ω0 (using the quadrature nodes and weights determined in step 4). Solve for weights and evaluate 6. Compute [AĎ ] for the preselected set of basis functions (step 1) and quadrature nodes (step 2) using say QR factorization [50]. 7. Using the {M} from Step 5 and [AĎ ] from step 6, compute the approximate weights for the arbitrary domain Ω from {W} = [AĎ ]{M} (Eq. (8)).  8. Finally, evaluate the approximation to the required integral by using these weights in Eq. (4) i.e. Ω f (x)dΩ ≈ n k=1 Wk f (xk ). It should be clear that the same algorithm may be used in most situations, with steps 4 and 5 depending strongly on type, dimension, and representations of the geometric domain Ω and its approximation Ω0 . Note that these two steps also dominate the computational complexity of the proposed approach, as we discuss below.

4 An n-point quadrature in 1D integrates a polynomial of degree up to 2n − 1 exactly and optimally. As of now, in higher dimension, there are no criterion available for determining the optimal number of quadrature nodes required to integrate any chosen set of (polynomial) basis functions [13]. However, an n-point tensor-product rule in nD, although not optimal, would still integrate a polynomial of degree up to 2n − 1 exactly.

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

1691

5.2. Algorithm analysis The main operations that are needed to implement the AW integration procedure are: construction of simplified approximation of Ω (step 3), distance to boundary computation (for design velocity computation), normal computations for polygonal edges/faces (for computing normal component of design velocity), and ability to sample points on a 2D/3D domain/boundary (for quadrature node generation). All these operations are supported by any standard CAD software and hence are readily available. For a more careful asymptotic analysis, let 2k − 1 be the degree of the domain boundary (Γ ), o be the order of integration, n∗ be the number of edges/faces in the approximating polygon/polyhedron Ω0 , and m(o) = O(od ) be the number of d-variate polynomial basis functions in d-dimensional space. We will assume that the cost of all arithmetic operations is constant. With this, we can estimate the running time of the algorithm. The initialization steps 1–2 are performed only once and do not affect the worst-case running time; step 3 is a construction of piecewise linear domain Ω0 , a task that is outside the scope of this paper. We note, however, that such approximations are often used by many CAD systems and may be available at no additional cost. Steps 6–8 are dominated by computation of [AĎ ] (step 6) and matrix multiplication (step 7). The cost of computing [AĎ ] using QR factorization can easily be proved to be O(o3d ) [50]. Also, the matrix multiplication costs O(o3d ). This leaves steps 4 and 5 that are at the very heart of the AW approach. Typically, the simplified domain Ω0 is bounded by triangles or convex quadrilaterals, step 4 is a trivial constant time operation. Close examination of Eq. (20) reveals that evaluation of moments in step 5 requires:

• • • •

O(n∗ ) of normal computations; O(n∗ (o + k)d−1 ) of distance computations; O(n∗ (o + k)d−1 od ) of multiplications/additions; and O(n∗ (o + k)d−1 od ) of basis function evaluations.

Assuming that all of these operations are constant, the running time of step 5 is O(n∗ (o + k)d−1 od ). Generally speaking, computing distance between two arbitrary domains is not a constant time operation; however, when correspondence between the faces in Γ and Γ0 is known, the distance computations are localized, and may be considered constant time for any practical purpose. We conclude that in most practical situations, the total cost of AW integration algorithm is O(o3d + n∗ (o + k)d−1 od ). For a fixed order of integration o, degree of bounding surfaces k and dimension of space d, the procedure is linear in the number of bounding faces/edges n∗ of the simplified domain Ω0 . This conclusion underscores one of the most appealing properties of the AW procedure: it does not assume any particular representation of the original domain Ω (such as voxels, mesh, splines, implicit surfaces, etc.), and can be used with any such representation, as long as it supports efficient distance computations and can be approximated by a polygonal domain Ω0 . In other words, k can be viewed as a measure of frequency of oscillations in the boundary of Γ relative to the simplified boundary Γ0 . Large k leads to the increase in the number of required integration points and distance computations which could undermine the advantages of AW integration procedure. On the other hand, AW procedure can be also used with arbitrarily complex approximate domains Ω0 . When a 3D Ω0 is bounded by more general polygonal faces, such as those found in [51–53], step 4 of the procedure requires constructing a quadrature for each face using the moment fitting equations as described in Section 3, with moments computed from Eq. (14). In principle, this could increase the theoretical complexity of the algorithm to O(n∗ (o + k)3(d−1) ) but could also reduce the overall cost of integration by allowing a much larger class of approximate domains in AW integration. 6. Experimental results In this section, in the first three examples, we compare the computational properties of Adaptively Weighted Numerical Integration (AW) method to four other methods and a reference symbolic integration. The four selected methods are: SS: a direct application of shape sensitivity analysis as explained in Section 2.4; C: the Cartesian product or tensor-product rule scaled to fit inside the geometry of the cell as recommended in [19]; GA: the Geometrically Adaptive integration method described in [18]; P: the quadrature rule obtained by solving the linear moment fitting equations (in least square sense) over a polygonized boundary (without correction). In the fourth example we just compare AW with Min and Gibou [34]. All algorithms were implemented in MATLAB 7.10 (on a Intel Core i7 CPU with 2.8 GHz speed and 8 GB memory), which was also used for symbolic integration, and tested on four simple 2D domains: quadrant of a circle in Section 6.1, notched domain in Section 6.2, ‘wavy’ domain quadtree decomposition in Section 6.3 and an elliptical domain in Section 6.4. We will denote the value of integral obtained from MATLAB’s symbolic integrator as IA and use I∗ to stand for the value of integral computed by method (*), where (*) could be one of the above methods. With this notation, we will compare the accuracy of integration for each domain in terms of the following error measures:

• Relative Error (%) = | IA I−A I∗ | ∗ 100; • Predicted Error-SS = ISS − IC ;

1692

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

Fig. 7. Approximate polygon (Ω0 ) construction by ray casting. The approximate polygon vertices are generated by casting rays parallel to integration rays and intersecting it with the original domain (Ω ).

• Predicted Error-AW = IAW − IC ; • Actual Error= IA − IC . For the first three examples, we will perform a series of tests aimed to establish how the accuracy of the algorithms compare in terms of their ability to approximate integrals of polynomial basis functions and integrands, adapt with better polygonal approximation, and improve with increased order of integration. Specifically, for each domain Ω , we perform the following tests: (a) Integration of basis functions The integrand for this test would be the bivariate polynomial basis functions of order up to 5. The approximating polygon Ω0 will have 7 edges (i.e. ne = 7), and 3-pt tensor-product rule (i.e. n = 3 × 3) will be used to generate the quadrature nodes. (b) Effect of polygonization Integral of f = x2 y2 + x2 y3 + x3 + 100x + 10y + 2 with 3-pt tensor-product quadrature rule (i.e. n = 3 × 3) is evaluated, as the number of edges in the approximating polygon Ω0 is varied. It is worth noting that increasing the number of edges does not always improve the polygonal approximation to Ω as we will see in one of the examples (Example 3). In other words, increasing the number of edges (finer polygonization) does not necessarily reduce the distance between Ω and Ω0 (in Hausdorff sense). (c) Dependence on quadrature rule Integral of f = x2 y2 + x2 y3 + x3 + 100x + 10y + 2 is evaluated over polygonal approximation with seven edges (i.e. ne = 7), while the quadrature rule is varied. For the third example in Section 6.3, which is based on quadtree decomposition, we perform one additional experiment showing how the accuracy of the three methods (AW, SS, GA) improves with increased depth of the approximating quadtree, while fixing ne = 7 and n = 3 × 3 in the leaf cells for the integrand f as defined above. For the fourth example, we only study the order of convergence of AW with respect to mesh refinement. The approximating polygon (Ω0 ) could be constructed in a number of ways including (i) ray casting [54], (ii) coarse quadtree decomposition [16] of Ω , and (iii) coarse surface meshing [55] of Ω . For the sake of simplicity, we employed method (i) to construct Ω0 . Specifically, we casted rays parallel to integration rays and intersected it with the original domain Ω to generate the polygon vertices (see Fig. 7). 6.1. Example 1—quadrant of a circle We will consider the domain Ω to be a quadrant of a circle (Fig. 9(a)) with the approximating polygon Ω0 as shown in Fig. 9(b). This example would demonstrate that, in an average sense, AW is at least as good or better than SS, C and P methods. a. Integration of basis functions The results are shown in Fig. 9(c). The relative errors of AW, SS, and C are in the range 0.012%–1.718%, 0.012%–1.821%, and 0%–8.270% respectively. The absolute value of the integrals and the relative errors are listed in Table 3 (Appendix B). Clearly, AW and SS are better than C in most cases, except for some lower order basis functions. Note that AW and SS correlate well with each other. b. Effect of polygonization The results are shown in Fig. 9(d). We observe a significant decrease in error for the P method, as expected. However, the results also indicate that for the P method to be as accurate as the AW or SS methods, it needs a finer

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

1693

Fig. 8. Cases in marching squares.

polygonal approximation. Fig. 9(e) shows that the error predicted using AW (and SS) correlates well with the actual error. This suggests that one could use AW and C methods to derive an error measure that could be used to drive an adaptive grid refinement process when integrating an arbitrary function over an arbitrary domain using a quadtree/octree. c. Dependence on quadrature rule The results are shown in Fig. 9(f). It is clear from Fig. 9(f) that only for ‘6-pt’ quadrature rule, the error in C is smaller than that of AW. Note that in this example, one of the quadrature nodes lies outside the approximating polygon Ω0 . The presence of nodes outside the polygon usually decreases the accuracy of the integral computed. Hence, whenever possible, it is best to choose a polygonal approximation that contains all the quadrature nodes well within the polygonal boundary. In this particular example, the presence of outside node for 3-pt rule does not have a serious impact on the results computed by AW method, as the quadrature weight corresponding to the node lying outside the polygon is close to zero. However, for higher order quadrature (4-pt and higher), outside nodes do have an effect on the results as can be seen from Fig. 9(f). However, the error (in AW and SS) stabilizes and becomes a constant with increase in quadrature rule. This is because, unlike the C method, increasing the quadrature rule beyond a limit does not usually improve the accuracy of AW. This is consistent with the fact that AW does not rely on the distribution of quadrature points to capture the geometry of the integration domain. 6.2. Example 2—notched domain We will consider the notched domain Ω (Fig. 10(a)) with the approximating polygon Ω0 as shown in Fig. 10(b). This example demonstrates that AW produces accurate results where C fails in spite of using a higher-order quadrature. a. Integration of basis functions The results are shown in Fig. 10(c). The relative errors of AW, SS, and C are in the range 0.010%–1.301%, 0.008%–1.550%, and 0.510%–14.271% respectively. The absolute value of the integrals and the relative errors are listed in Table 4 (Appendix B). Clearly, AW and SS are better than C for all basis functions as the C method does not detect the notch. As before, AW and SS are found to correlate well with each other. b. Effect of polygonization The results are shown in Fig. 10(d). Again, there is a significant decrease in error for the P method as expected. And, once again, for the P method to be as accurate as the AW or SS methods, it needs to have a finer polygonal approximation. In addition, Fig. 10(e) shows that the error predicted using AW (and SS) correlates well with the actual error. This further supports our previous observation that AW and C methods could be used to derive error measures that drive the adaptive grid refinement process. c. Dependence on quadrature rule The results are shown in Fig. 10(f) and demonstrate that ‘C’ cannot compete with AW and SS even for 6-pt quadrature because it fails to detect the presence of the notch. 6.3. Example 3—wavy domain with quadtrees In this example, we will demonstrate how AW method can be used in a quadtree based integration. For comparison purposes, AW, SS and GA methods are applied to integration over the leaf cells of the quadtree that decomposes the wavy domain5 (Fig. 11(a)). The approximating polygons and the quadrature nodes were constructed in the leaf cells following the basic rules of the GA method [18], by accounting for the type of leaf cell that arises in marching squares (see Fig. 8). The interior cells (case 15) were integrated using the regular Cartesian product rule while the exterior cells (case 0) were ignored. The quadrature nodes for the leaf cells were allocated based on Cartesian coordinates (with appropriate scaling) for cases 3, 6, 7, 9, 11, 12, 13 & 14 and polar coordinates for cases 1, 2, 4 & 8. The ambiguous cases (case 5 and 10) were handled simply by recursively dividing the cell until the ambiguity is resolved. a. Integration of basis functions The results are shown in Fig. 11(c). The quadtree of depth one (i.e. D = 1) along with the approximating polygons (ne = 7) used in the leaf cells are shown in Fig. 11(b). The absolute value of the integrals and the relative errors are listed in Table 5 (Appendix B). The relative errors of AW, SS and GA are in the range 0.004%–0.557%, 0.009%–0.557% and 1.975%–12.180% respectively. It is apparent that AW and SS are more accurate than GA for all integrands. As in previous examples, AW and SS correlate well with each other.

5 The wavy surface is defined by the equation y(x) = 3.2 + 0.09 sin(10.5π x).

1694

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

a

b

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

c

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

8

Relative Error (%)

7 6 5 C SS AW

4 3 2 1 0

d

14

Relative Error (%)

12 10 8

C

6

SS AW

4 P 2

e

0

0

5

5

10

7

15

20

9

11

13 15 # of polygon Edges

f

25

-0.1

Error

-0.2 Actual Error

-0.3

Predicted Error-SS Predicted Error-MW

-0.4

19

21

4 3 2.5 C SS AW

2 1.5 1 0.5

-0.5 -0.6

17

3.5 Relative Error (%)

0

0 2 # of polygon Edges

3

4 5 Quadrature Rule

6

Fig. 9. (a) Quadrant of a circle (Ω ) (b) Approximating polygon (Ω0 ) (c) Rel. error (%) vs basis functions (for bij = xi yj , n = 3 × 3, ne = 7) (d) Rel. error (%) vs no. of polygon edges (for integrating f with n = 3 × 3) (e) error prediction and (f) Rel. error (%) vs quadrature rule (for integrating f with ne = 7).

b. Effect of polygonization The number of edges in polygons approximating the leaf cells is varied, while all other parameters remain the same: D = 1, n = 3 × 3 (Fig. 12(a)–(d)) and the integrand f over the wavy domain. For this test only, we also include integration results with P method applied to the leaf nodes of the quadtree. The relative errors for AW, SS, P and GA methods are in the range 0.003%–1.569%, 0.005%–1.424%, 0.05%–2.455% and 2.91%–2.91% respectively. However, GA method is not adapting to finer polygonization as one would expect, while P method is improving but somewhat inconsistently. This is because, in this particular case, finer polygonization does not necessarily mean better approximation

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

a

b

1

1

0.9

0.9

0.8

0.8 0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

c

1695

0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

14

Relative Error (%)

12 10 8

C SS AW

6 4 2 0

d

3.5 3

Relative Error (%)

2.5 C

2

SS

1.5

AW

1

P 0.5 0 6

8

10

12

14

16

18

20

22

# of polygon Edges

f

1.85

Error

3.5 3

1.8 1.75 Actual Error

1.7

Predicted Error-SS Predicted Error-AW

1.65 1.6 0

5

10 15 20 # of polygon Edges

25

Relative Error (%)

e

2.5 2

C SS AW

1.5 1 0.5 0

2

3

4 5 Quadrature Rule

6

Fig. 10. (a) Notched domain (Ω ) (b) Approximating polygon (Ω0 ) (c) Rel. error (%) vs basis functions (for bij = xi yj , n = 3 × 3, ne = 7) (d) Rel. error (%) vs no. of polygon edges (for integrating f with n = 3 × 3) (e) Error prediction and (f) Rel. error (%) vs quadrature rule (for integrating f with ne = 7).

of the original domain Ω . In other words, increasing the number of polygon edges does not necessarily reduce the distance between Ω0 and Ω (in Hausdorff sense). This also explains why the error in AW oscillates with finer polygonization in contrast to the previous two examples (Figs. 9(d) and 10(d)). The number of edges required for the relative error to fall below 0.1% are 9 for AW, 9 for SS and 14 for P method. This suggests that P method requires finer and better polygonization

1696

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

Fig. 11. Quadtree integration of bivariate polynomial basis functions (bij = xi yj ) of order up to five over the wavy domain for ne = 7 and n = 3 × 3. (a) Wavy domain (Ω ) (b) Quadtree with D = 1 and (c) Rel. error (%) vs basis functions. Table 1 Relative errors, order and computational time for integral of f over the wavy domain. Depth

h

1 2 3 4

1.929 0.965 0.482 0.241

Relative error (%)

Order

GA

SS

AW

2.9106 1.0998 0.0633 0.0007

0.2268 0.0640 0.0255 0.0004

0.2357 0.0643 0.0255 0.0004

GA 1.40 4.12 6.41

Time (ms) SS 1.83 1.33 5.94

AW

GA

SS

AW

1.87 1.34 5.94

47 94 187 562

125 172 343 655

78 172 343 780

to produce accurate integrals over domains with rapidly varying boundaries and likely to be more expensive than AW and SS methods. c. Dependence on quadrature rule Number of quadrature nodes are varied in the leaf cells with D = 1 and ne = 7 (Fig. 13(a)–(c)). Again, even with ‘6-pt’ rule in the leaf cells, the error in GA could not be made smaller than that of SS and AW. d. Accuracy vs depth of quadtree Depth of the quadtree is varied with ne = 7 and n = 3 × 3 (Fig. 14(a)–(d)). From the results, it is clear that AW and SS are able to produce significantly more accurate results even at coarser resolutions (say D = 1 or 2) when compared to GA method. In fact, from Table 1 it is clear that for the relative error to fall below 0.05%, AW (and SS) requires 343 ms as against GA’s 562 ms. Also, we find that, in general, the error decreases with refinement for all the 3 methods. However, the order of convergence varies. This is explained by the fact that the order of accuracy depends not only on the level of refinement but also on the kind of polygonal approximation at any given resolution. Since, in this example, the degree of polygonal approximation varies non-uniformly from one level to another we do not observe a consistent order of convergence. 6.4. Example 4—area of an ellipse For the sake of completeness, in this example we will compute the area of an ellipse represented by the zero level set x2 1.52

+

y2 0.752

− 1 = 0 using AW for various mesh sizes (h) and compare it with Min and Gibou [34]. For AW, we employed 6 or 7 edges per boundary cell (ne = 6 or 7) and 2-pt rule (n = 2 × 2) in each of the cells.

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

1697

Fig. 12. Quadtree integration of f over the wavy domain for various polygonal approximations of the leaf cells with D = 1 and n = 3 × 3. (a) ne = 4 (b) ne = 6 (c) ne = 11 and (d) Rel. error (%) vs no. of polygon edges (in leaf cells).

Table 2 Relative errors and order of convergence in computing the area of an ellipse. h

0.80 0.40 0.20 0.10 0.05 0.025 0.0125 0.00625

AW

Min and Gibou [34]

Relative error (%)

Order

Relative error (%)

Order

1.05E−03 1.64E−04 3.43E−05 6.35E−06 1.65E−06 3.77E−07 9.55E−08 2.33E−08

– 2.68 2.25 2.43 1.94 2.13 1.98 2.04

– – 1.59E−02 3.76E−03 9.46E−04 2.25E−04 5.78E−05 1.46E−05

– – – 2.08 1.99 2.07 1.97 1.99

As can be seen from Table 2, for this example, we have quadratic convergence similar to that of Min and Gibou [34]. However, for any given resolution, the method is found to be much more accurate than Min and Gibou [34]. In fact, the accuracy of AW with one of the coarsest meshes (say h = 0.10) cannot even be matched by the finest mesh (h = 0.00625) of Min and Gibou [34]. 7. Conclusions We presented a method to accurately and efficiently integrate arbitrary 2D/3D domains based on moment fitting equations [11–13], divergence theorem [8] and shape sensitivity analysis [15]. We demonstrated the use of the method in integrating arbitrary polynomial functions over arbitrary 2D domains and compared it with four other methods: scaled Cartesian product rule as recommended in [19], geometric adaptive (GA) integration method proposed in [18], polygonal

1698

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

Fig. 13. Quadtree integration of f over the wavy domain for various quadrature rules with D = 1 and ne = 7. (a) n = 5 × 5 (b) n = 6 × 6 and (c) Rel. error (%) vs quadrature rule.

(P) approximation method, and shape sensitivity (SS) method. The AW method was shown to be superior to the GA and to P methods in terms of accuracy and comparable to SS method in most cases. Recall from Section 2.4 that we abandoned the SS approach due to its unrealistic requirements such as function extension outside the domain, domain velocity computation inside the domain etc. Yet, experimental results show that SS and AW correlate well with each other when applied to functions that can be readily defined even outside the original domain of integration. In fact this is not a coincidence. To understand this, consider an alternate way to formulate AW. To begin with, let us rewrite the boundary integral in Eq. (2) in domain form as:

 Γ0

f (X)VN (X)dΓ0 =

 Ω0

 = Ω0

∇ · (f (X)V(X))dΩ0 f (X)Ψ (X)dΩ0

(22)

where

  ∇ · (f (X)V(X)) Ψ (X) = f (X)  ∇ · (f (X)V(X))

:X∈A : X ∈ Ω0 − A

where A = {X ∈ Ω0 |f (X) ̸= 0} is the set of points on the domain where the function f is non-zero. It can be easily proved that f (X)Ψ (X) is integrable and that the equality holds in Eq. (22). Substituting Eq. (22) in Eq. (2) we get:

 I ≈ Ω0

f (X)dΩ0 +

 Ω0

f (X)Ψ (X)dΩ0 =

 Ω0

f (X)(1 + Ψ (X))dΩ0 .

(23)

Setting W (X) = (1 + Ψ (X)), we get:

 I ≈ Ω0

W (X)f (X)dΩ0 .

(24)

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

1699

Fig. 14. Quadtree integration of f over the wavy domain for various quadtree depths with ne = 7 and n = 3 × 3. (a) D = 2 (b) D = 3 (c) D = 4 and (d) Rel. error (%) vs depth of quadtree.

Now, we have Mi ≈ Ω W (X)bi (X)dΩ0 with W (X) = (1 + Ψ (X)). Further, substituting for W (X) and simplifying 0 we obtain exactly the same expression for Mi as before i.e. Eq. (20). In other words, whether we first approximate the integral of the function f (over Ω by applying SSA over Ω0 ) and then derive Mi or approximate Mi directly using SSA, we obtain the same weights. Equivalently, this implies that AW is as good as SS but without the limitations of SS as outlined in Section 2.4. We anticipate that the main application of AW will be in non-conforming cell based integration. We observed that AW, when used in the leaf cells of quadtrees, required fewer subdivisions/quadrature nodes to resolve the geometry of the integration domain when compared to other methods. Based on our experiments, 8–10 edges per leaf cell seems to be a good compromise in achieving high accuracy at reasonable computational cost. More generally, the AW method satisfies many of the characteristics of an ideal numerical integrator postulated in Section 2. In particular, we note that the method is formulated irrespective of how the original domain Ω is represented, and it can be used with or without integration mesh. The latter observation suggests that AW is suitable for variety of meshless and meshfree applications such as Scan & Solve [56,57], Finite Cell Method [38] and other immersed boundary methods [36]. For the purposes of this paper, we have assumed that the integrand is continuous. However, AW method can be extended to integrate discontinuous and weakly discontinuous functions by splitting the domain into two subregions along the line of discontinuity and developing two different quadrature formulas for each of the subregions separately as in [2]. Thus, future promising research directions include verifying the efficacy of the method in handling discontinuous functions and 3D domains. In addition, appropriate and efficient selection of integration points or generation of approximate Fekete points in arbitrary 3D domains is an open problem that warrants further research.



Acknowledgments This research work was supported by the National Science Foundation grants CMMI-0856778 and CMMI-1029553. We would like to thank Dr. S.E. Mousavi and Dr. N. Sukumar of UC Davis for helpful discussions. We would also like to thank

1700

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

Table 3 Integral of basis functions over the quadrant. Basis function

b00 b10 b01 b20 b11 b02 b30 b21 b12 b03 b40 b31 b22 b13 b04 b50 b41 b32 b23 b14 b05

Integral value

= x0 y0 = x1 y0 = x0 y1 = x2 y0 = x1 y1 = x0 y2 = x3 y0 = x2 y1 = x1 y2 = x0 y3 = x4 y0 = x3 y1 = x2 y2 = x1 y3 = x0 y4 = x5 y0 = x4 y1 = x3 y2 = x2 y3 = x1 y4 = x0 y5

Relative error (%)

C

SS

AW

Analytical

C

SS

AW

0.7890 0.3372 0.3333 0.2006 0.1250 0.1961 0.1380 0.0667 0.0664 0.1333 0.1035 0.0417 0.0324 0.0417 0.0983 0.0825 0.0287 0.0185 0.0190 0.0287 0.0763

0.7850 0.3334 0.3332 0.1944 0.1247 0.1962 0.1309 0.0662 0.0665 0.1332 0.0965 0.0411 0.0326 0.0416 0.0981 0.0760 0.0282 0.0188 0.0191 0.0285 0.0762

0.7850 0.3333 0.3331 0.1945 0.1248 0.1963 0.1310 0.0662 0.0665 0.1333 0.0965 0.0411 0.0326 0.0416 0.0981 0.0759 0.0282 0.0188 0.0190 0.0285 0.0761

0.7854 0.3333 0.3333 0.1963 0.1250 0.1963 0.1333 0.0667 0.0667 0.1333 0.0982 0.0417 0.0327 0.0417 0.0982 0.0762 0.0286 0.0190 0.0190 0.0286 0.0762

0.4609 1.1685 0.0000 2.1605 0.0000 0.1056 3.5009 0.0000 0.3864 0.0000 5.4327 0.0000 1.1118 0.0000 0.0956 8.2698 0.6250 2.8577 0.4687 0.6021 0.0781

0.0497 0.0120 0.0487 1.0001 0.2166 0.0676 1.8206 0.6758 0.2352 0.1014 1.7039 1.2972 0.3209 0.2289 0.0649 0.1982 1.3631 1.2174 0.1990 0.1579 0.0422

0.0460 0.0122 0.0703 0.9365 0.1893 0.0280 1.7183 0.6660 0.2082 0.0414 1.6840 1.2858 0.4688 0.2513 0.0660 0.4042 1.2967 1.5349 0.1494 0.2575 0.1128

the anonymous reviewers for their helpful comments/suggestions that have helped improve the paper. Responsibility for errors and omissions lies solely with the authors. Appendix A. SSA of moments Recall that Mi =



Ω

bi (x)dΩ . It is required to compute the shape sensitivity of the moments Mi i.e.

dMi dt

   

. In order

t =0

to evaluate this, it is first required to transform the moment integral over the reference domain Ω0 . This can be done by making use of the Jacobian matrix of the mapping T(X, t ). It can be easily established that dΩ = |J |dΩ0 [15], where |J | is the determinant of the Jacobian matrix. Thus, we have the following:

     = d  b ( x ) | J | d Ω i 0   dt t =0 dt Ω0 t =0

dMi 

now, taking the derivative inside the integral and then applying the product and chain rule of differentiation we get: d



dt

Ω0

bi (x) |J |dΩ0

   



d

= Ω0

t =0

= Ω0

ˆ (x) = where V d



dt

Ω0

dx dt

dt





  (bi (x) |J |) dΩ0 

t =0

∇ bi (x) · Vˆ (x)|J | + bi (x)

d|J | dt



 

dΩ0 

is the design velocity vector at any point x ∈ Ω . Now, using the fact that

  bi (x) |J |dΩ0 

t =0

     ˆ ˆ = ∇ bi (x) · V(x)|J | + bi (x)|J |∇ · V(x) dΩ0  Ω0 t =0      = ∇ bi (x) · Vˆ (x) + bi (x)∇ · Vˆ (x) |J |dΩ0  Ω0

 = Ω0

t =0

  ∇.(bi (x)Vˆ (x))|J |dΩ0 

t =0

evaluating the RHS of the above at t = 0, we obtain the desired sensitivity as:

   = ∇.(bi (X)Vˆ (X))dΩ0 dt t =0 Ω0

dMi 

where X ∈ Ω0 .

t =0 d|J | dt

= |J |∇ · Vˆ (x) [15], we get:

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702

1701

Table 4 Integral of basis functions over the notched domain. Basis functions b00 b10 b01 b20 b11 b02 b30 b21 b12 b03 b40 b31 b22 b13 b04 b50 b41 b32 b23 b14 b05

= x0 y 0 = x1 y 0 = x0 y 1 = x2 y 0 = x1 y 1 = x0 y2 = x3 y0 = x2 y1 = x1 y2 = x0 y3 = x4 y0 = x3 y1 = x2 y2 = x1 y3 = x0 y4 = x5 y0 = x4 y1 = x3 y2 = x2 y3 = x1 y4 = x0 y5

Integral value C 0.9556 0.4778 0.4578 0.3222 0.2289 0.2932 0.2444 0.1561 0.1466 0.2118 0.1972 0.1197 0.1011 0.1059 0.1636 0.1653 0.0974 0.0783 0.0738 0.0818 0.1320

Relative error (%) SS 0.9845 0.4922 0.4843 0.3294 0.2429 0.3198 0.2481 0.1627 0.1594 0.2372 0.1990 0.1225 0.1083 0.1174 0.1868 0.1662 0.0983 0.0827 0.0810 0.0919 0.1526

AW 0.9845 0.4922 0.4845 0.3294 0.2429 0.3197 0.2481 0.1627 0.1594 0.2371 0.1990 0.1226 0.1081 0.1175 0.1868 0.1662 0.0984 0.0825 0.0808 0.0922 0.1528

Analytical 0.9843 0.4921 0.4850 0.3294 0.2425 0.3189 0.2480 0.1629 0.1595 0.2362 0.1990 0.1231 0.1075 0.1181 0.1867 0.1661 0.0990 0.0815 0.0798 0.0934 0.1539

C 2.9195 2.9195 5.6048 2.1693 5.6048 8.0692 1.4248 4.1489 8.0692 10.3262 0.8730 2.7146 5.9522 10.3262 12.3890 0.5107 1.6582 3.8806 7.5919 12.3890 14.2709

SS 0.0256 0.0082 0.1274 0.0182 0.1856 0.2904 0.0315 0.1151 0.0311 0.4531 0.0401 0.4730 0.7364 0.5875 0.0544 0.0436 0.7489 1.5504 1.4125 1.5496 0.8892

AW 0.0244 0.0096 0.1024 0.0181 0.1595 0.2474 0.0295 0.0922 0.0218 0.3820 0.0355 0.3917 0.6204 0.4889 0.0458 0.0357 0.6254 1.3015 1.1856 1.2967 0.7470

Table 5 Integral of basis functions over the wavy domain. Basis functions

b00 b10 b01 b20 b11 b02 b30 b21 b12 b03 b40 b31 b22 b13 b04 b50 b41 b32 b23 b14 b05

= x0 y 0 = x1 y 0 = x0 y 1 = x2 y0 = x1 y1 = x0 y2 = x3 y0 = x2 y1 = x1 y2 = x0 y3 = x4 y0 = x3 y1 = x2 y2 = x1 y3 = x0 y4 = x5 y0 = x4 y1 = x3 y2 = x2 y3 = x1 y4 = x0 y5

Integral value

Relative error (%)

GA

SS

AW

Analytical

GA

SS

AW

6.2743 6.2743 9.8417 8.3715 9.8417 20.5837 12.5659 13.1405 20.5837 48.4319 20.1172 19.7380 27.5020 48.4319 121.5551 33.5444 31.6175 41.3386 64.7551 121.5551 317.7952

6.4135 6.4115 10.2718 8.5511 10.2704 21.9308 12.8300 13.7040 21.9315 52.6733 20.5350 20.5698 29.2759 52.6771 134.9193 34.2398 32.9364 43.9604 70.3402 134.9168 359.8569

6.4238 6.4159 10.2772 8.5525 10.2716 21.9259 12.8296 13.7030 21.9295 52.6618 20.5339 20.5689 29.2746 52.6786 134.9275 34.2407 32.9373 43.9614 70.3403 134.9165 359.8555

6.4055 6.4055 10.2615 8.5442 10.2615 21.9271 12.8218 13.6936 21.9271 52.7321 20.5235 20.5578 29.2733 52.7321 135.3215 34.2202 32.9203 43.9657 70.4284 135.3215 361.8728

2.0479 2.0479 4.0907 2.0217 4.0907 6.1268 1.9955 4.0392 6.1268 8.1548 1.9799 3.9877 6.0509 8.1548 10.1731 1.9750 3.9572 5.9753 8.0554 10.1731 12.1804

0.1255 0.0951 0.1002 0.0800 0.0864 0.0167 0.0640 0.0761 0.0199 0.1116 0.0560 0.0582 0.0089 0.1043 0.2972 0.0573 0.0490 0.0120 0.1252 0.2991 0.5571

0.2866 0.1636 0.1524 0.0969 0.0985 0.0055 0.0605 0.0682 0.0106 0.1333 0.0503 0.0538 0.0042 0.1015 0.2912 0.0599 0.0518 0.0098 0.1251 0.2993 0.5574

Appendix B. Tables See Tables 3–5. References [1] A. Requicha, Representations of rigid solid objects, in: Computer Aided Design Modelling, Systems Engineering, CAD-Systems, 1980, pp. 1–78. [2] S.E. Mousavi, N. Sukumar, Numerical integration of polynomials and discontinuous functions on irregular convex polygons and polyhedrons, Comput. Mech. 47 (2011) 535–554. [3] S.E. Mousavi, N. Sukumar, Generalized Gaussian quadrature rules for discontinuities and crack singularities in the extended finite element method, Comput. Methods Appl. Mech. Eng. 199 (2010) 3237–3249. [4] Y.T. Lee, A. Requicha, Algorithms for computing the volume and other integral properties of solids. I. Known methods and open issues, Commun. ACM 25 (1982) 635–641. [5] C. Cattani, A. Paoluzzi, Boundary integration over linear polyhedra, Comput. Aided Des. 22 (1990) 130–135.

1702 [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57]

V. Thiagarajan, V. Shapiro / Computers and Mathematics with Applications 67 (2014) 1682–1702 F. Bernardini, Integration of polynomials over n-dimensional polyhedra, Comput. Aided Des. 23 (1991) 51–58. B. Mirtich, Fast and accurate computation of global geometric properties of solid objects, J. Graph. GPU Game Tools 1 (1996) 31–50. Gautam Dasgupta, Integration within polygonal finite elements, J. Aerosp. Eng. 16 (2003) 9–18. P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, Vol. 1, Academic Press, New York, 1975. A.R. Krommer, C.W. Ueberhuber, Computational Integration, Vol. 53, Society for Industrial Mathematics, 1998. J. Lyness, D. Jespersen, Moderate degree symmetric quadrature rules for the triangle, IMA J. Appl. Math. 15 (1975) 19–32. S. Wandzura, H. Xiao, Symmetric quadrature rules on a triangle, Comput. Math. Appl. 45 (2003) 1829–1840. H. Xiao, Z. Gimbutas, A numerical algorithm for the construction of efficient quadrature rules in two and higher dimensions, Comput. Math. Appl. 59 (2010) 663–676. A. Sommariva, M. Vianello, Gauss–Green cubature and moment computation over arbitrary geometries, J. Comput. Appl. Math. 231 (2009) 886–896. K.K. Choi, N.H. Kim, Structural Sensitivity Analysis and Optimization I: Linear Systems, Springer, New York, 2005. J.J. Laguardia, E. Cueto, M. Doblare, A natural neighbor Galerkin method with quadtree structure, Internat. J. Numer. Methods Engrg. 63 (2005) 789–812. O. Klaas, M. Shephard, Automatic generation of octree based three dimenisonal discretizations for partition of unity methods, Comput. Mech. 25 (2000) 296–304. B. Luft, V. Shapiro, I. Tsukanov, Geometrically adaptive numerical integration, in: 2008 ACM Symposium on Solid and Physical Modeling, NY, 2008, pp. 147–157. William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery, Numerical Recipes in C: The Art of Scientific Computing, second ed., Cambridge University Press, 1992. Dalton R. Hunkins, Cubatures of precision 2k and 2k + 1 for hyperrectangles, Math. Comp. 29 (1975) 1098–1104. A.H. Stroud, Approximate Calculation of Multiple Integrals, Prentice Hall, Englewood Cliffs, NJ, 1971. R. Cools, P. Rabinowitz, Monomial cubature rules since stroud: a compilation, J. Comput. Appl. Math. 49 (1993) 309–326. R. Cools, An encyclopaedia of cubature formulas, J. Complexity 19 (2003) 445–453. R. Cools, E. Goor, Encyclopaedia of cubature formulas, http://nines.cs.kuleuven.be/ecf/. A. Sommariva, M. Vianello, Product Gauss cubature over polygons based on Green’s integration formula, BIT 47 (2007) 441–453. S.E. Mousavi, H. Xiao, N. Sukumar, Generalized Gaussian quadrature rules on arbitrary polygons, Int. J. Numer. Methods Eng. 82 (2010) 99–113. J.B. Lasserre, Integration on a Convex Polytope, vol. 126, The American Mathematical Society, 1998, pp. 2433–2441. J.B. Lasserre, Integration and Homogeneous Functions, vol. 127, The American Mathematical Society, 1999, pp. 813–818. R. Sarraga, Computation of surface areas in GMSolid, IEEE Comput. Graph. Appl. 2 (1982) 65–70. K. Höllig, Finite Element Methods with B-Splines, in: Frontiers in Applied Mathematics, SIAM, 2003. L.F. Shampine, MATLAB program for quadrature in 2D, Appl. Math. Comput. 202 (2008) 266–274. V. Rvachev, A. Shevchenko, V. Veretelnik, Numerical integration software for projection and projection-grid methods, Cybernet. Systems Anal. 30 (1994) 154–158. I. Tsukanov, V. Shapiro, The architecture of SAGE—a meshfree system based on RFM, Eng. Comput. 18 (2002) 295–311. Chohong Min, Frédeŕic Gibou, Geometric integration over irregular domains with application to level-set methods, J. Comput. Phys. 226 (2007) 1432–1443. Peter Smerka, The numerical approximation of a delta function with application to level set methods, J. Comput. Phys. 211 (2006) 77–90. R. Mittal, G. Iaccarino, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 239–261. A. Abedian, J. Parvizian, A. Düster, H. Khademyzadeh, E. Rank, Performance of different integration schemes in facing discontinuities in the finite cell method, Int. J. Comput. Methods 10 (2013). J. Parvizian, A. Düster, E. Rank, Finite cell method, Comput. Mech. 41 (2007) 121–133. G.R. Liu, Meshfree Methods: Moving Beyond the Finite Element Method, CRC, 2009. G.R. Liu, T.T. Nguyen, Smoothed Finite Element Methods, CRC, 2010. Y. Abdelaziz, A. Hamouine, A survey of the extended finite element, Comput. Struct. 86 (2008) 1141–1151. T. Rabczuk, P.M.A. Areias, T. Belytschko, A meshfree thin shell method for nonlinear dynamic fracture, Internat. J. Numer. Methods Engrg. 72 (2007) 524–548. S. Osher, R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Springer-Verlag, 2003. A. Kumar, J. Lee, Step function representation of solid models and application to mesh free engineering analysis, Trans. ASME, J. Mech. Des. 128 (2006) 46–56. G. Ventura, On the elimination of quadrature subcells for discontinuous functions in the eXtended finite-element method, Internat. J. Numer. Methods Engrg. 66 (2006) 761–795. Ming Li, Shuming Gao, Estimating defeaturing-induced engineering analysis errors for arbitrary 3d features, Comput. Aided Des. 43 (2011) 1587–1597. K.K. Choi, N.H. Kim, Structural Sensitivity Analysis and Optimization II: Non-Linear Systems, Springer, New York, 2005. A. Sommariva, M. Vianello, Computing approximate Fekete points by QR factorizations of vandermonde matrices, Comput. Math. Appl. 57 (2009) 1324–1336. G. Strang, The fundamental theorem of linear algebra, Amer. Math. Monthly 100 (1993) 848–855. N.T. Lloyd, B. David, Numerical Linear Algebra, SIAM, 1997. M. Wicke, M. Botsch, M. Gross, A finite element method on convex polyhedra, Comput. Graph. Forum 26 (2007) 355–364. S. Martin, P. Kaufmann, M. Botsch, M. Wicke, M. Gross, Polyhedral finite elements using harmonic basis functions, in: Eurographics Symposium on Geometry Processing (Copenhagen, Denmark, July 2–4, 2008), Computer Graphics Forum, 2008, pp. 1521–1529. P. Kaufmann, S. Martin, M. Botsch, M. Gross, Flexible simulation of deformable models using discontinuous Galerkin FEM, in: 2008 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, 2008, pp. 105–115. D.R. Scott, Ray casting for modeling solids, Comput. Graph. Image Process. 18 (1982) 109–144. Steven J. Owen, A survey of unstructured mesh generation technology, in: 7th International Meshing Roundtable, Sandia National Lab, 1998, pp. 239–267. M. Freytag, Scan and solve: acquiring the physics of artifacts, in: ASME 2007 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, Las Vegas, USA, 2007. M. Freytag, V. Shapiro, I. Tsukanov, Finite element analysis in situ, Finite Elem. Anal. Des. 47 (2011) 957–972.