Spline‐ and hp‐basis functions of higher differentiability in the finite cell method

In this paper, the use of hp‐basis functions with higher differentiability properties is discussed in the context of the finite cell method and numerical simulations on complex geometries. For this purpose, Ck hp‐basis functions based on classical B‐splines and a new approach for the construction of C1 hp‐basis functions with minimal local support are introduced. Both approaches allow for hanging nodes, whereas the new C1 approach also includes varying polynomial degrees. The properties of the hp‐basis functions are studied in several numerical experiments, in which a linear elastic problem with some singularities is discretized with adaptive refinements. Furthermore, the application of the Ck hp‐basis functions based on B‐splines is investigated in the context of nonlinear material models, namely hyperelasticity and elastoplasicity with finite strains.

functions. The IGA approach was originally proposed by Hughes et al. [14] in 2005 and is typically based on the use of B-splines or NURBS functions. It has shown a high accuracy per-degree-of-freedom in many applications ranging from solids and structures [15][16][17][18] to fluids [19] and fluid-structure interaction, [20] also allowing for discretizations of higher-order PDEs in primal form. [21,22] IGA has also been combined with the FCM. [5,23,24] Very similar technologies have emerged recently such as immersogeometric analysis [25] in the context of fluid-structure interaction and isogeometric B-Rep analysis for shell structures. [26] In the present work, we discuss the use of ℎ -basis functions with higher differentiability properties in the context of the FCM and numerical simulations on complex geometries. To this end, we introduce ℎ -basis functions which are based on the classical B-splines approach and present a new approach for the construction of 1 ℎ -basis functions with minimal local support. Both approaches allow for hanging nodes, whereas the new 1 approach also includes varying polynomial degrees.
An essential drawback of classical B-splines approaches is that the enforcement of higher differentiability properties requires supports that consist of large patches of mesh elements. Typically, this leads to a strong coupling of the degrees of freedom and, thus, to a dense structure of the stiffness matrix. To overcome this drawback and to construct basis functions with a more local support, approaches based on truncated B-splines have been developed. [27] In this paper, we introduce an alternative approach for 1 ℎ -basis functions, in which the construction is directly based on a Hermite interpolation and, thus, causes minimal supports of the basis functions. In particular, the supports have the same sizes as the supports of the usual 0 basis functions (for instance, four elements for nodal modes, two mesh elements for edge modes, and one element for inner modes in the case of uniform quadrilateral meshes). Therefore, assembling procedures as well as ℎ -adaptivity as used in the context of 0 basis functions can be applied straightforwardly. The construction approach of the 1 ℎ -basis functions is, in principle, transferable to (which is, however, beyond the scope of this paper and will be a future work). Currently, paraxial meshes are assumed for their construction, but this is not a restriction with respect to the FCM as the use of paraxial meshes is typical for this method.
The properties of the ℎ -basis functions based on classical B-splines as well as the new 1 ℎ -basis functions are studied in several numerical experiments, in which these basis functions are also applied to the FCM. Here, we consider a linear elastic problem with some singularities which is discretized with adaptive refinements. Using the 1 ℎ -basis functions with minimal support ℎ -adaptive refinements can be carried out in the usual elementwise way. The underlying idea of adaptivity with the ℎ -basis functions based on B-splines is to consider a sequence of refined hierarchical levels, each defining a set of basis functions. A suitable subset of all these functions is then chosen according to some refinement criteria. Since the ℎ -basis functions based on B-splines enable adaptivity with respect to local differentiability properties, we also study this type of adaptivity in some numerical experiments, where we observe a significant improvement of the convergence when low differentiability of the basis functions is chosen near to singularities and high differentiability away from them. For both ℎ -basis functions we obtain very similar convergence results. However, we observe essentially larger condition numbers in the case of the 1 ℎ -basis functions with minimal support. Finally, we also discuss the application of the ℎ -basis functions based on B-splines in the context of more involved nonlinear material models. Here, we consider hyperelasticity and elastoplasticity with finite strains. We investigate the usage of reduced differentiability close to the plastic front in order to capture the lower regularities in the stress field while maintaining maximum smoothness in regions far from the plastic zones.
The structure of the paper is as follows: Sections 2.1 and 2.2 introduce the ℎ -basis functions based on B-splines as commonly presented in the literature together with some fundamental concepts. Section 2.3 discusses the construction of 1 ℎ -basis functions with minimal support. Finally, Section 4 presents numerical experiments in the context of linear elasticity and nonlinear material models.

BASES WITH 1 AND HIGHER CONTINUITY
In the following we will introduce different basis functions with higher continuity. Both communities, IGA and ℎ -fem use the term continuity and differentiability. At times different interpretations are associated to this term. Thus, in the following, we will use the expressions continuity and differentiability not as given strictly by their mathematical definition but rather as synonyms describing a degree of smoothness of a given function.

B-splines
We herein briefly introduce the basic definitions and notations about B-splines. For further details, readers are referred to Piegl and Tiller [28] and Cottrell et al., [29] and references therein. A B-spline basis function of degree is generated starting from a nondecreasing sequence of real numbers referred to as knot vector where is the number of basis functions (equal to the number of the associated control points). A univariate B-spline basis function , ( ) can then be constructed using the following Cox-de Boor recursion formula: starting from = 0, where 0 otherwise the basis functions for > 0 are obtained from where the convention 0∕0 = 0 is assumed. Given the multiplicity of a knot, the smoothness of the B-spline basis is − at that location, while it is ∞ everywhere else. In so-called open knot vectors, the first and the last knots have multiplicity = + 1 and the basis is interpolatory at the ends of its domain. A B-spline curve can then be constructed as the linear combination of the basis functions where the coefficients P ∈ R of the linear combination are the so-called control points, being the dimension of the physical space. Multivariate B-splines are generated through the tensor product of univariate B-splines. If denotes the dimension of the parametric space, univariate knot vectors are needed: is the polynomial degree in the parametric direction , and is the associated number of basis functions. Denoting the univariate basis functions in each parametric direction by , , the multivariate basis functions i,p ( ) are obtained as } denotes the position in the tensor product structure, p = { 1 , … , } indicates the polynomial degrees, and = { 1 , … , } is the vector of the parametric coordinates in each parametric direction . B-spline surfaces and solids are obtained, for = 2 and = 3, respectively, from a linear combination of multivariate B-spline basis functions and control points as follows where the summation is extended to all combinations of the multi-index i.
The hierarchical B-spline basis  are defined as The domains Ω are chosen by the user to define the refinement areas. See Figure 1A for an example of hierarchical basis refined towards the right part of the domain. The colored functions compose the hierarchical B-spline basis, while the light gray functions are not active. For further details, we refer to Giannelli et al. [31,32] and D'Angella et al. [33] , and the references therein. Note that the nestedness of 0 ⊂ 1 ⊂ · · · ⊂ allows for different continuities across the hierarchy of functions. For example, see Figure 1B where the continuity is reduced only on the finest level. A further possibility (case 3, Figure 1C) is to use repeated knot insertion to generate 0 -continuous B-splines also on  1 and on  0 at the knot shared with the finer level close to the expected irregularity in the solution, but to keep maximum continuity away from it.

-basis functions with minimal support
When B-splines are used, the construction of the global basis functions as introduced in the previous sections is patch-wise: the domain is decomposed into the union of subdomains called patches, which are the images of multidimensional intervals under sufficiently smooth transformations, on which basis functions are defined. Thereby, the requested continuity is attained on each patch. Global -continuity is then obtained by enforcing the -continuity between the patches.
In contrast, the finite-element view on the construction of global basis functions is rather node-wise. By the nodes of a subdivision we denote its vertices, edges, faces, and elements. Typically, a set of global 0 basis functions is associated with

F I G U R E 2 Visualization of a patch consisting of elements of the refinement history. (A) A subdivision with multilevel hanging nodes. (B)
The patch ℘( ) ∶= { 1 , … , 4 } for the marked vertex each node, and the support of such basis functions is minimal in the sense that it spans exactly the elements in which the node lies. In the following, we call such basis functions nodal.
In this section, we explain how nodal 1 -ℎ basis functions for arbitrary-level hanging nodes can be defined. The construction is very similar to the well-known 0 multilevel space [34] and provides a more classical view on higher differentiability. In contrast to the patch-wise construction based on B-splines where the polynomial degree is typically fixed on each patch, varying polynomial degrees are supported in the same way as in standard ℎ finite elements.
Recall that, in subdivisions without hanging nodes, the construction of globally continuous basis functions is simple: here, local basis functions with matching values along common facets are combined to form continuous global basis functions for each node of the subdivision. This process is also called gluing together of local basis functions. [35] Typically, these local basis functions are found on the set of elements  ( ) containing the node , which is commonly called a patch for the node. The central idea of the multilevel space is to replicate this simple construction in the presence of hanging nodes by, again, determining a patch of elements ℘( ) sharing the node on which basis functions can be combined to form a globally continuous function. However, the elements of ℘( ) are not required to be contained in the subdivision  . Instead, each of these patch elements is formed by taking the union of certain elements of  .
We remark that global basis functions are only constructed for nonhanging nodes, which are called the regular nodes of the subdivision. Therefore, a method to decide whether a given node in  is regular or hanging needs to be available. For arbitrary-level hanging nodes, this decision can become computationally expensive. We refer to Byfut and Schröder [36] for possible ways that overcome this issue.
For determining the patch ℘( ) of a regular node , we assume that  =  is part of a sequence of subdivisions  1 , … ,  resulting from refinement out of an initial subdivision  1 of Ω. Moreover, we assume that the patch ℘( ) for each regular node of  can be retrieved from the refinement history  of  , that is, the union  = ⋃ =1  . Since ℘( ) is designed as a replacement of the usual patch  ( ), it has to fulfill the following properties: contains an Ω-neighborhood of , that is, Ω contains the intersection of Ω with a ball around . 3. ℘( ) is a regular subdivision of Ω , that is, ℘( ) has no hanging nodes. Some comments on these properties are in order. The first property ensures that is indeed a node of each ∈ ℘( ), so that local basis functions defined on each such can be glued together as expected. By the second property, it is ensured that there is an element contained in ℘( ) in each direction around . Finally, the third property secures that the construction of global, linearly independent functions is simple: because ℘( ) is a regular subdivision of Ω , a finite-element space can be defined on ℘( ) easily, for example, by constructing the usual nodal global basis functions for on Ω . In order to express these functions as linear combinations of basis functions defined on elements in  , constrained approximation as described [37][38][39] can be employed. As these functions vanish on the boundary of Ω , they can be extended by zero onto Ω to form global basis functions. Patches ℘( ) fulfilling these properties can be easily chosen from  if certain restrictions are imposed on the refinement patterns: in any case, the symmetric refinement of a quad resp. hex element into 4 resp. 8 child elements is supported. [34] In Figure 2, a patch for a vertex is found among the elements of the refinement history. Note that the set  ( ) does not fulfill the third requirement for a patch.
While the construction of the global 0 functions is well known, the same concept can be extended to the 1 case. Here, one starts with nodal functions spanning the polynomial space ([−1, 1]) suitable for 1 -continuity: the functions need to satisfy four constraints (one constraint on the value and one on the value of the derivative at the two endpoints), which implies that ≥ 3. Assuming = 3, for each of the two endpoints and each of the two derivative levels = 0, 1, there is exactly one nodal function satisfying ( ) ( ) = 1 and the other three functionŝfulfill̂( ) ( ) = 0. The resulting functions of degree = 3 are displayed in Figure 3A. For higher degrees, ( − 3) local inner functions of degree = 4, … , are introduced whose value The set −3∕2 , = 4, … , 7 is visualized in Figure 3B.
The tensor product of the presented functions yields a basis of 1 , 2 ([−1, 1] 2 ). A typical difference to the tensor product functions in the 0 case is the increased number of functions per boundary node: assuming 1 = 2 = , there are four local functions per vertex, 2( − 3) local functions per edge, and ( − 3) 2 inner local functions per element. For dimension = 3, the distribution of the ( 1 + 1)( 2 + 1)( 3 + 1) local functions to the vertices, edges, faces, and the element of [−1, 1] 3 is as follows: While the values of the functions remain unchanged after transformation, the values of the derivatives do not. Therefore, one needs to take into account the derivative D of the element mapping. In order not to be forced to solve an interpolation problem explicitly, we assume that D is a scaling and, thus, each element is a multidimensional interval. We point out that this is the same reasoning as applied to classical 1 finite elements, such as the Bogner-Fox-Schmit element, [40] where only elements of rectangular shape are permitted. Another reason for the restriction to intervals is the computational efficiency: for more general geometries of the patch, one could still solve an interpolation problem to obtain linearly independent, global 1 basis functions as linear combinations of local basis functions of patch elements. However, the geometry then dictates that these linear combinations have a large number of nonzero coefficients, and the computation of these coefficients requires the solution of local linear systems which leads to additional numerical errors. In contrast, if all elements are multidimensional intervals, each global basis function corresponds to exactly one local basis function per element and, thus, the resulting linear combination consists of only one nonzero coefficient.
We remark that, in the context of fictitious domain methods such as the FCM, using only rectangles or rectangular cuboids as elements is no restriction since the geometry of the domain is resolved by the use of appropriate quadrature rules.

FEATURES AND DIFFERENCES
In this section we compare the different bases introduced previously. We commence with a boundary conforming setting in Section 3.1 before inspecting their features in the context of FCM in Section 3.2.

Poisson's problem on a square
With the purpose of showing the difference on a simple example, we consider the Poisson's problem on a square domain Ω visualized in Figure 4 along with the governing equations. The data , ℎ, of the problem are obtained from the manufactured solution ( , ) = 6 6 . The problem is solved numerically on a Cartesian mesh of 2 elements per side, for = 0, … , 7, and different polynomial degrees = 3, … , 5. The continuity is fixed to 1 for both bases. Figure 5A shows the error in energy ( , ) − ( ℎ , ℎ ) of the numerical solution ℎ , where Figure 5B plots the condition number of the stiffness matrix. While the bases yield virtually the same accuracy in the solution, the conditioning differs considerably. In particular, the condition numbers of the nodal 1 basis increase with 2 , where is the number of degrees of freedom. By contrast, the condition of the stiffness matrix resulting from applying the IGA basis grows with 1 .

Cook's membrane
In this section, we compare the IGA basis with the nodal 1 basis solving the well-known Cook's membrane problem in the linear elastic case using the FCM. The geometry of the domain is displayed in Figure 6A. We choose Young's modulus as = 210 000 MPa and Poisson's ratio as = 0.3. The loading is given by a constant traction of 0 = 1.0 MPa on the upper right edge of the domain. The reference solution is provided by a standard -FEM approximation based on a fixed mesh graded towards the four corners of the domain as visualized in Figure 6B. The reference solution yields a strain energy of 1.341698174 × 10 −2 N m.
In the FCM, the computational domain is embedded into a larger domainΩ which can be discretized easily into a subdivision  . Then, an approximation space ℎ on  can be constructed as usual by, for example, a finite-element or a multilevel space. To obtain a unique solution to the discrete problem, the space is reduced by considering only those basis functions that have support on Ω. While, in theory, the problem then admits a unique solution, the method might require additional stabilization if some basis functions have very little support in Ω. To this end, it is common to introduce an indicator function ∶Ω → R Hereby, the contributions of basis functions with little support on Ω are increased. For the comparison of the 1 -continuous bases in the FCM regime, we embed the domain into the rectangle [0, 48] mm × [0, 60] mm which is initially discretized using a 2 × 2 mesh consisting of rectangles, see Figure 7. These are chosen such that the Dirichlet part of the boundary is an edge of the mesh, in order for the Dirichlet condition to be enforced in a strong manner. An initial polynomial degree distribution of = 3 is used, which is the minimum degree possible for the nodal 1 basis. We compare the behavior of the two 1 -continuous bases for a refinement strategy starting on this initial configuration. The mesh is refined locally towards the four corners six times by equally spaced bisections into four subelements. Then, the polynomial degree is increased globally from = 3 to = 9. The stabilization parameter is chosen to be = 10 −8 . The resulting energy errors and stiffness condition estimates are displayed in Figure 8A,B, respectively.
The energy error is depicted in Figure 8A while the condition estimate is depicted in Figure 8B. Both exhibit the same behavior as in the nonembedded case, that is, both strategies deliver a similar energy error; however, the nodal 1 basis exhibits larger condition estimates. Clearly, it is important to reduce the condition number when using the 1 basis in particular in the FCM setting. We aim to address this issue by implementing an algebraic Gram-Schmidt procedure and an additive Schwarz preconditioner as applied by Jomoet al. [41] and de Prenter et al. [42] These methods have proven to be effective in particular for large-scale FCM computations.
In contrast to Poisson's example, Cook's membrane exhibits four singularities towards which the mesh is refined. It is known from classical -FEM that a better accuracy per degree of freedom may be obtained by choosing a low polynomial degree towards

F I G U R E 9
Cook's membrane. Error in energy norm for = 3 and different continuities of hierarchical B-splines the singularity and successively increasing it away from it, see, for example, Byfut and Schröder [36] for a recent overview with extensions to unsymmetrical choices of polynomial degrees and refinements driven by error estimates. Although this strategy leads to an optimal choice of degrees of freedom, it is not yet readily available neither for the B-spline basis nor for the 1 continuous basis presented in Section 2.3.
Instead, we would like to point out that not only does the choice of ℎand -refinements influence the accuracy of the approximation but so does the continuity. All curves in Figure 9 depict the convergence of the error for a B-spline discretization of = 3 for a base mesh of 2 elements per parametric direction, = 1, 2, …, on which five local hierarchical refinements are carried out towards the corners. This is analogous to Figure 7, where now instead of a 2 × 2, we now use a 2 × 2 base level.
The curve labeled " 2 everywhere" in Figure 9 is obtained by the classical hierarchic overlays all of which exhibit a 2 -continuity at the knots and ∞ everywhere else. For the purpose of clarity, this choice is illustrated in a one-dimensional setting in Figure 1A for = 2. In addition to this classic, hierarchical discretization, all possible combinations of introducing a 0 continuity at the levels of the hierarchic refinements are investigated. To this end, first, only the continuity of the knots of the finest level is reduced to 0 . The corresponding curve is labeled " 0 only on finest level" in Figure 9. In a one dimensional setting, this choice corresponds to the discretization depicted in Figure 1B. The introduction of a 0 -continuity only on the finest overlay leads to an error which is initially higher per degree of freedom, but the asymptotic performance is the same as if the continuity were not reduced at all.
Next, the continuity is reduced to 0 for all levels, including the base level. The corresponding curve is presented in Figure 9, labeled " 0 everywhere" and marked by the brown circles. This type of reduction of continuity drastically improves the performance of the discretization by about 1.5 orders of magnitude. A further possible choice is to reduce the continuity on all overlays, that is, all levels of enrichments but not on the base level. The corresponding curve is depicted in Figure 9, labeled " 0 on refined levels," and marked by the black stars. This choice exhibits practically the same behavior as the previous choice of reducing the continuity on all levels, including the base level.
Interestingly, the best approximation w.r.t. the necessary number of degrees of freedom is obtained by a discretization which is also nonstandard: 0 on all overlays while the maximum 2 continuity is maintained on most of the base level except along the line defined by the knots which are closest to the singularity. It is only at this line that the continuity is also reduced to 0 at the base level. The results obtained by this choice of basis are marked by the blue diamonds in Figure 9 and is additionally labeled "case 3." This case 3 is a combination of a higher order ℎ-refinement with lower continuity closer to the singularity but a higher continuity elsewhere. For clarity, this case 3 is illustrated in Figure 1C for = 2 in a one dimensional setting. This kind of discretization delivers results which are approximately two orders of magnitude more accurate than the classic choice of maintaining 2 (case 1 in Figure 1A).
This indicates that a proper choice of continuity on the base level combined with 0 -refinement of higher order leads to a much more efficient discretization than the classic choice of a refinement which maintains maximum continuity everywhere.
It is worth pointing out that, at least in this example, the standard choice of keeping the highest continuity everywhere delivers the worst possible discretization among the investigated choices.

NONLINEAR NUMERICAL EXAMPLES
In this section, we utilize the B-spline basis for finite strain elasto-plasticity. After the introduction of the basics of finite 2 plasticity in Section 4.1 in the context of the FCM, we investigate the performance of different refinements on a benchmark example for hyperelasticity and finite strain elasto-plasticity in Section 4.1.1. In Section 4.1.2 we present the computation of a representative body produced by selective beam melting. The body has a complex shape and is computed under both compression and tension.

The FCM for finite 2 elastoplasticity
In the following, we will briefly describe the governing equations of the underlying material model which is based on the 2 flow theory of elastoplasticity for finite strains accounting for nonlinear isotropic hardening. Please refer to Wriggers and Hudobivnik, [43] Korelc and Stupkiewicz, [44] Simo and Miehe, [45] and Igelbüscher et al. [46] for a more detailed overview of the theory.
For the elastic part of the deformation, we assume an isotropic compressible neo-Hookean material behavior using the following strain energy function Ψ = 4 In (1), and are the Lamé constants and 1 = tr and 3 = det (2) denote the first and third invariant of the elastic left Cauchy-Green tensor . Based on the strain energy function the Cauchy and the Kirchhoff stress tensors can be computed: The von Mises yield criterion used to restrict the elastic region of the material model reads In Equation (4) defines the deviatoric part of the Kirchhoff stress tensor and (̄) is a function describing the hardening. Thereby,̄denotes the hardening variable that is often referred to as the equivalent plastic strain. For this material model (̄) is composed of a linear and an exponential function Here, 0 denotes the initial yield stress, the linear hardening parameter, ∞ the saturation stress, and the hardening exponent. According to Simo and Miehe, [45] the associative plastic flow rule and the evolution of the hardening variable read where  denotes the Lie derivative of the elastic left Cauchy-Green tensor and ≥ 0 is the nonnegative plastic multiplier.
T A B L E 1 Elastoplastic material properties [47] Lamé's first parameter ( ) 110,743 MPa Since the FCM can suffer from ill-conditioning, due to badly broken finite cells, we follow the same approach presented by Taghipour et al. [11] and Hubrich and Düster [47] that has been applied for stabilization of the FCM for problems in finite strain elastoplasticity. In doing so, we distribute additional integration points within the fictitious domain of broken finite cells. Then, for the material behavior of those stabilization points, we assign the same material model used for points within the physical domain. However, the Young's modulus of the stabilization points is scaled by a factor of = 10 − . Moreover, we set the initial yield stress to infinity 0 = ∞, thus, the material behavior of the stabilization points is only described by the hyperelastic part of the material model.

Plate with a cylindrical hole
In this example, we consider a plate with a cylindrical hole. The motivation is to compare different basis functions, however a comparison with the basis introduced in Section 2.3 is not yet available. Instead, we compare the B-spline basis introduced in Section 2.1 to the somewhat more classic integrated Legendre basis with 0 -continuity. For the material response, we utilize the elastoplastic constitutive model described in Section 4.1. For the hyperelastic simulations, the initial yield stress is set to infinity which leads to the neo-Hookean model. All elastoplastic material parameters are listed in Table 1.
The geometry of the model as well as the boundary conditions (BCs) can be seen in Figure 10A. Symmetry BCs are applied, by fixing the bottom side in -direction, the right-hand side in -direction, and the back side in -direction. On the top face, a displacement of 12 mm is applied which leads to a total deformation of 12%.

Hyperelasticity
In the first part, we investigate the different bases using the neo-Hookean model. In doing so, we first generate a reference solution using conforming high-order finite elements. The curved boundary is accounted for by applying the blending functions method. Thereby, the geometry is discretized with 3420 hexahedral elements as can be seen in Figure 10B. For the FCM, a plate of 100 × 100 × 10 mm 3 size is created which is cut by a cylinder that is defined using the following level set function The plate is then discretized with 78 cells as can be seen in Figure 11A. For capturing the physical domain, an adaptive octree is used for integration with a tree depth of = 3. Figure 12 shows the load-displacement curves for = 3, using B-spline and the Legendre basis. The results are in excellent agreement with the reference solution.

Elastoplasticity
In the second part, we investigate the elastoplastic behavior employing the same meshes as before. By applying a displacement of 12 mm at the top face, the load-displacement curves are plotted in Figure 13. We observe a slightly stiffer behavior for the basis with higher continuities. Subsequently, we apply one ℎ-refinement for cells at the base using the B-spline basis to generate a high resolution at the plastic area as depicted in Figure 11B. The load-displacement curves are plotted in Figure 14A for = 2 to 4. The B-spline discretization converges well to the reference solution using 0 for the refined elements and −1 -continuity elsewhere. Figure 14B depicts for = 2 the load-displacement curve for different choices of the continuities in the construction of the overlays (see Figure 1 for a one dimensional representation of those different constructions).
We also plot the convergence of the load versus the number of degrees of freedom at the last load step as depicted in Figure 15. Using B-splines with maximum continuity and no refinements yields the least number of DOFs, nevertheless the force is too high. On the other hand, utilizing a B-spline basis, with 0 -continuity for cells in the plastic zone and −1 -continuity elsewhere, yields the same load as the Legendre basis with less DOFs for the B-spline discretization. Finally, using B-splines with one ℎ-refinement and 0 for the refined elements gives the best investigated load convergence.

Swiss cheese domain
In the last example, we study the "swiss cheese domain," a more complex geometry from additive manufacturing using the multilevel Bezier extraction.
The domain consists of eight cheese blocks of different geometries with a size of 3 × 3 × 3 mm 3 for one block, and a total size of 6 × 6 × 6 mm 3 . The whole model is depicted in Figure 16.
The geometry of a single cheese block is defined using the level set function [47] ( where the center coordinates , , and as well as the inner radius , the outer radius , and the parameter d of each of the cheese blocks are listed in Table 2. Using the FCM, a Cartesian grid of 20 × 20 × 20 cells is created with BCs as described in Figure 16B. Cells which are not intersected by the physical domain are removed. We apply one ℎ-refinement towards the plastic zone which leads to a total number of 8552 cells. The bottom face is fixed in -direction and a displacement at the top face in tension or compression is applied. The load-displacement curve is plotted for a displacement of z = 0.24 mm (tension) and z = −0.24 mm (compression) for = 2 and = 10 −5 as can be seen in Figure 17, obtaining a maximum load of 3.7 kN for tension, and −5.2 kN in compression. We again observe a difference in the load-displacement behavior for the different continuities of the discretization given by case 1 and case 3 (see also Figure 1 for a one dimensional visualization of the type of discretization). Figures 18 and 19 show the contour plots for the von Mises stress vM as well as the equivalent plastic strain̄for the final load step in tension and compression. Both load cases feature a strong localization of strain at the fillets where necking occurs. In tension, higher plastic values are located at the top four blocks whereby in compression, plastic zones can be seen on both the top and the bottom blocks.
In this example, we also noted effects on stability of the discretization. Different continuities lead to different computability in the sense that some combinations enabled larger displacements than others. These effects still need further investigation. Furthermore, the computational results need to be validated.

CONCLUSIONS
This contribution investigates several basis functions for the FCM in linear and nonlinear settings whereby we focus on the class of discretizations providing higher-order continuity. The hierarchic B-spline basis is its most common representative, having a considerably large support. As an alternative we present a newly developed 1 ℎ -basis with minimal support. For the investigated, linear examples both bases lead to the same convergence in the energy norm but the basis delivering the small support constructs a stiffness matrix with a higher condition number. The new 1 ℎ -basis still needs to be tested for nonlinear examples.
Classically, the hierarchic B-spline refinements keep maximum continuity in all of the patch. We investigate hierarchical refinements whose continuities may be reduced to 0 , for example, only on the finest refinement level, only on the base level, or all levels of refinement. We observe that it is possible to construct refinements which deliver an accuracy of two orders of magnitude higher than the classical refinement. Our observation is that a lower continuity towards and in the zone of an irregular solution and a higher continuity far away from it is a good choice. However, this new aspect also needs further investigation, especially because no mathematical proof is available which suggests guidance on how to choose continuity in hierarchic overlays.
In the last part of the article we attempt a comparison of the B-spline basis to the standard integrated Legendre basis used in -FEM. While direct performance comparisons are difficult, a more practical observation is that clearly both bases converge towards the same solution, also for nonlinear examples of finite-strain. In the strongly nonlinear setting where, additionally, elastoplasticity was added to the finite-strain kinematics, local refinements using reduced continuity towards the plastic zone helped to capture local behavior.