    # Poly-Spline Finite Element Method

We introduce an integrated meshing and finite element method pipeline enabling black-box solution of partial differential equations in the volume enclosed by a boundary representation. We construct a hybrid hexahedral-dominant mesh, which contains a small number of star-shaped polyhedra, and build a set of high-order basis on its elements, combining triquadratic B-splines, triquadratic hexahedra (27 degrees of freedom), and harmonic elements. We demonstrate that our approach converges cubically under refinement, while requiring around 50 similarly dense hexahedral mesh composed of triquadratic hexahedra. We validate our approach solving Poisson's equation on a large collection of models, which are automatically processed by our algorithm, only requiring the user to provide boundary conditions on their surface.

## Authors

##### 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

The numerical solution of partial differential equations is ubiquitous in computer graphics and engineering applications, ranging from the computation of UV maps and skinning weights, to the simulation of elastic deformations, fluids, and light scattering.

The finite element method (FEM) is the most commonly used discretization of PDEs, especially in the context of structural and thermal analysis, due to its generality and rich selection of off-the-shelf commercial implementations. Ideally, a PDE solver should be a “black box”: the user provides as input the domain boundary, boundary conditions, and the governing equations, and the code returns an evaluator that can compute the value of the solution at any point of the input domain. This is surprisingly far from being the case for all existing open-source or commercial software, despite the research efforts in this direction and the large academic and industrial interest.

To a large extent, this is due to treating meshing and FEM basis construction as two disjoint problems. The FEM basis construction may make a seemingly innocuous assumption (e.g., on the geometry of elements), that lead to exceedingly difficult requirements for meshing software. For example, commonly used bases for tetrahedra are sensitive to the tetrahedron shape, so tetrahedral mesh generators have to guarantee good element shape everywhere: a difficult task which, for some surfaces, does not have a fully satisfactory solution. Alternatively, if few assumptions are made on mesh generation (e.g., one can use elements that work on arbitrary polyhedral domains), the basis and stiffness matrix constructions can become very expensive.

This state of matters presents a fundamental problem for applications that require fully automatic, robust processing of large collections of meshes of varying sizes, an increasingly common situation as large collections of geometric data become available. Most importantly, this situation arises in the context of machine learning on geometric and physical data, where a neural network could be trained using large numbers of simulations and used to compute efficiently an approximated solution

[Chen et al., 2018; Kostrikov et al., 2018]. Similarly, shape optimization problems often require solving PDEs in the inner optimization loop on a constantly changing domain [Panetta et al., 2015].

#### Overview.

We propose an integrated pipeline, considering meshing and element design as a single challenge: we make the tradeoff between mesh quality and element complexity/cost local, instead of making an a priori decision for the whole pipeline. We generate high quality, simple, and regularly arranged elements for most of the volume of the shape, with more complex and poor quality polyhedral shapes filling the remaining gaps [Sokolov et al., 2016; Gao et al., 2017a]. Our idea is to match each element to a basis construction, with well-shaped elements getting the simplest and most efficient basis functions and with complex polyhedral element formulations used only when necessary to handle the transitions between regular regions, which are the ones that are topologically and geometrically more challenging.

A spline basis on a regular lattice has major advantages over traditional FEM elements, since it has the potential to be both accurate and efficient: it has a single degree of freedom per element, except at the boundary, yet, it has full approximation power corresponding to the degree of the spline. This observation is one of the foundations of isogeometric analysis in 3D [Hughes et al., 2005; Cottrell et al., 2009]. Unfortunately, it is easy to define and implement only for fully regular grids, which is not practical for most input geometries. The next best thing are spline bases on pure hexahedral meshes: while smooth constructions for polar configurations exist [Toshniwal et al., 2017]

, a solution applicable to general hexahedral meshes whose interior singular curves meet is still elusive, restricting this construction to simple shapes. Padded hexahedral-meshes

[Maréchal, 2009a] are necessary to ensure a good boundary approximation for both regular and polycube [Tarini et al., 2004] hexahedral meshing methods, but they unfortunately cannot be used by these constructions since their interior curve singularities meet in the padding layer.

We propose a hybrid construction that sidesteps these limitations: we use spline elements only on fully regular regions, and fill the elements that are touching singular edges, or that are not hexahedra, with local constructions (harmonic elements for polyhedra, triquadratic polynomial elements for hexahedra). This construction further relaxes requirements for meshing, since it works on general hexahedral meshes (without any restriction on their singularity structure) but also directly supports hex-dominant meshes, which can be robustly generated with modern field-aligned methods [Gao et al., 2017a; Sokolov et al., 2016]. These meshes consist mostly of well-shaped hexahedra with locally regular mesh structure, but also contain other general polyhedra. Our construction takes advantage of this high regularity, adding a negligible overhead over the spline FEM basis only for the sparse set of non-regular elements.

We demonstrate that our proposed Poly-Spline FEM retains, to a large extent, both the approximation and performance benefits of splines, at the cost of the increasing basis construction complexity, and at the same time, works for a class of meshes that can be robustly generated for most shapes with existing meshing algorithms.

Our method exhibits cubic convergence on a large data set, for a degree of freedom budget comparable to trilinear hexahedral elements, which have only quadratic convergence. To the best of our knowledge, this paper is the first FEM method exploiting the advantages of spline basis that has been validated on a large collection of complex geometries.

## 2. Related Work

When numerically solving PDEs using the finite element method, one has to discretize the spatial domain into finite elements and define shape functions on these elements. Since shape functions, element types, and mesh generation are closely related, we discuss the relevant approaches in tandem.

For complex spatial domains, the discretization is frequently based on the Delaunay triangulation [Shewchuk, 1996] or Delaunay tetrahedrization [Si, 2015], respectively, since those tessellations can be computed in a robust and automatic manner. Due to their simplicity and efficiency, linear shape functions over triangular or tetrahedral elements are often the default choice for graphics applications [Hughes, 2000], although they are known to suffer from locking for stiff PDEs, such as nearly incompressible elastic materials [Hughes, 2000].

This locking problem can be avoided by using bilinear quadrangular or trilinear hexahedral elements ( elements), which have the additional advantage of yielding a higher accuracy for a given number of elements [Cifuentes and Kalbag, 1992; Benzley et al., 1995]. Triquadratic hexahedral elements () provide even higher accuracy and faster convergence under mesh refinement (cubic converge in -norm for vs. quadratic converge for ), but their larger number of degrees of freedom (8 vs. 27) leads to high memory consumption and computational cost.

The main idea of isogeometric analysis (IGA) [Hughes et al., 2005; Cottrell et al., 2009; Engvall and Evans, 2017] is to employ the same spline basis for defining the CAD geometry as well as for performing numerical analysis. Using quadratic splines on hexahedral elements results in the same cubic convergence order as elements, but at the much lower cost of one degree of freedom per element (comparable to elements). This efficiency, however, comes at the price of a very complex implementation for non-regular hexahedral meshes. Moreover, generating IGA-compatible meshes from a given general boundary surface is still an open problem [Aigner et al., 2009; Martin and Cohen, 2010; Li et al., 2013].

Concurrent work [Wei et al., 2018]

introduces a construction that can handle irregular pure hex meshes, with tensor-product cubic splines used on regular parts. However, we focus on handling general polygonal meshes and we use quadratic splines (note that, our approach can be easily extended to cubic polynomials if desired).

A standard method for volumetric mesh generation is through hierarchical subdivision of an initial regular hexahedral mesh, leading to so-called octree meshes [Maréchal, 2009b; Ito et al., 2009; Zhang et al., 2013]. The T-junctions resulting from adaptive subdivision can be handled by using T-splines [Sederberg et al., 2004; da Veiga et al., 2011] as shape functions. While this meshing approach is very robust, it has problems representing geometric features that are not aligned with the principal axes.

Even when giving up splines or T-splines for standard / elements, the automatic generation of the required hexahedral meshes is problematic. Despite the progress made in this field over the last decade, automatically generating pure hexahedral meshes that (i) have sufficient element quality, (ii) are not too dense, and (iii) align to geometric features is still unsolved. Early methods based on paving or sweeping [Owen and Saigal, 2000; Yamakawa and Shimada, 2003; Staten et al., 2005; Shepherd and Johnson, 2008] require complicated handling of special cases and generate too many singularities. Polycube methods [Gregson et al., 2011; Li et al., 2013; Livesu et al., 2013; Huang et al., 2014; Fang et al., 2016; Fu et al., 2016], field-aligned methods [Nieser et al., 2011; Huang et al., 2011; Li et al., 2012; Jiang et al., 2014], and the surface foliation method [Lei et al., 2017] are interesting research venues, but they are currently not robust enough and often fail to produce a valid mesh.

However, if the strict requirement of producing hexahedral elements only is relaxed, field-aligned methods [Sokolov et al., 2016; Gao et al., 2017a] can robustly and automatically create hex-dominant polyhedral meshes, that is, meshes consisting of mostly, but not exclusively, of hexahedral elements. The idea is to build local volumetric parameterizations aligned with a specified directional field, and constructing the mesh from traced isolines of that parameterization, inserting general polyhedra if necessary. Their drawback is that the resulting hex-dominant meshes, are not directly supported by most FEM codes.

One option is to split these general polyhedra into standard elements, leading to a mixed FEM formulation. For instance, the field-aligned meshing of Sokolov et al.  extract meshes that are composed of hexahedra, tetrahedra, prisms, and pyramids. However, the quality of those split elements is hard to control in general. An interesting alternative is to avoid the splitting of polyhedra and instead incorporate them into the simulation, for instance though mimetic finite differences [Lipnikov et al., 2014], the virtual element method [Beirão Da Veiga et al., 2013], or polyhedral finite elements [Manzini et al., 2014]. The latter employ generalized barycentric coordinates as shape functions, such as mean value coordinates [Floater et al., 2005; Ju et al., 2005], harmonic coordinates [Joshi et al., 2007], or minimum entropy coordinates [Hormann and Sukumar, 2008]. From those options, harmonic coordinates seem most suitable since they generalize both linear tetrahedra and trilinear hexahedra to general non-convex polyhedra [Martin et al., 2008; Bishop, 2014]. While avoiding splitting or remeshing hex-dominant meshes, the major drawback of polyhedral elements is the high cost for computing and integrating their shape functions.

In the above methods the meshing stage either severely restricts admissible shape functions, or the element type puts (too) strong requirements on the meshing. In contrast, we use the most efficient elements where possible and the most flexible elements where required, which enables the use of robust and automatic hex-dominant mesh generation.

## 3. Algorithm Overview

In this section, we introduce the main definitions we use in our algorithm description, and outline the structure of the algorithm. We refer to Appendix A for a brief introduction to the finite element method and the setup of our mathematical notation.

#### Input complex and subcomplexes.

The input to our algorithm is a 3D polyhedral complex , with vertices , , consisting of polyhedral cells , , most of which are hexahedra. Figure 2 shows a two-dimensional example of such complex. The edges, faces, and cells of the mesh are defined combinatorially, that is, edges are defined by pairs of vertices, faces by sequences of edges, and cells by closed surface meshes formed by faces. We assume that 3D positions of vertices are also provided as input and that is three-manifold, i.e., that there is a way to identify vertices, edges, faces, and cells with points, curves, surface patches and simple volumes, such that their union is a three-manifold subset of .

We assume that for any hexahedron there is at most one non-hexahedral cell sharing one of its faces or edges, which can be achieved by refinement. We also assume that no two polyhedral cells are adjacent, and that no polyhedron touches the boundary, which can also be achieved by merging polyhedral cells and/or refinement. This preprocessing step (i.e., one step of uniform refinement) is discussed in Section 6. As a consequence of our refinement, all faces of are quadrilateral.

One of the difficulties of using general polyhedral meshes for basis constructions is that, unlike the case of, for example, pure tetrahedral meshes, there is no natural way to realize all elements of the mesh in 3D just from vertex positions (e.g., for a tetrahedral mesh, linear interpolation for faces and cells is natural). This requires constructing bases on an explicitly defined parametric domain associated with the input complex.

For this purpose, we define a certain number of complexes related to the original complex (Figure 2). There are two goals for introducing these: defining the parametric domain for the basis, and defining the geometric map, which specifies how the complex is realized in three-dimensional physical space. Figure 2. Complexes involved in our construction. In green we show S, in red Q, and in blue P.
• [leftmargin=*]

• is the hexahedral part of , consisting of hexahedra .

• is the non-hexahedral part of , consisting of polyhedra .

• is the complex consisting of spline-compatible hexahedra defined in Section 4.1.

• is the spline-incompatible pure-hexahedral part of .

Note that the sub-complexes of are nested: .

In the context of finite elements, the distinction between parametric space and physical space is critical: the bases on the hexahedral part of the mesh are defined in terms of parametric space coordinates, where all hexahedra are unit cubes; this makes it possible to define simple, accurate, and efficient bases. However, the derivatives in the PDE are taken with respect to physical space variables, and the unknown functions are naturally defined on the physical space. Remapping these functions to the parametric space is necessary to discretize the PDE using our basis. We define parametric domains , , , and corresponding to , , , and , respectively. consists of unit cubes , one per hexahedron with corresponding faces identified, and and are its subcomplexes. The complete parametric space is obtained by adding a set of polyhedra for , defined using the geometric map as described below. For polyhedra, physical and parametric space coincide.

#### Geometric map and complex embedding.

The input complex, as it is typical for mesh representations, does not define a complete geometric realization of the complex: rather it only includes vertex positions and element connectivity. We define a complete geometric realization as the geometric map , from the parametric domain to the physical space. We use for points in the parametric domain, and for points in the physical space, and denote the image of the geometric map by (Figure 3).

The definition requires bootstrapping: is first defined on . For example, the simplest geometric map on can be obtained by trilinear interpolation: restricted to a unit cube is a trilinear interpolation of the positions of the vertices of its associated hexahedron . We make the following assumption about : the map is bijective on the faces of , corresponding to the boundary of any polyhedral cell , and the union of the images of these faces does not self-intersect and encloses a volume . Section 6 explains how this is ensured. Then we complete by adding the volume as the parametric domain for . We add this volume to the parametric domain , identifying corresponding faces with faces in , and defining the geometric map to be the identity on these domains.

The simplest trilinear map is adequate for elements of the mesh outside the regular part , but is insufficient for accuracy on the regular part, as discussed below. We consider a more complex definition of ensuring smoothness across interior edges and faces of , described in Section 5, after we describe our basis construction. Our construction is isoparametric, that is, it uses the same basis for the geometric map as for the solution.

In other words, on we use the standard tri-quadratic geometric map that maps each reference cube to the actual hex-element in the mesh. On we use a spline mapping, explained in Section 5. On the polyhedral part, the geometric map is the identity, thus all quantities are defined directly on the physical domain.

#### Overview of the basis and discretization construction.

Given an input complex , we construct a set of bases , , such that:

• [leftmargin=*]

• the restriction of basis function to spline compatible hexahedral domains is a spline basis function;

• the restriction to hexahedra is a standard triquadratic () element function;

• the restriction to polyhedra (or ) is a harmonic-based nonconformal, third-order accurate basis function.

The degrees of freedom (dofs) corresponding to basis functions are associated with:

• [leftmargin=*]

• each hexahedron either in or adjacent to a spline-compatible one (spline cell dofs);

• each boundary vertex, edge, or face of (spline boundary dofs); these are needed to have correct approximation on the boundary;

• each vertex, edge, face, and cell of (triquadratic element degrees of freedom).

The total number of degrees of freedom is denoted by . While most of the construction is independent of the choice of PDE (we assume it to be second-order), with a notable exception of the consistency condition for polyhedral elements, we use the Poisson equation to be more specific.

Note that hexahedra adjacent to , but not in (i.e., hexahedra in ) get both spline dofs and triquadratic element dofs: such a cell may have dofs instead of 27.

Polyhedral cells are not assigned separate degrees of freedom: the basis functions with support overlapping polyhedra are those associated with dofs at incident hexahedra.

We assemble the standard stiffness matrix for an elliptic PDE, element-by-element, performing integration on the hexahedra of and polyhedra . The entry of the stiffness matrix for the Poisson equation is computed as follows:

 (1) Kij=∑ˆC∈ˆM∫g(ˆC)∇ϕi(x)⋅∇ϕj(x)\ddx,

where . The actual integration is performed on the elements in the parametric domain , using a change of variables for every element:

 (2) Kij=∑ˆC∈ˆM∫ˆC∇ˆϕi(ˆx)TA(ˆx)∇ˆϕj(ˆx)\absDg\ddˆx

where is the metric tensor of the geometric map at , given by , with being the Jacobian of .

In the next sections, we describe the construction of the basis on each element type, the geometric map, and the stiffness matrix construction.

## 4. Basis construction

We seek to construct a basis on that has the following properties:

1. it is everywhere on , at regular edges and vertices, and within each and (polynomials on hexahedra).

2. it has approximation order 3 on each and . Figure 4. Spline local grid (shown in dark), for an internal and a boundary quadrilateral. The color codes are as defined in Figure 2. Figure 5. Spline hex degrees of freedom for a central element and a corner one.

The unknown function on the domain is approximated by , where are the basis functions. The support of each basis function is a union of a set of the images under of cells in .

The actual representation of the basis, which allows us to perform per-element construction of the stiffness matrix, consists of three parts. The first two parts are local: we define a local set of dofs and a local basis. For hexahedral elements, there are several types of local polynomial bases, each coming with its set of local dofs, associated with a local control mesh for the element. These basis functions are encoded as sets of polynomial coefficients. For polyhedral elements, all local basis functions are weighted combinations of harmonic kernel functions and a triquadratic polynomial, so these are encoded as kernel centers, weights and polynomial coefficients.

The third part is the local-to-global linear map that represents local dofs in terms of the global ones. Importantly, unlike most standard FEM formulations, our local-to-global maps are not necessarily simply identifying local dofs the global ones: some local dofs are linear combinations of global ones. These maps are formally represented by matrices, where is a small number of local dofs, and is the total number of global dofs. However, as the elements local dofs depend only on nearby global dofs, these matrices have a small number of nonzeros and can be encoded in a compact form.

In the following, we consider the construction of these three elements (set of local basis functions, set of local dofs, local-to-global map) for each of our three element types. But before we can construct the basis for each element, hexahedral elements need to be classified into

(spline-compatible) and .

### 4.1. Spline-compatible hexahedral elements

We define a hexahedron to be spline-compatible, if its one-ring cell neighborhood is a regular grid, possibly cut on one or more sides if is on the boundary, see Figure 4.

The local dofs of this element type form a grid (for interior elements), with the element in the center (Figure 5 left); for boundary elements, there are still 27 dofs, ensuring a full triquadratic polynomial reproduction. If a single layer with 9 dofs is missing, we add an extra degree of freedom for each face of the local grid corresponding to the boundary. Other cases are handled in a similar manner; e.g. the configuration for a regular corner is shown in Figure 5, right.

The basis functions

in this case are just the standard triquadratic uniform spline basis functions for interior hexahedra. For the boundary case, we use the knot vector

in the direction perpendicular to the boundary. Figure 6 shows an example of the bases in 2D, for an internal node on the left and for a boundary node on the right. Finally, the local-to-global map simply identifies local basis dofs with corresponding global ones. Figure 6. Plot of the spline bases for a regular 2D grid.

Compared to a standard element, the ratio of degrees of freedom to the number of elements is much lower (a single degree of freedom per element for splines), although the approximation order is the same.

### 4.2. Q2 hexahedral elements

This element is used for all remaining hexahedra. It is a standard element, widely used in finite element codes. Local dofs for this element are associated with the element vertices, edge midpoints, face centers, and cell centers (Figure 26).

The local basis functions for the element are obtained as the tensor product of the interpolating quadratic bases on the interval , consisting of , , and (Appendix C). The only complicated part in the case of elements is the definition of the local-to-global map. For the two-dimensional setting, it is illustrated in Figure 7. Figure 7. Local-to-global map for a Q2 element (gray) adjacent to a single spline element (green).

The difficulty in the construction of this map is due to the interface between spline elements and elements, and the need to ensure continuity between the two. In the two-dimensional case, suppose that a element shares an edge with exactly one quad spline element . Let , , be the global dofs of the spline element, and let , , be the degrees of freedom of the element, as shown in the picture.

In this case, we ensure continuity of the basis by expressing the values of the polynomials on at the shared boundary points in terms of global degrees of freedom. Since both the and the spline basis restricted to an edge are quadratic polynomials, they only need to be equal on three distinct points of the edge to ensure continuity. By noticing that the basis is interpolatory at the nodes, it is enough to evaluate the spline basis at these edge nodes.

For the two-dimensional example in Figure 7, the local-to-global map for the local dofs , and along the edge (in blue) that the element shares with the spline is obtained as follows:

 (3) q31=14(u11+u12+u21+u22),q32=38(u12+u22)+116(u11+u21+u13+u23),q33=14(u12+u13+u22+u23).

In 3D, the construction is similar. We first identify all spline bases overlapping with a local dof on the boundary of a element (i.e., a vertex, edge, or face dof). To determine the weights of the local-to-global map, we evaluate each spline basis on the local dof and set it as weight.

The remaining degrees of freedom of the element are identified with global degrees of freedom at the same locations. We note once again that at the center of cells in with neighboring cells in , there are two dofs, one spline dof and one dof. Figure 8 shows an example of two basis functions on the transition from the regular part on the left to the “irregular” part on the right. We clearly see that on the regular part the bases are splines and on the irregular one are the standard basis function: on the interface the functions are only . Figure 8. Plot of the bases on a junction between a regular (green) and an irregular (red) part for a regular 2D grid.

### 4.3. Basis construction on polyhedral cells

The construction of the basis on the polyhedral cells is quite different from the construction of the basis on hexahedra. For hexahedra, the basis functions are defined on the parametric domain , and are remapped to via the geometric map. For polyhedra, we construct the basis directly in physical space.

On possible option to construct the basis on polyhedral cells is to split each polyhedral cells into tetrahedra. This approach has two main disadvantages: (i) it requires the use of pyramids to ensure conformity to the neighboring hexahedra, (ii) it is difficult to guarantee a sufficient element quality after subdivision. Instead, we follow the general approach of [Martin et al., 2008] with two important alterations designed to ensure third-order convergence.

Recall that all polyhedron faces are quadrilateral, and all polyhedra are surrounded by hexahedra, specifically hexahedra as their neighborhood is not regular. Moreover, since we always perform an initial refinement step, there are no polyhedral cells touching each other. We use the degrees of freedom on the faces of these elements as degrees of freedom for the polyhedra, therefore the local-to-global map in this case is trivial.

Each dof is already associated with a basis function defined on the hexahedra adjacent to the polyhedron. We construct the extension of to the polyhedron from harmonic kernels centered at positions outside the polyhedron and quadratic monomials , , as

 (4) ϕj∣∣∣∣P(x) =k∑i=1wjiψi(x)+10∑d=1ajdqd(x) =wj⋅ψ(x)+aj⋅q(x),

where , , , and . The coefficients and are and matrices, respectively, with (scalar PDEs) or (vector PDEs). Following , the weights are determined using a least squares fit to match the values of the basis evaluated on a set of points sampled on the boundary of the polyhedron .

In [Martin, 2011], it is shown that this construction automatically guarantees reproduction of linear polynomials if are linear; the quadratic case is fully analogous. However, this condition is insufficient for high-order convergence, because our basis is non-conforming, that is non . In the context of the second-order PDEs we are considering, it means that it lacks continuity on the boundary of the polyhedron. For this type of elements, additional consistency conditions are required to ensure high-order convergence. These conditions depend on the PDE that we need to solve. Figure 9. The local basis for a polygon consists of the set of triquadratic polynomials qd and harmonic kernels ψi centered at shown locations zi.

#### FEM theory detour.

To achieve higher order convergence three conditions need to be satisfied: (1) polynomial reproduction; (2) consistency, which we discuss in more detail below and (3) quadrature accuracy. We refer to standard FEM texts such as [Braess, 2007] for details, as well as to virtual element method literature (e.g., [de Dios et al., 2016] is closely related).

To satisfy the third condition, we use high-order quadrature on the polyhedron: we decompose it into tetrahedra and use Gaussian quadrature points in each tetrahedron (the decomposition is detailed in Section 6.1). The first condition, polynomial reproduction, is ensured by construction of the basis above.

The second constraint, consistency, requires further elaboration. We first derive it for the Poisson equation, and then summarize the general form. We leave as future work the complete proof of the convergence properties of our method (cf. [de Dios et al., 2016]), which requires, in particular, a proper stability analysis. Nevertheless, in Section 7 we provide numerical evidence that our method does converge at the expected rate, and that its conditioning is not affected in a significant way by the presence of nonconforming polyhedral elements.

The standard way to find the solution of a PDE for a finite element system is to consider its weak form. For the Poisson equation, find such that

 (5) ∫ΩΔuv=−∫Ω∇u⋅∇v=∫Ωfv,∀v

Remark. We omit, for readability, the integration variable . In the remaining formulas we use integration over the physical space exclusively, in practice carried over to the parametric space by adding the Jacobian of the geometric map.

Then, is approximated by , and is taken to be in the space spanned by the basis functions . The stiffness matrix entries are obtained as , where the integral is computed per element , leading to the discrete system (Appendix A).

For general non-conforming elements, however, we cannot rely on this standard approach. For example, if we consider piecewise-constant elements for the Poisson equation, the stiffness matrix would be all zeros.

However, for a given PDE, one can construct converging non-conforming elements. One condition that is typically used, is that the discrete matrix, constructed per element as above, gives us exact values of the weak-form integral for all polynomials reproduced by the basis (cf. -consistency property in [de Dios et al., 2016]).

As our basis reproduces triquadratic monomials (i.e., they are in the span of bases ), we have . To ensure consistency, we require that any nonconforming basis function satisfies

 (6) −∫g(ˆM)Δqdϕj=∑iKijqid

To convert this equation to an equation for the unknown coefficients and , we observe that

 (7) ∑iKijqid=∫g(ˆM)(∑iqid∇ϕi)⋅∇ϕj=∫g(ˆM)∇qd⋅∇ϕj

due to the polynomial reproduction property. Separating the integral into the part over the hexahedra and over the polyhedron , we write

 (8) ∑iKijqid

where

 CH=∑ˆC∈ˆM∖P∫g(ˆC)∇qd⋅∇ϕj,b=∫P∇qd⋅∇ψ,c=∫P∇qd⋅∇q.

Similarly, the left-hand side of Equation 6 is reduced to a linear combination of and . This forms a set of additional constraints for the coefficients of the basis functions on the polyhedron. To enforce them on each polyhedron, we solve a constrained least squares system for each nonconforming basis function and store the obtained coefficients.

Importantly, the addition of constraints to the least squares system does not violate the polynomial reproduction property on the polyhedron. This can be seen as follows. Let be the linear combination of basis functions overlapping that yields a triquadratic mononomial when restricted to . Then is continuous on : the samples at the points of the boundary are from a quadratic function, therefore, match exactly the quadratic continuation to adjacent hexahedra.

The consistency condition (Equation 6) applied to simply states that it satisfies the integration by parts formula, which it does as it is at the element boundaries, and smooth on the elements:

 −∑C∫CΔqdvh=∑C∫C∇qd⋅∇vh.

We conclude that is in the space defined by the consistency constraint, and imposing this constraint preserves polynomial reproduction. See Appendix D for the complete list of constraints for the Poisson equation.

More generally, for a linear PDE and for any polynomial (for vector PDEs, e.g., elasticity, this means that all components are polynomial) we require

 a(qd,vh)=ah(qd,vh),wherea(u,v)=∫ΩF(x,u,∇u,Δu)v

where is a linear function of its arguments depending on , and is defined as a sum of integral over after formal integration by parts of , to eliminate the second-order derivatives. For a conforming basis, this condition automatically follows from the integration by parts formulas, which are applicable. We now split the two bilinear forms as and where and contains the integral over the hexahedral known part, and and the integral over the polyhedral unknown part. Thus, for a basis we obtain the following set of constraints

 aH(qd,ϕj)−aHh(qd,ϕj) =aH(qd,wj⋅ψ(x)+aj⋅q(x)) −aHh(qd,wj⋅ψ(x)+aj⋅q(x)).

For a scalar-valued PDE, we have the same number of constraints (5 in 2D and 9 in 3D) as monomials , thus we are guaranteed to have a solution that respects the constraints for any . For vector PDEs (e.g., elasticity), we impose the additional constraints such that the coefficients are the same for all dimensions , or , which simplifies the implementation, but increases the number of required centers , so that all constraints can be satisfied. More explicitly, for vector PDEs we require that the constraints

for are satisfied, with and denoting scalar polynomials and scalar basis functions respectively, defined as in (4) for dimension 1, and is the unit vector for axis . For dimensions 2 and 3, the number of monomials is and respectively. The number of constraints is given by , and thus we will need at least 15 in 2D and 72 in 3D to ensure that the constraints are respected.

### 4.4. Imposing boundary conditions

We consider two standard types of boundary conditions: Dirichlet (fixed function values on the boundary) and Neumann (fixed normal derivatives at the boundary). Neumann (also known as natural) boundary conditions are handled in the context of the variational formulation of the problem as extra integral terms, in the case of inhomogeneous conditions. Homogeneous conditions do not require any special treatment and are imposed automatically in the weak formulation.

We assume that the Dirichlet conditions are given as a continuous function defined on the boundary of the domain. For all boundary dofs, we sample the boundary condition on the faces of the domain and perform a least-squares fit to retrieve the nodal values.

## 5. Geometric map construction

The geometric map is a map from to , defined per element. Its primary purpose is to allow us to construct basis functions on reference domains (i.e., the elements of that are unit cubes), and then to remap them to the physical space as . As the local basis on the polyhedral elements is constructed directly in the physical space, is the identity on these elements.

The requirements for the geometric map are distinct for the spline and elements, and are matched by using spline basis itself for and trilinear interpolation for elements.

Because of the geometric mapping , for the quadratic spline, the basis does not reproduce polynomials in the physical space; nevertheless, the approximation properties of the basis are retained [Bazilevs et al., 2006].

For elements, Arnold et al.  shows that bilinear maps are sufficient, and in fact allow to retain reproduction of triquadratic polynomials in the physical space. This is very important for the basis construction on polyhedral elements, as polynomial reproduction on these elements depends on reproduction of polynomials on the polyhedron boundary.

#### Computing the geometry map.

If we assume that the input only has vertex positions for , we solve the equations , which is a linear system of equations in terms of coefficients of in the basis we choose. In the trilinear basis, the system is trivial, as the coefficients coincide with the values at , and these are simply set to . For the triquadratic basis, this is not the case, and a linear system needs to be solved. If the system is under-determined, we find the least-norm solution.

## 6. Mesh preprocessing and refinement

Without loss of generality, we restrict the meshing discussion to 2D, as the algorithm introduced in this section extends naturally to 3D.

For the sake of simplicity, in this discussion the term polygon refers to non-quadrilateral elements. As previously mentioned, our method can be applied to hybrid meshes without two adjacent polygons and without polygons touching the boundary, which we ensure with one step of refinement. While our construction could be extended to support these configurations, we favored refinement due to its simplicity. Refining polygonal meshes is an interesting problem on its own: while there is a canonical way to refine quads, there are multiple ways to refine a polygon. We propose the use of polar refinement (Section 6.2), which has the added benefit of allowing us to resample large polygons to obtain a uniform element size. However, to avoid self-intersections between edges during the refinement, we impose each polygon be star-shaped. This condition is often, but not always, satisfied by existing hybrid meshers: we thus introduce a simple merging and splitting procedure to convert hybrid meshes into star-shaped polyhedral meshes (Section 6.1), and then detail our refinement strategy (Section 6.2).

Another advantage of restricting ourselves to star-shaped polygons is that partitioning it into triangles (respectively tetrahedra in 3D) is trivial, by introducing a point in the kernel and connecting it to all the boundary faces. This step is required to generate quadrature points for the numerical integration (Section 4.3): the quality of the partitioning is usually low, but this is irrelevant for this purpose.

### 6.1. Mesh preprocessing Figure 10. Our algorithm iteratively merges polygons (gray polygon in the first image), until the barycenter of the merged polygon is inside its kernel (gray polygon in the second).

We propose a simple and effective algorithm to convert polygonal meshes into star-shaped polygonal meshes, by combining existing polygons until they are star-shaped (and eventually splitting them if they contain a concave part of the boundary).

For every non-star-shaped polygon, we compute its barycenter and connect it to all its vertices (Figure 10, left). This procedure generates a set of intersecting segments (red in Figure 10), which we use to grow the polygon by merging it with the faces incident to each intersecting segment. The procedure is repeated until no more intersections are found, which usually happens in one or two iterations in our experiments. If we reach a concave boundary during the growing procedure, it might be impossible to obtain a star-shaped polyhedron by merging alone: In these cases, we triangulate the polygon, and merge the resulting triangles in star-shaped polygons if possible.

### 6.2. Polar refinement

Each star-shaped polygon is refined by finding a point in its kernel (Figure 11, a), connecting it to all its vertices (b), splitting each edge with mid-point subdivision and connecting them to the point in the kernel (c), and finally adding rings of quadrilaterals around the boundary (d). Figure 12 shows an example of polar refinement in two and three dimensions. The more splits are performed in the edge, the more elements are added. This is a useful feature to homogenize the element size in case the polygons were expanded too much during the mesh preprocessing stage. In our implementation, we split the edges evenly, ensuring that the shortest segment has a length as close as possible to the average edge length of the input mesh. Figure 11. Polar refinement for polygons. Figure 12. Example of polar refinement for a polygon and a polyhedron. The bottom view is a cut-through of the actual 3D mesh.

## 7. Evaluation

We demonstrate the robustness of our method by solving the Poisson equation on a dataset of pure hex and hybrid meshes, consisting of 205 star-shaped polygonal meshes in 2D, 165 pure hexehedral meshes in 3D, and 29 star-shaped polyhedral meshes in 3D. The dataset can be found at https://cims.nyu.edu/gcl/papers/2019-Polyspline-Dataset.zip. All those meshes were automatically generated using [Gao et al., 2017a, b]. We show a selection of meshes from our dataset in Figures 1 and 13.

We evaluated the performance, memory consumption, and running time of our proposed spline construction compared with standard and elements. For our experiments, we compute the approximation error on a standard Franke’s test function [Franke, 1979] in 2D and 3D (Appendix B). Note that in all these experiments, we enforced the consistency constraints on the bases spanning the polyhedral elements, to ensure the proper convergence order.

The 2D experiments were run on a PC with an Intel® Core™ i7-5930K CPU @ 3.50GHz with 64 GB, while the 3D dataset was run on a HPC cluster with a memory limit of 64 GB.

#### Absolute Errors.

Figure 14 shows a scatter plot of the and errors on both 2D and 3D datasets, with respect to the number of bases created by each type of elements (, , Splines), after one step of polar refinement. The plot shows that in 2D both the and errors are about orders of magnitude lower for our splines compared to , while keeping a similar number of dofs. In comparison, has lower error, but requires a much larger number of dofs. In 3D the spread of both errors is much larger, and the gain in is less visible, but still present, compared to . Figure 15. Peak memory for the direct solver as reported by Pardiso. Left: 2D results. right: 3D results.

#### Memory.

A histogram of the memory consumption of the solver is presented in Figure 15. The figure shows the peak memory usage as reported by Pardiso [Petra et al., 2014b, a] when solving the linear system arising from the FEM. Out of the 159 pure hexahedral models we tested, 33 went out of memory when solving using elements, while only 2 are too big to solve with our spline bases. On the star-shaped hybrid meshes, one model is too big to solve for both and our spline construction. More detailed statistics are reported in Table 1. We remark that the error for our method is higher than because our method has less dofs (50% less in average) since both meshes have the same number of vertices.

#### Time.

Figure 18 shows the assembly time and solve time for solving a Poisson problem on an unit square (cube) under refinement in two (three) dimensions. Note that both steps (assembly and solve) are performed in parallel. For the 2D experiment we used a 3.1 GHz Intel Core i7-7700HQ with 8 threads, while in 3D we used a 3.5 GHz Intel Core i7-5930K with 12 threads (both machines use hyper-threading). In Table 1 we summarize the timings for the large dataset using a 2.6 GHz Intel Xeon E5-2690v4 with 8 threads. In all cases, the total time is dominated by the solving time.

#### Convergence.

Figures 16 and 19 show the convergence of spline elements vs and for the , , and norms, in the ideal case of a uniform grid, both in 2D and 3D. This is in a sense the best-case scenario that can be expected for our spline construction: every element is regular and has a or neighborhood. In this situation, splines exhibit a superior convergence under both , , and norms.

On a 2D test mesh with mixing polygons and splines (model shown in Figure 12 top), we achieved a convergence rate of 2.8 in , and 3.1 in (Figure 17, left). Figure 17 also displays the convergence we obtained on a very simple hybrid 3D mesh, starting from a cube marked as a polyhedron, to which we applied our polar refinement described in Section 6. On this particular mesh, the splines exhibited a convergence similar to , albeit producing an error that is somewhat larger.

#### Consistency Constraints.

Figure 20 shows the effect of our consistency constraint on the convergence of a polygonal mesh under refinement (the one shown in Figure 12, top), with elements used on the quadrilateral part. Without imposing any constraint on the bases overlapping the polygon, one can hope at best for a convergence of , whereas pure elements should have a convergence rate of . With a constraint ensuring linear reproduction for the bases defined on polyhedra, the convergence rate is still only . Finally, with the constraints we describe in Section 4.3 to ensure the bases reproduce triquadratic polynomials, we reach the expected convergence rate of .

#### Polyhedral Basis Resilience.

Our polyhedral bases are less susceptible to badly shaped elements than . We computed the and interpolation errors for the gradients of the Franke function for 14 badly shaped hexahedra, Figure 21 shows some of them. The and maximum and average errors are 3 times smaller with our polygonal basis.

#### Conditioning and Stability.

An important aspect of our new FE method is the conditioning of the resulting stiffness matrix: this quantity relates to both the stability of the method, and to its performances when an iterative linear solver is used (important only for large problems where direct solvers cannot be used due to their memory requirements). We compute the condition number of the Poisson stiffness matrix on a regular and perturbed grid (Figure 22). In both cases, our discretization has a good conditioning number, slightly higher than pure linear elements, but lower than pure quadratic elements (while sharing the same cubic convergence property). To evaluate the conditioning of the polyhedral bases we started from a base mesh of good quality, marked of the quads as polygons, and pushed one of the vertices inwards. Even for this extreme distortion of polyhedral elements, the conditioning remained similar to the case when no polyhedral elements are used on the same mesh. Figure 21. Low-quality polyhedra used to evaluate the interpolation errors. Figure 22. Evolution of the condition number of the stiffness matrix for the Poisson problem under refinement. For each level of refinement we artificially marked 5% of the quads as polyhedra and move one random vertex on the diagonal between 20% 40%, as shown in the insets figures in blue. Note that some of the curves coincide, that is Q1 with Q1 poly and Q2 with Q2 poly.

#### Elasticity.

While most of our testing was done for the Poisson equation, we have performed some testing of linear elasticity problems. Figure 23 top shows the solution of a linear elasticity problem on a pure hexahedral mesh. The outer loops of the knots are pulled outside of the figure, deforming the knot. The color in the figure represents the magnitude of the displacement vectors. On the bottom we show the result for a Young’s modulus of 2e5.

Figure 24 shows a plot for the linear elasticity PDE with Young’s modulus 200 and Poisson’s ratio 0.35 on a regular grid, and similar results are obtained an hybrid mesh, Figure 25. The convergence plots for and are obtained by mixing regular / bases with the polyhedral construction (Section 4.3). Figure 23. Displacements computed solving linear elasticity on a pure hexahedral 3D model, using spline bases. Top a complicated model with λ=1 and μ=1, bottom a bended bar ν=0.35 and large young modulus E=2e5.

## 8. Limitations and concluding remarks

We introduced Poly-Spline FEM, an integrated meshing and finite element method designed to take advantage of recent developments in hexahedral-dominant meshing, opening the doors to black box analysis with an high-order basis and cubic convergence under refinement. Our approach is to use the best possible basis for each element of the mesh and is amenable to continuous improvement, as the mesh generation methods and basis constructions improve. For instance, in this setting, one can avoid costly manual mesh repair and improvement, at the expense of modest increases in solution time, by switching to more expensive, but much less shape-sensitive elements when a hexahedron is badly shaped.

While our basis construction is resilient to bad element quality, the geometric map between the splines and the elements might introduce distortion (and even inversions in pathological cases), lowering convergence rate. These effects could be ameliorated by optimizing the positions of the control points of the geometric map, which is an interesting avenue for future work.

Our current construction always requires an initial refinement step to avoid having polyhedra adjacent to other polyhedra or to the boundary. This limitation could be lifted by generalizing our basis construction, and would allow our method to process very large datasets, that cannot be refined due to memory considerations. Another limitation of our method is that the consistency constraints in our basis construction (Section 4.3) are PDE-dependent, and they thus require additional efforts to be used with a user-provided PDE: a small and reusable investment compared to the cost of manually meshing with hexahedra every surface that one wishes to analyze using elements. The code can be found at https://polyfem.github.io/ and provides an automatic way to generate such constraints relying on both the local assembler and automatic differentiation.

Poly-Spline FEM is a practical construction in-between unstructured and fully-structured pure splines: it requires a smaller number of dofs than (thanks to the spline elements) while preserving cubic convergence rate. We believe that our construction will stimulate additional research in the development of heterogeneous FEM methods that exploit the regularity of spline basis and combine it with the flexibility offered by traditional FEM elements. To allow other researchers and practitioners to immediately build upon our construction, we will release our entire software framework as an open-source project.

## Acknowledgements

We are grateful to the NYU HPC staff for providing computing cluster service. This work was partially supported by the Sponsor NSF CAREER Rl award Grant #3, the NSF grant Grant #3, the NSF grant Grant #3, the NSF grant Grant #3, the SNSF grant Grant #3, a gift from Adobe Research, and a gift from nTopology.

## References

• 
• Aigner et al.  M. Aigner, C. Heinrich, B. Jüttler, E. Pilgerstorfer, B. Simeon, and Vuong. 2009. Swept volume parameterization for isogeometric analysis.
• Arnold et al.  D. Arnold, D. Boffi, and R. Falk. 2002. Approximation by quadrilateral finite elements. Math. Comput. (2002).
• Bazilevs et al.  Y. Bazilevs, L. Beirao da Veiga, J. A. Cottrell, T. J. Hughes, and G. Sangalli. 2006.

Isogeometric analysis: approximation, stability and error estimates for h-refined meshes.

Math. Meth. Appl. Sci. (2006).
• Beirão Da Veiga et al.  L. Beirão Da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo. 2013. Basic Principles Of Virtual Element Methods. Math. Meth. Appl. Sci. (2013).
• Benzley et al.  S. E. Benzley, E. Perry, K. Merkley, B. Clark, and G. Sjaardema. 1995. A comparison of all hexagonal and all tetrahedral finite element meshes for elastic and elasto-plastic analysis. In Proceedings of the 4th International Meshing Roundtable.
• Bishop  J. Bishop. 2014. A displacement-based finite element formulation for general polyhedra using harmonic shape functions. Int. J. Numer. Methods Eng. (2014).
• Braess  D. Braess. 2007. Finite elements: Theory, fast solvers, and applications in solid mechanics.
• Chen et al.  R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud. 2018. ArXiv e-prints (2018).
• Cifuentes and Kalbag  A. O. Cifuentes and A. Kalbag. 1992. A performance study of tetrahedral and hexahedral elements in 3-D finite element structural analysis. Finite Elements in Analysis and Design (1992).
• Cottrell et al.  J. A. Cottrell, T. J. R. Hughes, and Y. Bazilevs. 2009. Isogeometric Analysis: Toward Integration of CAD and FEA.
• da Veiga et al.  L. B. da Veiga, A. Buffa, D. Cho, and G. Sangalli. 2011. Isogeometric analysis using T-splines on two-patch geometries. Comput. Meth. Appl. Mech. Eng. (2011).
• de Dios et al.  B. A. de Dios, K. Lipnikov, and G. Manzini. 2016. The nonconforming virtual element method. ESAIM: Mathematical Modelling and Numerical Analysis (2016).
• Engvall and Evans  L. Engvall and J. A. Evans. 2017. Isogeometric unstructured tetrahedral and mixed-element Bernstein-Bezier discretizations. Comput. Meth. Appl. Mech. Eng. (2017).
• Fang et al.  X. Fang, W. Xu, H. Bao, and J. Huang. 2016. All-hex meshing using closed-form induced polycube. ACM Trans. Graph. (2016).
• Floater et al.  M. S. Floater, G. Kós, and M. Reimers. 2005. Mean Value Coordinates in 3D. Comput. Aided Geom. Des. (2005).
• Franke  R. Franke. 1979. A Critical Comparison of Some Methods for Interpolation of Scattered Data. (1979).
• Fu et al.  X. Fu, C. Bai, and Y. Liu. 2016. Efficient Volumetric PolyCube-Map Construction. Comput. Graph. Forum (2016).
• Gao et al. [2017a] X. Gao, W. Jakob, M. Tarini, and D. Panozzo. 2017a. Robust Hex-dominant Mesh Generation Using Field-guided Polyhedral Agglomeration. ACM Trans. Graph. (2017).
• Gao et al. [2017b] X. Gao, D. Panozzo, W. Wang, Z. Deng, and G. Chen. 2017b. Robust Structure Simplification for Hex Re-meshing. ACM Trans. Graph. (2017).
• Gregson et al.  J. Gregson, A. Sheffer, and E. Zhang. 2011. All-Hex Mesh Generation via Volumetric PolyCube Deformation. Comput. Graph. Forum (2011).
• Hormann and Sukumar  K. Hormann and N. Sukumar. 2008. Maximum Entropy Coordinates for Arbitrary Polytopes. Comput. Graph. Forum (2008).
• Huang et al.  J. Huang, T. Jiang, Z. Shi, Y. Tong, H. Bao, and M. Desbrun. 2014. L1-based Construction of Polycube Maps from Complex Shapes. ACM Trans. Graph. (2014).
• Huang et al.  J. Huang, Y. Tong, H. Wei, and H. Bao. 2011. Boundary aligned smooth 3D cross-frame field. ACM Trans. Graph. (2011).
• Hughes et al.  T. J. Hughes, J. A. Cottrell, and Y. Bazilevs. 2005. Isogeometric analysis: cad, finite elements, NURBS, exact geometry and mesh refinement. Comput. Meth. Appl. Mech. Eng. (2005).
• Hughes  T. J. R. Hughes. 2000. The Finite Element Method. Linear Static and Dynamic Finite Element Analysis.
• Ito et al.  Y. Ito, A. M. Shih, and B. K. Soni. 2009. Octree-based reasonable-quality hexahedral mesh generation using a new set of refinement templates. Int. J. Numer. Methods Eng. (2009).
• Jiang et al.  T. Jiang, J. Huang, Y. T. Yuanzhen Wang, and H. Bao. 2014. Frame Field Singularity Correction for Automatic Hexahedralization. IEEE Trans. Vis. Comput. Graph. (2014).
• Joshi et al.  P. Joshi, M. Meyer, T. DeRose, B. Green, and T. Sanocki. 2007. Harmonic Coordinates for Character Articulation. ACM Trans. Graph. (2007).
• Ju et al.  T. Ju, S. Schaefer, and J. Warren. 2005. Mean Value Coordinates for Closed Triangular Meshes. ACM Trans. Graph. (2005).
• Kostrikov et al.  I. Kostrikov, Z. Jiang, D. Panozzo, D. Zorin, and J. Bruna. 2018. Surface Networks. CVPR (2018).
• Lei et al.  N. Lei, X. Zheng, J. Jiang, Y.-Y. Lin, and D. X. Gu. 2017. Quadrilateral and hexahedral mesh generation based on surface foliation theory. Comput. Meth. Appl. Mech. Eng. (2017).
• Li et al.  B. Li, X. Li, K. Wang, and H. Qin. 2013. Surface mesh to volumetric spline conversion with generalized polycubes. IEEE Trans. Vis. Comput. Graph. (2013).
• Li et al.  Y. Li, Y. Liu, W. Xu, W. Wang, and B. Guo. 2012. All-hex meshing using singularity-restricted field. ACM Trans. Graph. (2012).
• Lipnikov et al.  K. Lipnikov, G. Manzini, and M. Shashkov. 2014. Mimetic finite difference method. J. Comput. Phys. (2014).
• Livesu et al.  M. Livesu, N. Vining, A. Sheffer, J. Gregson, and R. Scateni. 2013. PolyCut: monotone graph-cuts for PolyCube base-complex construction. ACM Trans. Graph. (2013).
• Manzini et al.  G. Manzini, A. Russo, and N. Sukumar. 2014. New perspectives on polygonal and polyhedral finite element methods. Math. Meth. Appl. Sci. (2014).
• Maréchal [2009a] L. Maréchal. 2009a. Advances in Octree-Based All-Hexahedral Mesh Generation: Handling Sharp Features.
• Maréchal [2009b] L. Maréchal. 2009b. Advances in octree-based all-hexahedral mesh generation: handling sharp features. In proceedings of the 18th International Meshing Roundtable.
• Martin  S. Martin. 2011. Flexible, unified and directable methods for simulating deformable objects. (2011). Phd thesis, ETH Zurich.
• Martin et al.  S. Martin, P. Kaufmann, M. Botsch, M. Wicke, and M. Gross. 2008. Polyhedral Finite Elements Using Harmonic Basis Functions. In Proceedings of the Symposium on Geometry Processing.
• Martin and Cohen  T. Martin and E. Cohen. 2010. Volumetric parameterization of complex objects by respecting multiple materials. Comput. Graph. (2010).
• Nieser et al.  M. Nieser, U. Reitebuch, and K. Polthier. 2011. CubeCover - Parameterization of 3D Volumes. Comput. Graph. Forum (2011).
• Owen and Saigal  S. J. Owen and S. Saigal. 2000. H-Morph: an indirect approach to advancing front hex meshing. Int. J. Numer. Methods Eng. (2000).
• Panetta et al.  J. Panetta, Q. Zhou, L. Malomo, N. Pietroni, P. Cignoni, and D. Zorin. 2015. Elastic Textures for Additive Fabrication. ACM Trans. Graph. (2015).
• Petra et al. [2014a] C. G. Petra, O. Schenk, and M. Anitescu. 2014a. Real-time stochastic optimization of complex energy systems on high-performance computers. IEEE Computing in Science & Engineering (2014).
• Petra et al. [2014b] C. G. Petra, O. Schenk, M. Lubin, and K. Gärtner. 2014b. An augmented incomplete factorization approach for computing the Schur complement in stochastic optimization. SIAM J. Sci. Comput. (2014).
• Sederberg et al.  T. W. Sederberg, D. L. Cardon, G. T. Finnigan, N. S. North, J. Zheng, and T. Lyche. 2004. T-spline Simplification and Local Refinement. ACM Trans. Graph. (2004).
• Shepherd and Johnson  J. F. Shepherd and C. R. Johnson. 2008. Hexahedral Mesh Generation Constraints. Eng. Comput. (2008).
• Shewchuk  J. R. Shewchuk. 1996. Triangle: Engineering a 2D quality mesh generator and Delaunay triangulator.
• Si  H. Si. 2015. TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator. ACM Trans. Math. Softw. (2015).
• Sokolov et al.  D. Sokolov, N. Ray, L. Untereiner, and B. Lévy. 2016. Hexahedral-Dominant Meshing. ACM Trans. Graph. (2016).
• Staten et al.  M. L. Staten, S. J. Owen, and T. D. Blacker. 2005. Unconstrained Paving & Plastering: A New Idea for All Hexahedral Mesh Generation.
• Tarini et al.  M. Tarini, K. Hormann, P. Cignoni, and C. Montani. 2004. PolyCube-Maps. ACM Trans. Graph. (2004).
• Toshniwal et al.  D. Toshniwal, H. Speleers, R. R. Hiemstra, and T. J. Hughes. 2017. Multi-degree smooth polar splines: A framework for geometric modeling and isogeometric analysis. Comput. Meth. Appl. Mech. Eng. (2017).
• Wei et al.  X. Wei, Y. J. Zhang, D. Toshniwal, H. Speleers, X. Li, C. Manni, J. A. Evans, and T. J. Hughes. 2018. Blended B-spline construction on unstructured quadrilateral and hexahedral meshes with optimal convergence rates in isogeometric analysis. Comput. Meth. Appl. Mech. Eng. (2018).
• Yamakawa and Shimada  S. Yamakawa and K. Shimada. 2003. Fully-automated hex-dominant mesh generation with directionality control via packing rectangular solid cells. Int. J. Numer. Methods Eng. (2003).
• Zhang et al.  Y. J. Zhang, X. Liang, and G. Xu. 2013. A robust 2-refinement algorithm in octree or rhombic dodecahedral tree based all-hexahedral mesh generation. Comput. Meth. Appl. Mech. Eng. (2013).

## Appendix A Brief Finite Element Introduction

Many common elliptic partial differential equations have the general form

 F(x,u,∇u,Δu)=f(x),x∈Ω,

subject to

where is the surface normal, is the Dirichlet boundary where the function is constrained (e.g., positional constraints) and is the Neumann boundary where the gradient of the function is constrained. The most common PDE in this class is the Poisson equation .

#### Weak Form

The first step in a finite element analysis consists of introducing the weak form of the PDE: find such that

 ∫ΩF(x,u,∇u,Δu)v(x)\ddx=∫Ωf(x)v(x)\ddx,

holds for any test function vanishing on the boundary. This reformulation has two advantages: (1) it simplifies the problem, and (2) it weakens the requirement on the function . For instance, in case of the Poisson equation, the strong form is well defined only if is twice differentiable, which is a difficult condition to enforce on a discrete tesselation. However, the weak form requires only that the second derivatives of are integrable, allowing discontinuous jumps. Using integration by parts it can be further relaxed to

 ∫Ω∇u(x)⋅∇v(x)\ddx=∫Ωf(x)v(x)\ddx,

where only the gradient of needs to be integrable, that is , and can thus be represented using piecewise-linear basis functions.

#### Basis Functions

The key idea of a finite element discretization is to approximate the solution space via a finite number of basis functions , , which are independent from the PDE we are interested in. The number of nodes (and basis functions) per element and their position is directly correlated to the order of the basis, see Figure 26. We note that the nodes coincide the mesh vertices only for linear basis functions. Instead of solving the PDE, the goal becomes finding the coefficients , of the discrete function that approximates the unknown function . For a linear PDE this results in a linear system , where is the stiffness matrix, captures the boundary conditions, and is the vector of unknown coefficients . For instance, for the Laplace equation the entries of the stiffness matrix are

 Kij=∫Ω∇ϕi(x)⋅∇ϕj(x)\ddx.

#### Local Support

Commonly used basis functions are locally supported. As a result, most of the pairwise intgerals are zero, leading to a sparse stiffness matrix. The pairwise integrals can be written as a sum of integrals over the elements (e.g., quads or hexes) on which both functions do not vanish. This representation enables so-called per-element assembly: for a given element, a local stiffness matrix is assembled.

For instance, if and element has four non-zero basis functions , , , (this is the case for linear quad) the local stiffness matrix for the Poisson equation is

 KLo,p=∫C∇ϕn(x)⋅∇ϕm(x)\ddx,

where and . By using the mapping of local indices to global indices , the local stiffness matrix entries are summed to yield the global stiffness matrix entries.

#### Geometric Mapping

The final piece of a finite element discretization is the geometric mapping . The local integrals need to be computed on every element. The element stiffness matrix entries are computed as integrals over a reference element (e.g., a regular unit square/cube) through change of variables

 ∫C∇ϕn(x)⋅∇ϕm(x)\ddx=∫^C(Dg−T∇^ϕn(x))⋅(Dg−T∇^ϕm(x))\absDg\ddx,

where is the Jacobian matrix of the geometric mapping , and are the bases defined on the reference element . While usually is expressed by linear combination of , leading to isoparametric elements, the choice of is independent from the basis.

All integrals are computed numerically by means of quadrature points and weights, which translates the integrals into weighted sums. Although there are many strategies to generate quadrature data (e.g., Gaussian quadrature), all of them integrate exactly polynomials up to a given degree to ensure an appropriate approximation order. For instance, if we use one quadrature point in the element’s center with weight 1, we can integrate exactly constant functions.

#### Right-hand Side

The setup of the right-hand side is done in a similar manner: its entries are .

Dirichlet boundary conditions are treated as constrained degrees of freedom. The Neumann boundary conditions are imposed by setting

 bj=∫∂ΩNϕj(x)⋅n(x)\ddx

for any node in .

As for the stiffness matrix assembly the basis and node construction for the right-hand side is performed locally.

## Appendix B Test Functions

In our experiments, we use the test functions proposed in [Franke, 1979]. For 2D:

 f2D(x1,x2) =34e−(9x1−2)2+(9x2−2)24+34e−(9x1+1)249−9x2+110 +12e−(9x1−7)2+(9x2−3)24−15e−(9x1−4)2−(9x2−7)2,

and 3D:

 f3D (x1,x2,x3)= 34e−(9x1−2)2+(9x2−2)2+(9x3−2)24+34e−(9x1+1)249−9x2+110−9x3+110 +12e−(9x1−7)2+(9x2−3)2+(9x3−5)24−15e−(9x1−4)2−(9x2−7)2−(9x3−5)2.

## Appendix C Basis functions

We use Lagrange tensor product function to interpolate between the nodes in quadrilateral and hexahedral elements. We provide the explicit formulation for and both in 2D, the 3D formulation follows. The four linear bases are constructed from the 1D linear bases

 α1(t)