Finite Element Quadrature of Regularized Discontinuous and Singular Level Set Functions in 3d Problems

Regularized Heaviside and Dirac delta function are used in several fields of computational physics and mechanics. Hence the issue of the quadrature of integrals of discontinuous and singular functions arises. In order to avoid ad-hoc quadrature procedures, regularization of the discontinuous and the singular fields is often carried out. In particular, weight functions of the signed distance with respect to the discontinuity interface are exploited. 19: 527–552) proved that the use of compact support weight function is not suitable because it leads to errors that do not vanish for decreasing mesh size. They proposed the adoption of non-compact support weight functions. In the present contribution, the relationship between the Fourier transform of the weight functions and the accuracy of the regularization procedure is exploited. The proposed regularized approach was implemented in the eXtended Finite Element Method. As a three-dimensional example, we study a slender solid characterized by an inclined interface across which the displacement is discontinuous. The accuracy is evaluated for varying position of the discontinuity interfaces with respect to the underlying mesh. A procedure for the choice of the regularization parameters is proposed.


Introduction
Regularization of singular and highly localized fields is a widely diffuse practice in several computational mechanics fields, such as thermal analysis, advection-diffusion problems [1], in Boundary Element Method applications [2], electronic structure calculations [3], for applications to modeling delamination in multilayered composites [4], and for capturing the behaviour of plastic hinges in Timoshenko beam elements [5].Among other methods [6,7], the eXtended Finite Element Method (XFEM) [6,[8][9][10] is a powerful technique for modelling engineering problems characterized by discontinuous and singular functions.XFEM is based on incorporating Heaviside and highly localized functions reproducing the expected shape of the interpolated fields.In this context, quadrature of terms containing discontinuous and singular functions is performed by subdividing the elements crossed by discontinuities in quadrature sub-domains [8], or, alternatively, by making use of adaptive quadrature [3,11].Subdivision into quadrature sub-domains can be cumbersome in the multidimensional case.Ventura [12] proposed an alternative technique for the integration of discontinuous integrands based on the use of equivalent polynomials.
Alternatively, regularization makes it possible to use standard Gauß quadrature and to avoid sub-domain quadrature.Regularized XFEM is also an effective alternative to nonlocal and cohesive crack-models [13][14][15][16].The key concept underlying regularization of discontinuous and highly localized functions is simple: the Dirac delta function is replaced by a weight function whose support is governed by a length parameter ρ, while the standard Heaviside function is constructed by integration of such a weight function.For vanishing ρ, the Heaviside and delta functions are recovered [17].Regularization obviously introduces approximations in the computation.Tornberg et al. [18,19] proved that the total error of integration can be split into the sum of the analytical and the numerical (quadrature) error.Specifically, the analytical error is governed by the approximation of the delta function, while the numerical error is governed by the differentiability of the regularized delta function besides the order of the quadrature rule.
The level set method [20] is a highly successful computational technique for tracking the spatial evolution of three-dimensional interfaces.An attractive feature is that interfaces are implicitly described and do not need to be parameterized.The use of the level set method [20] leads to a straightforward extension of regularized Heaviside and delta functions to the multidimensional case.Albeit regularization can be carried out by means of compact support weight functions, in [1,3,19,21] it has been proved that three-dimensional extensions based on compact support weight functions can lead to inconsistent results, as errors do not reduce with discretization refinement.The problem is that the moment conditions no longer hold for any shift (dilation) of the weight function with respect to the quadrature grid.Hence non-compact support weight functions should be used.
The structure of the paper is briefly summarized.In Section 2, we define a displacement field discontinuous across an interface.The discontinuity is expressed by means of the Heaviside function.The strain obtained by differentiation displays a delta function dependence, i.e., it is singular.The necessity of numerically integrating these functions arises.In order to perform quadrature in a more straightforward way, we replace the delta function with a regularized version.Regularization is achieved by introducing a smooth weight function whose support is governed by a regularization length (Section 3).Following [22], we exploit the relationship between Fourier transforms of the weight functions and accuracy of the regularization procedure in Section 3.1.As compact support functions were shown to be non-satisfying, we introduce non-compact support weight functions in Section 3.3.The truncation error is discussed in Section 3.2.Section 4 contains a brief presentation of the regularization procedure based on XFEM.A three-dimensional bar characterized by a discontinuous interface is discussed as a benchmark test in Section 5.In particular, the case of inclined discontinuities with respect to the underlying mesh is shown to be sensitive to the truncated support width.The influence of the truncation of the non-compact support is illustrated.A procedure for choosing the regularization parameter and the truncated support width is proposed.

Regularization in the Level Set Context
Let the volume V be divided into V 1 and V 2 by the interface Γ.At any point x ∈ Γ, the interface is characterized by the normal vector n(x).The position of Γ is implicitly defined by means of the signed distance function with respect to Γ [20] We introduce the Heaviside function ≤ 0, and focus on discontinuous functions of the type where v(x) is the regular part of the field u, while j(x) denotes the amplitude of the jump of u(x) across Γ. Differentiation of the discontinuous vector function (2) leads to a singular function of the type where ∇ denotes differentiation.The delta function δ(d(x)) is introduced through the generalized gradient Our purpose is to numerically integrate terms containing either discontinuous functions such as Equation (2), or singular functions such as Equation (3).In order to avoid quadrature subdomains, we want to replace the Heaviside function H by a regularized Heaviside function H ρ and the delta function δ by a regularized version δ ρ .In the one-dimensional case, denoting with t the scalar distance from the discontinuity, regularization can be carried out by replacing the delta function δ(d(x)) with [18,23] where the weight function φ is a smooth function depending on a length ρ such that The normalization parameter V ρ depends on the type of the weight function and the dimension of the space [17].Examples of weight functions are shown in Table 1; other examples can be found in [24].The one-dimensional regularized Heaviside function is then obtained as A straightforward extension to the multi-dimensional case can be carried out by assuming the regularized delta function to be a function of the signed distance function d (1) In this paper, the level set definition ( 7) is adopted.
Table 1.Examples of weight functions.
weight function

Error Analysis for Regularized Delta and Heaviside Functions
Let F(x) be a discontinuous or singular function, and consider the problem of the exact integration of the product where v(x) is a sufficiently regular scalar function for the Equation ( 8) to make sense [18,23].Assume that the function F ρ is the regularized version of F, namely F ρ is either equal to the regularized Heaviside function H ρ or the regularized Delta function δ ρ .Our aim is to compute the integral Equation ( 8) through a quadrature procedure applied to its regularized version.Let be the analytical error when replacing the Heaviside function by its regularized version and be the numerical error related to the numerical quadrature.Tornberg [18] showed that the total error associated with the replacement of the exact integral Equation ( 8) with its approximated quadrature can be expressed as Following [18,25], analytical and numerical errors are competing source of error as • the analytical error E ρ,v (F ρ ) typically increases with increasing ρ ; • the numerical error E quad,v (F ρ ) V is large for ρ small compared with the mesh size and decreases for increasing ρ.
More precisely, the q−order of the moment condition influences the analytical error of the regularized delta approximation [18].Moreover, while the analytical error depends on the number of satisfied moment conditions, differentiability influences the numerical error.
In particular, consider an error analysis based on the progressive mesh refinement for fixed ρ.If the mesh size h is rather large compared with the regularization parameter ρ, the numerical error E quad,v (F ρ ) is expected to compete with the analytical error E ρ,v (F ρ ).Subsequently, when h/ρ decreases, the analytical error will be more evident while the numerical error decreases.
In the following sections we focus on the analytical and the numerical errors corresponding to regularizing the Delta function, namely F ρ = δ ρ .

Fourier Transforms
Let us consider a finite element spatial discretization of the geometrical domain where the representative mesh size is equal to h.In the case of regularized functions with compact support ℓ ρ we can write ℓ ρ = m h where m is a positive integer or non-integer number.Tornberg and Engquist [19] proved that, whenever m is a non-integer number in the one-dimensional case, or, alternatively, the two-dimensional and the three-dimensional cases are considered, the error Equation (11) does not decrease for decreasing mesh size.This drawback was discussed by taking into consideration the discrete form of the first moment condition-also denoted as mass condition-whose satisfaction is necessary in order for the regularized formulation to preserve homogeneous fields.Zahedi and Tornberg [22] showed that, even in the one-dimensional case, simple compact support functions such as the hat function or the cosine function lead to a satisfying convergence in the one-dimensional case when m = 1 and half the support is taken equal to m h.However, compact support functions fail to satisfy the mass condition for any non-integer m.In the two-dimensional and the three-dimensional case, this limitation becomes more evident [19].This is because the number of quadrature points within the support is not constant, but it changes depending on the discretization and the geometry of the interface with respect to the mesh.It should be emphasized that, though effective, we are not focusing on mesh adaptive schemes [3].We consider instead meshes that do not necessarily conform to the interface position, in line with the spirit of the XFEM [24].If we choose indeed compact support functions, the width of the support will depend on the position of the interface with respect to the mesh.Consequently, the mesh has to conform to the weight function.
In order to discuss the accuracy of the quadrature procedure, we exploit the results recently obtained by [22] based on the Fourier transforms φ of φ Conversely, we have that Under the hypothesis that the space is covered by a regular grid {x j }, with j ∈ Z, of quadrature points, Zahedi and Tornberg [22] showed that the discrete version of the moment conditions ( 12) is related to the derivatives of the Fourier transform φ as follows Now, suppose that φ(κ) ≈ 0 for κ ≥ β; consequently, the last term on the r.h.s. of Equation ( 15) turns out to vanish for any m ≥ λ.Hence if the conditions are satisfied for r = 0, 1, . . ., q − 1, then δ ρ satisfies q discrete moment conditions.In particular, if is satisfied [19].Moreover, the derivatives of the Fourier transform and the moment conditions on φ are related by the relationship hold for any r > 1.Note that the Fourier transform of a compact support weight function is generally non-compact [22].

Truncation Error
Considering that the support of a non-compact support function will be truncated at a certain distance from the interface Γ, truncation will introduce a further source of approximation.However, it can be shown that the analytical error decreases for decreasing mesh size h while this does not generally happen in the case of compact support functions.However, the analytical error mainly depends on the truncation error, which is introduced in neglecting the contributions of quadrature points placed beyond a certain distance [22].Hence we can fix ρ independently of the mesh size and increase the support until the numerical error is below an acceptable value.
To recapitulate, we have two possibilities: 1. we choose weight functions with compact support equal to the finite element size multiplied by a suitable integer number m that depends on h, namely ℓ ρ ∝ m h [1,3].This means that whenever the mesh changes, the weight function has to be changed accordingly, otherwise non-vanishing errors for decreasing mesh size are obtained.
2. we choose non-compact support weight functions whose Fourier transform φ satisfies differentiability conditions; the length ℓ ρ will denote the truncated support of φ.This makes it possible to have ℓ ρ independent of h when ℓ ρ ≥ γ φ h, where γ φ is the width of the truncated support of the Fourier transform φ.
We want to investigate option 2. This requires evaluating the error due to truncation of the weight function.Denote the truncated delta function δ ω ρ as where w = ℓ ρ /2, ℓ ρ being the width of the truncated support of the weight function.Tornberg and Engquist [19] showed that the analytical error E ω associated with truncation can be neglected provided that Consequently, Therefore Equation (23) implies that the truncated weight function satisfies the moment condition (21b).
In order to prove condition (21a), we estimate For instance taking φ = e − |t| ρ leads to I 1 = I 2 = ρ e − ω ρ .Hence the analytical error decreases for increasing support width ℓ ρ = 2ω.The maximum value of the error, equal to ρ, is reached for ω = 0.In the following we will omit ω and use δ ρ to indicate the truncated regularized delta.

Example
In Table 2, we compare the Fourier transforms and the corresponding derivatives of modified versions φL and φG of the exponential and the Gauß weight functions, respectively.The Fourier transform was calculated using the forward definition (13) [26].The weight function φL (t) = (1/2)e −|ξ| where ξ = t ρ and its Fourier transform φL (κ) are compared in Figure 1.We can note that φL (κ) ≤ 10 −4 for κ ≥ 20, which suggests that we can reach at maximum the accuracy 10 −4 if we truncate the weight function φ L assuming ρ ≥ 20 h.The value 20 h provides an upper bound for ρ, any smaller value of ρ should be discarded.
Table 2. Fourier transforms of two weight functions.According to the previous results, φL and φG satisfy the moment conditions up to order 2. We have checked that the results pertaining to the width of the support of the simplified weight functions φL and φG can be extended to the weight functions φ L and φ G .In the numerical analysis, φ L = e −|ξ| is exploited, the discussion of the case φ G is part of a research in progress.

Basic Statements of the Regularized XFEM Formulation
We consider a regularized XFEM developed for strain localization and cohesive interfaces [13,17,27].The displacement field is approximated by enriching the space of the shape functions with additional functions that fulfill the partition of unity property of the standard finite element method.The enrichment functions consist of the product of the Heaviside function by the standard shape functions N u(x) = N(x)V + H(d(x))N(x)J (25) In Equation ( 25), V is the vector containing the degrees of freedom of the standard part of the displacement field v, while J collects the enriched degrees of freedom related to the displacement jump across the interface.The strain field is computed via compatibility from Equation (25) as where B is the compatibility matrix, N is related to the term in Equation ( 3) containing the (generalized) gradient ∇H = δn.Following the methodology proposed in [13,17] regularization is carried out.The regularized displacement and strain fields become where Let the material occupying volume V be elastic with Young's modulus E. The regularized interface is characterized by the Young's modulus by unit length Ē/δ ρ = E/(βδ ρ ), where β < 1 and β > 1 mean that the regularized interface is weaker or stronger, respectively, than the bulk, while for β = 1 the Young's moduli of the interface and the bulk coincide.The dimensional consistency and the motivations for assuming Ē in this way were thoroughly discussed in [17].Let us omit the spatial dependence.It was proved in [17] that the internal virtual work has to be assumed as where the constitutive laws are adopted.Equation (29) can be written in a more compact form as where the stiffness matrix collects the following submatrices The model has been implemented in X3D c ⃝ , a Fortran three-dimensional code based on the eXtended Finite Element method.All the submatrices (33) were computed by means of standard Gauß quadrature.More details on the quadrature procedure are given in Section 5.

Numerical Results
The accuracy of the proposed model is evaluated by considering the slender three-dimensional solid shown in Figure 2 with length L = 2 cm, cross section A = H × B, H = B = L/10.The element is clamped at one end and subjected to load control at the free end.An interface at L/2 is placed.Meshes of hexahedra are considered consisting of at least 5 × 5 elements in the cross section and increasing number of elements along the loading direction (see Figure 3).Enriched hexahedra possess a quadrature grid of 8 3 Gauß points, while non-enriched hexahedra contain 2 3 Gauß points.The analytical solution predicts a uniform state of stress along the loading direction.The uniform norm of the error, namely the maximal error on the stress component is evaluated [23].The parameter β characterizing the ratio between the Young's modulus of the bulk and that of the interface is set equal to 10 −4 .Tables 3 and 4 illustrate the errors obtained in the case of vertical interface with α = 0 • for fixed ρ and fixed h, respectively.In both cases, the support width ℓ ρ = 20ρ was adopted.The ratio ρ/h is incremented, in the former case by decreasing h, and in the latter case by increasing ρ.The same analysis is repeated for one case of inclined interface, namely α = 30 • .The corresponding results are reported in Tables 5 and 6.In Table 5, the mesh of 255 elements along the length was refined in the cross section in order to reduce mesh distortion [28].In particular, in Table 6, the influence on the numerical error of the support width is illustrated.While here, unlike in the case α = 0 • , the support width ℓ ρ = 20ρ leads to non-satisfying results, the minimal error obtained by adopting a support ℓ ρ = 40ρ is of the order 10 −3 , while the quality of the convergence is non-satisfying as E oscillates.Conversely, adoption of ℓ ρ = 60ρ improves convergence.Note that in Tables 5 and 6, a smaller value of the ratio ρ/h is used than in Table 5, because larger ρ may lead to ℓ ρ larger than the length L. The map of the stress error in percentage for varying truncated support width are displayed in Figures 5 and 6, corresponding to ℓ ρ = 40ρ and ℓ ρ = 60ρ, respectively.In the former figure, the total error is dominated by the truncation error, because it is localized at the enriched volume boundary.In the latter figure, the truncation error disappears, while the relevant error is localized along the inclined interface.In both figures, the ratio ρ/h = 5 is used.Consequently, the enrichment involves several layers of elements.This reflects the fact that the displacement jump across the discontinuity is smeared over a width larger than the representative element size, as shown in Figure 7.We have obtained a level-set dependent delta function approximation δ ρ .The function δ ρ satisfies three moment conditions up to a precision C when φ(κ) ≤ C for any κ ≥ C.This circumstance occurs if the ratio ρ/h is greater than 20.Moreover, the width of the support width ℓ ρ should be chosen at least equal to 40ρ in order to contain the truncation error.Hence, we can conclude that, generally, the support width should be larger than 40ρ.

Algorithm
An algorithmic scheme taking into account the previous analysis is the following: Step 1 Fix the tolerance parameters C, C E .
Step 4 Perform the analysis and evaluate E.
end while

Conclusions
In this paper, Gauß quadrature was used in a three-dimensional XFEM analysis of regularized interfaces.The relationship between Fourier transforms of the weight function and the analytical error was exploited.The main advantages are the simplicity of the approach and the fact that adaptive remeshing and quadrature are avoided.The numerical accuracy of regularized non-compact support was discussed in a benchmark test consisting of a 3D slender solid subjected to a tensile loading.The truncated support width ℓ ρ and the ratio ρ/h are the main parameters to be considered for controlling the approximation.It was shown that non-compact support weight functions make it possible to get a reasonable accuracy of the results.In particular: • truncation of the support width ℓ ρ influences the analytical error, as the moment conditions (12) are not exactly satisfied; • the computed error E decreases for increasing values of the ratio ρ/h, that is either increasing ρ for fixed h or decreasing h for fixed ρ; • discontinuities and singularities along interfaces that are non-aligned with the mesh are prone to severe errors when the support width ℓ ρ < 40ρ.
The results shown in this contribution are part of an ongoing research on general quadrature procedures for discontinuous and singular functions in three-dimensional structural problems.

Figure 2 .
Figure 2. Geometry of the uniaxial test.
Figure 4 displays the error in percentage for the case of ρ = L/100, ρ/h = 5 and the support width ℓ ρ = 20ρ.