The interoperability limitation of Computer Aided Design (CAD) and Finite Element Analysis (FEA) systems has become one of the major bottlenecks in simulation-based design. CAD and FEA are inherently incompatible because they use for historical reasons different mathematical representations. As advocated in isogeometric analysis the use of identical basis functions for CAD and FEA can facilitate their integration. Today most of the research on isogeometric analysis focuses on NURBS Hughes et al. (2005); Cottrell et al. (2009) and the related t-splines Bazilevs et al. (2010) and subdivision basis functions Cirak et al. (2002). The inherent tensor-product structure of NURBS means that additional techniques are required for geometries that are composed out of several NURBS patches. Specifically, around extraordinary (or irregular) points where the number of patches that join together is different than four, i.e. , alternative techniques are necessary to maintain smoothness. One prevalent approach in geometric design is to introduce additional higher order patches around the extraordinary point and to ensure that all patches match up continuously at their boundaries. refers to the notion of geometric continuity and, for instance, implying tangent plane continuity. As recently pointed out by Groisser et al. Groisser and Peters (2015) in isogeometric analysis leads to
continuity because the geometry and field variables are interpolated with the same basis functions. The utility of classicalconstructions in isogeometric analysis has recently been demonstrated in a number of papers Nguyen et al. (2014, 2016). constructions have also been explored in the context of isogeometric analysis with t-splines Scott et al. (2013). A different approach for dealing with extraordinary points is provided by subdivision surfaces. The neighbourhood of the extraordinary point is replaced by a sequence of nested continuous patches which join continuously at the point itself Peters and Reif (2008); Zorin and Schröder (2000). Subdivision basis functions for finite element analysis have originally been proposed in Cirak et al. (2000) and have been more intensely studied in a number of recent papers Cirak and Long (2011); Jüttler et al. (2016); Wei et al. (2015).
We introduce in this paper an isogeometric analysis technique that builds on manifold-based basis functions for geometric modelling and analysis. As known from differential geometry, manifolds provide a rigorous framework for describing and analysing surfaces with arbitrary topology; see do Carmo (1976); Schutz (1980). Informally, with manifolds a surface in Euclidean space is constructed from mapping and blending together planar patches from . Manifold techniques for mesh-based construction of smooth continuous surfaces were first considered in Grimm et al. Grimm and Hughes (1995). Other mesh-based manifold constructions have later been proposed, e.g., in Navau and Garcia (2000); Ying and Zorin (2004); Gu et al. (2006); Della Vecchia et al. (2008). In a manner similar to the description of splines, a continuous surface is described with a quadrilateral or triangular control mesh and each vertex has a corresponding basis function with a local support, see Figure 1. In contrast to the aforementioned constructions which rely on matching up separate patches, in manifolds a continuous surface is created by smoothly blending of overlapping patches. The idea of blending surfaces from overlapping patches is a common theme in geometric modelling and has been used, for instance, for increasing the smoothness of subdivision surfaces around the extraordinary vertices Levin (2006); Zorin (2006); Antonelli et al. (2013) or (meshfree) point-based surface processing Levin (2004); Pauly et al. (2002). In Millan et al. Millán et al. (2010, 2013) point-based surface blending techniques have also been used for meshfree thin-shell analysis.
In the present work we follow Ying and Zorin Ying and Zorin (2004) and construct smooth basis functions by combining manifold techniques with conformal parameterisations and partition of unity method. The control mesh consists of quadrilateral elements with some extraordinary vertices (i.e. for some non-boundary vertices) and the construction gives one basis function for each vertex. The first step is to assign each vertex of the control mesh and its adjacent elements a planar sub-mesh with the same connectivity. The sub-meshes serve as control meshes for planar surface patches, which can be understood as parameter spaces for basis functions. For continuous basis functions the planar patches have to have a smooth parameterisation. Although other choices are conceivable, the patches are parameterised using conformal (angle-preserving) maps. Since each surface point is represented on several patches transition functions composed of conformal maps are used to navigate between adjacent patches. In the second step of the construction, on each planar patch the conventional partition of unity method (PUM) of Melenk et al. Melenk (1995); Melenk and Babuska (1996) is used for constructing basis functions. According to PUM, the basis functions are the product of a partition of unity function and a patch specific polynomial approximant. In computer graphics literature the partition of unity function and the patch specific polynomial basis are usually referred to as the blending function and the embedding function, respectively. We use as partition of unity functions b-splines that have zero value and zero derivatives at the patch boundaries. In order to enforce partition of unity the b-splines on different patches overlapping the same point on the manifold surface are first identified with transition functions and subsequently normalised as in usual PUM. The last step in the basis function construction is to express the local PUM polynomial approximant in dependence of vertex values using a least-squares approximation. In this mesh-based approach the degree of the polynomial approximant and the number of vertices in a sub-mesh are correlated. In order to increase the polynomial degree the sub-meshes are enlarged with mesh refinement by quadrisectioning. The basis functions depend only on the connectivity of the control mesh but not its geometry so that they can be precomputed and tabulated for different valences .
The outline of this paper is as follows. Section 2 reviews the relevant manifold concepts from differential geometry and the partition of unity method. The mesh-based manifold basis functions are introduced in Section 3. First, one-dimensional polygonal control meshes are considered, even though it is straightforward to combine one-dimensional manifolds with the partition of unity method by using simple transition functions. Subsequently, two-dimensional quadrilateral meshes are considered for which conformal maps are used as transition functions. In both dimensions, it is shown how the polynomial coefficients in the partition of unity approximation can be expressed as a linear combination of vertex coefficients. In Section 4 the derived mesh-based basis functions are applied to a number of Poisson and thin-shell problems. Numerical convergence with increasing mesh size on meshes with and without extraordinary vertices is investigated.
2.1 Review of manifold concepts
In the following we provide an informal introduction to differentiable manifolds with the aim to introduce the necessary terminology. For clarity, our discussion is restricted to surfaces, i.e. two-dimensional manifolds, embedded in the three-dimensional Euclidean space . A similar introduction, but oriented more towards geometric modelling, can be found in Grimm and Zorin (2006). The manifold concept is much more general than needed in this paper. A more rigorous discussion is found in most standard differential geometry textbooks, see, e.g., do Carmo (1976); Schutz (1980).
A regular two-manifold, or surface, is defined as the set of points in , which can be locally continuously one-to-one mapped onto a set of points in . This definition naturally excludes surfaces with t-joints or isolated points. In applications this is, however, not a major restriction since a geometry, for instance, with a t-joint can be represented as a collection of several manifolds.
By definition on a regular surface around each point there is an open region that can be mapped to an open planar region , see Figure 2. In line with the partition of unity method terminology, in the following will be referred to as a patch and as a planar patch. Furthermore, we denote the function for mapping between two patches with . On the surface there are many overlapping patches such that
and each point lies at least on one patch. The pair consisting of is called a chart. The set of all charts is referred to as an atlas for representing the surface .
As illustrated in Figure 2, each planar patch has its own coordinate system . The same point in the intersection between the two patches has the coordinates in and in . In order to compute the underlying coordinate transformations we introduce the transition functions
that are composed out of the mappings . The transition functions are symmetric and satisfy the cocyle condition when the preimage of lies in three planar patches , and . For a surface to be continuous the transition maps must be continuous. Evidently, for a differentiable surface has to be equal or larger than one.
2.2 Review of the Partition of Unity Method (PUM)
The first use of PUM for creating finite element basis functions goes back to the seminal work of Babuska et al. Babuška et al. (1994) and was subsequently further developed, for instance, in Melenk et al. Melenk and Babuska (1996) and Duarte et al. Duarte and Oden (1996). In PUM a given domain in the Euclidean space , with , is partitioned into overlapping patches such that
As opposed to the manifolds introduced in previous Section 2.1, there is only one single coordinate system in the Euclidean space and a point has the same coordinates on and all patches . Hence, the transition functions between the different patches are identity maps.
Next, a blending function is defined on each patch . By definition the sum of the blending functions over all the regions is
The set of blending functions is also referred to as the partition of unity subordinate to the set of open patches . As will become clear, in order to obtain smooth PUM basis functions the blending functions have to be smooth. In addition, on the region boundaries the function value and derivatives of have to be zero. After choosing on each patch an arbitrary functions that have the prerequisite properties they can be normalised to yield a blending function
In our applications are usually b-spline basis functions.
On each patch , in addition to the blending functions a local polynomial approximant is considered
is vector containing a complete polynomial basis,is the corresponding vector of the coefficients and the dot represents their scalar product. For instance, for one-dimensional domains and a monomial basis the two vectors are of the form
The global approximant is the sum of the local approximants and their multiplication with the blending functions , that is,
The smoothness of this function depends on the smoothness of and . Since the polynomials in are infinitely smooth, the smoothness of is exclusively controlled by the blending functions
. According to the convergence estimates given inMelenk and Babuska (1996), the convergence rates for the approximant depend on the degree of the polynomial basis and the constants on the blending functions and the layout of the overlaps .
The illustrative one-dimensional example in Figure 3 showcases the construction of a smooth function using PUM. On each patch one single cubic b-spline basis functions is used as for defining the (normalised) blending function according to (5). The support size of a cubic b-spline is equal to the size of its corresponding patch . The cubic b-spline is a continuous function and at the boundaries of its value and first and second derivatives are zero Farin (2002). In addition to the blending functions , in Figure 3 the local quadratic polynomials are also shown. The smooth function constructed according to (9) is shown in Figure 3(b). It is evident that this curve is by construction continuous.
For the considered one-dimensional example, it is straightforward to use a manifold-based approach for constructing the partition of unity function. As shown in Figure 4, the coordinate systems on each patch can be chosen differently. For switching between the different coordinate systems the transition functions are used. In this example one of the chosen transition functions, i.e. , is quadratic and the other, i.e. ,is linear. It is however possible to chose any other monotone function. Choosing other transition functions will lead to a change in the shape of the constructed function .
3 Mesh-based manifold basis functions
We are now in a position to introduce the construction of manifold-based basis functions on one- and two-dimensional meshes. The idea of using manifolds for smooth interpolation on meshes has been originally introduced in computer graphics by Grimm et al. Grimm and Hughes (1995). The one-dimensional case is straightforward and is only discussed in order to provide some intuition for the two-dimensional case. The essential difficulty in two-dimensions lies in defining suitable transition functions. We use the conformal maps as introduced in Ying and Zorin Ying and Zorin (2004) for defining the transition functions. Alternative definitions have been provided in Grimm and Hughes (1995); Navau and Garcia (2000); Antonelli et al. (2013).
3.1 One-dimensional meshes
First we aim to construct a smooth curve, i.e. a one-dimensional smooth manifold for a given coarse control polygon in the Euclidean space . We begin with defining charts for each vertex, see Figure 5. The planar patch is formed from two segments and the attached three vertices. This is usually referred to as the one-ring of the centre vertex. It is also possible to increase the size of to a two-ring or even larger. The chosen size of the patches influences the number of overlapping patches at each point. As typical for manifolds, each planar patch has its own coordinate system. The scalar transition functions enable to navigate between the patches and . The transition functions are chosen as linear maps.
The coordinates of a point on the smooth curve is now determined with the partition of unity method. According to (9) we can write, for instance, for the coordinate of the point with a preimage on the planar patch and coordinate
The summation is over all patches and in order to evaluate the sum it is necessary to use the transition maps. Although we give here and in the following only the expression for , the other two coordinates and are expressed similarly. On each one-ring there are three vertices, which motivates the choice of a quadratic basis for . For a quadratic Lagrangian basis the three coefficients are simply the coordinates of the vertices in the one-ring:
Both vectors and have three components and the entries of are the coordinates of the three control polygon vertices. This equation is illustrated in Figure 5. As can be seen, the smooth curve passes exactly through the vertices of the control polygon.
Next, we rewrite equation (11) in index notation to define basis functions that can be used for finite element analysis:
where are the three basis functions corresponding to the three vertices in the patch . In Figure 6 the non-zero basis functions and their derivatives in one patch are shown. In the underlying construction the blending functions are normalised cubic b-splines, local polynomials are quadratic and the transition function are linear. Note that the support size of one basis function is two elements. Due to the overlaps there are four non-zero basis functions in one element. The resulting basis functions are continuous.
3.2 Two-dimensional quadrilateral meshes
We now consider the construction of a smooth surface, i.e. smooth two-manifold, for a given coarse control mesh. Our approach follows closely the construction originally introduced in Ying and Zorin Ying and Zorin (2004). Although only quadrilateral meshes are considered, it is straightforward to extend the technique to triangular meshes.
In addition to the smoothness properties of blending functions and local polynomials, the smoothness of transition functions is central in generating smooth surfaces. For implementation purposes, it is also important that the transition functions and their inverses are readily computable. Similar to the one-dimensional construction, we define charts for each vertex of the mesh. The planar patch is chosen for now as the one-ring of elements around a vertex. The number of elements in the one-ring of a vertex is referred to as the valence of the vertex. On structured meshes all vertices inside the domain have valence and on unstructured meshes it can be arbitrary. Hence, in the unstructured case the overlapping one-ring patches can have different valences, which makes the computation of a smooth transition function challenging. Ying et al. Ying and Zorin (2004) proposed conformal maps as smooth and easy computable transition functions. Recall here that conformal mapping is an angle preserving transformation. The generation of a conformal parameterisation for a one-ring of elements proceeds in several steps. In Figure 7 the procedure for a vertex with valence is illustrated. The smooth parameterisation is obtained by conformally mapping, rotating and combining unit squares. The points of the unit square have the coordinates and are expressed as a complex number . The conformal transformation maps the square to a wedge. In computing the mapping recall the following standard relations:
After the mapping the coordinates of a point become according to:
The wedge-shaped image of the conformally mapped unit square forms one sector of the one-ring patch. This wedge is copied and suitably rotated to form a smoothly parameterised one-ring forming patch . The angle preservation property of the conformal map ensures that the parameter lines on are smoothly connected across element edges, see Figure 7. For the sequence of transformations from the unit square to the smoothly parameterised patch we abstractly write
where the non-linear mapping is comprised of the conformal transformation and a rotation. It it worth emphasising that it is straightforward to compute the inverse and derivatives of the mapping (15). Moreover, depends only on the valence of the considered patch, but not on the vertex coordinates. See Figure 8 for conformal parameterisation of one-ring patches with valences .
The functions for computing the normalised blending functions are also assembled from smooth functions defined on unit squares. In our computations is one quarter of a tensor-product b-spline and covers the entire unit square. Figure 9 shows the procedure for constructing the blending function on a valence five patch. It can be seen that the b-spline has its maximum at the corner which maps to the central vertex of the one-ring patch . To normalise the functions , for a given point on patch the corresponding point on an overlapping patch is computed with transition functions, that is, , cf. (15). The blending function is then first mapped onto the conformal wedge and then appropriately rotated to construct the blending function on the patch. In Figure 10 the normalised blending functions for one-ring patches with are shown.
In Ying and Zorin (2004) the weight functions are chosen such that they are constant in a small neighbourhood of width close to the unit square boundaries. This is motivated by the need to circumvent the singularity of the conformal map at the central vertex of the one-ring patch. Our numerical experiments indicate that the finite element solutions are insensitive to the choice of so that is chosen. Note that for evaluating the finite element integrals the surface is only evaluated at quadrature points, which are usually away from the vertices. Alternatively, the one-ring patches can also be parameterised with the characteristic map of Catmull-Clark subdivision surfaces, which is smooth and does not have a singularity at the extraordinary vertex, see Antonelli et al. (2013); Levin (2006) for details.
We now proceed to the construction of a smooth surface, i.e. a two-dimensional smooth manifold, for a given coarse control mesh. The overall approach is very similar to the one-dimensional case introduced in the previous section. We consider again a point on the manifold with a preimage with the coordinates on the planar patch . Each of the coordinate components of are interpolated with the partition of unity method. We write, for instance, for the coordinate
where is a vector containing the components of a polynomial basis and are the unknown coefficients. Next, the coefficients are expressed in dependence of the vertex coordinates of the control mesh. The dimension of the polynomial basis has to be equal or smaller than the number of vertices in a patch . Because on unstructured meshes the valence of vertices is not fixed, we choose a polynomial basis based on the one-ring with the smallest valence in the mesh. If a higher degree polynomial is desired, two or more rings of elements can be considered as patches . When the dimension of the polynomial basis is smaller than the number of vertices in the patch , the polynomial coefficients
are determined with a local least-squares projection on each patch. Neglecting for the moment the patch index, and denoting the parametric coordinates of thevertices on patch as and the corresponding nodal control mesh coordinates as , the least-squares fit on reads
Note for a specific polynomial basis the matrix on the left hand side depends only on the valence of the one-ring and can be precomputed and stored. We abbreviate the least-squares projection (17) with
where is the projection matrix and the vector contains the coordinates of all the vertices in the one-ring, i.e. on the planar patch.
This gives rise to the definition of basis functions , where is the vertex id on patch . In conventional finite element implementations usually system matrices and vectors are evaluated by iterating over the elements in the control mesh. Moreover, during numerical integration the basis function values at pretabulated points in the integration element are needed. In manifold finite elements the unit square in the coordinate system is chosen as the integration element. Hence, the basis functions have to be evaluated in a given element and integration point .
We briefly consider the hand geometry shown in Figure 1 for illustrating the process of evaluating the basis functions; see A for more details. On the control mesh in Figure 1 (left) one element is highlighted and the four one-rings belonging to its four vertices are indicated. One of the element’s vertices has valence and the other three have valence . The interpolation within the highlighted element depends on the eighteen vertices in the union of the four one-rings. In Figure 11 the four patches used for partition of unity construction are shown. The unit square in the centre represents the integration element, equivalent to the parent element in isoparametric finite elements. The conformal mapping of the unit square to the four planar one-ring patches depends on the valence of the respective one-ring. For the three patches with valence four, the mapping is (up to some rotations) essentially an identity map. With the transition maps implied by the mappings it is straightforward to compute the basis functions defined in (19) for a given integration point . In this specific example the minimum number of vertices in a patch is nine so that the polynomial basis can be chosen either as a bilinear or biquadratic Lagrangian basis. Note that on patches with valence and a biquadratic Lagrangian basis the least-squares projection matrix
is an identity matrix. Moreover, as also can be deduced from Figure11 each basis function has a support size consisting of two rings of elements around its associated vertex. In Figure 12 basis functions for vertices with valence are shown.
We consider second and fourth-order partial differential equations to demonstrate the accuracy and convergence of the introduced manifold basis functions when used in finite element analysis. For manifold construction, around each vertex patches consisting of either one or two-rings of elements are considered. As blending functions we use either normalised linear, quadratic or cubic b-splines. The element integrals are evaluated withGauss integration points in all examples. This high number of integration points has been chosen in order to minimise the effect of integration errors. In convergence studies the control meshes are refined with the Catmull-Clark subdivision scheme Catmull and Clark (1978). The number of extraordinary vertices in a mesh remains constant because the new vertices introduced during the refinement are all ordinary.
Surfaces with boundaries require modified charts for manifold constructions, such as those introduced in Tosun and Zorin (2011). The specialised treatment of elements close to the boundary can be avoided by introducing ghost elements just outside the domain. This is achieved by reflecting sufficient number of internal elements and vertices along the boundary. Furthermore, the manifold basis functions are non-interpolating at the boundaries. Therefore, we use the penalty method for applying Dirichlet boundary conditions.
4.1 Two-dimensional Poisson problems
4.1.1 Square domain with a structured mesh
As an introductory example, we solve the Poisson-Dirichlet problem on the domain , discretised with a Cartesian grid, see Figure 13. The loading is chosen such that the analytical solution is
where the variables represent coordinates.
In the partition of unity construction we use a biquadratic Lagrangian basis as the local polynomial basis and as blending functions we consider normalised linear, quadratic and cubic b-splines. To begin with, the number of Gauss points required for adequate integration is determined, see Figure 14. In Figure 14 only normalised cubic b-splines are used. Since the basis functions are rational, a large number of Gauss points is unavoidable, especially for finer meshes. Furthermore, our findings indicate that blending functions with higher continuity, in general, require more Gauss points. According to Figure 14, the chosen Gauss points for all examples in the paper seems to provide a good trade-off between accuracy and efficiency.
Figure 15 shows the norm and semi-norm of error as the mesh is uniformly refined with the Catmull-Clark scheme. For this structured mesh Catmull-Clark refinement is equivalent to refinement by bisection. The error norms for three different blending functions, namely normalised linear, quadratic and cubic splines, are shown. It can be inferred from these convergence plots that the convergence rates are optimal and are unaffected by the blending functions. Interestingly, the constants in the convergence plots increase when the smoothness of the blending function is increased.
4.1.2 Square domain with an unstructured mesh
This example underlines the performance of manifold basis functions on unstructured meshes, and studies how the convergence rates are influenced in the presence of extraordinary vertices. To this end, the Poisson-Dirichlet problem with the analytical solution
is considered. Figure 16 shows the unstructured mesh used in the computations and a representative finite element solution. The mesh lines on the displaced solution in Figure 16 represent the edges of the elements on the exact surface. The mesh has eight extraordinary vertices, with four vertices of valence and the other four of valence .
In the convergence studies we use as local polynomials bilinear and biquadratic Lagrangian polynomials. As patches one- and two-rings of elements are considered. However, in one-ring patches with valence there are only seven vertices so that instead of a biquadratic Lagrangian polynomial locally a complete quadratic polynomial has to be used. In all cases the blending functions are normalised cubic b-splines.
Figure 17 shows the norm and semi-norm of the error as the mesh is successively refined with the Catmull-Clark scheme. It can be seen that for all cases, the convergence rates are close to optimal. In the norm and for quadratic polynomials () the convergence rates for one- and two-ring patches are approximately and respectively. In the semi-norm, the corresponding convergence rates are and respectively. We believe that the reduction of convergence rates with increasing patch size is primarily due to the suboptimal integration of rational polynomials.
Next, we study the pointwise convergence of the solution at selected vertices. For this part of the studies we use one- and two-ring patches and consider only biquadratic Lagrangian polynomials. As mentioned before for one-ring patches with valence instead of the biquadratic Lagrangian a complete quadratic has to be used. The blending functions are normalised cubic b-splines. Figure 18 (left) shows the convergence of the error at three selected vertices with valences when one-ring patches are used. Figure 18 (right) shows the corresponding plots when two-ring patches are used. For both types of patches, the convergence rate at extraordinary vertices is . It is interesting to note that two-ring patches, in general, have smaller errors than the one-ring patches.
4.1.3 Circular domain
In this example the Poisson-Dirichlet problem with the same solution as in (21) is solved on a circular domain with radius , see Figure 19. The aim is to illustrate the treatment of curved boundaries when manifold basis functions are used. As previously mentioned, in our present implementation, ghost nodes are introduced at the boundaries, which circumvents the use of modified charts close to the boundaries. With specialised boundary charts the treatment of boundaries would be different.
Figure 19 (left) shows the coarse mesh containing four extraordinary vertices with valence . The domain boundary is approximated by least-squares fitting an approximate circle described by the manifold basis functions to the exact circle. In the least-squares problem the unknowns are the positions of the vertices close to the boundary. This is performed as a preprocessing step each time after the mesh is refined with Catmull-Clark subdivision.
In the convergence study shown in Figure 20 only one-ring patches are used. The local polynomials are either biquadratic Lagrangian polynomials for patches with valence or complete quadratic polynomials for patches with . The blending functions are normalised cubic b-splines. As can be inferred from Figure 20, optimal convergence rates for the norm and the semi-norm error are achieved.
4.2 Thin-plate and thin-shell problems
The linear Kirchhoff-Love model is used for computing the thin-plate and thin-shell problems introduced in the following. The corresponding weak form depends on the metric and curvature tensors of the mid-surface in the reference and deformed configurations. Due to the presence of the curvature tensor, the basis functions have to be smooth, or more precisely in space . The presented computations have been performed by replacing the subdivision basis functions in our software Cirak and Long (2011); Long et al. (2012) with manifold-basis functions. Due to the algorithmic similarities between subdivision and manifold basis functions it is straightforward to replace one with the other. Out-of-plane shear deformations have been neglected although it would be possible to take them into account as proposed in Long et al. (2012).
4.2.1 Simply supported square plate
We consider the deformation of a simply supported square plate of unit length subjected to an applied pressure loading . The thickness of the plate is , the Young’s modulus is , and the Poisson ratio is . Its analytical solution according to Timoshenko and Woinowsky-Krieger (1964) is
where is the flexural rigidity.
The unstructured mesh shown in Figure 16, previously used for the Poisson-Dirichlet problem, is reused for this plate problem. Recall that the mesh has eight extraordinary vertices, namely four with and the other four with . In the convergence study, one-rings and two-rings of elements are considered. Except on valence one-ring patches, the local polynomials are biquadratic Lagrangian polynomials. As previously mentioned, on one-ring patches with only a complete quadratic polynomial can be used, because there are only seven vertices in a patch. In all cases the blending function is a normalised cubic b-spline. In Figure 21 (left), the deflected shape of the simply-supported plate for a relatively coarse mesh is shown. Figure 21 (right) illustrates the convergence of the norm error as the mesh is successively refined with the Catmull-Clark scheme. The constant in the convergence plot decreases with increasing patch size. For both patch sizes, the convergence rate is approximately , which is slightly lower than the optimum of . One possible reason for this is the inadequate integration of the rational polynomials.
4.2.2 Pinched cylinder
The pinched cylinder is one of the benchmark examples for shell finite elements suggested in Belytschko et al. Belytschko et al. (1985). The unstructured coarse control mesh, material properties and the general problem setup are shown in Figure 22 (left). The ends of the cylinder are unconstrained and the two diametrically opposite forces are applied within the middle section of the cylinder.
The deflected pinched cylinder is shown in Figure 22 (right). In Figure 23 the convergence of the normalised maximum change in diameter of the Catmull-Clark and manifold finite element solutions are compared, both structured and unstructured meshes are considered. The change in diameter is normalised with the exact solution for a membrane shell given in Timoshenko and Woinowsky-Krieger (1964). It can be seen that the Catmull-Clark and manifold solutions converge to very similar values. The convergence of manifold functions with cubic polynomials () on two-ring patches is comparable with subdivision basis functions. Note that for structured meshes Catmull-Clark subdivision basis functions are identical to tensor-product cubic b-splines Peters and Reif (2008). The slower convergence of the manifold basis functions with quadratic polynomials () on one-ring patches is as expected.
4.2.3 Pinched hemisphere
Our last example is the pinched hemisphere, which has also been suggested in Belytschko et al. Belytschko et al. (1985) as a benchmark for shell finite elements. The coarse control mesh, material properties and the general problem setup are shown in Figure 24 (left). The edge of the hemisphere is unconstrained and the four radial forces have alternating signs. The sum of the applied forces is zero. In the control mesh the valence of the vertices range between 3 and 5.
The deformed surface of the pinched hemisphere is shown in Figure 24 (right). Figure 25 shows the convergence of the normalised maximum radial displacement. The displacements are normalised by , given in Belytschko et al. Belytschko et al. (1985). The same plot also includes the convergence of the finite element solution when Catmull-Clark basis functions are used. As in case of pinched cylinder example, the manifold basis functions are constructed with quadratic local polynomials () on one-rings and cubic polynomials () on two-ring patches. In both cases, normalised cubic b-splines were used as blending functions. Figure 25 illustrates that the manifold basis functions constructed using cubic polynomial patches converges slightly faster than the Catmull-Clark subdivision basis functions.
We introduced an isogeometric analysis technique that uses manifold smooth basis functions on quadrilateral meshes. Manifold techniques have a long history in computer graphics and computer aided design and several variants have been proposed over the years. Our implementation closely follows Ying and Zorin Ying and Zorin (2004) and combines manifold techniques with conformal parameterisations and the partition of unity method. The smoothness of the basis functions is determined by the smoothness of the blending, or partition-of-unity, functions. The approximation properties of the basis functions is mainly determined by the polynomial degree used in each patch. In the presented computations, the blending function was chosen either as a normalised linear, quadratic or cubic b-spline leading to , or
continuous basis functions, respectively. As patch sizes for the manifold construction, one- or two-ring layers of elements around a vertex were considered. The number of vertices in a patch determines the maximum degree of the local polynomial that can be used in the partition-of-unity interpolation. The polynomial coefficients in the partition of unity interpolation are expressed as vertex coefficients using a least-squares procedure. The finally obtained basis functions are smooth, rational, locally supported and are associated to vertices in the mesh (similar to b-splines of odd degree). The near optimal convergence of the introduced basis on meshes with extraordinary vertices could be numerically confirmed.
For future research the combination of manifold basis functions with b-splines and the related NURBS or subdivision surfaces appears especially promising. B-splines have several compelling properties, including refinability, that make them ideal for geometric modelling and numerical analysis on meshes with no extraordinary vertices. However, most b-spline techniques for dealing with extraordinary vertices, including subdivision and many constructions, do not lead to optimally convergent finite elements Jüttler et al. (2016); Nguyen et al. (2014, 2016). In contrast, as numerically demonstrated manifold basis functions yield nearly optimally convergent finite elements independent of the connectivity of the mesh. This suggests to use b-splines in most of the domain and to introduce manifold basis functions only around extraordinary vertices. In subdivision surfaces manifold techniques have already been used to obtain continuity around extraordinary vertices Levin (2006); Zorin (2006); Antonelli et al. (2013). In these three papers, instead of the conformal map the characteristic map of subdivision surfaces is used to parameterise the planar patches. The advantage of the characteristic map, in comparison to the conformal map, is that it provides a more uniform parameterisation and does not have a singularity at the extraordinary vertex. Additional directions for future research include the mathematical convergence analysis and the proof of the linear independence of the introduced manifold basis function. To this end, the large number of results for partition of unity methods provide a good starting point.
Appendix A Implementation
In the following we discuss the implementation of manifold basis functions focusing on data structures and algorithms. For clarity, the discussion is restricted to the case of patches consisting of one-ring elements. This section should be read in conjunction with Section 3.2 and specifically Figure 11. In line with conventional finite element implementations it is assumed that the element matrices and vectors are assembled by iterating over elements. For computing the element matrices and vectors the basis function values and derivatives at integration points are required. The basis functions depend only on the connectivity of the mesh but not on the coordinates of the vertices. Hence, they can be precomputed as part of a preprocessing step and stored for later use.
The construction of the manifold basis functions proceeds in several steps. First, for each (non-boundary) vertex in the mesh the elements and vertices in its one-ring are identified. Recall that we introduced one layer of ghost elements just outside the domain of interest and that the charts which belong to the boundary vertices do not intersect the domain. The elements and vertices in an one-ring can efficiently be identified, for instance, using a half-edge data structure.
Following the assembly of the one-rings, we endow each patch with a blending function and a local polynomial basis . The degree of is chosen such that it is equal or less than the number of the vertices in the corresponding one-ring, i.e., , where is the valence of the centre vertex. Subsequently, the basis functions and their derivatives for each quadrature (or, Gauss) point in an element are precomputed. Each four-noded element in the mesh lies within the overlap of four patches, hence, with each receiving four contributions. The contribution of one patch to belonging to the vertex is computed as follows:
Determine the conformal coordinate by applying the conformal map .
Identify the image of the considered element in the patch by comparing vertex ID’s. Based on that determine the coordinate in by applying the rotation
where the rotation is a multiple of .
Compute the contribution of patch to basis function according to (19)
where the summation is over the components in the polynomial basis . We choose for either a monomial or a Lagrangian basis. The matrix , which is the inverse of the least-squares matrix, does not depend on the coordinates and can be precomputed for all possible valences and stored. For computing the contribution of patch to the derivatives the above equation is differentiated, i.e.,
After summing up the contributions of the four overlapping charts, the basis functions, their derivatives and corresponding vertices are stored in maps. The number of (non-zero) basis functions in an element is equal to the number of unique vertices in the union of the element’s four charts.
- Hughes et al. (2005) T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering 194 (2005) 4135–4195.
- Cottrell et al. (2009) J. Cottrell, T. Hughes, Y. Bazilevs, Isogeometric Analysis: Toward Integration of CAD and FEA, John Wiley & Sons Ltd., 2009.
- Bazilevs et al. (2010) Y. Bazilevs, V. M. Calo, J. A. Cottrell, J. A. Evans, T. J. R. Hughes, S. Lipton, M. A. Scott, T. W. Sederberg, Isogeometric analysis using T-splines, Computer Methods in Applied Mechanics and Engineering (2010) 229–263.
- Cirak et al. (2002) F. Cirak, M. Scott, E. Antonsson, M. Ortiz, P. Schröder, Integrated modeling, finite-element analysis, and engineering design for thin-shell structures using subdivision, Computer-Aided Design 34 (2002) 137–148.
- Groisser and Peters (2015) D. Groisser, J. Peters, Matched -constructions always yield -continuous isogeometric elements, Computer Aided Geometric Design 34 (2015) 67–72.
- Nguyen et al. (2014) T. Nguyen, K. Karčiauskas, J. Peters, A comparative study of several classical, discrete differential and isogeometric methods for solving Poisson’s equation on the disk, Axioms 3 (2) (2014) 280–299.
- Nguyen et al. (2016) T. Nguyen, K. Karčiauskas, J. Peters, finite elements on non-tensor-product 2d and 3d manifolds, Applied Mathematics and Computation 272 (2016) 148–158.
- Scott et al. (2013) M. A. Scott, R. N. Simpson, J. A. Evans, S. Lipton, S. P. A. Bordas, T. J. R. Hughes, T. W. Sederberg, Isogeometric boundary element analysis using unstructured T-splines, Computer Methods in Applied Mechanics and Engineering (2013) 197–221.
- Peters and Reif (2008) J. Peters, U. Reif, Subdivision Surfaces, Springer Series in Geometry and Computing, Springer, 2008.
- Zorin and Schröder (2000) D. Zorin, P. Schröder, Subdivision for Modeling and Animation, SIGGRAPH 2000 Course Notes, 2000.
- Cirak et al. (2000) F. Cirak, M. Ortiz, P. Schröder, Subdivision surfaces: A new paradigm for thin-shell finite-element analysis, International Journal for Numerical Methods in Engineering 47 (2000) 2039–2072.
- Cirak and Long (2011) F. Cirak, Q. Long, Subdivision shells with exact boundary control and non-manifold geometry, International Journal for Numerical Methods in Engineering 88 (2011) 897–923.
- Jüttler et al. (2016) B. Jüttler, A. Mantzaflaris, R. Perl, M. Rumpf, On numerical integration in isogeometric subdivision methods for PDEs on surfaces, Computer Methods in Applied Mechanics and Engineering 302 (2016) 131–146.
- Wei et al. (2015) X. Wei, Y. Zhang, T. J. R. Hughes, M. A. Scott, Truncated hierarchical Catmull–Clark subdivision with local refinement, Computer Methods in Applied Mechanics and Engineering 291 (2015) 1–20.
- do Carmo (1976) M. P. do Carmo, Differential geometry of curves and surfaces, Prentice-Hall, Englewood Cliffs, NJ, 1976.
- Schutz (1980) B. F. Schutz, Geometrical methods of mathematical physics, Cambridge University Press, Cambridge, UK, 1980.
- Grimm and Hughes (1995) C. M. Grimm, J. F. Hughes, Modeling surfaces of arbitrary topology using manifolds, in: SIGGRAPH 1995 Conference Proceedings, 359–368, 1995.
- Navau and Garcia (2000) J. C. Navau, N. P. Garcia, Modeling surfaces from meshes of arbitrary topology, Computer Aided Geometric Design 17 (2000) 643–671.
- Ying and Zorin (2004) L. Ying, D. Zorin, A simple manifold-based construction of surfaces of arbitrary smoothness, in: SIGGRAPH 2004 Conference Proceedings, 271–275, 2004.
- Gu et al. (2006) X. Gu, Y. He, H. Qin, Manifold splines, Graphical Models 68 (2006) 237–254.
- Della Vecchia et al. (2008) G. Della Vecchia, B. Jüttler, M.-S. Kim, A construction of rational manifold surfaces of arbitrary topology and smoothness from triangular meshes, Computer Aided Geometric Design 25 (2008) 801–815.
- Levin (2006) A. Levin, Modified subdivision surfaces with continuous curvature, in: SIGGRAPH 2006 Conference Proceedings, 1035–1040, 2006.
- Zorin (2006) D. Zorin, Constructing curvature-continuous surfaces by blending, in: Proceedings of the fourth Eurographics symposium on Geometry processing, 31–40, 2006.
- Antonelli et al. (2013) M. Antonelli, C. V. Beccari, G. Casciola, R. Ciarloni, S. Morigi, Subdivision surfaces integrated in a CAD system, Computer-Aided Design 45 (2013) 1294–1305.
- Levin (2004) D. Levin, Mesh-independent surface interpolation, in: Geometric Modeling for Scientific Visualization, Springer, 37–49, 2004.
- Pauly et al. (2002) M. Pauly, M. Gross, L. P. Kobbelt, Efficient simplification of point-sampled surfaces, in: Proceedings of the conference on Visualization’02, 163–170, 2002.
- Millán et al. (2010) D. Millán, A. Rosolen, M. Arroyo, Thin shell analysis from scattered points with maximum-entropy approximants, International Journal for Numerical Methods in Engineering 85 (2010) 723–751.
- Millán et al. (2013) D. Millán, A. Rosolen, M. Arroyo, Nonlinear manifold learning for meshfree finite deformation thin-shell analysis, International Journal for Numerical Methods in Engineering 93 (2013) 685–713.
- Melenk (1995) J. M. Melenk, On generalized finite element methods, Ph.D. thesis, The University of Maryland, 1995.
- Melenk and Babuska (1996) J. M. Melenk, I. Babuska, The partition of unity finite element method: Basic theory and applications, Computer Methods in Applied Mechanics and Engineering 139 (1996) 289–314.
- Grimm and Zorin (2006) C. Grimm, D. Zorin, Surface modeling and parameterization with manifolds, in: SIGGRAPH 2006 Course Notes, 2006.
- Babuška et al. (1994) I. Babuška, G. Caloz, J. E. Osborn, Special finite element methods for a class of second order elliptic problems with rough coefficients, SIAM Journal on Numerical Analysis 31 (1994) 945–981.
- Duarte and Oden (1996) C. A. Duarte, J. T. Oden, An h-p adaptive method using clouds, Computer Methods in Applied Mechanics and Engineering 139 (1996) 237–262.
- Farin (2002) G. Farin, Curves and Surfaces for CAGD, Academic Press, 2002.
- Catmull and Clark (1978) E. Catmull, J. Clark, Recursively generated B-spline surfaces on arbitrary topological meshes, Computer-Aided Design 10 (1978) 350–355.
- Tosun and Zorin (2011) E. Tosun, D. Zorin, Manifold-based surfaces with boundaries, Computer Aided Geometric Design 28 (2011) 1–22.
- Long et al. (2012) Q. Long, P. Bornemann, F. Cirak, Shear-flexible subdivision shells, International Journal for Numerical Methods in Engineering 90 (2012) 1549–1577.
- Timoshenko and Woinowsky-Krieger (1964) S. Timoshenko, S. Woinowsky-Krieger, Theory of plates and shells, McGraw-Hill Book Company, 2nd edn., 1964.
- Belytschko et al. (1985) T. Belytschko, H. Stolarski, W. Liu, N. Carpenter, J.-J. Ong, Stress projection for membrane and shear locking in shell finite elements, Computer Methods in Applied Mechanics and Engineering 51 (1985) 221–258.