Isogeometric Analysis on V-reps: first results

TL;DR

This paper introduces a new isogeometric analysis method for solving elliptic PDEs on trimmed geometries represented as V-reps, with theoretical validation and practical testing in 2D and 3D.

Contribution

It develops a novel approach for isogeometric analysis on V-reps, including approximation tools and re-parametrization techniques for trimmed geometries.

Findings

Validated on 2D and 3D diffusion problems

Successfully applied to linear elasticity problems

Theoretical framework supports algorithmic choices

Abstract

Inspired by the introduction of Volumetric Modeling via volumetric representations (V-reps) by Massarwi and Elber in 2016, in this paper we present a novel approach for the construction of isogeometric numerical methods for elliptic PDEs on trimmed geometries, seen as a special class of more general V-reps. We develop tools for approximation and local re-parametrization of trimmed elements for three dimensional problems, and we provide a theoretical framework that fully justify our algorithmic choices. We validate our approach both on two and three dimensional problems, for diffusion and linear elasticity.

Equations

$$\Omega ={\Omega }_0\setminus \bigcup _{i=1}^N{\bar{\Omega }}_i.$$

$$V=\{u\in H^1(\Omega {)}^k:u{|}_{{\Gamma }_D}=0\},$$

$$\text{Find }u\in V:a(u,v)={\int }_{\Omega }fv+{\int }_{{\Gamma }_N}gv\forall v\in V.$$

$${\mathbb{B}}_{0,j}:=\{B_{0,j,k}={\hat{B}}_{0,j,k}\circ {\mathbf{F}}_{0,j}^{-1},k\in {\mathcal{I}}_{0,j}\};$$

$${\mathcal{T}}_h(\Omega ):=\{Q:\forall Q\in {\mathcal{T}}_h({\Omega }_0):Q\cap \Omega \ne \varnothing \}.$$

$$\dim ({\mathbb{V}}_h(\Omega ))=\#{\mathcal{I}}_0^{act}.$$

$${\mathbb{V}}_h(\Omega ,{\Gamma }_D)=\{v\in {\mathbb{V}}_h(\Omega ):v{|}_{{\Gamma }_D}=0\};$$

$$\Vert v-{\Pi }_h(v){\Vert }_{H^s(\Omega )}$$

$$\Vert v-{\Pi }_h^{(D)}(v){\Vert }_{H^s(\Omega )}$$

$$\Vert \tilde{v}{\Vert }_{H^t({\mathbb{R}}^3)}\le C\Vert v{\Vert }_{H^t(\Omega )}t\ge 0,$$

$$\text{Find }{u}_h\in V_h(\Omega ):a(u_h,v_h)={\int }_{\Omega }f{v}_h+{\int }_{{\Gamma }_N}g{v}_h\forall v_h\in V_h(\Omega ).$$

$${\gamma }_i=\{(X_i,\psi (X_i));X_i\in {\Delta }_i\}\subseteq \gamma ,{\gamma }_{h,i}=\{(X_i,{\psi }_h(X_i));X_i\in {\Delta }_i\}\subseteq {\gamma }_h.$$

$$\bigcup _{i=1}^{n_Q}{\gamma }_i=\gamma ;\bigcup _{i=1}^{n_Q}{\gamma }_{h,i}={\gamma }_h;\sum _{i=1}^{n_Q}\mathrm{meas}({\gamma }_i)\le C\mathrm{meas}(\gamma );\sum _{i=1}^{n_Q}\mathrm{meas}({\gamma }_{h,i})\le C\mathrm{meas}({\gamma }_h).$$

$$\Vert \psi -{\psi }_h{\Vert }_{W^{m,\infty }({\Delta }_i)}\le C{h}^{r+1-m}m=0,1,$$

$$\Vert \tilde{u}{\Vert }_{H^s({\Omega }_h)}\le C\Vert u{\Vert }_{H^s(\Omega )}\forall s>0,$$

$$\text{Find }\widetilde{u}_{h}\in V_{h}(\Omega_{h})\ :\ a_{h}(\widetilde{u}_{h},v_{h})=L_{h}(v_{h})=\int_{\Omega_{h}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{f}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}}\ v_{h}+\int_{\Gamma_{h,N}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{g}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}}\,v_{h}\qquad\forall v_{h}\in V_{h}(\Omega_{h})\,.$$

$$\Vert \tilde{u}-{\tilde{u}}_h{\Vert }_{H^1({\Omega }_h)}\le C{h}^t\Vert u{\Vert }_{H^{t+1}(\Omega )}+C{h}^r\Vert \tilde{u}{\Vert }_{H^{n/2+2+\alpha }(\Omega )}2\le t\le p,0\le r\le p,$$

$$\Vert \tilde{u}-{\tilde{u}}_h{\Vert }_{L^2({\Omega }_h)}\le C{h}^{t+1}\Vert u{\Vert }_{H^{t+1}(\Omega )}+C{h}^{r+1}\Vert \tilde{u}{\Vert }_{H^{n/2+2+\alpha }(\Omega )}2\le t\le p,0\le r\le p.$$

$$\begin{array}{rl}\alpha \Vert {\tilde{u}}_h-w_h{\Vert }_{H^1({\Omega }_h)}^2& \le a_h({\tilde{u}}_h-w_h,{\tilde{u}}_h-w_h)=a_h({\tilde{u}}_h,{\chi }_h)-a_h(w_h,{\chi }_h)=L_h({\chi }_h)-a_h(w_h,{\chi }_h)\\ & \le a_h(\tilde{u}-w_h,{\chi }_h)+\Vert {\chi }_h{\Vert }_{H^1({\Omega }_h)}\underset{{\xi }_h\in V_h({\Omega }_h)}{\sup }\frac{L_h({\xi }_h)-a_h(\tilde{u},{\xi }_h)}{\Vert {\xi }_h{\Vert }_{H^1({\Omega }_h)}}.\end{array}$$

$$a_h(\tilde{u},{\xi }_h)={\int }_{{\Omega }_h}A\tilde{u}{\xi }_h+{\int }_{{\Gamma }_{h,N}}(\nabla \tilde{u}\cdot n_h){\xi }_h,L_h({\xi }_h)={\int }_{{\Omega }_h}\tilde{f}{\xi }_h+{\int }_{{\Gamma }_{h,N}}\tilde{g}{\xi }_h,$$

$$L_h({\xi }_h)-a_h(\tilde{u},{\xi }_h)={\int }_{{\Omega }_h\setminus \Omega }(\tilde{f}-A\tilde{u}){\xi }_h+{\int }_{{\Gamma }_{h,N}}(\tilde{g}-\nabla \tilde{u}\cdot n_h){\xi }_h=I+II,$$

$$\Vert w{\Vert }_{L^2({\text{Q}}_h\setminus \text{Q})}\le C(h^{r+1}|w{|}_{H^1(\text{Q}\setminus {\text{Q}}_h)}+h^{\frac{1}{2}(r+1)}\Vert w{\Vert }_{L^2(\partial {\text{Q}}_h)}),$$

$$h^{r+1}|w|_{H^{1}(\textsf{Q}_{h}\setminus\textsf{Q})}+h^{\frac{1}{2}(r+1)}\|w\|_{L^{2}(\partial\textsf{Q}_{h})}\leq Ch^{\frac{1}{2}(r+1)}\|{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}}\|_{H^{1}(Q)}\,.$$

$${\Vert \nabla w\cdot n_h-\tilde{g}\Vert }_{L^2({\gamma }_h)}\le C{h}^r\Vert w{\Vert }_{W^{2,\infty }(Q)}.$$

$$I$$

$$\displaystyle\leq C\sum_{Q\in{\cal T}_{h}^{\Gamma}}\big{[}h^{r+1}|\widetilde{f}-A\widetilde{u}|_{H^{1}(\textsf{Q}_{h}\setminus\textsf{Q})}+h^{\frac{1}{2}(r+1)}\|\widetilde{f}-A\widetilde{u}\|_{L^{2}(\partial\textsf{Q}_{h})}\big{]}\cdot$$

$$\displaystyle\qquad\qquad\big{[}h^{r+1}|\xi_{h}|_{H^{1}(\textsf{Q}_{h}\setminus\textsf{Q})}+h^{\frac{1}{2}(r+1)}\|\xi_{h}\|_{L^{2}(\partial\textsf{Q}_{h})}\big{]}\,.$$

$$I$$

$$\le C(h^{r+1}\Vert {\xi }_h{\Vert }_{L^2(\partial {\Omega }_h)}+h^{\frac{3}{2}(r+1)}\Vert {\xi }_h{\Vert }_{H^1({\Omega }_h)})\Vert u{\Vert }_{H^3(\Omega )}$$

$$\le C{h}^{r+1}\Vert {\xi }_h{\Vert }_{H^1({\Omega }_h)}\Vert u{\Vert }_{H^3(\Omega )},$$

$$II$$

$$\le C{h}^r\Vert \tilde{u}{\Vert }_{W^{2,\infty }(\Omega )}\Vert {\xi }_h{\Vert }_{H^1({\Omega }_h)}\le C{h}^r\Vert \tilde{u}{\Vert }_{H^{n/2+2+\alpha }({\Omega }_h)}\Vert {\xi }_h{\Vert }_{H^1({\Omega }_h)}.$$

$$|L_h({\xi }_h)-a_h(\tilde{u},{\xi }_h)|\le C{h}^r\Vert {\xi }_h{\Vert }_{H^1({\Omega }_h)}\Vert u{\Vert }_{H^{n/2+2+\alpha }(\Omega )},$$

Full text

Keywords: numerical methods for PDEs, isogeometric methods, trimmed geometries, V-reps.

1 Introduction

Back in 2005, T. Hughes and co-authors introduced isogeometric analysis (IGA) [29] with the promise to alleviate the burden of conversion between computer-aided design (CAD) geometries and finite element (FE) computational domains. IGA has been one of the most successful ideas in the last decades in computational mechanics, and there are really a huge amount of contributions to the field, also because splines have been proved to be a very powerful tool for the approximation of solutions of partial differential equations (PDEs), [6]. Indeed, IGA has been a tremendous success since 2005 with a wide range of applications (see for instance the special issue [30]), and it is becoming a mature method: its mathematical analysis is now well understood [6], fast assembly and solvers exist today [1, 52] and strategies for adaptive refinement with a sound mathematical theory are now available, see [7] and the references therein.

Nevertheless, the construction of geometric representations suitable to IGA, from CAD geometries, remains the major unsolved issue and IGA is indeed a very efficient strategy on “simple” geometries (i.e., on 2D or $2.5$D geometries [10, 4, 27]) but not for general three dimensional problems. For 3D objects the conversion of a CAD geometry into a geometry on which isogeometric methods are defined (IGA-ready) is challenging and has many similarities with the conversion to a FE domain: its automation is far from being at the state of the art. The main reason why this problem is still open is that CAD geometric descriptions are based two main ingredients: i) they represent only the boundary of geometries and not their interior, i.e., B-reps [8, 48] ; ii) they allow for boolean operations among spline surfaces, primitives and so on. The result is a complex boundary representation that by no means corresponds to a valid mesh. For two dimensional problems, this issue was successfully addressed for the first time in [44], and then by many other authors. Interest readers may refer to the review [40] and the many references therein.

As a matter of fact, this has limited the use of isogeometric techniques for complex three dimensional domains. The only viable approach has been, until now, the use of spline basis on unstructured hexahedral (or tetrahedral) meshes, which has been object of several contributions [38, 59, 61, 63]. In this context, the presence of extraordinary points and edges make the construction of regular B-spline functions preserving accuracy extremely challenging (see [62, 33] and the references therein).

On the other hand, a new paradigm has recently been introduced by Elber et al. in [42], according to which geometries are described as V-reps, i.e., via volumetric representations. V-reps are defined as a generalization of B-reps as boundaries of V-reps are always B-reps. Geometries are described thorough the geometric representation of the volume they occupy and the basic building blocks for V-reps are trivariate (as opposed to bivariate) representations. Boolean operations are handled at the level of volumes instead of surfaces. Loosely speaking, V-reps are a class of volumes in which volumetric boolean operations are allowed and produce V-reps.

As a matter of fact, V-reps are not IGA-ready geometries as the handling of boolean operations at the level of simulation of PDEs is not at the state of the art, and is very much linked with the numerical treatment of trimming. In this paper, we initiate the construction of accurate isogeometric methods on V-reps. In particular, we focus on the design of a numerical methodology for the simulation of elliptic problems on computational domains described in terms of V-reps, and the main rule we adopt in our design is the following: we wish to construct numerical methods that uses the V-rep representation without asking for meshing or global re-parameterization, while we will make use of local high-order re-parameterizations.

Our work is related to previous efforts in the construction of numerical methods on trimmed domains and of immersed finite elements approaches. In particular, we should mention the Finite Cell approach, which was successfully applied in the isogeometric framework in [53, 51, 35, 60] and the references therein; the CutFEM and CutIGA methods [15, 24]; as well as other high-order unfitted Galerkin discretizations, e.g. [37]. Subdivision techniques are used in [49] and boundary corrections are proposed in [50].

These families of immersed methods present numerous similiarities with the approach presented in this and previous works (see, e.g. [44, 40], for the case of 2D geometries). Indeed, those unfitted methods define the solution discretization by means of a background mesh that is completely disconnected from geometry boundary; while, following an isoparametric paradigm, the approach that we present here considers the underlying V-rep’s spline space as a basis for the solution discretization. Nonetheless, the results discussed in this work, namely, the presented error estimate and the construction of high-order local re-parameterizations, can be used directly also in the case of finite element and isogeometric immersed methods.

We say that a volume $\Omega$ is a non-conforming multipatch trivariate volume (nCMTV), if it is constructed as a collection of trivariates (patches), $T_{1},\ldots T_{k}$ such that:

• •

For each $T_{\ell}$, there is a spline diffeomorphism ${\bf F}_{\ell}:\widehat{T}\to T_{\ell}$ , $\ell=1,\ldots,k$ ($\widehat{T}$ stands for the parameter space);

• •

$\bar{T}_{i}\cap\bar{T}_{j}$ is either empty or is the image of a full face $\widehat{f}_{i}$ of $\widehat{T}$ for $T_{i}$ and a full face $\widehat{f}_{j}$ of $\widehat{T}$ for $T_{j}$;

• •

the parameterizations $F_{i}|_{\widehat{f}_{i}}$ and $F_{j}|_{\widehat{f}_{j}}$ may differ.

If the parameterizations $F_{i}|_{\widehat{f}_{i}}$ and $F_{j}|_{\widehat{f}_{j}}$ are identical, then we say that $\Omega$ is a conforming multipatch trivariate volume (CMTV). Following [42], volumetric representations are obtained as the closure of the nCMTV volumes with respect to the boolean operations of union, intersection and their combinations.

In this first contribution we consider only domains generated by subtraction (or intersection). We assume that we are given a collection of nCMTVs $\Omega_{0},\Omega_{1}...\Omega_{N}$, and we consider the computational domain $\Omega$ defined as follows:

[TABLE]

Finally, for the sake of simplicity we will assume that $\Omega_{0}$ has indeed conforming meshes. If this is not the case, non conforming multipatch interfaces should be handled via mortaring techniques [11]. The aim of this paper is to propose algorithms to tackle the isogeometric simulation of linear elliptic problems on such a computational domain, without any re-parameterization and within the isoparametric paradigm. As we concentrate of Neumann problems only, after reminding how the main approximation properties of splines ensure approximability of solutions on $\Omega$, we concentrate on the main issues: the definition of an integration formula on trimmed elements, and the study of the conditioning and pre-conditioning of the underlying linear system. We discuss the approximation of the trimmed boundaries and provide approximation estimates which are validated against numerical benchmarks.

Our contribution is organized as follows: in Section 2 we defined our model problem and corresponding isogeometric discretization, together with our assembly algorithms and a brief discussion about the linear system conditioning. In Section 3 we will provide numerical testing, we discuss both academic tests for validation of our algorithms, and also one rather complex application in linear elasticity. We draw our conclusions in Section 4.

2 V-reps as computational domain for PDEs

Given the domain $\Omega$ defined in (1), we think of it as the domain of definition of elliptic partial differential equations. Indeed, we consider two model problems: the Laplace problem as model of diffusion phenomena and compressible linear elasticity as a simplified model of the deformation of solids.

The boundary of the domain $\Omega$ defined in (1) naturally splits in two parts: a non-trimmed part of boundary $\partial\Omega\cap\partial\Omega_{0}$, and a trimmed boundary $\partial\Omega\setminus\partial\Omega_{0}$.

We denote by $\Gamma_{D}$ a connected subset (with Lipschitz boundary) of $\partial\Omega$. For the time being, we assume that $\Gamma_{D}\subset\partial\Omega_{0}\cap\partial\Omega$: this means that we allow for essential boundary conditions only on the non-trimmed part of the boundary. For the sake of simplicity, we assume that $\Gamma_{D}\neq\emptyset$. If this is not the case, all what follows remain valid, once the definition of spaces are adapted and suitable compatibility conditions of data are satisfied. As usual, we set $\Gamma_{N}=\partial\Omega\setminus\bar{\Gamma}_{D}$.

We define:

[TABLE]

where $k$ is either $1$ for diffusion problems or $3$ for elasticity problems. For both problems, we assume to have a bilinear form $a(\cdot,\cdot)$ acting on $V\times V$ and verifying:

• •

Continuity: there exists a $M>0$ such that $a(u,v)\leq M\|u\|_{V}\|v\|_{V}$;

• •

Coercivity: there exists a $\alpha>0$ such that $a(u,u)\geq\alpha\|u\|^{2}_{V}$.

The continuous problem to solve is formulated as follows and admits a unique solution:

[TABLE]

If needed, we will denote by $A$ the differential operator in strong form. I.e., $\langle Au,v\rangle_{V^{\prime}\times V}=a(u,v)$.

Remark 2.1**.**

Essential boundary conditions can be imposed only on the non-trimmed part of boundary because in the discrete setting described below, the imposition of essential boundary condition on the trimmed part of the boundary ask for special care [16, 17, 13]. The design of robust techniques to deal with this issue in the framework of trimmed isogeometric analysis is beyond the scope of the present paper. Contributions in this direction can be found in [23, 28, 41].

2.1 The isogeometric method on V-reps

Thanks to the construction above, the master volume $\Omega_{0}$ is described by $n_{0}$ trivariates $T_{0,1},\ldots T_{0,n_{0}}$ and $n_{0}$ parameterizations ${\bf F}_{0,1},\ldots,{\bf F}_{0,n_{0}}$. These parameterizations are constructed on open knots vectors that are defined in the parameter trivariate space $\widehat{T}$. We denote by $h_{j}$ the diameter of the largest knot span in $T_{0,j}$.

On each patch composing $\Omega_{0}$, we have a collection of B-splines defined on such open knot vectors:

[TABLE]

where $\mathcal{I}_{0,j}$ denotes the running indices on the B-spline functions on the patch $T_{0,j}$.

We denote by $\mathbb{B}^{\textrm{disc}}(\Omega_{0}):=\bigcup_{j=0}^{n_{0}}\mathbb{B}_{0,j}$ and by $\mathbb{V}_{h}^{\textrm{disc}}(\Omega_{0})=\mathrm{span}\{\mathbb{B}^{\textrm{disc}}(\Omega_{0})\}$, where we have set $h=\max_{j}h_{j}$. From now on, and for the sake of simplicity, we assume that $\Omega_{0}$ is a CMTV. We also construct the space $\mathbf{v}_{h}(\Omega)=\mathbb{V}_{h}^{\textrm{disc}}(\Omega_{0})\cap H^{1}(\Omega_{0})$. It is well known that a basis for this space can be constructed starting from $\mathbb{B}^{\textrm{disc}}(\Omega_{0})$ and identifying the indices of functions corresponding to coincident control points at the interfaces. In what follows, we suppose to have a set of indices $\mathcal{I}_{0}:=\{0,1,\ldots,I_{0}\}$ and a collection of B-splines functions $\mathbb{B}(\Omega_{0})$, indexed by $\ell\in\mathcal{I}_{0}$, that span $\mathbb{V}_{h}(\Omega)$. The generic basis function is denoted as ${\mathcal{B}}_{0,\ell}\,,\;\ell\in\mathcal{I}_{0}$. Finally we denote by $\mathbb{V}_{h}(\Omega_{0})$ the space spanned by these functions.

Mapping knot surfaces through the mappings ${\bf F}_{0,j}$, we obtain a mesh, $\mathcal{T}_{h}({\Omega_{0}})$ which is, indeed, a curvilinear hexahedral partition of $\Omega_{0}$. As $\Omega_{0}$ is a CMTV, this partition is a conforming decomposition of the domain $\Omega_{0}$.

Refinement on such a mesh can be performed by knot insertion, and for the sake of our subsequent analysis, we suppose we have a family of refined meshes (obtained by dyadic refinement of all non empty knot spans) and corresponding refined basis functions. Such a refinement will not change the geometry of the volume $\Omega_{0}$ but only its representation. In this way, we will assume to have a family of meshes $\mathcal{T}_{h}({\Omega_{0}})$, of spaces $\mathbb{V}_{h}(\Omega_{0})$ indexed with $h$. The dependence on $h$ will be omitted in what follows when it does not matter.

As our computational domain is $\Omega$ and not $\Omega_{0}$, we define $\mathcal{T}_{h}(\Omega)$ as follows:

[TABLE]

Note that, by an abuse of notation, we could also write: $\mathcal{T}_{h}(\Omega)=\mathcal{T}_{h}({\Omega_{0}})\cap\Omega$.

For further use, among the elements in $\mathcal{T}_{h}(\Omega)$, we distinguish two families:

• •

$Q\in\mathcal{T}_{h}(\Omega)$ such that ${Q}\cap\Omega=Q$, and we denote their collection as $\mathcal{T}^{int}_{h}(\Omega)$;

• •

$Q\in\mathcal{T}_{h}(\Omega)$ such that ${Q}\cap\Omega\neq Q$, and we denote their collection as $\mathcal{T}^{\Gamma}_{h}(\Omega)$; when $Q\in\mathcal{T}^{\Gamma}_{h}(\Omega)$, we will also make use of $\mathsf{Q}=Q\cap\Omega$.

The indices $\ell\in\mathcal{I}_{0}$, such that $\operatorname{supp}\{{\mathcal{B}}_{0,\ell}\}\cap\Omega\neq\emptyset$, are denoted by $\mathcal{I}_{0}^{act}$, whereas $\#\mathcal{I}_{0}^{act}$ stands for its cardinality. We have the following:

Lemma 2.2**.**

It holds:

[TABLE]

Proof.

This is an immediate consequence of the local linear independence of B-splines: i.e., their restriction to any element $Q$ is a collection of $(p+1)^{3}$ linear independent functions representing polynomials of degree $p$ in each coordinate direction. ∎

Thanks to our assumption on $\Gamma_{D}$, i.e. $\Gamma_{D}\subset\partial\Omega_{0}\cap\partial\Omega$, we can define the space:

[TABLE]

For both spaces $\mathbb{V}_{h}({\Omega,\Gamma_{D}})$ and $\mathbb{V}_{h}(\Omega)$, we also have the following approximation property:

Lemma 2.3**.**

There exist continuous interpolation operators $\Pi_{h}:L^{2}(\Omega)\to\mathbb{V}_{h}(\Omega)$ and $\Pi_{h}^{(D)}:L^{2}(\Omega)\to\mathbb{V}_{h}(\Omega,\Gamma_{D})$ such that:

[TABLE]

for all $s=0,1$ and $1\leq t\leq p+1.$

Proof.

Given a $v\in H^{s+1}(\Omega)$, we denote by $\widetilde{\cdot}$ the Sobolev extension operator as defined in [55, theorem 5, page 181] that continuously extends any function defined in $\Omega$ to functions defined in $\mathbb{R}^{3}$ and verifies:

[TABLE]

where the constant $C$ does not depend on $t$. Now, given a $v\in H^{t}(\Omega)$, the local quasi-interpolant defined in [6], or the one defined in [12], can be applied to $\widetilde{v}$ and the result follows. A similar idea applies to the case with Dirichlet boundary conditions, thanks to the assumption $\Gamma_{D}\subset\partial\Omega_{0}\cap\partial\Omega$. ∎

We are now in the position to write our first discrete problem. We set $V_{h}(\Omega)=\{\mathbb{V}_{h}(\Omega,\Gamma_{D})^{k}\}$ (as before, $k=1$ or $3$), and we wish to solve:

[TABLE]

As $V_{h}(\Omega)\subset V$, both continuity and coercivity hold true, which implies that the discrete problem is well posed and that the solution depends continuously on the data.

As a matter of fact, two types of problems may appear when trying to construct the linear system for the problem (8):

We need to be able to compute integrals as $\int_{Q\cap\Omega}\xi$ where $\xi$ is a regular function (say a multiplication of splines, regular coefficients, Jacobians of transformations) but the integration domain is not known analytically for all $Q\in\mathcal{T}^{\Gamma}_{h}(\Omega)$. 2. 2.

The contribution of all those basis functions whose support intersects the computational domain only for a small portion may generate bad conditioning of the stiffness matrices.

Clearly the first observation is crucial since, as it is, the assembly of the stiffness matrix for the problem (8) is not possible. The next two subsections are devoted to the design of an algorithm to compute integrals accurately. The latter problem is discussed in Section 2.4.

2.2 Numerical integration error estimate

In order to compute each contribution to the stiffness in an accurate way, we propose to:

• a)

approach each boundary $\gamma=Q\cap\partial\Omega$, for $Q\in\mathcal{T}_{h}^{\Gamma}(\Omega)$, with a suitable piecewise polynomial approximation surface (or line in 2D), that we call $\gamma_{h}$.

• b)

Construct a re-parameterization of each approximated element which allows for the definition of a suitable integration formula.

In the sequel of this section we prove that if the approximation at a) verifies certain assumptions, then the accuracy of the isogeometric method is preserved. In the next subsection instead, we provide algorithms to construct such (piecewise) polynomial approximations together with a strategy to construct local accurate integration formulae.

Given a $Q\in\mathcal{T}_{h}^{\Gamma}(\Omega)$, we denote by $\gamma=Q\cap\partial\Omega$ and by $\gamma_{h}$ its approximation of degree $r$ (when needed, a subindex $Q$ will be used to identify the element to which $\gamma$ and $\gamma_{h}$ belong) . In what follows we assume that, for all elements $Q\in\mathcal{T}_{h}^{\Gamma}(\Omega)$, the construction of $\gamma_{h}$ fulfills the following Assumption, with constants that do not depend on the specific $Q$.

Assumption 2.4**.**

We assume that

$\gamma_{h}\subset\bar{Q}$** 2. 2.

there exists a finite collection of local charts so that both $\gamma$ and $\gamma_{h}$ can be described as local graphs. More precisely, there exists a finite number $n_{Q}$ of local systems of coordinates $(X_{i},Y_{i})$, of corresponding neighborhoods $\Delta_{i}\subseteq\mathbb{R}^{n-1}$ and, for each of those, two functions $\psi_{i}:\Delta_{i}\to\mathbb{R}$ and $\psi_{h,i}:\Delta_{i}\to\mathbb{R}$:

[TABLE]

Moreover, there exists a constant C such that

[TABLE] 3. 3.

the following approximation property holds:

[TABLE]

where $r$ is the degree of the approximation $\gamma_{h}$ of $\gamma$.

Reminding that $\mathsf{Q}=Q\cap\Omega$, we denote now by $\mathsf{Q}_{h}$ the element obtained from $\mathsf{Q}$ by replacing its boundary $\gamma\subset\partial\Omega$ by the corresponding $\gamma_{h}$.

Similarly, we denote by $\Omega_{h}$ the domain obtained by replacing each $\mathsf{Q}$, associated with $Q\in\mathcal{T}_{h}^{\Gamma}(\Omega)$ with the corresponding $\mathsf{Q}_{h}$; by $\Gamma_{h,N}$ the approximation of the Neumann boundary; and by $a_{h}(\cdot,\cdot)$ and $L_{h}(\cdot{})$ the bilinear and right hand side forms for the problem (8) when $\Omega$ is replaced by $\Omega_{h}$. Moreover, we assume that the data of the problem $f$ and $g$ are restrictions of $\widetilde{f}$ and $\widetilde{g}$ that are defined in $\mathbb{R}^{3}$ and belong to the corresponding Sobolev space (i.e., for $f\in H^{s}(\Omega)$ then $\widetilde{f}\in H^{s}(\mathbb{R}^{3})$, and if $g\in H^{s-1/2}(\Gamma_{N})$ then $\widetilde{g}\in H^{s}(\mathbb{R}^{3})$). As for the solution $u$, by a little abuse of notation, we denote by $\widetilde{u}$ its Sobolev extension, as defined in the well known [55, Theorem 5, p.181]. We will make use of the following

[TABLE]

where the constant $C$ does not depend on $s$. Thus, instead of (8), we solve the following problem:

[TABLE]

In the sequel we prove that, under the Assumption 2.4, the problem (10) is an optimal approximation of (2). Note that, for the sake of simplicity, we assume that on each approximated element, the integration is performed exactly. Of course, this is not the case in practice but, instead, the quadrature is chosen in order for the integration error to verify the classical assumptions of the finite element theory. Taking into account this error here would make our proofs unreasonably tedious, without adding much. We refer the interested reader to [18].

Our approach is largely inspired by the papers [2, 3] and uses the technical results obtained in [58]. We basically adapt their approach to our problem.

We are going to prove the following approximation result.

Theorem 2.5**.**

Let $\Omega\subset\mathbb{R}^{n}$, $n=2,3$ be a computational domain defined as in (1). Let $u$ be the solution of (2). We assume that $u\in H^{t+1}(\Omega)$, with $t+1>\frac{n}{2}+2$, and we denote by $\widetilde{u}$ its Sobolev extension. Let $\widetilde{u}_{h}$ be the solution of (10). Then, under the Assumption 2.4, it holds:

[TABLE]

where $p$ denotes the degree of the spline approximation, $r$ the degree of the local re-parameterization, $n=2,3$ is the dimension of the ambient space, and $\alpha>0$. Moreover, it holds:

[TABLE]

Remark 2.6**.**

Indeed, the error due to the geometry approximation depends upon the $W^{2,\infty}(\Omega)$ norm of the (extension of) the solution $u$, and we have used here the continuous Sobolev embedding $H^{n/2+\alpha}(\Omega)\subset L^{\infty}(\Omega)$.

Proof.

By application of the Strang Lemma, see e.g., [18], and using coercivity and continuity, we can control the distance between $\widetilde{u}_{h}$ and any function $w_{h}\in V_{h}(\Omega_{h})$ (sometimes below we also use $\chi_{h}=\widetilde{u}_{h}-w_{h}$) as follows:

[TABLE]

We need to estimate the supremum in the right hand side.

[TABLE]

where $n_{h}$ denotes the unit vector normal to $\Gamma_{h,N}$, and then

[TABLE]

where we have used that $A\widetilde{u}=f$ pointwise for almost every $x\in\Omega$. The two terms $I$ and $II$ are studied separately, and we will make use of a few important estimates which hold true under the Assumption (2.4). We use the estimates shown in, e.g., [3, Lemma 4.1], but consider their local versions. For all $Q\in{\cal T}_{h}^{\Gamma}$ that for all $w\in H^{1}(Q)$, letting $\textsf{Q}=Q\cap\Omega$ and $\textsf{Q}_{h}$ is its approximation, it holds:

[TABLE]

where the constant $C$ does not depend on the specific $Q$, neither on the way $\gamma$ or $\gamma_{h}$ cuts the element $Q$. Keeping in mind that the trace operator is continuous from $L^{2}(\partial\textsf{Q}_{h})$ to $H^{1}(Q)$, with a continuity constant that does not depend on the way $\gamma_{h}$ cuts the element $Q$, sometimes we make use of the following estimate of the right hand side:

[TABLE]

We also use then the estimates proved by Čermák in [58, Lemma 3.2 and Lemma 3.3], and we apply them again locally to a Q. For all $w\in W^{2,\infty}(Q)$, verifying $\nabla w\cdot n-g=0$ on $\gamma=Q\cap\partial\Omega$, it holds:

[TABLE]

We start by estimating the term $I$ in the r.h.s. of (15), using (16):

[TABLE]

By (17), we can bound the r.h.s. as:

[TABLE]

where we have used the continuity of the trace operator from $L^{2}(\partial\Omega_{h})$ to $H^{1}(\partial\Omega_{h})$, and that its continuity constant does not depend on $h$, but may depend on the constants in the Assumption 2.4.

We should now estimate the second term in the r.h.s. of (15), using (18) and the Sobolev embedding theorem. We have, for $\alpha>0$:

[TABLE]

Thus, using (19), (20) and the continuity of the Sobolev extension, we obtain,

[TABLE]

which, together with (13), proves (11). We are left now with the proof of (12) and we are going to use a simple and classical duality argument. Our procedure is very similar to [9] and [3].

It holds

[TABLE]

We solve now,

[TABLE]

As a consequence of the regularity assumption we have done on the solution $u$, $z\in H^{2}(\Omega)$, and it holds:

[TABLE]

Now, we compute, for the right hand side of (21)

[TABLE]

where we have used that $u-u_{h}=0$ on $\Gamma_{D}$, and $\nabla z\cdot n=0\text{ on }\Gamma_{N}.$ Now, we estimate the two terms $I$ and $II$ separately. First, the term $II$ is simple to estimate by (16) and (17):

[TABLE]

as $r\geq 1$. Now, we concentrate on the term $I$. It holds:

[TABLE]

where the second term can easily be bounded with (16) and (17) as

[TABLE]

Now, the first term can be written as, for $z_{h}\in\mathbb{V}(\Omega_{h})$ being the best fit of $z$:

[TABLE]

The second term can be estimated, similarly as above, and with arguments similar to the ones used in [3, Lemma 4.2], as:

[TABLE]

Thus, the estimate (12) holds true. ∎

Remark 2.7**.**

Going through the various steps in the proof, it is clear that the estimates (11) and (12) may be improved, with the respect to the degree $r$ of the parameterization, as in few places we use suboptimal Sobolev inclusions, or estimates. Indeed, as we will see in the next section, numerical validation seems to suggest that the convergence is indeed slightly (less than half a order) faster for the examples we discuss.

2.3 Element-wise re-parameterizations for trimmed domains

As it was stated above, the use of boolean operations in the definition of the computational domain $\Omega$ makes the computation of integrals over $\Omega$ far from trivial.

Integrals are computed over every single element and their contributions are added to the global result, as \int_{\Omega}\xi=\sum_{Q\in\Omega}\int_{Q{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\cap\Omega\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}}}\xi. Thus, we will distinguish the integrals over the two different families of elements $\mathcal{T}^{int}_{h}(\Omega)$ and $\mathcal{T}^{\Gamma}_{h}(\Omega)$:

[TABLE]

The integrals over the Bézier elements $Q\in\mathcal{T}^{int}_{h}(\Omega)$ (including the integrals on their boundaries) are computed, in a standard way, using a tensor-product Gauss-Legendre quadrature rule. On the other hand, the integrals over the cut Bézier elements $Q\in\mathcal{T}^{\Gamma}_{h}(\Omega)$ require a special procedure.

In this section we present a novel technique for computing integrals over the active part $\mathsf{Q}$ of cut Bézier elements. We describe our approach in 3D, but it works verbatim in 2D. Our technique has a lot in common with the integration strategy used in the Finite Cell methods (FCM), [60, 36]. Even if FCM is normally used as an immersed method, i.e., with $\mathbf{F}$ being the identity map, it has been generalized to non-trivial $\mathbf{F}$ in [47], for thin walled structures.

While a comparison of the two approaches is beyond the scope of this paper, we anticipate that our method is likely to have a lower complexity, i.e. requires less evaluations, for coarse meshes. Indeed, when the mesh is coarse, and the trimmed Bézier elements have complex shapes, we do not perform any dyadic splitting which, in principle, generates several quadrature points. On the other hand, we expect the two techniques to be comparable for fine enough meshes.

An alternative approach for a high-order re-parameterization of the trimmed Bézier elements has been recently proposed in [43] for the case of 3D V-reps in which $\mathbf{F}$ is a non-trivial map.

In order to compute the integrals over cut Bézier elements the first step is to identify them and to obtain a V-rep of each of them. This operation, denominated slicing, is detailed in Section 2.3.1. In a second stage, a high-order re-parameterization is created for every cut Bézier element. This re-parameterization will help us in building quadrature rules adapted to the active part of every element. This operation is described in Section 2.3.2.

Even if the techniques presented in this section are valid for both 2D and 3D domains, we restrict our discussion to the 3D case. In addition, and for the sake of simplicity, hereinafter we assume that the domain $\Omega_{0}$ is composed by a single trivariate $T$, whose associated spline diffeomorphism is ${\bf F}:\widehat{T}\to T$. This simplification does not constitute a limitation and does not prevent the methodology described below to be applied to cases in which $\Omega_{0}$ is a multipatch trivariate.

2.3.1 Slicing

The slicing process consists in splitting the domain $\Omega$ into Bézier elements. The goal of this operation is to identify the two families of elements $\mathcal{T}^{int}_{h}(\Omega)$ and $\mathcal{T}^{\mathsf{Q}}_{h}(\Omega):=\{\mathsf{Q}\ :\ \mathsf{Q}=Q\cap\Omega\ :\ \forall Q\in\mathcal{T}^{\Gamma}_{h}(\Omega)\}$.

The splitting is performed by recursively bisecting $\Omega$ along all the knot surfaces in the parametric domain $\widehat{T}$. We denote knot surface as the composition with ${\bf F}$ of a knot plane of the parametric domain $\widehat{T}$. A knot plane is defined as a plane perpendicular to a parametric direction $d$ at a fixed knot value $d_{k}$.

The slicing procedure followed in this work is closely related to the one described in [43, Section 3.1], being its main steps summarized in Algorithm 1. Additionally, the procedure workflow is also illustrated in Figure 1, where the knot surfaces are represented in yellow color.

We initially consider the V-rep $\omega$ ($\omega\equiv\Omega$ if it corresponds to the full domain) that has an associated B-spline tensor-product trivariate $T_{\omega}$ and a boundary representation $\gamma_{\omega}$ (a set of trimmed surfaces). $\omega$ is then split according to a knot surface perperdicular to the longest parametric direction of the domain $\widehat{T}_{\omega}$. As a result of this splitting operation, denoted as SubdvWithKnotSrf in Algorithm 1, two new V-reps are obtained, $\omega_{1}$ and $\omega_{2}$, one on each side of the surface. The new V-reps are subsequently sliced by applying the same algorithm recursively.

The slicing procedure finishes when all the Bézier elements have been identified as members of $\mathcal{T}^{int}_{h}(\omega)$ or $\mathcal{T}^{\mathsf{Q}}_{h}(\omega)$. As mentioned above, the goal of the slicing operation is to identify cut and non-cut Bézier elements. Therefore, if at a certain step the algorithm finds a V-rep that is not trimmed, i.e., $\gamma_{\omega}$ coincides with the boundary of its associated trivariate $\partial T_{\omega}$, the slicing procedure is stopped for that particular V-rep. All the Bézier trivariates contained in that particular V-rep belong to $\mathcal{T}^{int}_{h}(\omega)$ (steps 14-16 of Algorithm 1).

In this work, this procedure has been implemented using the geometric kernel Open CASCADE [45] for the splitting of the B-reps and IRIT [22] for managing trivariates.

Remark 2.8**.**

When $\mathbf{F}$ is the identity map (as in many immersed methods), knot surfaces are planes, and the slicing is a simple task. Instead, when $\mathbf{F}$ is not trivial, the slicing involve several geometric operations. In this work, all these geometrical operations, e.g. computations of surface to surface intersections, points inversion, curves approximation, etc., are performed using the geometric kernels IRIT [22] and Open CASCADE [45]. In both tools, and probably in most of the geometric kernels available, these operations are executed with limited precisions (around $10^{-8}$ at best), that are far from the maximum precision of 64-bits floating point numbers (around $10^{-16}$).

Even if these precisions are small enough for most of the industrial applications they target, they impact the asymptotic behavior of our methods, for fine enough meshes. Some of the results presented in Section 3 reflect this issue.

2.3.2 Local re-parameterizations of cut Béziers

After the slicing procedure described above we obtained a V-rep associated to the active part $\mathsf{Q}$ of every cut element. In order to be able to compute numerical integrals over the cut elements, we re-parameterize each of those V-reps with the purpose of building suitable integration formulae.

The element re-parameterization we propose in this work consists in creating a high-order hybrid polynomial mesh $\mathsf{Q}_{h}$, composed of tetrahedra and / or hexahedra. The created mesh $\mathsf{Q}_{h}$ is the union of $n_{\mathcal{L}}$ Lagrangian elements $\mathsf{Q}_{h}=\cup^{n_{\mathcal{L}}}_{i=1}\mathcal{L}_{i}$, where $\mathcal{L}_{i}\subset\mathbb{R}^{3}$ denotes a single Lagrangian element. Roughly, each element $\mathcal{L}_{i}$ is defined by a set of basis functions (Lagrange polynomials) and the nodes associated to those functions. This operation is sketched in the first step of Figure 2.

In this work, the mesh $\mathsf{Q}_{h}$ is created using Gmsh software [25]. Gmsh is able to create high-order (up to degree 10) tetrahedral meshes from geometries defined by means of B-reps. In addition, Gmsh can use internally Open CASCADE as geometric kernel for the treatment of B-rep geometries, what makes the integration of both tools seamless. Thus, in our implementation, we directly use the B-rep of every cut Bézier element created with Open CASCADE and generate a mesh $\mathsf{Q}_{h}$ with Gmsh.

In order to construct a 3D mesh, Gmsh first constructs meshes with linear elements for the edges and faces of the B-rep. The nodes of the 1D and 2D elements are created by sampling points on the surfaces and edges of the B-rep model, thus, it is guaranteed that those points lay on the B-rep surfaces. Afterwards, the mesh that fills the interior of the B-rep is created.

In a second step, Gmsh performs a degree elevation of the tetrahedral mesh elements by introducing new nodes. As before, the new nodes created in this process, belonging to external faces of the mesh, are sampled in the spline surfaces of the B-rep. Hence, the boundary $\gamma$ is approximated with a high-order polynomial representation $\gamma_{h}$ (yellow surfaces in Figure 2). Finally, as it is described in [31, 57], the mesh quality is optimized, with an algorithm that guarantees that all vertices lying on the B-rep surfaces will still belong to those surfaces. As our mesh is only used to place quadrature points, our lonely quality criteria is the absence of negative Jacobians: in practice, these relaxed quality requirements makes the task of creating high-order meshes much simpler that in finite element methods. In addition, for the sake of computational efficiency, our aim is to create the coarsest possible valid mesh.

The approximation $\mathsf{Q}_{h}$ obtained with this approach complies with all the requirements stated in Assumption 2.4, upon which are obtained the estimate results presented in Section 2.2. Indeed, by construction, the boundary mesh $\gamma_{h}=\partial\mathsf{Q}_{h}$ is a high-order approximation of $\gamma=\partial\mathsf{Q}$. It may happen though that $\gamma_{h}\not\subset Q$: in this case, a dyadic subdivision of $\mathsf{Q}$ is performed and the meshing process is restarted on each concerned sub-element.

In a standard way, for calculating the operators involved in the problem (10), their integrals are transformed from the physical to the parametric domain as:

[TABLE]

Thus, once a valid mesh is created with Gmsh, we have a hybrid tetrahedral / hexahedral mesh of degree $r$, but defined in the physical space (green mesh in Figure 2). In order to use that mesh for building quadrature rules, we need to pull it back to the parametric domain $\widehat{T}$.

Through a point inversion algorithm, all nodes of the mesh are pulled back to the parametric domain $\widehat{T}$ and the corresponding mesh $\widehat{\mathsf{Q}}_{h}$, pre-image of $\mathsf{Q}_{h}$, is generated. If needed, the Gmsh optimization referenced above is also applied to $\widehat{\mathsf{Q}}_{h}$ in order to avoid negative Jacobians.

Once the mesh $\widehat{\mathsf{Q}}_{h}$ is created, in order to compute integrals over a single cut Bézier element we proceed as follows: we place a suitable quadrature rule (Gauss Legendre rule for hexahedra and collapsed Gauss Legendre for tetrahedra [39]) defined over the reference element (e.g., the reference tetrahedron in Figure 3); the quadrature coordinates and weights are transformed by pushing them forward to all the elements in $\widehat{\mathsf{Q}}_{h}$ in $\widehat{T}$ (right part of Figure 3); using those transformed quadratures, the contributions of all the elements are added to the global result.

Remark 2.9**.**

In order to avoid the creation of the mesh $\widehat{\mathsf{Q}}_{h}$ from $\mathsf{Q}_{h}$, an alternative procedure would be to set the quadrature points in the elements of the physical mesh $\mathsf{Q}_{h}$ and then transform their coordinates and weights back to the parametric domain $\widehat{T}$ by applying the inverse transformation ${\bf F}^{-1}$. This procedure presents the advantage of bypassing the creation of the mesh $\widehat{\mathsf{Q}}_{h}$, however, it is dependent of the integration formula chosen and must be repeated for different quadrature rules.

Remark 2.10**.**

The choice of the quadrature rule we made is likely not optimal, and, more generally, the optimisation of the quadrature is not the main focus of this contribution and deserves further studies.

Remark 2.11**.**

Even if the procedure presented above is applicable also to the 2D domains, in that case is far from optimal. The first operation, the slicing, is common for both 2D and 3D cases, however, the way in which local re-parameterization of the trimmed elements is created can be greatly simplified in the 2D case. See, for instance, [40, 35, 44]. The procedure used in this work for the 2D experiments presented in Section 3 is based on the use of patterns as described in [41, 54, 34]. When a cut Bézier element does not correspond to any pattern, we use the methodology presented above.

2.4 Linear system preconditioning

The error estimates obtained in Section 2.2 state that, under certain assumptions, guaranteed by the constructions presented in Section 2.3, the accuracy and optimality of the isogeometric method are preserved.

Nevertheless, due to the fact that the support of some basis functions is cut, the linear system can be potentially ill-conditioned. This is due to the fact that by reducing the active support of some basis functions, we are also reducing the minimum eigenvalues of the linear system matrix, while the maximum eigenvalue is not affected.

A detailed study of conditioning problems that derive from the use of trimmed domains is beyond the scope of this work, but we include here a brief discussion of the preconditioner used in the numerical examples presented in Section 3.

In this work we use a diagonal scaling preconditioner. Thus, the original linear system can be written as $\bm{\mathsf{A}}\bm{\mathsf{x}}=\bm{\mathsf{b}}$, where $\bm{\mathsf{A}}\in\mathbb{R}^{n\times n}$ is the stiffness matrix, $\bm{\mathsf{x}},\,\bm{\mathsf{b}}\in\mathbb{R}^{n}$ are the solution and right-hand-side vectors, respectively, and $n$ is the system size. By symmetrically preconditioning the system it becomes:

[TABLE]

where $\bm{\mathsf{D}}$ is the diagonal scaling preconditioner defined as

[TABLE]

where $A_{i,i}$, with $i=1,\dots,n$, are the diagonal entries of the matrix $\bm{\mathsf{A}}$. By applying this algebraic preconditioner, all the basis functions are renormalized respect to their contribution to the matrix $\bm{\mathsf{A}}$.

This preconditioner’s effectiveness is illustrated, for a series of elliptic problems, in the numerical experiments gathered in Section 3. Although in the context of the Finite Cell Method [21], this simple choice seems insufficient to restore a stable behavior of the condition number, i.e., a condition number independent of the way elements are cut, in our numerical tests we consistently experience a very good behavior of the preconditioned system. We believe this is due to the use, in our tests, of maximum continuity splines and we defer the reader to [13] for further discussions.

Possible alternatives to the diagonal scaling preconditioning are the ghost penalty techniques [14, 15]; the use of additive Schwarz preconditioners as proposed in [19, 32] for finite cell and isogeometric discretizations, and extended in [20] to case of immersed IGA combined with multigrid methods; or the stable removal of basis functions, as proposed in [24] in the context of CutIGA methods.

3 Numerical experiments

In this section we perform a series of numerical experiments that aim at illustrating the methodology presented in previous sections.

As a first example, we present in Section 3.1 a Poisson problem in a 2D circular domain which will help us in supporting the theoretical results presented in Section 2. Secondly, the well known plate with hole test for 2D linear elasticity is presented (Section 3.2). A Poisson problem is used again in Section 3.3 to validate our theoretical findings, but this time on a three-dimensional computational domain: a sphere, which is trimmed from a cube. Finally, the presented methodology is applied to study the elastic behavior of a 3D mold with a trimmed interior cooling channel, which has a complex geometry.

The results shown in this section have been obtained using an implementation of the methodology proposed above in an in-house code built upon the isogeometric analysis library Igatools [46], the solid modeling kernels IRIT [22] and OpenCASCADE [45] and the 3D mesh generator Gmsh [25].

3.1 Poisson equation in a 2D domain

In this Section we study a Poisson problem in the interior of a circular domain $\Omega$:

[TABLE]

The solution to the problem is

[TABLE]

and the associated source and Neumann terms are:

[TABLE]

No Dirichlet condition is imposed, alternatively we impose weakly the condition

[TABLE]

in order to remove the constant part of $u$.

3.1.1 Poisson equation in a 2D non-distorted domain

The computational domain is built as the boolean intersection of a $L$-length square domain and a circle, with radius $R$, centered at the origin (c.f. Figure 4(a)). The squared domain carries a Cartesian mesh. As a matter of example, in Figure 4(b) the re-parameterization of trimmed elements (in red) is shown together with the non-trimmed ones for a domain with $8\times 8$ elements.

For this re-parameterization, as well as for all the numerical results gathered in this section, $R=1$ and $L=2\,R/0.7$ were used. As the image of each knot line is straight, the re-parameterization of the trimmed elements is performed using a numerical precision $10^{-15}$, it is of the same degree as the one used for solution, and is the coarsest possible.

As a first result, the solution error in the $L^{2}$ and $H^{1}$ norms for different mesh sizes and degrees is shown in Figures 5(a) and 5(b), respectively, with the choice $p=r$, in agreement to Theoreom 2.5. For $p=4,5$ one quadrature point per direction is added in each element of the re-parameterization, and for $p=6$, two points per direction are added in order to recover the expected convergence rates.

Moreover, 5(c) and 5(d) show the error on the computation of area and perimeter which are bounded from below by the chosen geometric precision. Note that the limited accuracy of the re-parameterization for $p=6$ and $h\leq L/64$ (c.f. Figure 5(d)) has an effect in the computed $L^{2}$ error norm for $p=6$ and $h=L/128$. Finally, from 5(c) and 5(d), the area and perimeter errors converge as $h^{p+1}$ and $h^{p+2}$ for odd and even degrees, respectively. This fact is due to the use of equidistant interpolation points for the approximation of the boundary, and the known fact that quadrature formulae on equidistant odd points are always even degree.

In order to illustrate the necessity of high-order re-parameterizations for the trimmed elements, in Figure 6 we show the effect of using re-parameterization degrees $p_{t}$ lower than the discretization degree $p$. In this figure we split the contribution to the $H^{1}$ error norm in two parts: the contribution of the trimmed elements $\mathcal{T}^{\Gamma}_{\Omega}$ (solid lines) and of the non-trimmed ones $\mathcal{T}^{int}_{\Omega}$ (dashed lines).

As it can be seen, when $p_{t}<p$ and $h$ becomes small the consistency error starts to be dominant respect to the discretization error and consequently the error optimality is lost. For higher $p_{t}$, this phenomenon alleviates and the sub-optimality shows up on finer and finer meshes.

On the other hand, in Figure 7, despite the solution degree $p$, the trimmed elements are re-parameterized with a low-order approximation ($p_{t}=1$). Nevertheless, instead of creating the coarsest possible re-parameterization, the trimmed elements are approximated using more tiles. $h_{t}$ is the size of the tiles and $h^{-1}_{t}$ measures the number of tiles used along each side of every trimmed element. The use of a higher number of tiles has the effect of reducing the consistency error, thus its sub-optimality affects the total error for finer meshes.

As anticipated in Remark 2.7, the convergence rates we see in Figures 6 and 7 are slightly better than the ones proved true in Theorem 2.5.

Finally, in Figure 8 the stiffness matrices’ condition numbers are reported for the preconditioned and non-preconditioned cases. As it can be seen, the non-preconditioned case (Figure 8(a)) presents very high condition numbers for all degrees. This is controlled, as it can be appreciated in Figure 8(b), using the diagonal preconditioning described in Section 2.4. Additionally, the expected $h^{-2}$ asymptotic behavior of the conditioning is recovered in the latter case.

3.1.2 Poisson equation in a 2D distorted domain

The same Poisson problem is solved but using a different parameterization of the computational domain. The computational domain is still the interior of a circle, however the parameterization does not longer correspond to a Cartesian mesh, but to a distorted one. In Figure 9 the trimmed parametric domain is shown in the left side and the physical distorted domain in the right side. As described in Section 2.3, this time the geometric operations are performed approximately with a numerical precision of $10^{-8}$ that is far from the maximum precision of 64-bits floating point numbers used in the calculations (around $10^{-16}$).

The aim of this section is to show the impact of this geometry approximation on the accuracy of the numerical method.

In Figure 10 the solution error in the $L^{2}$ and $H^{1}$ norms and the absolute error of the area and external perimeter are shown for the distorted domain case. It can be seen that both the area and perimeter errors, that have a noisy behavior due to the distortion of the parameterization, are bounded from below by the geometric precision. On the other hand, the solution $L^{2}$ and $H^{1}$ errors (c.f. Figures 10(a) and 10(b)) show optimal converges properties until they reach the geometric precision (which bounds them from below).

In order to confirm the interpretation of our numerical results, for this specific (simple) case, we were able to compute an approximation of the curve $\partial\widehat{\Omega}$ with a precision under $10^{-12}$, by enforcing manually the precision of the geometric algorithms. The obtained results are gathered in Figure 11. All the errors present the same behavior as in the previous case, but solution errors are not longer bounded by the geometric precision.

3.2 Infinite plate with circular hole in a 2D domain

In this section the infinite plate with a circular hole under tension is studied, a well-known case in the isogeometric literature (e.g. [29, 5]).

Considering plain strain conditions for the problem, the analytic solution, expressed in polar coordinates, is:

[TABLE]

The involved quantities and the problem boundary conditions are described in Figure 12(a). Dirichlet boundary conditions are on the non-trimmed part of the boundary and are imposed strongly.

We refer the reader to [56] for further details. The computational domain is built as the boolean difference of a $L$-length square domain minus a circle with radius $R$.

In the Figure 13 the solution errors in the $H^{1}$ norm for different mesh sizes and degrees are shown, together with the absolute errors of the area computation. In these results, the trimmed elements are re-parameterized with the same degree of the solution discretization and using the coarsest possible re-parameterization (see Figure 12(b) as a matter of example). Thanks to the underlying Cartesian mesh, the geometric precision is down to $10^{-15}$. The solution errors converge optimally for all the degrees, and it is also the case for the area error, that presents the same odd/even behavior described in previous section.

Figure 14 reports the condition number of the stiffness matrices for the preconditioned (using the diagonal scaling preconditioner presented in Section 2.4) and non-preconditioned cases. As it was already seen in the Poisson example presented in Section  3.1, the diagonal scaling preconditioner keeps the matrix conditioning under control for all considered degrees.

3.3 Poisson equation in a 3D domain

A Poisson problem in the interior of a 3D sphere is studied in this section:

[TABLE]

The problem analytic solution, in spherical coordinates, is:

[TABLE]

where $\theta$ is the azimuth and $\varphi$ the co-latitude. The laplacian of $u$ is

[TABLE]

where the Neumann condition on the sphere boundary is $\nabla u\cdot\bm{n}=g_{n}=0$. As in Section 3.1, we impose weakly the condition

[TABLE]

The computational domain is built as the boolean intersection of a unit cube domain and a sphere with radius $R=1$.

As it can be seen in Figure 15, optimal approximation properties are achieved both in $L^{2}$ and $H^{1}$ norms for different mesh sizes and degrees. The computation of the domain volume is also optimal and presents the same odd/even behavior observed before.

Finally, and as in previous numerical experiments (Sections 3.2 and 3.1), the matrix conditioning, with and without preconditioning, is analyzed in Figure 16. As before, the diagonal scaling preconditioner reduced drastically the matrix’s condition number, allowing to recover the expected asymptotic behavior $h^{-2}$.

3.4 A 3D mold with a cooling channel

Finally, in order to illustrate the capabilities of the presented method for general geometries in this last numerical example we study the elastic response of a 3D printed mold. The mold presents an interior cooling channel whose mission is to reduce its temperature by means of a circulating cold fluid.

As it can be seen in Figure 17, the cooling channel follows a complex path in the interior of the mold in order to successfully refrigerate it: in enters from the interior cylindrical face, travels throw the interior and exterior faces of the mold, and exists from the exterior cylindrical face. In addition, some small cylindrical holes are also present in the top face of the mold (2 in the interior wall and 4 in the exterior one).

The mold geometry is described with a single spline trivariate whereas the cooling channel has been created as the extrusion of a circle along the channel path. The final geometry is built by means of a boolean operation: the mold minus the cooling tube. In Figure 18 the non-trimmed elements (in green) together with the trimmed ones (in blue) are shown.

The trimmed Bézier elements are re-parameterized with cubic tetrahedra. As a matter of example, some re-parameterized elements are shown in Figure 18(b).

The trivariate has degree 3 along the longitudinal direction and degree 2 in the other two directions ($56\times 6\times 36$ elements). The mechanical properties of the material mold are $E=10\,\text{GPa}$ and $\nu=0.2$. These values are similar to the material properties of a metal specimen manufactured with an additive manufacturing procedure: they were inferred from [26, Table 3].

The performed simulation consists in a linear elasticity analysis in which a traction force (Neumann condition) is applied on the curved face of the mold. The internal top face is fully fixed encastre while the exterior top face is fixed just vertically. The value of the traction is $1.0\,\text{MPa}$ along $z$-direction (pointing downwards) and $0.5\,\text{MPa}$ along $x$-direction.

The deformed mold together with the stress magnitude can be seen in Figure 19(a). The discontinuity that can be appreciated in the stress field near the corners is due to the fact that geometry parameterization is only $C^{0}$ continuous at those points. Additionally, a view of the stress in the interior of the mold is shown in Figure 19(b). As it can be seen, despite the fact that the mesh size is similar to the channel diameter, the stress distribution is affected by the presence of the channel.

4 Conclusions

We have presented a novel approach for the construction of numerical methods for elliptic PDEs on trimmed geometries, seen as a special class of more general V-reps. More specifically, this method is applicable to the case of trimmed non-conforming multipatch trivariate volumes.

Our approach is based on the local re-parameterization of Bézier elements that are cut by the trimming, in contrast to the creation of a global re-parameterization of the full domain. Thus, by performing a re-parameterization of the trimmed elements with the same degree as the discrete solution, we proved that our method guarantees optimal convergence rates when possibly non-homogeneous natural (Neumann) boundary conditions are imposed on the part of the boundary affected by trimming.

This theoretical result is supported by numerical experiments that illustrate how the use of low-order re-parameterizations lead inevitably to a lack of optimal convergence. The use of finer, but still low-order element re-parameterizations, mitigates this problem for coarse solution mesh sizes $h$. But, for sufficiently small values of $h$, fine low-order element re-parameterizations yield a sub-optimal behavior.

Our numerical experiments also illustrated how the limited precision values inherent to the available geometric modeling tools limit the potential accuracy achievable by our numerical methods.

In all the considered experiments the use of a diagonal scaling preconditioner proved itself efficient in controlling the ill-conditioning derived from the trimming of the basis functions’ support.

Regarding the re-parameterization of the trimmed elements, in this work we proposed a novel methodology for building adapted quadrature schemes based on the creation of coarse high-order tetrahedral and / or hexahedral finite element meshes. The only requirement to these meshes is the absence of sign change in tetrahedra / hexahedra Jacobians. The potential of this methodology was in the local re-parameterization of a complex V-rep.

Finally, we would like to remark that the imposition of essential (Dirichlet) boundary conditions on trimmed boundaries ask for stabilization: in this case our approach should be combined with the stabilization proposed and analyzed in [13] and this is object of further studies. The extension of the presented method to the case of domains constructed as the union of V-reps will be addressed in the future.

Acknowledgements

Pablo Antolin and Annalisa Buffa gratefully acknowledge the support of the European Research Council, through the ERC AdG n. 694515 - CHANGE. Massimiliano Martinelli has been supported by the European Union’s Horizon 2020 research and innovation programme under grant agreement n. 680448 - CAxMan. The authors also acknowledge R. Vázquez for his helpful insights and suggestions, and would like to thank G. Elber and F. Massarwi for the many inspiring discussions on volumetric modeling.