## 1 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. [PERSON] [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 Gauss 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 \(\rho\), while the standard Heaviside function is constructed by integration of such a weight function. For vanishing \(\rho\), the Heaviside and delta functions are recovered [17]. Regularization obviously introduces approximations in the computation. [PERSON] _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.