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 monographsFarin (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-splinesDe 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.
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.
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 vectorrepresents 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. .
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 .
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.
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.
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.
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.
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
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).
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.
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.
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.
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 Figure10. 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
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.
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.
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.
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.
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).
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
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
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).
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:
Deduce the set of cell edges from and .
Determine the intersection points between edges and the boundary with a bisectioning algorithm.
Introduce new vertices at the points determined in step 2.
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:
Deduce the set of cell edges from and .
Determine the intersection points between the edges and the six half planes in turn while keeping track of the vertices inside the half-spaces.
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 .
- Sobolev spaces. Elsevier Science. Cited by: Appendix A, Appendix A, §1.
- Voronoi diagrams — a survey of a fundamental geometric data structure. ACM Comput. Surv. 23 (3), pp. 345–405. Cited by: Appendix B.
- A meshfree unification: reproducing kernel peridynamics. Computational Mechanics 53, pp. 1251–1264. Cited by: §1.
- 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.
- 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.
- Higher order local accuracy by averaging in the finite element method. Mathematics of Computation 31, pp. 94–111. Cited by: §1.
- Meshfree methods: progress made after 20 years. Journal of Engineering Mechanics 143 (4), pp. 04017001: 1–37. Cited by: §1, §1.
- 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.
- Higher-order xfem for curved strong and weak discontinuities. International Journal for Numerical Methods in Engineering 82, pp. 564–590. Cited by: §4.
- Numerical integration of homogeneous functions on convex and nonconvex polygons and polyhedra. Computational Mechanics 56, pp. 967–981. Cited by: §1.
- 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.
- Computational geometry: algorithms and applications. 3rd edition, Springer. Cited by: Appendix B, Appendix D, §1.
- B(asic)-spline basics.. Technical report Madison mathematics research center, Winsconsin University. Cited by: §1, §2.
- 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.
- Centroidal voronoi tessellations: applications and algorithms. SIAM Rev. 41 (4), pp. 637–676. Cited by: Appendix B.
- Theory and practice of finite elements. Applied Mathematical Sciences, Springer New York. Cited by: Appendix A, Appendix A.
- 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.
- Partial differential equations. Graduate Studies in Mathematics, American Mathematical Society. Cited by: §1.
- Curves and surfaces for cagd. Academic Press. Cited by: §1.
- Smoothed particle hydrodynamics: theory and application to non-spherical stars. Monthly notices of the royal astronomical society 181, pp. 375–389. Cited by: §1.
- 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.
- 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.
- Variationally consistent domain integration for isogeometric analysis. Computer Methods in Applied Mechanics and Engineering 284, pp. 521 – 540. Cited by: §4.
- Meshfree methods. Encyclopedia of Computational Mechanics Second Edition, pp. 1–38. Cited by: §1.
- Enrichment and coupling of the finite element and meshless methods. International Journal for Numerical Methods in Engineering 48, pp. 1615–1636. Cited by: §6.
- 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.
- Reproducing kernel triangular b-spline-based fem for solving pdes. Computer Methods in Applied Mechanics and Engineering 267, pp. 342–358. Cited by: §1.
- Construction of analysis-suitable planar multi-patch parameterizations. Computer-Aided Design 97, pp. 41–55. Cited by: §1.
- Efficient collision detection using bounding volume hierarchies of k-DOPs. IEEE Transactions on Visualization & Computer Graphics 4, pp. 21–36. Cited by: Appendix C.
- Efficient and accurate numerical quadrature for immersed boundary methods. Advanced Modeling and Simulation in Engineering Sciences 2. Cited by: §4.
- A gradient-based, parameter-free approach to shape optimization. Computer Methods in Applied Mechanics and Engineering 200 (9), pp. 985–996. Cited by: §1.
- Meshfree and particle methods and their applications. Applied Mechanics Reviews 55, pp. 1–34. Cited by: §1.
- Reproducing kernel element method. part i: theoretical formulation. Computer Methods in Applied Mechanics and Engineering 193, pp. 933–951. Cited by: §1.
- Reproducing kernel particle methods. International Journal for Numerical Methods in Fluids 20, pp. 1081–1106. Cited by: §1.
- A numerical approach to the testing of the fission hypothesis. The Astronomical Journal 82, pp. 1013–1024. Cited by: §1.
- Isogeometric analysis using manifold-based smooth basis functions. Computer Methods in Applied Mechanics and Engineering 316, pp. 547–567. Cited by: §1.
- Fast and accurate computation of polyhedral mass properties. Journal of Graphics Tools 1, pp. 31–50. Cited by: §3.2.
- 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.
- Finding the intersection of two convex polyhedra. Theoretical Computer Science 7, pp. 217–236. Cited by: Appendix C.
- Ü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.
- A discontinuous hp finite element method for diffusion problems. Journal of Computational Physics 146 (2), pp. 491–519. Cited by: §1, §4.
- Subdivision surfaces. Springer Series in Geometry and Computing, Springer. Cited by: §1.
- Blending isogeometric analysis and local maximum entropy meshfree approximants. Computer Methods in Applied Mechanics and Engineering 264, pp. 95–107. Cited by: §6.
- 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.
- Voro++: a three-dimensional voronoi cell library in c++. Technical report Lawrence Berkeley National Lab, Berkeley, CA. Cited by: Appendix B.
- Analysis and design of univariate subdivision schemes. Springer. Cited by: §2.
- 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.
- Isogeometric boundary element analysis using unstructured T-splines. Computer Methods in Applied Mechanics and Engineering, pp. 197–221. Cited by: §1.
- 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.
- 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.
- High order local approximations to derivatives in the finite element method. Mathematics of Computation 31, pp. 652–660. Cited by: §1.
- Theory of elasticity. 3rd edition, McGraw-Hill Higher Education. Cited by: §5.2.2.
- 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.
- 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.
- An immersed discontinuous galerkin method for compressible navier-stokes equations on unstructured meshes. arXiv preprint arXiv:1902.10232. Cited by: §4.
- Subdivision surfaces with isogeometric analysis adapted refinement weights. Computer-Aided Design 102, pp. 104–114. Cited by: §1.