An immersed-isogeometric model: Application to linear elasticity and implementation with THBox-splines

Journal of Computational and Applied Mathematics 349 (2019) 410–423

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

An immersed-isogeometric model: Application to linear elasticity and implementation with THBox-splines ∗

Carlotta Giannelli a , Tadej Kanduč b , , Francesca Pelosi c , Hendrik Speleers c a

Dipartimento di Matematica e Informatica ‘‘U. Dini’’, Università degli Studi di Firenze, Italy Istituto Nazionale di Alta Matematica, Unità di Ricerca di Firenze c/o DiMaI ‘‘U. Dini’’, Università degli Studi di Firenze, Italy c Dipartimento di Matematica, Università degli Studi di Roma ‘‘Tor Vergata’’, Italy b

article

info

Article history: Received 16 March 2018 Received in revised form 13 September 2018 Keywords: Linear elasticity Immersed boundary method Isogeometric analysis Box splines Truncated hierarchical splines Local refinement Weak boundary conditions

a b s t r a c t We investigate the application of immersed boundary approaches in isogeometric analysis for the treatment of flexible domains by suitably incorporating trimming operations and geometry mappings. The considered immersed-isogeometric model is framed in the context of an automatic adaptive scheme to solve linear elasticity problems. The proposed method leads to a symmetric system of linear equations, and it is essentially free of user-defined penalty and stabilization parameters. Adaptivity is achieved by employing hierarchically nested spline spaces. In particular, we focus on truncated hierarchical box splines (THBox-splines) defined over regular triangulations. Several numerical examples demonstrate the optimal convergence of the adaptive scheme. © 2018 Elsevier B.V. All rights reserved.

1. Introduction An accurate representation of computational domains in numerical simulation is one of the key features in the design of isogeometric methods. The isogeometric approach extends finite element methods by considering spline constructions and related representations in line with computer-aided design systems [1]. The tensor-product structure of standard B-splines and their rational version, however, precludes a proper treatment of complex geometries, possibly trimmed, and lacks local refinement capability needed for the development of adaptive schemes. In order to achieve flexible isogeometric methods, we suitably combine immersed boundary methods with a geometry mapping and hierarchical box spline constructions. Tensor-product B-splines with equally spaced knots form a special case of the more general box spline model, the multivariate generalization of cardinal B-splines [2,3]. Box splines are nonnegative piecewise polynomial functions with compact support and enjoy convolution properties similar to (cardinal) B-splines. They also have a similar recurrence relation and are completely determined by certain direction vectors which, in turn, determine the properties of a single box spline and characterize the space generated by the translates of a fixed box spline. For example, focusing on the bivariate case, box splines on two-directional meshes are the standard tensor-product B-splines on uniform knot configurations with rectangular support, while box splines on three-directional meshes are characterized by hexagonal supports. The possibility of introducing additional directions in the mesh structure opens the path to more flexible geometry representations. Immersed boundary methods facilitate a proper treatment of non-trivial geometries by extending the mesh of the physical domain outside the domain boundary and usually require a weak imposition of the boundary conditions. The most common techniques to weakly enforce boundary conditions include the penalty method and Nitsche’s method; see, e.g., [4] ∗ Corresponding author. E-mail addresses: [email protected] (C. Giannelli), [email protected] (T. Kanduč), [email protected] (F. Pelosi), [email protected] (H. Speleers). https://doi.org/10.1016/j.cam.2018.09.027 0377-0427/© 2018 Elsevier B.V. All rights reserved.

C. Giannelli et al. / Journal of Computational and Applied Mathematics 349 (2019) 410–423

411

and references therein in the context of the finite cell method. From the extensive literature on the topic, we just mention here a few references connected with our setting [5–7]. Specifically, in this paper we follow the formulation presented in [5]. The weakly imposed Dirichlet boundary conditions do not require any tuning of user-defined penalty or stabilization parameters, as it is typically needed with the penalty and Nitsche’s method. Moreover, the system matrices are symmetric for symmetric problems. Isogeometric analysis on trimmed geometries poses crucial challenges and is an extremely active area of research; see, e.g., [8–10]. Particular attention is needed for numerical integration on trimmed geometries. When the immersed boundary approach is suitably combined with a geometry mapping, the issue of constructing robust quadrature rules on mesh elements cut by the boundary can be alleviated. The advantage of applying box spline constructions in connection with immersed boundary methods has a twofold motivation. First, it eliminates the need of identifying special constructions for the boundary functions [11], and second, it can profit from the regular structure of the set of box spline elements in the matrix assembly [12]. Adaptivity can be achieved through a local spline refinement using the hierarchical approach and by constructing truncated hierarchical bases, originally proposed in [13,14] and generalized to generating systems in [15]. The application of truncated hierarchical box splines (THBox-splines) in isogeometric analysis was recently considered in [16] for solving the advection–diffusion problem. In particular, the construction of a suitable adaptive domain strip along the domain boundary was proposed for the weak imposition of boundary conditions. In this paper, we continue this line of research and extend it in the following directions:

• by following the main idea of the immersed method in [5] we introduce an immersed boundary model for linear elasticity that is essentially free of user-defined penalty and stabilization parameters, and we illustrate the numerical performance of the model with THBox-splines; • we show that combining a weak imposition of boundary conditions with a geometry map can outperform a standard (mapping-free) immersed method, leading to better conditioned formulations; we refer to the numerical example in Section 4.3 for a comparison between the two approaches; • we extend the framework, previously derived in [11,16], by validating a fully automatic adaptive isogeometric method with THBox-splines, and we show its potential with respect to different domain configurations on three-directional meshes. The structure of the paper is as follows. We start by presenting our immersed model for linear elasticity in Section 2. Then, in Section 3 we review the (truncated) hierarchical box spline construction and describe the considered adaptive scheme. A selection of numerical examples is presented in Section 4, and Section 5 concludes the paper. 2. Immersed model for linear elasticity In this section we describe an immersed model for linear elasticity incorporating boundary conditions in a weak form. The immersed model is inspired by the one for Poisson and Stokes problems proposed in [5]. For notational convenience, we just focus on the 2D case, but it can be quite easily generalized to higher dimensions. 2.1. Problem statement ◦

◦

Let Ω ⊂ R2 be a domain with boundary ∂ Ω = ΓD ∪ ΓN and Γ D ∩ Γ N = ∅ (Dirichlet and Neumann parts of the boundary). We study the linear elastic behavior of the displacement field u : Ω → R2 , with u := [u1 , u2 ]T , described by

− div(S(u)) = f , in Ω , u = g , on ΓD , S(u) n = k , on ΓN , [ ] where S(u) := Sij (u) i,j=1,2 stands for the stress tensor and n := [n1 , n2 ]T the outward unit normal. In addition, f := [f1 , f2 ]T , [ ] g := [g1 , g2 ]T , and k := [k1 , k2 ]T are given vector functions. The strain tensor is denoted by ε (u) := εij (u) i,j=1,2 , and ( ) ∂ uj 1 ∂ ui εij (u) := + , i, j = 1, 2. 2 ∂ xj ∂ xi For a small deformation of an isotropic homogeneous material, Hooke’s law gives a relationship between stress and strain, S(u) = λ tr(ε (u)) I2 +2µ ε (u) = λ div(u) I2 +2µ ε (u), where tr(·) is the matrix trace and I2 the 2 × 2 identity matrix. The two constants λ and µ, called the Lamé parameters, are given by

λ :=

Eν (1 + ν )(1 − 2ν )

,

µ :=

E 2(1 + ν )

,

where E denotes the Young’s modulus and ν the Poisson’s ratio.

412

C. Giannelli et al. / Journal of Computational and Applied Mathematics 349 (2019) 410–423

Fig. 1. Three-directional mesh and subdivision of Ω with boundaries ΓN (dotted line) and ΓD (solid line).

Given a mesh that partitions the entire plane R2 , we define Ωh as the smallest collection of mesh cells such that Ω ⊆ Ωh . The boundary Γ := ∂ Ω does not need to coincide with the boundary of Ωh . More precisely, we set

• Ωh = Ω in ∪ ΩΓ , where ΩΓ contains at least all cells cut by or contiguous to Γ , and Ω in is formed by a set of mesh cells completely in the interior of Ω ; • ΩΓ = ΩΓD ∪ ΩΓN , where ΩΓD and ΩΓN are disjoint sets corresponding to ΓD and ΓN , respectively; • ΩΓinD := ΩΓD ∩ Ω . Note that the inner part of ΩΓ considered in the last item is defined only for the part of the boundary related to ΓD , since the imposition of the Neumann conditions does not need the definition of the boundary strip on ΓN ; see Section 2.2. In Fig. 1 an example of the domain subdivisions is shown. We can then reformulate the linear elasticity problem by introducing an additional stress unknown σ :

− div(S(u)) = f ,

in Ω ,

(1)

σ = λ div(u) I2 +2µ ε(u), in ΩΓinD ,

(2)

u = g,

on ΓD ,

(3)

S(u) n = k ,

on ΓN .

(4)

2.2. Discrete formulation We are now ready to formulate an immersed discrete version of problem (1)–(4), whose exact solution is [u, σ]. We select two discretization spaces Vh ⊂ (H 1 (Ωh ))2 and Wh ⊂ (L2 (ΩΓD ))2×2 , defined over Ωh and ΩΓD , respectively, and we look for an approximate solution [uh , σ h ] ∈ [Vh , Wh ]. Let {·, ·}ω denote the L2 inner product for matrix functions over ω, (·, ·)ω denote the L2 inner product for vector functions over ω, and ⟨·, ·⟩ω denote the L2 inner product for scalar functions over ω. The weak formulation of (1) using the test function vh ∈ Vh , after applying Green’s formula, leads to

λ⟨div(u), div(vh )⟩Ω + 2µ{ε (u), ε(vh )}Ω − (σ n, vh )ΓD = (f , vh )Ω + (k , vh )ΓN .

(5)

The formula (5) incorporates the natural imposition of the Neumann conditions (4) and the stress (2). Let τ h ∈ Wh . By weakly enforcing (2), tested against 1

η(λ + 2µ)

(λ div(vh ) I2 +2µ ε(vh ))

and

−

1

η(λ + 2µ)

τh,

respectively, we obtain

λ 2µ ⟨tr(σ ), div(vh )⟩Ω in + {σ, ε(vh )}Ω in ΓD ΓD η(λ + 2µ) η(λ + 2µ) 2λ 2µ = ⟨div(u), div(vh )⟩Ω in + {2µ ε(u), ε(vh )}Ω in , ΓD ΓD η η(λ + 2µ) 1 λ 2µ − {σ, τ h }Ω in = − ⟨div(u), tr(τ h )⟩Ω in − {ε (u), τ h }Ω in . ΓD ΓD ΓD η(λ + 2µ) η(λ + 2µ) η(λ + 2µ)

(6) (7)

Finally, by weakly enforcing the Dirichlet boundary condition (3), tested against −τ h n, we get

−(u, τ h n)ΓD = −(g , τ h n)ΓD .

(8)

C. Giannelli et al. / Journal of Computational and Applied Mathematics 349 (2019) 410–423

413

By adding up (5)–(6) and (7)–(8), and replacing the continuous solution [u, σ] with an approximate one [uh , σ h ] ∈ [Vh , Wh ], we arrive at the following discrete variational formulation: ⎧ λ⟨div(uh ), div(vh )⟩Ω + 2µ{ε (uh ), ε(vh )}Ω − (σ h n, vh )ΓD ⎪ ⎪ ⎪ ⎪ 1 2µ λ ⎪ ⎪ ⎪ tr(σ h ) − 2 div(uh ), div(vh )⟩Ω in + {σ h − 2µ ε(uh ), ε(vh )}Ω in + ⟨ ⎪ ⎪ ΓD ΓD η λ + 2 µ η ( λ + 2µ) ⎨ = (f , vh )Ω + (k , vh )ΓN , ∀vh ∈ Vh , (9) ⎪ ⎪ ⎪ λ 1 ⎪ ⎪ ⟨div(uh ), tr(τ h )⟩Ω in + {2µ ε(uh ) − σ h , τ h }Ω in − (uh , τ h n)ΓD ⎪ ⎪ ΓD ΓD η(λ + 2µ) η(λ + 2µ) ⎪ ⎪ ⎩ = −(g , τ h n)ΓD , ∀τ h ∈ Wh . The parameter η > 1 can be fixed a priori. For a stability analysis one could follow a similar line of arguments as in [5]. In particular, it can be shown that the model is not sensitive to the choice of η. The parameter can influence the condition number of the system matrix but the problem remains well-posed for any η > 1. For that reason, the model is essentially free of user-defined parameters. In the numerical experiments the value is fixed to η = 2 for all the tests. 2.3. Implementation By choosing bases for Vh and Wh , the discrete weak formulation (9) can be rewritten in a discrete algebraic form. Let nh and mh be the dimensions of Vh and Wh , respectively. Let U ∈ Rnh and Σ ∈ Rmh be the unknown coefficient vectors describing uh and σ h , respectively. Then, we obtain the linear system

[

K uu + K in uu K σ u + Gσ u

K uσ + G uσ K σσ

][ ] U

Σ

]

[

f + k uk , = uf gσg

(10)

where notations of the different blocks are described below: K uu U ↔ λ⟨div(uh ), div(uh )⟩Ω + 2µ{ε (uh ), ε (vh )}Ω , 2λ

G uσ Σ ↔ −(σ h n, vh )ΓD ,

4µ2

{ε (uh ), ε(vh )}Ω in , ΓD η(λ + 2µ) λ 2µ K uσ Σ ↔ ⟨tr(σ h ), div(vh )⟩Ω in + {σ h , ε(vh )}Ω in , ΓD ΓD η(λ + 2µ) η(λ + 2µ) λ 2µ ⟨div(uh ), tr(τ h )⟩Ω in + {ε (uh ), τ h }Ω in , K σ uU ↔ ΓD ΓD η(λ + 2µ) η(λ + 2µ) K in uu U ↔ −

K σσ Σ ↔ −

η

⟨div(uh ), div(uh )⟩Ω in − ΓD

1

η(λ + 2µ)

(σ h , τ h )Ω in , ΓD

G σ u U ↔ −(uh , τ h n)ΓD , f uf

↔ (f , vh )Ω ,

gσg

↔ −(g , τ h n)ΓD ,

k uk

↔ (k , vh )ΓN .

The linear system (10) is clearly symmetric. When the constants λ and µ are large and of the same order, the condition number of the matrix in (10) could be quite large. This could be improved by considering a properly scaled problem. We propose 1

1

− div(S(u)) = f , in Ω , κ κ κσ = λ div(u) I2 +2µ ε(u), in ΩΓinD , u = g,

on ΓD ,

S(u) n = k ,

on ΓN ,

where κ is the introduced scaling parameter, for example κ ≈ λ + 2µ. Then, the linear system (10) changes into

[1 κ

K uu + κ1 K in uu K σ u + Gσ u

K uσ + G uσ κK σσ

][ ] U

Σ

[1 =

f κ uf

+ κ1 k uk

gσg

]

.

3. THBox-splines To test our immersed model for linear elasticity, we focus on (truncated) hierarchical box splines as basis functions for the spaces Vh and Wh . In this section we first summarize their construction and then describe our fully automatic adaptive scheme.

414

C. Giannelli et al. / Journal of Computational and Applied Mathematics 349 (2019) 410–423

3.1. Hierarchical box spline spaces A d-variate box spline MΞ is determined by the matrix Ξ := [v1 · · · vn ], whose columns are a collection of (possibly repeated) direction vectors vk ∈ Zd \ {0}, k = 1, . . . , n. For simplicity, we assume n ≥ d and the submatrix Ξ d := [v1 · · · vd ] is non-singular. Furthermore, following the notation in [2], the matrix Ξ \ v is obtained from Ξ by omitting the vector v once, and the set Ξ d [0, 1)n contains all points in Rd obtained by multiplying Ξ d with all the points in [0, 1)n . Definition 3.1. The box spline MΞ is defined by successive convolutions as follows: 1

∫

MΞ \vn (x − t vn ) dt ,

MΞ (x) :=

n > d,

(11)

0

starting from 1/|det(Ξ d )|, 0,

{ MΞ d (x) :=

if x ∈ Ξ d [0, 1)d , otherwise.

(12)

The construction of the box spline MΞ does not depend on the ordering of the direction vectors vi , i = 1, . . . , n. Any box spline is a piecewise polynomial that is non-negative, locally supported, and with a certain global smoothness depending on the direction vectors. Box splines and their derivatives can be evaluated through simple recurrence relations, and elegant expressions are known for their inner products [12]. There is a well-established theory for spaces spanned by integer translates of box splines, namely MΞ (· − i), i ∈ Zd .

{

}

(13)

We refer to [2] for further information on box splines. Example 3.2. A space of C 2 cubic splines defined over two-directional meshes (d = 2) and a space of C 2 quartic splines defined over three-directional meshes (d = 2) are obtained by considering integer translates of the box spline MΞ generated by the matrices

[

1 0

Ξ =

0 1

1 0

0 1

1 0

0 1

1 0

0 1

]

[ and

Ξ =

1 0

0 1

1 1

1 0

0 1

1 1

]

,

(14)

respectively. More details on the computation of C 2 quartic box splines can be found in [11]. In the following we recall the definition of hierarchical box spline spaces from [14,16]. We consider a nested sequence of ˆ 0 in Rd . Each of these spaces Vˆ ℓ is spanned ˆ 0 ⊂ Vˆ 1 ⊂ Vˆ 2 ⊂ · · · restricted to a given domain Ω d-variate box spline spaces V by scaled integer translates of a certain box spline MΞ using a proper scaling factor hℓ according to the level ℓ (we assume hℓ > hℓ+1 ), i.e., Bˆℓ :=

{ MΞ

( · hℓ

} ) − i , i ∈ Zd .

(15)

For simplicity of notation, we denote an element of Bˆℓ by βˆ ℓ . We assume that the chosen sequence of Bˆℓ , ℓ ≥ 0, possesses the following properties:

• all elements in Bˆℓ are linearly independent; • two-scale relations between Bˆℓ and Bˆℓ+1 have only non-negative coefficients, so ∑ βˆ ℓ = cβˆ ℓ+1 (βˆ ℓ ) βˆ ℓ+1 , cβˆ ℓ+1 (βˆ ℓ ) ≥ 0.

(16)

ˆ 1 βˆ ℓ+1 ∈Bℓ+

The above properties are satisfied for several popular families of box splines defined over nested sequences. In particular, they are fulfilled for box splines on two- or three-directional meshes constructed by dyadic refinement [2]. ˆ0 ⊇ Ω ˆ1 ⊇ Ω ˆ 2 ⊇ · · · be a nested sequence of subsets, where each Ω ˆ ℓ is chosen to be aligned with the box Let Ω ˆ N = ∅ for some N ∈ N, and we denote the corresponding finite sequence by spline mesh of level ℓ. We assume that Ω ˆ := {Ω ˆ 0, Ω ˆ 1, . . . , Ω ˆ N −1 }. This sequence represents the regions to be refined at different levels of resolution. Ω

ˆ ) of hierarchical box splines (HBox-splines) related to the hierarchy Ω ˆ is defined as ˆ (Ω Definition 3.3. The set H

{

}

ˆ ) := βˆ ℓ ∈ Bˆℓ : suppΩ (βˆ ℓ ) ⊆ Ω ˆ ℓ ∧ suppΩ (βˆ ℓ ) ̸⊆ Ω ˆ ℓ+1 , ℓ = 0, . . . , N − 1 , ˆ (Ω H ˆ0

ˆ0

ˆ. ˆ 0 , i.e., the largest set in Ω where suppΩ (βˆ ℓ ) stands for the intersection of the support of βˆ ℓ with Ω ˆ0

C. Giannelli et al. / Journal of Computational and Applied Mathematics 349 (2019) 410–423

415

Fig. 2. A sequence of locally refined two-directional (top) and three-directional (bottom) hierarchical meshes. The anchors of C 2 (T)HBox-splines are indicated with different markers for the different levels (a square for level 1, a circle for level 2, and a star for level 3).

Fig. 3. HBox-spline basis (left) and THBox-spline basis (right) defined on the three-directional hierarchical mesh in Fig. 2(c).

ˆ ℓ ⊂ Vˆ ℓ+1 expressed with respect By exploiting the two-scale relations (16), we define the truncation of a function sˆ ∈ V ˆ ℓ+1 , to Bˆℓ+1 as the linear combination of only those basis functions in Bˆℓ+1 whose support is not included in Ω truncℓ+1 sˆ :=

cβˆ ℓ+1 (sˆ) βˆ ℓ+1 .

∑ ˆ0 Ω

βˆ ℓ+1 ∈Bˆ ℓ+1 , supp

(17)

ˆ ℓ+1 (βˆ ℓ+1 )̸ ⊆Ω

By iterating the truncation over the spline hierarchy, we define the truncated basis for hierarchical box splines as follows.

ˆ ) of truncated hierarchical box splines (THBox-splines) related to the hierarchy Ω ˆ is defined as Definition 3.4. The set Tˆ (Ω

{

}

ˆ ) := Truncℓ+1 (βˆ ℓ ) : βˆ ℓ ∈ Bˆℓ ∩ Hˆ (Ω ˆ ), ℓ = 0, . . . , N − 1 , Tˆ (Ω where Truncℓ+1 (βˆ ℓ ) := truncN −1 (truncN −2 (· · · (truncℓ+1 (βˆ ℓ )) · · ·)).

ˆ ) and Tˆ (Ω ˆ ) span the same hierarchical space. Furthermore, the THBox-splines are linearly independent and ˆ (Ω The sets H form a convex partition of unity; see [14]. They also allow for a very simple construction of quasi-interpolation schemes [17]. Example 3.5. Fig. 2 shows a sequence of locally refined two- and three-directional hierarchical meshes, together with the anchors of the corresponding C 2 box splines (defined by (15) with Ξ in (14)) involved in the hierarchical construction. The anchors are located in the center of the box spline supports. The HBox-splines and THBox-splines defined on the final three-directional hierarchical mesh are depicted in Fig. 3.

416

C. Giannelli et al. / Journal of Computational and Applied Mathematics 349 (2019) 410–423

3.2. Immersing the boundary and the boundary strip The discretization method elaborated in Section 2 fits in the class of immersed methods. One of the main difficulties of immersed methods is an accurate numerical integration over possibly very small mesh cells in ΩΓinD cut by the boundary. This is a subtle source of ill-conditioning (see, e.g., Section 4.3). In order to avoid such problems, we propose to combine the method with a geometry map F describing the physical domain Ω , according to the isogeometric philosophy. This allows us to take Ωh = Ω and ΩΓD = ΩΓinD . In this way we combine the benefits offered by immersed boundary and isogeometric methods in the hierarchical box spline context. We now specify the approximation spaces Vh and Wh in terms of THBox-splines (of course, the use of HBox-splines is ˆ 0 , suppose that the domain Ω can be described by a sufficiently smooth geometry analogous). Given a parametric domain Ω ˆ 0 → Ω , which is invertible and satisfies F (∂ Ω ˆ 0 ) = ∂ Ω . Then, using the THBox-spline basis related to the function F : Ω 0 N − 1 ˆ ˆ ,...,Ω ˆ hierarchy Ω := {Ω }, we specify

{

}

ˆ) , Vh := β := βˆ ◦ F −1 : βˆ ∈ Tˆ (Ω

Vh := span(Vh ).

(18)

ˆ Γ := {Ω ˆ0∩Ω ˆ ΓD , . . . , Ω ˆ N −1 ∩ Ω ˆ ΓD } where Ω ˆ ΓD := F −1 (ΩΓD ), we Similarly, using the hierarchy of restricted subsets Ω D specify

{

}

ˆΓ ) , Wh := β := βˆ ◦ F −1 : βˆ ∈ Tˆ (Ω D

Wh := span(Wh )2 .

(19)

The hierarchical construction in Definition 3.4 guarantees that both sets Vh and Wh form a basis, and that span(Vh ) is a subspace of span(Wh ) if we restrict the domains to ΩΓD . We refer to [16, Section 3.1] for a more detailed explanation of this phenomenon. ˆ Γ , and in particular Ω ˆ ΓD . For this, we can use one of the It still remains to explain how to choose the boundary strip Ω strategies described in [16, Section 3.2]. In particular, we opt for the so-called THBox-spline boundary strip because it results in the smallest number of degrees of freedom.

ˆ 0 contained in the union of the supports Definition 3.6. Let TˆΓℓ be a strip of thickness hℓ along the part of the boundary ∂ Ω ℓ+1 ˆ ℓ ℓ ℓ ˆ ˆ ˆ of Trunc (β ) where β ∈ H Bˆ . The THBox-spline boundary strip is defined as the union of the substrips TˆΓℓ , for ⋃N(−Ω1 ) ∩ ℓ ˆ ˆ ℓ = 0, . . . , N − 1, as ΩΓ := ℓ=0 TΓ . Example 3.7. The boundary strips corresponding to the sequence of hierarchical meshes in Fig. 2 (Example 3.5) are depicted in Fig. 4.

Fig. 4. The THBox boundary strips corresponding to the hierarchical meshes in Fig. 2.

C. Giannelli et al. / Journal of Computational and Applied Mathematics 349 (2019) 410–423

417

Fig. 5. Refinement patch for two types of marked cells on three-directional meshes.

3.3. Adaptive scheme We now describe an adaptive scheme that iteratively computes a numerical solution for our model problem in terms of THBox-splines defined on locally refined three-directional meshes. At each iteration step, the construction of the next hierarchical mesh requires suitable (local) error estimates, a marking strategy, and a refinement procedure that enlarges the set of initially marked cells. In our case, we consider a standard residual-based error estimator; see, e.g., [18]. The a posteriori error analysis of adaptive isogeometric methods with hierarchical spline constructions was recently discussed in [19] for a simple model problem. As marking procedure, we use the quantile marking strategy. The cells are first arranged in an order {π1 , π2 , . . . πK }, such that ϵπ21 ≥ ϵπ22 ≥ · · · ≥ ϵπ2K , where ϵπ2 is the local error estimator on a cell π . Subsequently, a marking threshold parameter ρ ∈ [0, 1] is chosen to specify the relative amount of cells to be refined. The first k cells with the largest estimated error are marked, M := {π1 , π2 , . . . , πk }, with k = ⌈ρ K ⌉. Finally, we define the refinement patch associated with every marked cell π ∈ M as the collection of its neighboring cells, as shown in Fig. 5. Every cell in the refinement patch of any marked cell is also refined. The introduced refinement patches prevent a possible occurrence of mesh refinements of isolated cells that do not contribute to any refinement of the spline space. 4. Numerical examples In this section we present some numerical experiments illustrating the performance of C 2 quartic THBox splines on three-directional meshes, introduced in Section 3.1, in the context of linear elasticity. In the first example we validate the linear elasticity model described in Section 2.2 by considering a simple rectangular domain. In the second example the optimal convergence rate of the method is demonstrated in a problem with an analytical exact solution by applying a uniform refinement strategy for the approximate solution. In the third example we compare two techniques to construct the computational domain. The first is the mapping-free approach with a cell trimming step, while in the second case we construct a suitable geometry mapping and avoid any cell trimming. In the last two numerical examples we test the convergence rate of the approximate solutions when applying the adaptive refinement strategy from Section 3.3. The examples comprise of problems with singular exact solutions on an L-shape domain and on a more geometrically complex non-simply connected triangular support structure. In all the experiments we consider homogeneous conditions in the interior of the domain Ω , i.e., f ≡ [0, 0]T . The material elasticity constants are set to E = 2.1 · 106 , ν = 0.28, and the stabilization parameters to η = 2, κ = 106 . In the convergence plots we always show the L2 norm of the (displacement/stress) error versus degrees of freedom (DOF). When needed, for the cut triangles at the domain boundary we adopt the quadrature scheme presented in [11]. More precisely, the area of interest is locally split into a few smaller triangles, so that each triangle has at most one curved edge. Then, a standard triangular quadrature rule with a local remapping is applied. 4.1. Bending of a cantilever loaded at the end In the first benchmark example we verify the correctness of the derived model. We consider a cantilever of length L and thickness D that is fixed at the left edge; see Fig. 6(a). The vertical force P applied at the right edge causes parabolic traction and the cantilever to bend. The top and bottom edges are free from load. The analytical solution can be found in [20]; the exact displacement is given by u1 (x, y) =

2Py

(

E0 D3

u2 (x, y) = −

2P E0 D 3

(

(6L − 3x)x + (2 + ν0 ) y2 −

(

3ν0 y2 (L − x) +

D2 4

where E0 := E /(1 − ν 2 ),

ν0 := ν/(1 − ν ).

D2 4

))

,

(4 + 5ν0 )x + (3L − x)x2

)

,

418

C. Giannelli et al. / Journal of Computational and Applied Mathematics 349 (2019) 410–423

Fig. 6. Bending of a cantilever.

Fig. 7. Bending of a curved bar.

In the test we take L = 48 m, D = 12 m and P = 1000 N. Since the considered quartic box splines reproduce cubic polynomials, the numerical solution is accurate up to machine precision even on a coarse mesh; see Fig. 6(b). Stress components S11 and S12 are depicted in Fig. 6(c)–(d). 4.2. Bending of a curved bar by a force at the end In this example we simulate another plane strain problem presented in [20]. A curved bar with internal radius a and external radius b is fixed at θ = −π/2, and we apply a force P on the top of the bar (at θ = 0), as shown in Fig. 7(a). The remaining two edges are traction free. In the test we consider a = 13 m, b = 17 m and P = 1 N. The exact displacement in

C. Giannelli et al. / Journal of Computational and Applied Mathematics 349 (2019) 410–423

419

polar coordinates is equal to ur (r , θ ) = − uθ (r , θ ) = −

sin θ

(

E0 cos θ

D(1 − ν0 ) log r + A(1 − 3ν0 )r 2 +

(

E0

D(1 − ν0 ) log r − A(5 + ν0 )r 2 −

B(1 + ν0 ) r2

B(1 + ν0 ) r2

) +K

+

D cos θ E0

+ D(1 + ν0 ) + K

(2θ + π ),

) −

D sin θ E0

(2θ + π ),

for the following parameters M = a2 − b2 + (a2 + b2 ) log(b/a),

A = P /(2M),

B = −Pa b /(2M),

D = −P(a2 + b2 )/M ,

r0 = (a + b)/2,

K = −D(1 − ν0 ) log(r0 ) − A(1 − 3ν0 )r02 − B(1 + ν0 )/r02 .

2 2

The initial mesh in the parametric space is shown in Fig. 7(b) and the stress components of the exact solution are depicted in Fig. 7(c)–(d). The error plot in Fig. 7(e) reveals the optimal order of convergence for the displacement when uniform refinement is applied. Note that we expect order O(h4 ) with h the mesh size since the considered box spline space contains the polynomial space of maximal degree 3, and in 2D it holds O(h4 ) ≈ O(DOF−2 ). 4.3. Infinite plate with a circular hole In this example we model an infinite plate with a circular hole under constant in-plane tension in the x-direction (see, e.g., [1,21]). Due to the symmetry of the problem, only a quarter of a finite part of the infinite plate needs to be modeled (see Fig. 8(a)). The components of the exact stress tensor have the following form in polar coordinates (the Cartesian stress components S11 and S22 are depicted in Fig. 8(b)–(c)):

σrr (r , θ ) =

Tx

σθ θ (r , θ ) =

Tx

( 1−

2

( 1+

2

σr θ (r , θ ) = −

Tx 2

R2

) +

r2 R2

)

r2

( 1+2

− R2 r2

Tx

( 1−4

2 Tx 2

−3

( 1+3 R4 r4

)

R2 r2 R4

+3

)

r4

R4

)

r4

cos 2θ,

cos 2θ,

sin 2θ,

where Tx is the amount of stress applied in x-direction at the left and the right edge of the finite plate, R is the radius of the hole, and L is the length of the quarter of the plate in x- and y-direction. The exact stress is applied as Neumann boundary conditions on the plate (see Fig. 8(a)). Let R = 1 m, L = 4 m and Tx = 10 N. We consider two different techniques to construct the domain geometry. In the first approach, we immerse the domain Ω inside a bigger rectangular domain Ωh and apply the trimming of the mesh as depicted in Fig. 9(a). In the second approach, we achieve the mesh lines to be aligned with the domain boundary using a map F a , which depends on a parameter a. The flexibility of working with a three-directional mesh in the parameter space allows us to easily construct a regular map starting from a pentagonal parameter domain (the length of its bottom and right edge is equal to a; see Fig. 9(b)) to the desired physical domain Ω . In our example we consider two cases, a = 2,3; the resulting mapped meshes are shown in Fig. 9(c)–(d). Adaptive refinement is steered by the quantile marking strategy with threshold parameter ρ = 1/8 (see Section 3.3).

Fig. 8. Infinite plate with a circular hole: problem definition and exact stress components.

420

C. Giannelli et al. / Journal of Computational and Applied Mathematics 349 (2019) 410–423

Fig. 9. Infinite plate with a circular hole: initial meshes in the physical space for the immersed and domain map approach.

Fig. 10. Infinite plate with a circular hole: adaptive meshes (top) and boundary strips (bottom) obtained with trimming (left), map F 2 (center) and map F 3 (right).

Fig. 11. Infinite plate with a circular hole: errors plots related to the displacement (left) and the stress (right). For easier reading, the error estimator values are multiplied by the factor 10−5 .

Fig. 10 depicts the obtained meshes for Ω and ΩΓ after the final adaptive refinement (N = 6). Both the trimmedimmersed and mapped-immersed domain approaches give similar convergence rates to the exact solution (Fig. 11) but

C. Giannelli et al. / Journal of Computational and Applied Mathematics 349 (2019) 410–423

421

Table 1 Infinite plate with a circular hole: condition numbers with respect to the three cases of Fig. 10. Step

1

2

3

4

5

6

trimmed F2 F3

3.40 · 1015 1.59 · 108 1.99 · 108

1.41 · 1017 1.92 · 108 2.50 · 108

1.27 · 1017 9.22 · 108 4.43 · 109

1.11 · 1018 5.06 · 109 6.08 · 1010

1.26 · 1018 5.86 · 1010 8.05 · 1010

3.70 · 1019 5.19 · 1010 5.38 · 1011

Fig. 12. L-shape domain. For easier reading, the error estimators are multiplied by a factor 10−9 .

in the mapped approaches the condition numbers are several orders of magnitude smaller (Table 1). As expected, the local refinement strategy refines the mesh near the circular hole while maintaining the optimal order of convergence. 4.4. L-shape domain The following example is taken from [22]. We consider the L-shaped domain shown in Fig. 12(a), and the exact solution is given in polar coordinates by ur (r , θ ) = uθ (r , θ ) =

1 2µ 1 2µ

r α ((c2 − α − 1) c1 cos((α − 1)θ ) − (α + 1) cos((α + 1)θ )), r α ((c2 + α − 1) c1 sin((α − 1)θ ) + (α + 1) sin((α + 1)θ )),

where c1 = −

cos((α + 1) 34π ) cos((α −

1) 34π

)

,

c2 =

2(λ + 2µ)

λ+µ

.

The parameter α is the nontrivial solution of the equation sin(3π/2)α + sin(3π/2α ) = 0, hence α ≈ 0.5445. All boundary conditions are of Dirichlet type. The displacement is depicted in Fig. 12(b)–(c). Values of stress components have poles at the corner (0, 0). The starting mesh is depicted in Fig. 12(d). Due to the nature of the problem, the convergence order is reduced after a few iterations when we apply the uniform refinement strategy (Fig. 12(f)). Instead, the optimal convergence order is maintained

422

C. Giannelli et al. / Journal of Computational and Applied Mathematics 349 (2019) 410–423

Fig. 13. Triangular support structure.

for the adaptive refinement strategy and the quantile marking parameter ρ = 1/15. The hierarchical mesh after the final refinement step is shown in Fig. 12(e). 4.5. Loading of a triangular support structure In the final example we consider a problem of loading a triangular shaped support structure (Fig. 13(a)). The structure is fixed at the left edge and it is loaded on the top edge with a constant force P, causing a stress S22 = −1. All other remaining edges are free of traction. The exact solution is simulated by taking a numerical solution on a very fine mesh (hmin = 1/128). The displacement components are shown in Fig. 13(b)–(c). The exact solution has a reduced regularity at the three internal corners, and at the two top corners. The initial and the final mesh for the quantile marking adaptive strategy and the threshold parameter ρ = 1/5 are shown in Fig. 13(d)–(e). From Fig. 13(f) we can observe that the optimal convergence order is achieved with the adaptive refinement. 5. Conclusions We have presented the application of an immersed boundary method to an adaptive isogeometric linear elasticity model, based on truncated hierarchical box splines. In particular, we have shown that the possibility of combining a geometry mapping with an immersed approach can improve the condition number with respect to the standard mapping-free immersed case based on trimmed mesh cells. In addition, the adaptive nature of the method guarantees high accuracy with a reduced number of degrees of freedom, with respect to the uniform refinement case. Additional flexibility could be achieved by constructing parameterizations more suitable for immersing the domain in (possibly trimmed) geometries. Different requirements could then be taken into account for a better coupling of isogeometric and immersed methods, which, for example, does not necessarily require an exact description of the physical domain. The approach can be applied to general domains in 2D and 3D. Acknowledgments This work was partially supported by the MIUR ‘‘Futuro in Ricerca’’ programme through the project DREAMS (RBFR13FBI3) and by the Mission Sustainability Program of the Università degli Studi di Roma ‘‘Tor Vergata’’ through the project IDEAS (E81I18000060005). The authors are members of the INdAM Research group GNCS. The INdAM support through GNCS and Finanziamenti Premiali SUNRISE is gratefully acknowledged.

C. Giannelli et al. / Journal of Computational and Applied Mathematics 349 (2019) 410–423

423

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs, Isogeometric Analysis: Toward Integration of CAD and FEA, John Wiley & Sons, 2009. C. de Boor, K. Höllig, S. Riemenschneider, Box Splines, Springer–Verlag, 1993. L.L. Schumaker, Spline Functions: Basic Theory, third ed., Cambridge University Press, 2007. D. Schillinger, M. Ruess, N. Zander, Y. Bazilevs, A. Düster, E. Rank, Small and large deformation analysis with the p- and B-spline versions of the finite cell method, Comput. Mech. 50 (2012) 445–478. J. Baiges, R. Codina, F. Henke, S. Shahmiri, W.A. Wall, A symmetric method for weakly imposing Dirichlet boundary conditions in embedded finite element meshes, Int. J. Numer. Methods Engrg. 90 (2012) 636–658. S. Kollmannsberger, A. Özcan, J. Baiges, M. Ruess, E. Rank, A. Reali, Parameter-free, weak imposition of Dirichlet boundary conditions and coupling of trimmed and non-conforming patches, Int. J. Numer. Methods Engrg. 101 (2015) 670–699. R.A.K. Sanches, P.B. Bornemann, F. Cirak, Immersed b-spline (i-spline) finite element method for geometrically complex domains, Comput. Methods Appl. Mech. Engrg. 200 (2011) 1432–1445. M.C. Hsu, C. Wang, A.J. Xu, A. Krishnamurthy, Direct immersogeometric fluid flow analysis using B-rep CAD models, Comput. Aided Geom. Design 43 (2016) 143–158. B. Marussig, T.J.R. Hughes, A review of trimming in isogeometric analysis: challenges, data exchange and simulation aspects, Arch. Computat. Methods Eng. (2017). C. Wang, F. Xu, M.C. Hsu, A. Krishnamurthy, Rapid B-rep model preprocessing for immersogeometric analysis using analytic surfaces, Comput. Aided Geom. Design 52–53 (2017) 190–204. F. Pelosi, C. Giannelli, C. Manni, M.L. Sampoli, H. Speleers, Splines over regular triangulations in numerical simulation, Comput. Aided Design 82 (2017) 100–111. H. Speleers, Inner products of box splines and their derivatives, BIT Numer. Math. 55 (2015) 559–567. C. Giannelli, B. Jüttler, H. Speleers, THB-splines: The truncated basis for hierarchical splines, Comput. Aided Geom. Design 29 (2012) 485–498. C. Giannelli, B. Jüttler, H. Speleers, Strongly stable bases for adaptively refined multilevel spline spaces, Adv. Comput. Math. 40 (2014) 459–490. U. Zore, B. Jüttler, Adaptively refined multilevel spline spaces from generating systems, Comput. Aided Geom. Design 31 (2014) 545–566. T. Kanduč, C. Giannelli, F. Pelosi, H. Speleers, Adaptive isogeometric analysis with hierarchical box splines, Comput. Methods Appl. Mech. Engrg. 316 (2017) 817–838. H. Speleers, C. Manni, Effortless quasi-interpolation in hierarchical spaces, Numer. Math. 132 (2016) 155–184. R. Verfürth, A review of a posteriori error estimation techniques for elasticity problems, Comput. Methods Appl. Mech. Engrg. 176 (1999) 419–440. A. Buffa, C. Giannelli, Adaptive isogeometric methods with hierarchical splines: error estimator and convergence, Math. Models Methods Appl. Sci. 26 (2016) 1–25. S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, McGraw–Hill, 1970. P.L. Gould, Introduction to Linear Elasticity, Springer–Verlag, 1999. J. Alberty, C. Carstensen, S.A. Funken, R. Klose, Matlab implementation of the finite element method in elasticity, Computing 69 (2002) 239–263.