Mollified finite element approximants of arbitrary order and smoothness

by   Eky Febrianto, et al.
University of Cambridge

The approximation properties of the finite element method can often be substantially improved by choosing smooth high-order basis functions. It is extremely difficult to devise such basis functions for partitions consisting of arbitrarily shaped polytopes. We propose the mollified basis functions of arbitrary order and smoothness for partitions consisting of convex polytopes. On each polytope an independent local polynomial approximant of arbitrary order is assumed. The basis functions are defined as the convolutions of the local approximants with a mollifier. The mollifier is chosen to be smooth, to have a compact support and a unit volume. The approximation properties of the obtained basis functions are governed by the local polynomial approximation order and mollifier smoothness. The convolution integrals are evaluated numerically first by computing the boolean intersection between the mollifier and the polytope and then applying the divergence theorem to reduce the dimension of the integrals. The support of a basis function is given as the Minkowski sum of the respective polytope and the mollifier. The breakpoints of the basis functions, i.e. locations with non-infinite smoothness, are not necessarily aligned with polytope boundaries. Furthermore, the basis functions are not boundary interpolating so that we apply boundary conditions with the non-symmetric Nitsche method as in immersed/embedded finite elements. The presented numerical examples confirm the optimal convergence of the proposed approximation scheme for Poisson and elasticity problems.



There are no comments yet.


page 3

page 5

page 7

page 8

page 9

page 13

page 15

page 19


A simple third order compact finite element method for 1D BVP

A simple third order compact finite element method is proposed for one-d...

Trefftz Finite Elements on Curvilinear Polygons

We present a Trefftz-type finite element method on meshes consisting of ...

A numerical study of the dispersion and dissipation properties of virtual element methods for the Helmholtz problem

We study numerically the dispersion and dissipation properties of the pl...

De Rham compatible Deep Neural Networks

We construct several classes of neural networks with ReLU and BiSU (Bina...

The Hellan-Herrmann-Johnson method with curved elements

We study the finite element approximation of the Kirchhoff plate equatio...
This week in AI

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

1 Introduction

Smooth high-order finite element approximants are often more efficient and, in general, integrate better with prevalent computer-aided geometric design (CAGD) descriptions Cirak et al. (2000); Hughes et al. (2005); Evans et al. (2009)

. The construction of mesh-based smooth high-order approximants is currently an active area of research as partly motivated by recent academic and industrial interest in isogeometric analysis. In most mesh-based approaches such approximants are defined as the tensor-products of univariate approximants. Such constructions do not generalise to unstructured meshes and auxiliary techniques are needed in the vicinity of the so-called extraordinary vertices where the tensor-product structure breaks. In CAGD a range of ingenious constructions has been conceived to generate smooth high-order approximants around the extraordinary vertices, see the monographs 

Farin (2002); Peters and Reif (2008) for an overview. Unfortunately, most of these constructions, including Scott et al. (2013); Jia et al. (2013); Majeed and Cirak (2017); Toshniwal et al. (2017b, a); Zhang et al. (2018); Kapl et al. (2018), target bivariate manifolds and do not generalise to the arbitrary variate case. Indeed, there are currently no sufficiently flexible and intuitive non-tensor-product arbitrary-variate constructions that can yield smooth polynomial high-order basis functions. By contrast, the proposed mollified approximation scheme over polytopic partitions is easy to construct, is polynomial and can have arbitrary order and smoothness.

Convolutional techniques are widely used in the analysis and numerics of partial differential equations. The convolution of a function with a mollifier, i.e. a specific type of kernel, yields a function that is smoother than the mollifier and the original function. This smoothing property is, for instance, used to recursively define uniform B-splines 

De Boor (1986), to analyse non-smooth functions and partial differential equations Hilbert (1973); Evans (1998); Adams and Fournier (2003), to postprocess finite element solutions Thomée (1977); Bramble and Schatz (1977); Mirzaee et al. (2013) or to regularise optimisation problems Sigmund and Petersson (1998); Le et al. (2011); Bletzinger (2014). Indeed, some of the classical meshless methods, like the smooth particle dynamics (SPH) Lucy (1977); Gingold and Monaghan (1977) and the reproducing kernel particle method (RPKM) Liu et al. (1995), are defined via convolutions, see also the reviews Li and Liu (2002); Bessa et al. (2014); Chen et al. (2017); Huerta et al. (2018). Although these methods were conceived as meshless, it is possible to define their mesh-based cousins as well Liu et al. (2004). Different from our mollified approximation scheme, SPH is intrinsically restricted to low order approximants and the RPKM yields high-order and arbitrarily smooth approximants which are rational. The kernels used in RPKM depend on the local node distribution and are determined so that they can exactly reproduce a polynomial of a given order.

Figure 1: Illustrative three-dimensional finite element computation using mollified basis functions. The domain boundary is described with the triangular mesh in (a) and the domain is partitioned with the Voronoi tessellation in (b). The solution of a Poisson problem is shown in (c). Note that in (b) and (c) the cells intersected by the domain boundary have been clipped and others omitted for visualisation purposes.

In the proposed mollified approximation scheme each non-overlapping polytopic cell has an independent local polynomial approximant of prescribed degree . The local approximants are discontinuous across cell boundaries. The convolution of a local approximant with a mollifier yields a smoother  approximant. The chosen mollifiers are compactly supported symmetric polynomials and have a unit volume. It is clear that the convolution of a global polynomial  of degree  with the chosen mollifiers gives a polynomial . However, it is straightforward to find a polynomial  of degree  such that . This implies that the mollified approximants can exactly reproduce any global function  of degree . Mollified basis functions for finite element analysis are defined by convolving the local approximant of each cell individually. At a given evaluation point the basis functions are evaluated first by computing the intersection between the support of the mollifier and the cell. Subsequently, the convolution integral over the resulting intersection polytope is evaluated numerically, but exactly (up to round-off error). We apply the divergence theorem to reduce the dimension of the integrals, but other methods of integrating polynomials over polytopes can be used, see e.g. Chin et al. (2015); Sudhakar et al. (2014). The obtained basis functions consist of several  continuously joined polynomial pieces. It is worth emphasising that we use, in contrast to RPKM, a fixed kernel which does not depend on the local node distribution and its convolution with a polynomial  of degree  is not required to yield the same polynomial, i.e. .

The derived mollified basis functions can be used as usual in the finite element discretisation of partial differential equations. For ease and efficiency of implementation we assume that each of the polytopic cells representing a finite element is convex. The required geometric operations, like the intersection computations, are significantly simplified by the convexity assumption. We partition the problem domain into a set of convex polytopic cells using a Voronoi diagram, see Figure 1. The evaluation of the finite element integrals requires some care because the  continuous breakpoints of the basis functions are not aligned with cell boundaries. We use to this end the variationally consistent integration approach proposed in Chen et al. (2017), which significantly reduces the number of required integration points. The required support of the mollified basis functions is the Minkowski sum of the mollifier support with the respective cell de Berg et al. (2010). Furthermore, the present version of the mollified basis functions are non-boundary-interpolating so that the Dirichlet boundary conditions are applied weakly with the non-symmetric Nitsche method Nitsche (1971); Oden et al. (1998); Boiveau and Burman (2015); Schillinger et al. (2016) as in immersed/embedded finite element methods, see e.g. Rüberg and Cirak (2012).

The outline of this paper is as follows. In Section 2 we briefly review the convolution of univariate polynomials with a mollifier and characterise the properties of the resulting mollified polynomials. Subsequently we derive in Sections 3.1 and 3.2 first the univariate and then the multivariate mollified basis functions. The key difference between the two cases lies in the evaluation of the convolution integrals. In the univariate case the integrals are evaluated analytically and in the multivariate case numerically. The use of the derived mollified basis functions in finite element analysis, especially the integration and treatment of boundary conditions, is discussed in Section 4. Finally, in Section 5

we introduce several Poisson and elasticity examples to confirm the optimal convergence of the developed approach. The paper is supplemented by four appendices which provide convergence estimates and discuss implementation details.

2 Preliminaries

We consider the one-dimensional domain  partitioned into a set of  non-overlapping segments , referred to as cells, such that


On each cell  a compactly supported local polynomial is defined,


where the vector 

represents a polynomial basis of degree  and  are its polynomial coefficients. See Figure 2 for an illustrative example. We choose in each cell the same polynomial basis , although it is possible to change the type of basis and its polynomial degree. The sum of the local polynomials defined over the entire domain  is given by


Evidently, across the cell boundaries this function can be discontinuous, i.e. .

Figure 2: Mollification of piecewise discontinuous functions (black, dashed) with a linear mollifier (red, solid). The resulting mollified functions  (blue, solid) are continuous.

The smoothness of  is increased by convolving it with a mollifier . The mollifier is chosen such that it has the following properties


That is, the mollifier is positive, has a unit volume and a finite support of size . In addition we require that the mollifier is symmetric, i.e. , and that it has certain smoothness properties, as yet to be specified. The mollification of is defined with the convolution


The equality of both integrals can be shown by a simple substitution. We usually use the first integral expression in the following. Furthermore, we choose polynomial mollifiers  of degree  and  is, as stated above, of degree . The mollified function  has monomials up to degree .

If the derivative of the mollifier  exist, the derivative of the mollified function  is given by


The higher order derivatives are computed similarly. Considering that  is discontinuous and the mollifier is  we can deduce for the mollified smooth function  .

Finally, the cell-wise definition of  introduced in (5) and (6) yields




The function  is composed of infinitely smooth polynomial pieces that are smoothly connected, i.e. , across a finite number of breakpoints. In general the location of the breakpoints does not coincide with the cell boundaries. However, as is known from B-splines, on uniformly partitioned domains, it is possible to choose the support size of the mollifier such that the breakpoints fall on the cell boundaries De Boor (1986); Sabin (2010). In the illustrative example in Figure 2 the influence of the choice of the local polynomial  on the mollified function  is demonstrated. In this example the convolution integral (7) has been evaluated analytically.

A final remark concerns the reproduction of polynomials with mollified functions. It can be shown that it is possible to find for a given polynomial  of degree  a polynomial  of the same degree which yields after mollification . Specifically, the mollification of a polynomial


is given by



are the moments of the mollifier defined as


After rearranging the terms in (10) according to the powers of  it is easy to see how to choose a function  which is after mollification equal to the given function .

3 Smooth piecewise basis functions

We now use the mollification approach to derive basis functions on one- and multi-dimensional domains. Again the domain  is assumed to be partitioned into a set of non-overlapping convex polytopes , which are in the present paper obtained from a Voronoi diagram. In the following we refer to the polytopes as cells. The mollification approach yields a set of basis functions for each cell. The convolution integrals for obtaining the basis functions are evaluated analytically in the one-dimensional, i.e. univariate, case with and numerically in the multi-dimensional, i.e. multivariate, case with . Convergence estimates for the obtained mollified basis are provided in A.

3.1 Univariate basis functions

To set the stage for multivariate basis functions, we first consider the derivation of univariate basis functions. The mollified basis functions belonging to a cell  are defined according to the mollification (7) by


where the vector of mollified basis functions are defined with


Note that the local polynomial basis is outside the cell  zero and the integration domain is restricted to . As a local polynomial basis  different basis choices are possible, such as the monomial, Lagrange or Bernstein. While this choice has no influence on the approximation quality of the resulting mollified basis, it affects the interpretability of the coefficients  and the conditioning of the resulting finite element system matrices. In our examples we use in each cell  a scaled monomial basis


where  is its degree,  is the centre of the cell and  is the average length of all the cells  in the domain. The scaling by ensures that the obtained mollified basis functions have a similar maximum value. It is possible to apply a different scaling factor and to choose the scaling in each cell differently.

For any given point  in the domain the mollified basis functions  are evaluated by computing the convolution integral (13). Evidently, when the chosen local polynomial basis  is a monomial basis the mollified basis functions are simply the moments of the mollifier. The derivatives of the mollified basis functions are computed according to (6). In Figure 3 the mollified basis functions for a cubic monomial basis with  and a piecewise linear mollifier  with two different support sizes  are shown. In each case the support size of the mollified basis functions is  with the mollifier size  and the cell size . The basis functions consist of several polynomial pieces that are continuously connected at the breakpoints. The number of breakpoints in each cell depends on the number and arrangement of the breakpoints in the mollifier. Note, although not shown in the figure, the basis functions of the neighbouring cells are non-zero in the considered cell as well. The breakpoints of those neighbouring basis functions may not coincide with the breakpoints of the shown basis functions.

Figure 3: Univariate mollified basis functions with a bilinear mollifier () and a qubic scaled monomial basis () on a cell . The mollified basis functions in (b) are obtained with the narrow mollifier in (a) and the ones in (d) with the wide mollifier in (c). The dashed lines in (b) and (d) indicate the continuous breakpoints of the obtained basis functions.

3.2 Multivariate basis functions

Without loss of generality we focus in the following on trivariate basis functions. As in the univariate case the basis functions for a cell  are given by


The vector  contains the multivariate monomial basis of degree . In this paper, the mollifier


is composed of -continuous quartic splines


It is easy to verify that  and . Hence, the mollifier is continuous. The continuity of the mollifier can be increased by forcing more derivatives to be zero at  and , which can be achieved either by choosing a higher order polynomial, using a non-polynomial function or introducing more breakpoints. Obviously all approaches increase the cost of evaluating the convolution integrals.

Figure 4: -continuous quartic spline mollifier with the support size

The availability of efficient evaluation techniques for the multi-dimensional convolution integral (15) are vital to the proposed approach. We note that the integrand is a polynomial and that it is only non-zero on the intersection, i.e. boolean intersection, of the support of the mollifier and the considered cell, i.e.,


where  denotes the support of the mollifier centred at the evaluation point , see Figure 5. The intersection domain  is convex because both   and  are convex. Computing the intersection of polytopes and the integration of polynomials over polytopes are recurring tasks in computer graphics and many robust and efficient algorithms are available. In the C we introduce one such algorithm for determining the intersection between a cell and a mollifier. One possible approach to evaluate the integrals on  is first to tesselate it and then to integrate over the obtained simplices using Gaussian quadrature. Considering that the integrands are polynomials and the tessellation consists of affinely mapped simplices, the integration can be performed exactly using a sufficient number of quadrature points.

Figure 5: Evaluation of mollified basis functions  by numerically computing the convolution integral (15) at a point . (a) The support of the mollifier  centred at  (in red) and the cell (in blue). (b) The integration domain for the convolution integral is the boolean difference .

A more elegant approach is to reduce the domain integrals to line integrals by repeated application of the divergence theorem, see e.g. Mirtich (1996). We briefly sketch the conversion of volume integrals to surface integrals for completeness. The polytope  consists of a set of uniquely orientated faces , i.e. all the respective normals  point outside the domain. To integrate an arbitrary polynomial  it is first integrated, e.g. in the  direction,


The divergence theorem applied to this new function yields


where the surface normal  of the face  is constant. It is possible to stop at this point and to numerically evaluate the surface integrals after triangulating the faces . The number of quadrature points on each face is chosen exactly to integrate polynomials of degree . However, it is also possible to reapply the divergence theorem to reduce the surface integrals to line integrals, which can then be analytically evaluated. In our implementation we evaluate (20) numerically by triangulating the faces , as indicated in Figure LABEL:fig:basis2Db. Note that the change from volume to surface integrals already yields a significant reduction in the number of integration points.

As in the univariate case, the support of a mollified basis function  is larger than its respective cell. In finite element computations also the support of a basis function is required, which is given by the Minkowski addition


In D we introduce an algorithm for computing the Minkowski sum of two convex polyhedra.

Figure 6: A cutout from a Voronoi tesselation and the support of the mollifier. The mollified basis functions  of the highlighted cell  are supported on the larger polygonal domain . The polygonal domain is the Minkowski sum, i.e. , of the cell with the mollifier support .

4 Finite element discretisation with mollified basis functions

The smoothness and approximation properties of the mollified basis functions make them ideal for finite element analysis. In the following, we briefly outline the discretisation of a Poisson equation using the mollified basis functions. As in the previous sections we assume that a partitioning of the domain consisting of convex polytopes is given. We generate such a partitioning using a Voronoi diagram from a given implicit, i.e. level set, or parametric, i.e. surface mesh, description of the domain boundary. A set of points is placed within and outside the domain to generate the Voronoi diagram. The points outside the domain ensure that the mollifier’s support is fully covered by local polynomials when the mollifier is placed on the domain boundary. The Voronoi diagram and the respective mollified basis functions do not conform to the domain boundaries. Therefore, the boundary conditions are applied weakly on cells cut by the boundary, as in immersed or embedded finite elements.

Notwithstanding this, it appears to be possible to continuously shrink the mollifier support when approaching the boundaries from the inside. In the limit on the boundary the mollifier becomes a Dirac delta and the interpolation within the domain becomes independent of the outside of the domain. This idea has, however, not been further pursued in this paper in order to focus on other aspects of the method.

The Poisson equation on a domain  is given by

in (22)

where is the prescribed solution field on the Dirichlet boundary  and  is the prescribed flux on the Neumann boundary  with the outward normal . The weak formulation of the Poisson equation can be stated according to Nitsche Nitsche (1971) as: Find such that


The stabilisation parameter can be set to  when the sign of the last term  is reversed, as proposed in a number of papers Oden et al. (1998); Boiveau and Burman (2015); Schillinger et al. (2016). In our computations we use this so-called non-symmetric Nitsche method which has a non-symmetrical system matrix. The trial and test functions are discretised with the mollified basis functions


where  is the number of the polytopic cells in the mesh.

Introducing the interpolation equations (24) into the weak form (23) yields a linear system of equations with the unknowns , which are the coefficients of the local monomial basis in the cells. For instance, the bilinear form  becomes after discretisation


As usual, the domain integral is evaluated numerically after splitting it into cell contributions


We evaluate the integral over a cell  by first decomposing it into tetrahedra and then applying standard Gauss integration. A cell is tetrahedralised by introducing additional nodes at its centre and face centres. At each integration point the mollified basis functions are evaluated as described in Section 3. Although the sketched integration of the weak form (26) is straightforward, Gauss integration unfortunately requires too many quadrature points because of the breakpoints of the basis functions within the cells. Using too few quadrature points usually leads to suboptimal convergence rates.

The principal difficulties encountered in efficient and accurate integration of the weak form (26) are very similar to those encountered in meshless methods. In the variationally consistent integration techniques for meshless methods the test functions are modified to satisfy a consistency condition when integrated numerically, see e.g. Chen et al. (2013); Hillman et al. (2015). If the basis functions can reproduce a solution  of polynomial degree , the finite element scheme must be able exactly to solve such a problem even when the integrals are evaluated numerically. Assuming a problem with the solution of polynomial degree  and inserting it into the discretised weak form yields for each cell a consistency condition for integration. To satisfy this consistency condition in a cell  the gradient of the test functions is modified as


where  is a vector of all ones and the vector  contains a monomial basis of degree  with the yet to be determined matrix of coefficients . The support of the basis functions  is chosen to be same as of the basis functions . The matrix of coefficients  is determined by solving per basis function one small linear equation system so that variational consistency condition is satisfied. This equation system contains the integrals of the basis functions over their supports, which is given by the Minkowski sum of the mollifier support and the cell. See Chen et al. (2013) for further details.

In the cells cut by the domain boundary the element integrals are evaluated only over parts of the cell which lie inside the domain. The respective integration domains are obtained by clipping the cells, see B. The resulting polyhedron is tetrahedralised with the same approach used for a non-clipped cell. As in immersed, or embedded domain, methods the faces of the tetrahedra can be projected to the curved domain boundaries if higher-order boundary approximation is needed, see e.g. Cheng and Fries (2010); Kudela et al. (2015); Xiao et al. (2019).

5 Examples

We introduce in this section several examples of increasing complexity to experimentally verify the convergence of the proposed mollified finite element approach. In the one-dimensional problems both the convolution and finite element integrals are evaluated analytically, whereas in multi-dimensional problems both integrals are evaluated numerically. In multi-dimensional problems only the quartic spline mollifier (17) consisting of a single polynomial with no internal breakpoints is used. The mollified basis functions contain monomials of up to degree , although they are only complete up to degree , and are non-zero over several cells. Therefore, it is not obvious how many quadrature points to choose in each of the integration triangles used for evaluating the finite element integrals. In our present computations we determine a stable number of quadrature points by successively increasing the number of quadrature points until we have a stable solution. In two dimensional problems we choose for linear basis functions () three integration points for domain integrals and five for boundary integrals, and for quadratic () we choose four and five integration points respectively. In all problems the Dirichlet boundary conditions are enforced with the parameter-free non-symmetric Nitsche method.

5.1 One-dimensional Poisson problem

As a first example we consider the solution of the one-dimensional Poisson-Dirichlet problem on the domain . The source term  is chosen such that the solution is equal to


The initial coarse mesh consisting of  cells is chosen to be non-uniform. The cell sizes, starting from the left, are , , , , and

. In addition to these cells, each domain boundary is padded with an extra ghost cell to ensure that the obtained mollified basis functions have the same approximation properties over the entire domain. We obtain finer meshes by repeated bisectioning of all cells.

In the following set of experiments we study the influence of the choice of the local polynomial basis and the mollifier  on the convergence of the finite element solution. Firstly, we take different local polynomials of degrees  and a normalised linear B-spline, i.e. hat function, mollifier with a support width of


Note that the normalisation of the mollifier is essential to ensuring that it integrates as required to one. Figure 7 shows the optimal convergence of the mollified finite element approach in the norm and seminorm with and respectively. These convergence rates are in agreement with the analytic estimates provided in A.

(a) -norm error
(b) -seminorm error
Figure 7: One-dimensional Poisson problem. Convergence with normalised linear B-spline mollifier and local polynomial basis of degrees .

Next, we investigate the influence of the mollifier support width on the convergence order while keeping the normalised linear B-spline mollifier. The mollifier width is chosen according to


The increase in mollifier size leads to an increase in the support size of the mollified basis functions, which results in an increase of the number of non-zero basis functions in a cell. The obtained optimal convergence rates for  are shown in Figure 8. The increase in mollifier width leads to a somewhat decrease in the convergence constants, but the optimal support size appears to depend on the specific problem considered. The results for higher order local polynomials are similar and have been not included here.

(a) -norm error
(b) -seminorm error
Figure 8: One-dimensional Poisson problem. Convergence with local polynomial basis of degree and normalised linear B-spline mollifier with different support sizes of .

Finally, we study the effect of the mollifier smoothness on finite element convergence. The normalised B-spline mollifiers are of degree and the local polynomial is of degree . The mollifier width factor is chosen as . Note that the B-spline mollifiers are  continuous so that the obtained mollified basis functions are  continuous. Figure 9 shows that an increase of the kernel degree and smoothness does not have an effect on the optimal convergence rate, but leads to a significant decrease in convergence constants.

(a) -norm error
(b) -seminorm error
Figure 9: One-dimensional Poisson problem. Convergence with local polynomial basis of degree and normalised B-spline mollifiers of degrees .

5.2 Two-dimensional examples

5.2.1 Poisson problem on a square domain

We consider next the Poisson-Dirichlet problem on a square domain . The domain is partitioned into cells using the Voronoi diagram of 

non-uniformly distributed points, see Figure 

10. Starting from a set of uniformly distributed points we introduce non-uniformity by randomly perturbing their coordinates by with . Only the coordinates of points farther than a certain distance from the boundaries are perturbed.

As in the one-dimensional case the domain is padded with an extra layer of ghost cells (not shown in Figure 10) to ensure that the mollified basis functions are complete close to the boundaries. Depending on the number of Voronoi cells  the width of the continuous quartic spline mollifier is chosen with

Figure 10: Poisson problem on a square domain. Three different partitionings of the domain using the Voronoi diagram of  non-uniformly distributed points.

We firstly perform a patch test to verify that the mollified finite element method in combination with variationally consistent integration can exactly solve problems with polynomial degree . To this end, we consider on the mesh shown in Figure (b)b two problems with the exact solutions  and . Solving the linear problem using the linear mollified basis functions with  leads to an norm error of and an seminorm error of . The corresponding errors for the quadratic problem using quadratic mollified basis functions with  are and . This clearly confirms that the mollified finite element method satisfies the patch test.

With the consistency of the method confirmed, we proceed to establish its convergence under mesh refinement. The source term  is now chosen such that the solution is equal to


The used mollified basis functions are the continuous linear and quadratic basis functions with  and  respectively. Figure 11 shows the convergence of the errors in norm and seminorm as the mesh is refined. Note that the refined meshes are not nested so that some small kinks in the convergence curves may be expected. The average convergence rates are, however, close to optimal, as indicated by the dashed triangles in Figure 11.

(a) -norm error
(b) -seminorm error
Figure 11: Poisson problem on a square domain. Convergence with quartic spline mollifier and local polynomial basis of degree and .

5.2.2 Elastic plate with a hole

As a two-dimensional problem with a non-trivial geometry we compute the infinite elastic plate with a hole subjected to uniaxial tension. The radius of the hole is  and the applied uniaxial traction in the vertical direction is . The Young’s modulus of the material is and its Poisson’s ratio is . This problem has a known closed-form analytic solution Timoshenko (1970). Therefore, we discretise only the smaller plate of size shown in Figure 12 by applying Dirichlet boundary conditions over its entire boundary.

Figure 12: Geometry and boundary conditions of the elastic plate with a hole.

The initial mesh consists of a Voronoi diagram of 16 non-uniformly distributed points, see Figure 13. The refined meshes are obtained by subdividing the cells by introducing new vertices on the cell and edge centres. This refinement ensures that the meshes are nested. As can be inferred from Figure 13 along the circular boundary the mesh edges are not aligned with the boundary. In the respective cells cut by the boundary the finite element integrals are evaluated only over the cell areas inside the domain. The cut-cells for integration are obtained with the clipping process introduced in B. To achieve a higher order approximation the edges of the triangles used for integration are curved by introducing additional nodes on the faces. As in standard finite elements, to achieve an optimally convergent method the boundary geometry has to be approximated with the same polynomial order as the used mollified basis functions.

Figure 13: Elastic plate with a hole. Three different partitionings of the plate with a hole (highlighted in blue). The elements outside the domain are the ghost cells. The coarse mesh in (a) is a Voronoi diagram which is refined to obtain (b) and (c) by introducing new vertices on the cell and edge centres.

To analyse this problem we again use continuous linear and quadratic basis functions with and , respectively. Figure 14 shows the convergence of the errors in the energy norm. It is apparent that optimal convergence rates are achieved. A final comment concerns the possibly very small contributions of basis functions cut by the boundary to the system matrix. To this end, several approaches have been developed in immersed/embedded finite elements Rüberg and Cirak (2012); de Prenter et al. (2019); Gürkan and Massing (2019). A particularly simple approach is simply to scale the relevant basis functions according to their support size within the domain, i.e.

Figure 14: Elastic plate with a hole. Convergence of the relative energy norm error with quartic spline mollifier and local polynomial basis of degree  and .

5.3 Three-dimensional example

As an illustrative three-dimensional example with a complex boundary we consider the solution of a Poisson-Dirichlet problem on the domain contained in the Stanford bunny, see Figure 15. The geometry of the bunny is given as a triangle mesh with 66272 facets. The volume mesh shown in Figure (c)c is created in several steps. Firstly we introduce within the bounding box of the bunny a set of uniformly distributed points with each apart. The head is then refined by introducing additional points with each apart. Finally, one of the ears is refined by adding points with each apart. We then generate the Voronoi diagram of all the points and iteratively relax it to achieve a more even distribution of cell sizes. During the iterative relaxation the new position of each point is recomputed by convolving nodal coordinates with a box function. This relaxation is equivalent to standard Laplace smoothing of meshes. After the relaxation the Voronoi diagram is clipped with the technique described in B. The final mesh consists of 897 cells.

We solve on the generated polytopic mesh a Poisson-Dirichlet problem with a source term  such that the solution is equal to


The isocontours of the computed solution are shown in Figure (d)d. For the employed mollified basis functions we use a local polynomial basis with and a quartic spline mollifier with a support size of , which is twice the coarse cell size. The discretisation has in total  basis functions. The cut-cells are stabilised as described earlier by scaling the basis functions according to (33).

(a) Fine surface mesh
(b) Surface of the clipped Voronoi diagram
(c) Voronoi diagram
(d) Finite element solution
Figure 15: Poisson-Dirichlet problem on the domain contained within the Stanford bunny. The two boxes in (a) indicate the locations where additional nodes are introduced to refine the mesh. The mesh in (c), with some of the cells omitted, shows the unstructured Voronoi diagram used in the mollified finite element computation (d).

6 Conclusions

We introduced the mollified basis functions of arbitrary order and smoothness and verified their excellent finite element approximation properties with a selected set of examples. In the two- and three-dimensional examples we chose the Voronoi diagram of a given set of points as the partitioning of the domain. The mollified basis functions are obtained by convolving cell-wise defined local polynomial approximants with a compactly supported smooth mollifier with unit volume. We integrated the convolution integrals exactly (up to round-off errors) by first determining the geometry of the polytopic integration domain and then reducing the dimension of the integrals using the divergence theorem. Although we considered only one-, two- and three-dimensional domains, the proposed numerical integration approach to evaluating the convolution integrals applies to any dimension. Because the mollified basis functions are not boundary conforming we enforce boundary conditions with standard immersed/embedded finite element techniques. The obtained polynomial basis functions may have breakpoints, i.e. points or lines of reduced continuity, within the cells. Therefore, we evaluated the finite element integrals with a variationally consistent approach originally developed for meshless methods. As shown numerically and analytically the mollified basis functions in combination with the proposed finite element implementation can pass the patch test and achieve optimal rates of convergence.

There are several promising applications and extensions of the proposed mollified approximation scheme worth mentioning. Evidently, it is straightforward to apply h-, p- and hp-refinement. A given Voronoi diagram can be h-refined by incrementally adding new points and updating the Voronoi diagram. For p-refinement it is sufficient simply to choose in each cell the order of the polynomial approximant differently. A-priori and a-posteriori estimators are crucial to making efficient use of h-, p- and hp-refinement in applications. Furthermore, in our present implementation the mollifier support size is uniform throughout the domain. As our preliminary experiments indicate, it is possible to vary the mollifier support within a domain. This can be, for instance, used in creating boundary interpolating approximants by continuously shrinking the mollifier support size to zero while approaching the boundary. An alternative approach to easing the enforcement of boundary conditions is to blend mollified basis functions with finite elements similar to blending techniques developed for meshless methods Huerta and Fernández-Méndez (2000); Rosolen and Arroyo (2013). Lastly, returning to our original motivation in developing smooth basis functions for isogeometric analysis, it is appealing to develop mollified blending techniques for B-spline patches meeting at an extraordinary vertex. The convolutional definition of the B-splines can be used to derive mollified approximants which reduce to B-splines away from the mesh boundaries.

Appendix A Convergence estimates

We make use of convolution and polynomial approximation estimates to derive convergence estimates for the proposed mollified finite element approximation scheme Ern and Guermond (2004); Adams and Fournier (2003). In doing so, we adopt a multiindex notation and denote derivatives with . If the multiindex is an -tuple of nonnegative integers , then is a differential operator of order , with the convention that . We denote the norms of the Lebesgue and Sobolev spaces  and with


In this Appendix we assume that the mollifier is continuous. For instance, we may take


where the scalar is chosen so that  has unit volume. The scaled mollifier with the support size  is given by


so that


for any multiindex .

As shown, e.g. in Adams and Fournier (2003), the derivatives of the convolution satisfy for any pair of multiindices and the relation


and the Young’s inequality for convolutions reads


yielding the estimate


Moreover, according to the Bramble-Hilbert lemma, for a , , , there is a (polynomials of degree less than or equal to in one variable) such that


Let us now extend and by zero outside the domain and denote their mollifications with


Then, for , from (39), (41), (38) and (42) we have


where is a constant. The difference between a function and its mollification can be bounded with a standard approximation theorem for convolutions (Adams and Fournier (2003), Theorem 5.33) which reads


where is a constant. Combining the preceding two estimates we obtain


This bound is minimised by taking


which defines the optimal mollifier support size in dependence of the mesh size . Inserting the optimal  into the estimate (46) gives


where  denotes the mollification of the polynomial  with a mollifier with support sizel . We note that the estimate (48) provides control over -th order derivatives, whereas the initial estimate (42) does not.

As in standard finite element approximation theory, see e.g. Ern and Guermond (2004), applying estimate (48) cell-wise and considering their sum yields global convergence estimates. Subsequently, it is straightforward to confirm the optimal convergence of the proposed mollified finite elements as already suggested by our numerical experiments.

Appendix B Clipped Voronoi diagrams

We briefly review the properties of Voronoi diagrams and sketch the generation of clipped Voronoi diagrams that approximately fill a given domain . For a more detailed discussion see, e.g., Aurenhammer (1991); Du et al. (1999); de Berg et al. (2010). For a set of points in the Voronoi diagram is defined by a set of cells such that


As indicated in Figure 16, the cells are convex, are either bounded or unbounded and have planar faces. There are a number of efficient software libraries available for generating Voronoi diagrams, such as the Voro++ Rycroft (2009) library (for 3D) or Mathematica (for 2D) used in this work. To obtain a Voronoi diagram that approximately fills a given domain  the cells intersected by the boundary are clipped. To implement the clipping process we assume that the domain  is described implicitly with a signed distance function


where  is the boundary of the domain . Domains that are described with a parametric polygon mesh can first be converted to an implicit signed distance function representation using standard algorithms, see e.g. Rüberg and Cirak (2012).

(a) Polygon mesh
(b) Signed distance function
(c) Voronoi diagram
(d) Clipped Voronoi diagram
Figure 16: Clipped Voronoi diagram of a domain described by a polygonal mesh.

The minimal data structure for representing a Voronoi cell consists of its vertices and orientated faces , with all the face normals  pointing, e.g., outside the cell. The cells cut by the boundary are determined by evaluating the level set function at the vertices and checking whether


is satisfied. Each of the cut cells is clipped by performing the following steps:

  1. Deduce the set of cell edges from and .

  2. Determine the intersection points between edges and the boundary with a bisectioning algorithm.

  3. Introduce new vertices at the points determined in step 2.

  4. Generate a new clipped cell by determining the convex hull of vertices inside the domain and the vertices on the domain boundary.

In finite element computations the clipped cells represent the integration domain for element integrals for the cells touching the domain boundary. The clipping process introduced yields only convex clipped cells with planar boundaries. For domains with curved boundaries this limits the overall accuracy of finite element method to first order even when higher order mollified basis functions are used. In this setting, a standard approach to achieving higher order accuracy in immersed finite elements is to curve the planar faces by introducing additional vertices on the faces, which is clearly also applicable to mollified finite elements.

Appendix C Intersection of a convex cell with a box

The intersection between a convex cell and a box is required to evaluate the mollified basis functions. The box represents the support of the mollifier and is centred at the given evaluation point. The intersection of two convex solids is frequently required in computer graphics and a range of efficient and robust algorithms is available. Our specific approach is motivated by the more general algorithm presented in Muller and Preparata (1978). The key idea is that the intersection of a cell and a box can be determined by clipping the cell in turn by the six half-spaces defining the box. Each half-space is defined by a plane, or its respective normal, and has an inside and outside. Furthermore, recall that the intersection of two convex solids is always convex.

As mentioned in B, in our implementation, a cell  is represented by its vertices and orientated faces . With this in mind, the sequence of steps in computing the convex polytope representing the intersection domain is as follows:

  1. Deduce the set of cell edges from and .

  2. Determine the intersection points between the edges and the six half planes in turn while keeping track of the vertices inside the half-spaces.

  3. Generate a new polyhedron by determining the convex hull of vertices inside the domain and the intersection points on the edges.

In meshes with a large number of cells it is usually more efficient first to identify the small set of cells which are possibly intersected by a given box. Subsequent intersection computations have to be applied only to the few identified cells. The relevant cells can be efficiently identified with a standard hierarchical bounding volume tree, see e.g. Klosowski et al. (1998).

Appendix D Minkowski sum of two polytopes

The support of the mollified basis functions  corresponding to the cell  is obtained as the Minkowski sum of the cell with the support of the mollifier. For an in-depth introduction to Minkowski sums see e.g. de Berg et al. (2010). The Minkowski sum of two sets  is defined by


The domain  resulting from the Minkowski sum may be visualised as that obtained by sliding the centre of  along the boundaries of , see Figure 17. Recall here that the domain  is centred at the origin of the coordinate axis as implied by the subscript . It is easy to show that  is convex because both   and  are convex. We use this to devise a simple algorithm for computing the Minkowski sum. That is, we first generate a set of points by sliding the domain  along the boundaries of , which we subsequently combine with a convex hull algorithm to obtain . In generating the set of points it is sufficient to place  at the vertices  of the domain  and to take successively the union of the vertices of , see Figure 17 .

(a) and
Figure 17: Minkowski sum of the mollifier support  and the cell . The dashed  in (b) indicate the sliding of the centre of  along the boundaries of .


  • R.A. Adams and J.J.F. Fournier (2003) Sobolev spaces. Elsevier Science. Cited by: Appendix A, Appendix A, §1.
  • F. Aurenhammer (1991) Voronoi diagrams — a survey of a fundamental geometric data structure. ACM Comput. Surv. 23 (3), pp. 345–405. Cited by: Appendix B.
  • M. Bessa, J. T. Foster, T. Belytschko, and W. K. Liu (2014) A meshfree unification: reproducing kernel peridynamics. Computational Mechanics 53, pp. 1251–1264. Cited by: §1.
  • K.-U. Bletzinger (2014) A consistent frame for sensitivity filtering and the vertex assigned morphing of optimal shape. Structural and Multidisciplinary Optimization 49 (6), pp. 873–895. Cited by: §1.
  • T. Boiveau and E. Burman (2015) A penalty-free nitsche method for the weak imposition of boundary conditions in compressible and incompressible elasticity. IMA Journal of Numerical Analysis 36 (2), pp. 770–795. Cited by: §1, §4.
  • J. Bramble and A. Schatz (1977) Higher order local accuracy by averaging in the finite element method. Mathematics of Computation 31, pp. 94–111. Cited by: §1.
  • J. Chen, M. Hillman, and S. Chi (2017) Meshfree methods: progress made after 20 years. Journal of Engineering Mechanics 143 (4), pp. 04017001: 1–37. Cited by: §1, §1.
  • J. Chen, M. Hillman, and M. Rüter (2013) An arbitrary order variationally consistent integration for galerkin meshfree methods. International Journal for Numerical Methods in Engineering 95 (5), pp. 387–418. Cited by: §4.
  • K. W. Cheng and T. Fries (2010) Higher-order xfem for curved strong and weak discontinuities. International Journal for Numerical Methods in Engineering 82, pp. 564–590. Cited by: §4.
  • E. B. Chin, J. B. Lasserre, and N. Sukumar (2015) Numerical integration of homogeneous functions on convex and nonconvex polygons and polyhedra. Computational Mechanics 56, pp. 967–981. Cited by: §1.
  • F. Cirak, M. Ortiz, and P. Schröder (2000) Subdivision surfaces: A new paradigm for thin-shell finite-element analysis. International Journal for Numerical Methods in Engineering 47, pp. 2039–2072. Cited by: §1.
  • M. de Berg, O. Cheong, M. van Kreveld, and M. Overmars (2010) Computational geometry: algorithms and applications. 3rd edition, Springer. Cited by: Appendix B, Appendix D, §1.
  • C. De Boor (1986) B(asic)-spline basics.. Technical report Madison mathematics research center, Winsconsin University. Cited by: §1, §2.
  • F. de Prenter, C. V. Verhoosel, and E. H. van Brummelen (2019) Preconditioning immersed isogeometric finite element methods with application to flow problems. Computer Methods in Applied Mechanics and Engineering 348, pp. 604–631. Cited by: §5.2.2.
  • Q. Du, V. Faber, and M. Gunzburger (1999) Centroidal voronoi tessellations: applications and algorithms. SIAM Rev. 41 (4), pp. 637–676. Cited by: Appendix B.
  • A. Ern and J.L. Guermond (2004) Theory and practice of finite elements. Applied Mathematical Sciences, Springer New York. Cited by: Appendix A, Appendix A.
  • J. A. Evans, Y. Bazilevs, I. Babuska, and T. J. R. Hughes (2009) N-widths, sup-infs, and optimality ratios for the k-version of the isogeometric finite element method. Computer Methods in Applied Mechanics and Engineering 198, pp. 1726–1741. Cited by: §1.
  • L. C. Evans (1998) Partial differential equations. Graduate Studies in Mathematics, American Mathematical Society. Cited by: §1.
  • G. Farin (2002) Curves and surfaces for cagd. Academic Press. Cited by: §1.
  • R. A. Gingold and J. J. Monaghan (1977) Smoothed particle hydrodynamics: theory and application to non-spherical stars. Monthly notices of the royal astronomical society 181, pp. 375–389. Cited by: §1.
  • C. Gürkan and A. Massing (2019) A stabilized cut discontinuous galerkin framework for elliptic boundary value and interface problems. Computer Methods in Applied Mechanics and Engineering 348, pp. 466–499. Cited by: §5.2.2.
  • S. Hilbert (1973) A mollifier useful for approximations in sobolev spaces and some applications to approximating solutions of differential equations. Mathematics of Computation 27, pp. 81–89. Cited by: §1.
  • M. Hillman, J.S. Chen, and Y. Bazilevs (2015) Variationally consistent domain integration for isogeometric analysis. Computer Methods in Applied Mechanics and Engineering 284, pp. 521 – 540. Cited by: §4.
  • A. Huerta, T. Belytschko, S. Fernández-Méndez, T. Rabczuk, X. Zhuang, and M. Arroyo (2018) Meshfree methods. Encyclopedia of Computational Mechanics Second Edition, pp. 1–38. Cited by: §1.
  • A. Huerta and S. Fernández-Méndez (2000) Enrichment and coupling of the finite element and meshless methods. International Journal for Numerical Methods in Engineering 48, pp. 1615–1636. Cited by: §6.
  • T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs (2005) Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering 194, pp. 4135–4195. Cited by: §1.
  • Y. Jia, Y. Zhang, G. Xu, X. Zhuang, and T. Rabczuk (2013) Reproducing kernel triangular b-spline-based fem for solving pdes. Computer Methods in Applied Mechanics and Engineering 267, pp. 342–358. Cited by: §1.
  • M. Kapl, G. Sangalli, and T. Takacs (2018) Construction of analysis-suitable planar multi-patch parameterizations. Computer-Aided Design 97, pp. 41–55. Cited by: §1.
  • J. T. Klosowski, M. Held, J. S. Mitchell, H. Sowizral, and K. Zikan (1998) Efficient collision detection using bounding volume hierarchies of k-DOPs. IEEE Transactions on Visualization & Computer Graphics 4, pp. 21–36. Cited by: Appendix C.
  • L. Kudela, N. Zander, T. Bog, S. Kollmannsberger, and E. Rank (2015) Efficient and accurate numerical quadrature for immersed boundary methods. Advanced Modeling and Simulation in Engineering Sciences 2. Cited by: §4.
  • C. Le, T. Bruns, and D. Tortorelli (2011) A gradient-based, parameter-free approach to shape optimization. Computer Methods in Applied Mechanics and Engineering 200 (9), pp. 985–996. Cited by: §1.
  • S. Li and W. K. Liu (2002) Meshfree and particle methods and their applications. Applied Mechanics Reviews 55, pp. 1–34. Cited by: §1.
  • W. K. Liu, W. Han, H. Lu, S. Li, and J. Cao (2004) Reproducing kernel element method. part i: theoretical formulation. Computer Methods in Applied Mechanics and Engineering 193, pp. 933–951. Cited by: §1.
  • W. K. Liu, S. Jun, and Y. F. Zhang (1995) Reproducing kernel particle methods. International Journal for Numerical Methods in Fluids 20, pp. 1081–1106. Cited by: §1.
  • L. B. Lucy (1977) A numerical approach to the testing of the fission hypothesis. The Astronomical Journal 82, pp. 1013–1024. Cited by: §1.
  • M. Majeed and F. Cirak (2017) Isogeometric analysis using manifold-based smooth basis functions. Computer Methods in Applied Mechanics and Engineering 316, pp. 547–567. Cited by: §1.
  • B. Mirtich (1996) Fast and accurate computation of polyhedral mass properties. Journal of Graphics Tools 1, pp. 31–50. Cited by: §3.2.
  • H. Mirzaee, J. King, J. K. Ryan, and R. M. Kirby (2013) Smoothness-increasing accuracy-conserving filters for discontinuous galerkin solutions over unstructured triangular meshes. SIAM Journal on Scientific Computing 35 (1), pp. A212–A230. Cited by: §1.
  • D. E. Muller and F. P. Preparata (1978) Finding the intersection of two convex polyhedra. Theoretical Computer Science 7, pp. 217–236. Cited by: Appendix C.
  • J.A. Nitsche (1971) Über ein variationsprinzip zur lösung von Dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind. Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg 36, pp. 9–15. Cited by: §1, §4.
  • J. T. Oden, I. Babuŝka, and C. E. Baumann (1998) A discontinuous hp finite element method for diffusion problems. Journal of Computational Physics 146 (2), pp. 491–519. Cited by: §1, §4.
  • J. Peters and U. Reif (2008) Subdivision surfaces. Springer Series in Geometry and Computing, Springer. Cited by: §1.
  • A. Rosolen and M. Arroyo (2013) Blending isogeometric analysis and local maximum entropy meshfree approximants. Computer Methods in Applied Mechanics and Engineering 264, pp. 95–107. Cited by: §6.
  • T. Rüberg and F. Cirak (2012) Subdivision-stabilised immersed b-spline finite elements for moving boundary flows. Computer Methods in Applied Mechanics and Engineering 209–212, pp. 266–283. Cited by: Appendix B, §1, §5.2.2.
  • C. Rycroft (2009) Voro++: a three-dimensional voronoi cell library in c++. Technical report Lawrence Berkeley National Lab, Berkeley, CA. Cited by: Appendix B.
  • M. Sabin (2010) Analysis and design of univariate subdivision schemes. Springer. Cited by: §2.
  • D. Schillinger, I. Harari, M. Hsu, D. Kamensky, S. K. Stoter, Y. Yu, and Y. Zhao (2016) The non-symmetric nitsche method for the parameter-free imposition of weak boundary and coupling conditions in immersed finite elements. Computer Methods in Applied Mechanics and Engineering 309, pp. 625–652. Cited by: §1, §4.
  • M. A. Scott, R. N. Simpson, J. A. Evans, S. Lipton, S. P. A. Bordas, T. J. R. Hughes, and T. W. Sederberg (2013) Isogeometric boundary element analysis using unstructured T-splines. Computer Methods in Applied Mechanics and Engineering, pp. 197–221. Cited by: §1.
  • O. Sigmund and J. Petersson (1998) Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Structural optimization 16, pp. 68–75. Cited by: §1.
  • Y. Sudhakar, J. M. De Almeida, and W. A. Wall (2014) An accurate, robust, and easy-to-implement method for integration over arbitrary polyhedra: application to embedded interface methods. Journal of Computational Physics 273, pp. 393–415. Cited by: §1.
  • V. Thomée (1977) High order local approximations to derivatives in the finite element method. Mathematics of Computation 31, pp. 652–660. Cited by: §1.
  • S.P. Timoshenko (1970) Theory of elasticity. 3rd edition, McGraw-Hill Higher Education. Cited by: §5.2.2.
  • D. Toshniwal, H. Speleers, R. R. Hiemstra, and T. J. Hughes (2017a) Multi-degree smooth polar splines: a framework for geometric modeling and isogeometric analysis. Computer Methods in Applied Mechanics and Engineering 316, pp. 1005–1061. Cited by: §1.
  • D. Toshniwal, H. Speleers, and T. J. R. Hughes (2017b) Smooth cubic spline spaces on unstructured quadrilateral meshes with particular emphasis on extraordinary points: geometric design and isogeometric analysis considerations. Computer Methods in Applied Mechanics and Engineering 327, pp. 411–458. Cited by: §1.
  • H. Xiao, E. Febrianto, Q. Zhang, and F. Cirak (2019) An immersed discontinuous galerkin method for compressible navier-stokes equations on unstructured meshes. arXiv preprint arXiv:1902.10232. Cited by: §4.
  • Q. Zhang, M. Sabin, and F. Cirak (2018) Subdivision surfaces with isogeometric analysis adapted refinement weights. Computer-Aided Design 102, pp. 104–114. Cited by: §1.