IETI - Isogeometric Tearing and Interconnecting.

Finite Element Tearing and Interconnecting (FETI) methods are a powerful approach to designing solvers for large-scale problems in computational mechanics. The numerical simulation problem is subdivided into a number of independent sub-problems, which are then coupled in appropriate ways. NURBS- (Non-Uniform Rational B-spline) based isogeometric analysis (IGA) applied to complex geometries requires to represent the computational domain as a collection of several NURBS geometries. Since there is a natural decomposition of the computational domain into several subdomains, NURBS-based IGA is particularly well suited for using FETI methods. This paper proposes the new IsogEometric Tearing and Interconnecting (IETI) method, which combines the advanced solver design of FETI with the exact geometry representation of IGA. We describe the IETI framework for two classes of simple model problems (Poisson and linearized elasticity) and discuss the coupling of the subdomains along interfaces (both for matching interfaces and for interfaces with T-joints, i.e. hanging nodes). Special attention is paid to the construction of a suitable preconditioner for the iterative linear solver used for the interface problem. We report several computational experiments to demonstrate the performance of the proposed IETI method.

Introduction

Isogeometric analysis (IGA), introduced by Hughes et al. [22], is a promising concept that establishes a close link between the technologies of CAD (computer aided design) and numerical simulation via finite element analysis (FEA). In the IGA framework, the same function spaces, which are used for the geometric representation of the computational domain, are used for the approximation (e.g. Galerkin) of the problem unknowns. There are several computational geometry technologies that could serve as a basis for IGA. However, Non-Uniform Rational B-Splines (NURBS) are the most widely used and well established computational technique in CAD, which we shall pursue here. The main advantage consists in the fact that an exact representation of the computational domain is available at all times during the (possibly adaptive) simulation process. The need for geometric approximation processes, known as mesh generation, is thus eliminated.

Despite its short history, IGA has already attracted enormous interest which is documented by a substantial number of publications. IGA has been successfully applied to various simulation problems, including linear elasticity, simulation in electromagnetics, and flow simulations, e.g. [4,5,8,16,41]. It was shown that IGA not only eliminates the geometrical errors introduced by the approximation of the physical domain by a computational mesh, but also significantly reduces the number of unknowns (when compared to standard FEA) needed to achieve a certain accuracy in practice [10].

On the theoretical side, various issues such as error estimates, convergence rates, stability issues, and numerical quadrature rules have been analyzed thoroughly [2,6,23,50]. Many of these topics are already well understood and it was demonstrated that most of the theoretical background of FEA can be adapted and extended to cover the IGA framework. A notable exception, which is the subject of several on-going studies, is the search for generalizations of tensor-product spline spaces providing the possibility of local refinement. This is needed for adaptive simulation methods. Several concepts are currently being explored, including T-splines, hierarchical splines, THB-splines, and polynomial splines over hierarchical T-meshes [3,7,14,20,30,33,40,46,54].

Another new problem, which is located right at the crossroads between numerical simulation and geometric design, is the challenge of domain parameterization for IGA: given a boundary representation of the computational domain by a collection of NURBS patches, extend it to a volume parameterization that is suitable for IGA. Several methods for designing such parameterizations and for creating them from boundary data exist [1,9,21,39,55]. For the numerical simulation in real-world applications it will be essential to create and to use multi-patch parameterizations [38], where the computational domain is represented by several NURBS volumes. These structures, which combine a highly regular structure on subdomains with topological flexibility to describe complex geometries, open new perspectives for the design of solvers in numerical simulation, which will be explored in this paper.

To do so, we will extend the finite element tearing and interconnecting (FETI) method to the isogeometric context. FETI methods are powerful solvers for large-scale finite element systems in computational mechanics. They were introduced by Farhat and Roux [17] and belong to the class of iterative substructuring methods (also called non-overlapping domain decomposition methods), see e.g. [52]. The computational domain is subdivided into non-overlapping subdomains, and each subdomain gets its own set of equations derived from the global equation. To ensure the equivalence to the global equation, continuity conditions are introduced at the interfaces between subdomains using Lagrange multipliers.

In this paper, we combine the ideas of IGA and FETI and propose a new method called Isogeometric Tearing and Interconnecting. Following the nomenclature of FETI and BETI (Boundary Element Tearing and Interconnecting) methods, we abbreviate the proposed method as IETI (to be pronounced [′jεti], like Yeti). The FETI method, however, is not only a coupling method but provides a powerful solver design. By (carefully) eliminating the original variables from the resulting saddle point problem, one obtains a system in the Lagrange multipliers (i.e., only on the interface). The solution of the original problem can be easily computed from the solution of this interface problem. For generalizations to boundary element discretizations see e.g. [31,32]. For generalizations to spectral element discretizations see e.g. [25].

Since the number of Lagrange multipliers is typically large, the interface problem is usually solved iteratively by FETI preconditioned conjugate gradients [45]. Suitable preconditioners have been proposed in [17,27,28]. The analyses in [28,36] show that under suitable conditions the condition number of the preconditioned system is bounded by , where H is the subdomain diameter and h is the mesh size. This results in quasi-optimal complexity of the overall method. We mention that the classical FETI method involves the solution of a coarse system. An efficient alternative is the dual-primal FETI (FETI-DP) method, see [18,26,29,37], which is followed in this paper.

The remainder of this paper is organized as follows: in Section 2, we define the model problems considered herein and recall the definition of NURBS geometry mappings. In Section 3, we describe multi-patch NURBS mappings and multi-patch NURBS discretizations. In Section 4, we formulate the IETI method and investigate the solver and preconditioner aspects. Local refinement options which are introduced by the IETI method are described in Section 5. In Section 6, the performance of the IETI method is demonstrated in numerical examples. Finally, a conclusion is drawn in Section 7.

Preliminaries

We first briefly define the model problems considered in this paper, and then recall the definitions of B-spline basis functions and NURBS geometry mappings.

The model problems

Let be an open, bounded and connected Lipschitz domain with boundary . Let the boundary be composed of two disjoint sets, namely , where essential (Dirichlet or displacement) boundary conditions are prescribed, and , where natural (Neumann or traction) boundary conditions are prescribed. Let denote the outer unit normal vector to , and let and be given functions. We shall consider two model problems:(I) The scalar diffusion problem:

Find the scalar function such that(II) The linearized elasticity problem:

Find the displacement field such thatIn problem (I), denotes the diffusion coefficient. In problem (II), , where C is the fourth-order stiffness tensor and is the linearized strain tensor.

The variational forms of these problems can be written as follows:

Find such thatThe space contains test functions which vanish on . In problem (I), we have , and in problem (II), . The set consists of functions which fulfill the essential boundary conditions on .

We assume that the problem data is such that the bilinear form is bounded, symmetric and positive definite, and that is a bounded linear functional. Then, the problem (1) can be reformulated as the following minimization problem:

Find u such that

NURBS geometry mappings

To provide an overview, and to introduce notations, we briefly recall the definition of B-spline basis functions and NURBS mappings. Let p be a non-negative degree and let be a knot vector with for all k. We consider only open knot vectors, i.e. knot vectors s where the multiplicity of a knot is at most p, except for the first and last knot which we require to have multiplicity . For simplicity, we assume that and , which can be easily achieved by a suitable scaling. The univariate B-spline basis functions , are defined recursively as follows:Whenever a zero denominator appears in the definition above, the corresponding function is zero, and the whole term is considered to be zero. For open knot vectors, the first and last basis function are interpolatory at the first and the last knot, respectively.

Let and be two families of B-spline basis functions defined by the degrees p and q, and the open knot vectorsrespectively. We denote the set of all double-indices of NURBS basis functions byLet , be positive weights. We define the bivariate NURBS basis functions as follows:where the numerator and the denominator are given byGiven a control net of control points , where again , the two-dimensional NURBS-surface is defined bywhere . We refer to Q as the parameter domain and to as the physical domain.

For a detailed discussion of NURBS mappings and B-spline basis functions, and the properties of these functions we refer e.g. to [3,10,11,22] and the references therein. Note that NURBS basis functions of degree p are, in general, globally -continuous. However, in the presence of multiple knots, the continuity reduces according to the multiplicity. If a knot appears i times, the continuity of a NURBS basis function of degree p at that knot is .

For our purposes, we assume that the geometry mapping is continuous and bijective (i.e. not self-penetrating), which are natural assumptions for CAD-applications.

For better readability, we collapse the double index into one multi-index in the following. Hence, instead of (3), we write

Multi-patch geometry mappings

Single-patch NURBS discretization

In isogeometric analysis, the NURBS basis functions of the geometry mapping are also used to represent the discrete solution. The scalar -conforming finite element space in a single-patch setting is the following:whereFor problems of type (II), we need to use the corresponding vector-valued space . For simplicity, the following notation is only oriented on the scalar case.

We assume that the prescribed Dirichlet data is such that there exists (otherwise, we project the Dirichlet data to ). We defineA function is represented asThe real-valued coefficients are referred to as control variables or degrees of freedom (DOF).

Multi-patch NURBS discretization

Assume that the physical domain is represented by N single-patch NURBS geometry mappings , each of which maps the parameter domain to an open physical subdomainsuch thatWe use the superscript (i) to indicate that knot vectors, degrees, NURBS basis functions, index sets, DOF, etc. are associated with a mapping . For example, we denote the set of NURBS basis functions of the geometry mapping by . For each subdomain , the local function space (for ) is defined analogously to (4) and (5) aswhereWe denote the space of functions that are locally in byAs mentioned in Section 2.1, we have for problem (I), and for problem (II). Functions in are not necessarily continuous. We choose the following subsets of continuous functions in :A function is represented subdomain-wise by:We denote the interface of two subdomains and (see Fig. 1) byWe collect the index-tupels of all interfaces that are non-emptywith . For , we call a subdomain vertex (or vertex for brevity) if it consists of a single point, otherwise we call it an edge. For , we collect the indices of those basis functions in whose support intersects the interface :If , we say that the DOF is associated with the interface . In case of a subdomain vertex, we have , where indicates the cardinality.

Fig. 1

Fully matching subdomains and : All weights equal to 1, , , , .

Let be an edge. We say that the subdomains and are fully matching, if the following two conditions are fulfilled (see Fig. 1 for an illustration):For an illustration of two fully matching subdomains, see Fig. 1: The interface is the image of the entire northern edge of under the mapping , and the image of the entire western edge of under . Furthermore, . The knot vectors and are not equal, but due to the way they are mapped to , condition (ii) is fulfilled. Hence, and are fully matching.

The interface is the image of an entire edge of the respective parameter domains.

For each index , there must be a unique index , such thatThis is the case, if the two knot vectors are affinely related and the corresponding weights and degrees are equal.

The tensor-product structure of the NURBS basis functions is very convenient for collecting/identifying the DOF associated with an interface, i.e. the index-set . In particular, in combination with condition (i) for the fully matching case, one only has to know which side (north, east, south or west) of the parameter domain defines the interface to identify the associated DOF. For example, in Fig. 1, using the double-index notation, we haveAssume that the linear form can be assembled from contributions on , and let denote the restriction of to . Then, we can discretize the original minimization problem (2) as a sum of local contributions:

Find such thatCondition (ii) for the fully matching setting implies that each DOF , can be identified with a DOF , such that the corresponding basis functions match as in (12). If we identify the DOF corresponding to these matching basis functions, then, together with the remaining DOF, we get a representation of the space in (8). Employing a suitable numbering for the identified DOF, we can rewrite (13) in the formwhere is the coefficient vector corresponding to the DOF which are not on is the stiffness matrix and is the load vector. This system can be solved by standard sparse LU-factorization [12]. However, it is well known that for large problem size, the memory requirement and the runtime complexity of direct solvers are inefficient. Alternatively, one can use efficient iterative solvers such as conjugate gradient methods with appropriate preconditioners [45]. In the case of standard FEM discretizations, such preconditioners have been well studied in the literature, see e.g. [53] for geometric and algebraic multigrid methods, and see [52] for domain decomposition methods.

It is important to note that by assembling the system matrix from the subdomain contributions, the structural f(subdomain-wise) properties of the problem are lost, which are hard to regain from alone. To alleviate this difficulty, in the following section we present a solver (inspired by the FETI methods [15,17-19,42,52]) which inherently uses the local structure of (13). This approach is very suitable for parallelization, and moreover, since it mainly uses solvers on the local subdomains, their tensor product structure can be exploited (e.g. using wavelets or FFT).

Solver design

In the IETI method, we work in the spaceas defined in (7). Since these functions are, in general, discontinuous across subdomain interfaces, we need to impose the continuity conditions separately. In the following, let denote the ith component of a function .

Continuity constraints

Note that, in the index set , every interface is represented twice. Therefore, we define the setwhere each interface is represented only once. For each pair , we say that is the master subdomain and the slave subdomain. For all and a fixed index , we can rewrite (12) in the following form:where, for a fixed , all coefficients are zero, except for one coefficient which is 1. Note that the generalization of (12) given in (16) might seem superfluous in the fully matching setting, but it will be needed in Section 5 in the context of cases which are not fully matching. There, we will also adapt the definition of .

We can collect the coefficients in the permutation matrixIn the case of a subdomain vertex, where , the matrix has only one entry which is 1. Fig. 2(a) illustrates the coupling conditions between subdomains and at a subdomain vertex.

Fig. 2

Illustration of fully redundant coupling and all floating setting. Arrows indicate coupling conditions and point from master subdomain to slave subdomain.

To guarantee the continuity of across all interfaces, we impose the following condition on the DOF associated with interfaces, i.e. for all :The incorporation of essential boundary conditions is done in a similar way, see the all-floating BETI method [42] and the total FETI method [15]. We denote the interface between and bySimilar to (11), we collect the DOF of that are associated with , i.e. with the interface , in the following index set:Recall from Section 3.1 that we assume that there exists a function with . Hence, we can writewith real-valued coefficients . We incorporate essential boundary conditions by imposing the following constraints on the DOF that are associated with , i.e. for all i such that :These constraints can be thought of as continuity between the physical subdomains and a virtual neighbour subdomain (see Fig. 2(b) for an illustration).

Let J denote the total number of constraints of the form (18) and (19):We assume a fixed numbering of these constraints and introduce the following notation. For a vector denotes the component corresponding to the constraint (18), and denotes the component corresponding to the constraint (19). We define the jump operator B as follows:Hence, the conditions for -continuity (18) and the incorporation of essential boundary conditions (19) readwhere the entry of the vector is zero when corresponding to an interface condition (20), and it equals when corresponding to an essential boundary condition (21). For with , we have . Note that the linear operator B can be represented by a signed Boolean matrix. With (22), we obtain the following restricted minimization problem which is equivalent to (13):

Find u such that

Throughout this paper, we assume that the computational domain is represented as a collection of several patches which are joined with -smoothness along their interfaces. These patches simultaneously serve as the subdomains for the IETI method. The approach can be extended to unstructured meshes, such as T-spline representations [3,46,47], which are also frequently used in isogeometric analysis. Similar to the case of classical FETI, a computational domain represented by a T-spline mesh can be split into subdomains. The coupling of the DOF across the interfaces needs to take all test functions into account that do not vanish on the interface. Typically, these DOF form a strip whose width increases with the degree of smoothness of the T-spline representation. This is different both from the classical FETI method for piecewise linear finite elements and from the multi-patch NURBS discretization with -continuity, where these DOF are arranged along lines. Consequently, when extending the framework to unstructured meshes in IGA, a larger number of Lagrange multipliers will be needed as compared with the case. Moreover, for general T-spline meshes, it is no longer possible to benefit from the simple tensor-product structure of the subdomains.

Saddle point formulation

Using the local basis of each subdomain space , each function is uniquely represented by a vector . Correspondingly, each function has a representation as a vector of the formwhose blocks are the local vectors . Let denote the stiffness matrix corresponding to the local bilinear form and defineAnalogously, let denote the load vector whose blocks correspond to the local subdomain load vectors, and let be the matrix representation of the jump operator B.

The minimization problem (23) is then equivalent to the following saddle point problem:

Find with the vector representation as in (24) and Lagrange multipliers , such thatWe note that even though is only unique up to an element from , the solution is unique.

The common strategy of FETI-type methods is to reduce (25) to an equation that involves only . This is not straightforward, since in the case of our model problems (see Section 2.1) the local matrices are not invertible. More precisely, in the scalar elliptic case (I) the kernel is spanned by the constant functions, and for the two-dimensional linearized elasticity problem (II), the kernel is spanned by three rigid body modes.

In the classical FETI methods [17] and the total FETI method [15], additional unknowns are introduced that span the kernel. In this paper, however, we will follow the dual-primal FETI (FETI-DP) method, see [18] and [52, Section 6.4], thus obtaining IETI-DP.

Dual-primal formulation

Recall that we only consider NURBS geometry mappings with open knot vectors. Therefore, at every vertex of the parameter domain Q, there is exactly one index , such that at this vertex is 1, while all other basis functions are zero. Hence, we can distinguish DOF that are associated with the vertices of the parameter domain. Such DOF that are associated with the vertices of a parameter domain are called primal DOF. All other DOF are referred to as non-primal or remaining DOF.

We define the following subset of To achieve the continuity above, we identify all the primal DOF that are associated with a common point in the physical domain, and we fix a global numbering of these primal DOF. Then, a function can be uniquely represented by a vector of the formwhere the subscripts P and R refer to primal and remaining DOF, respectively. Note that in (24) and in (26) are different vector representations of the same function .

Let denote the jump operator on defined bywhere is the representation of u in the form of (26). Analogously to (26), we can write

Setting up the global system

Since the solution , we can replace the space in (25) by . In the following, we will derive an equivalent saddle point formulation.

By rearranging the local DOF in such a way that the primal DOF come first, the subdomain stiffness matrix and the load vector take the formFrom these local contributions, we obtain the global stiffness matrix and the global load vector Note that, due to the identification of the primal DOF, the components carrying a subscript P in (28) are assembled from local contributions. and have the formNote that is not block-diagonal.

The coupling conditions of the form (20) can be neglected at the primal DOF, but it has no effect on the algorithm if they are included. Depending on the implementation of the jump operator, one might decide to keep them or not.

With the same steps as before, we arrive at the following saddle point problem which is equivalent to (25):

Find , represented by as in (26), and Lagrange multipliers , such thatNote that the matrix in (29) is singular, because the space has no restrictions on for our model problems. In the next section, we will incorporate essential boundary conditions at those primal DOF that are associated with . This incorporation will lead to a non-singular matrix.

Essential boundary conditions

We distinguish between essential primal DOF associated with , and floating primal DOF that are in the interior of or associated with . We indicate this with the subscripts d and f, respectively. We assume for simplicity that in the vector the essential primal DOF are listed first, i.e.andLet be the vector whose entries are the values of at the essential primal DOF. Since , the saddle point problem (29) is equivalent to the following problem:

Find , represented by as in (30), and Lagrange multipliers , such thatwhereand is the identity matrix. We see that each entry of corresponding to a multiplier acting on an essential primal DOF vanishes. Also, these multipliers are superfluous as they do not influence the solution and can be left out completely.

As before, we denote the components of , , and in (32)–(34) that correspond to primal and remaining DOF with the subscripts P and R, respectively. Hence,Again, the coupling conditions of the form (20) and (21) can be neglected at all the essential primal DOF, but it has no effect on the algorithm if they remain.

Similar to the construction (32)–(35), we can also directly incorporate the essential boundary conditions at the remaining non-primal essential DOF. If this is done, the corresponding multipliers can again be left out. This approach would be closer to the original FETI-DP method as proposed in [18].

Dual problem

In our model problems the matrix in (32) is invertible, and the first line of (31) yieldsInserting this identity into the second line of (31), we obtain the dual problemwith and . To realize the application of , we use the block factorizationwhereRecall that is block diagonal. Hence, applying corresponds to solving local problems independently on each subdomain, e.g. by sparse LU factorization [12]. Note that appears three times in (37), but it has to be applied only twice, because two of the applications of are on the same vector.

The matrix can be assembled from local contributionsOne can show that is sparse and that it can be factorized using standard sparse LU factorization [12]. The size of is determined by the number of primal DOF. Since we only have primal DOF at the subdomain vertices, their number is bounded by . Typically, the number of subdomains, and therefore the size of , is much smaller than the size of , but even in the case of many subdomains is sparse.

We can solve the symmetric and positive definite system (36) for by a CG algorithm. Once we have obtained , we can calculateThe remaining local solutions are then given byIn [19], the unpreconditioned interface problem (36) is discussed for the classical FETI method, and it is shown that the condition number is of orderwhere H and h denote the characteristic subdomain size and the finite element mesh size, respectively. The numerical tests presented in Section 6 indicate that the IETI method behaves similarly. In the next section, following [18], we define a preconditioner for the interface problem which will be used for the numerical examples in Section 6.

Preconditioner

Our construction follows the scaled Dirichlet preconditioner that was introduced in [19,28,44] and extended to the dual-primal formulation in [18,29]. We indicate interior DOF with the subscript I, and the DOF associated with the boundary of a subdomain with the subscript B. Assume that the DOF are now numbered such that the interior DOF are listed first, then the local stiffness matrix takes the formThe dual-primal Dirichlet preconditioner is defined bywhereSince is the local stiffness matrix of with all boundary DOF fixed, it can be factorized as easily and cheaply as . The matrix in (42) is the restriction of to the interface conditions associated with . The matrix is a scaled diagonal matrix of size , where J is the number of Lagrange multipliers. Its entries arewhere is the number of subdomains which have interfaces associated with the Lagrange multiplier . In particular, takes the following values:These scalings can be found in e.g. [28,29], where the authors show that certain jumps in the diffusion coefficient (problem (I)) or the Lamé parameters (problem (II)) can be treated robustly.

if corresponds to an essential boundary condition.

if corresponds to a coupling condition that does not involve a subdomain vertex.

if corresponds to a coupling condition that involves a subdomain vertex.

In [29], it was shown that the condition number of the preconditioned FETI-DP interface problem behaves likewhere H and h are as defined at the end of Section 4.3. The numerical results presented in Section 6 show that a similar behaviour can be observed in the IETI method.

Isogeometric Tearing and Interconnecting algorithm

To summarize, the overall IETI-DP algorithm is as follows.

For each , locally on each subdomain (in parallel)

Assemble the local stiffness matrix and load vector using a fully local numbering of the DOF.

Partition and as in (28), factorize , calculate .

Partition as in (41), factorize .

Assemble and factorize , calculate .

Solve by PCG with preconditioner as in (42).

Calculate as in (38),

For each , obtain as in (39) (in parallel),

Refinement options

The tensor-product structure of NURBS basis functions is inconvenient for local refinement. The insertion of a knot affects the whole domain and may introduce superfluous DOF. Recent work on methods for local refinement in IGA includes analysis suitable T-splines (see, e.g. [3,34,35,47-49]), THB-splines (see [20]), and PHT-splines (see [13,40]), as well as the use of finite element-based strategies in the IGA framework (see [30]).

However, the IETI-approach introduces some possibilities for restricting the refinement to one or a few subdomains (of many), even when working with tensor-product NURBS basis functions and straightforward knot insertion. We will sketch two such methods: h-refinement on one subdomain and refinement by substructuring (see Fig. 3). Note that, in both cases, we assume that the initial setting is fully matching.

Fig. 3

Two options for refining the shaded area in (a).

h-refinement on one subdomain

Although a knot insertion affects the whole parameter domain due to the tensor-product structure of the NURBS basis functions, we can limit the refinement to a single subdomain, as depicted in Fig. 3(b). Such a local refinement procedure was proposed in [10] in the context of multi-patch NURBS discretizations.

We assume that on an interface, which is not fully matching, the knot vector on one subdomain (the fine side) is a refinement of the knot vector on the other subdomain (the coarse side), as in the example in Fig. 3(b). In Fig. 4, such a case is illustrated schematically: The knot vector is obtained from by one step of uniform h-refinement. In contrast to the fully matching case, the numbers of DOF of and on the interface are not equal, and condition (ii) for the fully matching case is not fulfilled. In reference to hanging nodes in finite element methods, we call this a setting with hanging knots (note that we still assume that the geometry is conforming). As a consequence, we cannot couple the DOF as simply as in (18). In particular, the matrix in (17) has to be modified accordingly.

Fig. 4

Interface with hanging knots: , is a refinement of .

The number of interface conditions on such an interface is determined by the fine side, which is chosen as the master subdomain. Hence, we adapt the definition (15) of as follows:

If is an interface with hanging knots, we defineIf is a fully matching interface, we follow the definition of in (15), i.e.,For example, in Fig. 4, the master subdomain is , i.e. .

Without loss of generality, we assume that is a refinement of and that the weights on the finer side are obtained by the knot insertion algorithm [43]. Hence, on the interface the coarse basis function can be represented exactly as a linear combination of fine basis functions . Therefore, for each , there exist coefficients , such thatThe coefficients can be obtained from well-known formulae for the refinement of B-Spline basis functions [43].

We require -continuity of across the interface , i.e. we requireBy comparing the coefficients, we obtaini.e. we obtain a continuity constraint in the same form as in (18) and (20), where the coefficients of the coupling matrix , see (17), are given by

In [24], the issue of coupling mortar discretizations in the FETI-DP context was addressed. The procedure for formulating the jump operator given therein would result in the same coefficients if it is applied to the setting considered here.

With the modified coupling matrices, one can then perform the same steps as in Section 4.3.2 and the sections thereafter. Here we would like to point out that there are Lagrange multipliers which connect primal essential DOF with a boundary condition as well as multipliers connecting one and the same primal DOF. Both types of multipliers are superfluous and can be left out. However, in the presence of an interface with hanging knots, there are also multipliers that connect primal and remaining DOF. These multipliers cannot be left out and the corresponding entries in the vector do not vanish in general.

Local refinement by substructuring

As an alternative, the number of DOF can be increased locally by subdividing one subdomain into smaller subdomains as illustrated in Fig. 3(c).

If one wanted to subdivide the NURBS geometry mapping, e.g. in order to be able to edit the geometry locally, one would split the mapping and construct new knot vectors, weights and control nets for the new geometry mappings. This is not necessary in the IETI context, since we are not interested in actually splitting the geometry representation. Instead, we split the parameter domain into smaller subdomains which are mapped to the physical domain by the original (unchanged) coarse geometry mapping.

A very simple method for substructuring the subdomain is the following: We split the parameter domain into four subdomainsWe refer to this substructuring method as cross insertion. The basis functions of the original parameter domain Q are pushed forward to the smaller subdomain by a linear mapping and then transformed to the physical domain by the original mapping (see Fig. 5 for an illustration). The basis functions on have the formThe domain decomposition obtained by substructuring is again a setting with hanging knots. The matrix in the interface condition (18) is adapted as described in Section 5.1.

Fig. 5

Embedding smaller subdomains into the original parameter domain .

When we refine by substructuring, we introduce situations where the vertex of one subdomain coincides with the edge of another subdomain. Such cases are illustrated in Fig. 3(c) and and Fig. 6(a). We call such a subdomain vertex a hanging subdomain vertex (or short hanging vertex). Note that not every T-shaped subdomain vertex is a hanging vertex, as illustrated in the example in Fig. 6(b).

Fig. 6

Examples for hanging and not hanging subdomain vertices.

The choice of primal DOF in substructured subdomains, where we have hanging vertices, is not as straightforward as in the fully matching case. In the example of a hanging vertex in Fig. 6(a), there is (in the scalar case) exactly one DOF on that is associated with the hanging vertex marked by the dashed circle (cf. the discussion at the beginning of Section 4.3). While the same applies to , this is not true on , where we have several NURBS basis functions which are nonzero at the marked hanging vertex. Instead of incorporating a special treatment of hanging vertices, we choose to omit primal DOF at hanging vertices and discuss under which conditions this is possible.

For the scalar elliptic problem (I), the kernel of the stiffness matrix of a floating subdomain is spanned by the constant function, i.e., the kernel has dimension one. In this case, it is sufficient to have at least one primal DOF on each subdomain. This is easily guaranteed, if we start from a fully matching setting, apply substructuring by cross-insertion as described in Section 5.2, and select primal DOF at all subdomain vertices which are not hanging. The example in Fig. 7(a) shows the positions of primal DOF after two cross insertions.

Fig. 7

Subdomains refined by substructuring. Positions of primal vertices marked by .

In a two-dimensional setting, we call a subdomain vertex a primal vertex, if the two DOF associated with this vertex are primal. For the two-dimensional linearized elasticity problem (II), where the kernel is spanned by the three rigid body modes, we need at least three primal DOF per subdomain. This is satisfied, if we choose at least two primal vertices per subdomain (i.e., four primal DOF in total). However, as illustrated in Fig. 7(a), this is not guaranteed if we apply substructuring by cross-insertion without additional considerations.

For linearized elasticity problems, we introduce refinement levels and we assign refinement level 0 to every subdomain in the initial setting. When a subdomain is split into four smaller subdomains by cross insertion, the levels of the new, smaller subdomains are increased by 1 (see Fig. 8 for an illustration). We call the refinement a 1-level substructuring, if the refinement levels of any two subdomains with an edge as their interface differ by at most 1. If we start from a fully matching setting, apply 1-level substructuring by cross-insertion, and choose all non-hanging vertices as primal vertices, then it is guaranteed that there are at least two primal vertices, i.e., at least four primal DOF, on each subdomain. The example in Fig. 7(b) illustrates the positions of primal vertices after two such 1-level substructuring steps. Note that, depending on the location of the refined area, 1-level substructuring can effect neighbouring subdomains. This disadvantage is accepted as a trade-off for avoiding an involved treatment of hanging vertices.

Fig. 8

Refinement levels of subdomains (initial setting as in Fig. 3(a)).

Note that the discretization is only -continuous along subdomain interfaces. By substructuring a subdomain, new interfaces are introduced, and thereby the discretization is changed.

Note that the refinement options discussed in Section 5.1 and Section 5.2 can be combined by applying them one after the other. Furthermore, the coupling methods described in this section can be combined with isogeometric local refinement methods such as those mentioned at the beginning of Section 5. If such local refinement methods are applied only in the interior of one or many subdomains, obviously, the coupling at the subdomain interfaces is not influenced. If the refined areas extend to subdomain boundaries, and one side is a refinement of the other, coupling can be done as described in this section.

Preconditioning in the presence of hanging knots

As mentioned in Section 5.1, when we have hanging knots, the coupling matrix is not a signed Boolean matrix any more. The Dirichlet preconditioner that was defined in (42) can still be applied in settings with hanging knots. However, as already mentioned in [24] in the context of mortar discretizations, while the asymptotic behaviour of the condition number remains the same, the condition number itself increases. This can also be observed in the numerical tests with the IETI-DP method (see Section 6 for the results).

We now adapt the preconditioner for settings with hanging knots by replacing the scaling matrix in (42) by a modified diagonal matrix . Its entries are defined as follows:The preconditioner for the case with hanging knots is defined analogous to (42) byAs it will be reported in Section 6, this preconditioner leads to lower condition numbers as compared to from Section 4.4 in settings with hanging knots. Note that we have in fully matching settings.

if corresponds to a fully matching interface. Here, is as defined in Section 4.4.

if corresponds to an interface with hanging knots and is the master subdomain.

otherwise.

Numerical examples

In this section, we present three numerical test examples for the IETI-DP method. We test the method and the refinement options presented in Section 5, and we study the performance of the proposed preconditioners. The conjugate gradient method was applied to solve the interface problem in (36), both without and with the discussed preconditioners. In the following tables, we display the numbers of iterations needed until the stopping criterionwas fulfilled, where is the residual in the kth iteration and denotes the Euclidean norm. The condition numbers in these tables were computed numerically using the Lanczos method.

Note that the aim of these examples (due to their sizes, they could be dealt with by using direct solvers) is to illustrate the method and its potential. The true computational advantage of the IETI method against direct solvers will become apparent in large problem sizes and/or three-dimensional problems.

Bracket under load

Our first example is a linearized elasticity problem (type (II) in Section 2.1). We consider the two geometries displayed in Fig. 9(a) and and Fig. 10(a), where the first one, which is taken from an illustration in [11], has a rounded reentrant corner, and the second one has a sharp reentrant corner. We refer to these geometries as case (A) and case (B), respectively.

Fig. 9

Case (A), bracket with rounded reentrant corner.

Fig. 10

Case (B), bracket with sharp reentrant corner.

In Fig. 9(a) and 10(a), we show the subdomain decomposition and indicate the boundary conditions. We fix the two lower holes by applying homogenous essential boundary conditions, while a constant downward pointing traction with magnitude 1000 N is applied at the walls of the rightmost hole. The remaining boundaries are free of traction. The material parameters are set to and . Note that the circular holes contained in the domains are represented exactly by NURBS geometry mappings of degree 2.

In Fig. 9(b) and 10(b), the calculated stress component is depicted for cases (A) and (B), respectively. Note that the scales are different, and that the scale in Fig. 10(b) has been cutoff below for a better visibility of the stress distribution. The results illustrate that the IETI-DP method can be applied to non-trivial geometries including holes and consisting of numerous subdomains.

The condition numbers and the (P)CG iteration numbers for these fully matching settings are given in the table in Fig. 11. The column labeled n shows the number of knot spans in the direction indicated by the small arrows at the bottom in Fig. 9(a) and 10(a). The column shows the number of Lagrange multipliers, i.e., the size of the interface problem (36). The columns labeled and display the condition numbers of the interface problems and the (P)CG iteration numbers without preconditioner, and with preconditioner as defined in Section 4.4, respectively. The entry “>200” indicates that the desired accuracy was not reached after 200 iterations. The results show the expected, moderate growth in the preconditioned case.

Fig. 11

Condition numbers and (P)CG iterations for cases (A) and (B) of a bracket under load.

In Fig. 10(b), the peak stress near the reentrant corner in case (B) is clearly visible. To obtain a better resolution of the peak stress, we introduce case (C), which is indicated in Fig. 12(a). Here, the subdomains near the corner have a finer discretization than the subdomains which are far from the corner, and we have interfaces with hanging knots. The number of knot spans on the finest and coarsest subdomain discretizations are denoted by and , respectively. These numbers are measured in the directions indicated by the small arrows in Fig. 12(a). The ratio is the same for all chosen meshes. The condition numbers and (P)CG iteration numbers for this setting with hanging knots are presented in Fig. 12(b). Clearly, the preconditioner defined in Section 5.3 performs better than from Section 4.4. The stress component is not plotted for case (C), since it is the same as in case (B), Fig. 10(b). The energy norm of the numerical solutions in case (B) and case (C) is compared in Fig. 13. It shows that for given DOF a faster convergence can be achieved by local h-refinement.

Fig. 12

Case (C), bracket with sharp reentrant corner and local h-refinement near the corner. The ratio nfine/ncoarse = 4 is the same for all chosen meshes. The stress component as in Fig. 10(b).

Fig. 13

Comparison of the energy norms of the discrete solutions in cases (B) and (C).

Bending of a cantilever

We now consider a linearized elasticity problem (type (II) in Section 2.1) on a cantilever of length L and thickness D. It is fixed at and subject to a parabolic traction at with resultant P as illustrated in Fig. 14(a). We choose the parameters as follows: , and .

Fig. 14

Bending of a cantilever, problem setting and discussed cases.

An analytical solution for the displacement field can be found, e.g. in [30,40,51]:where is the moment of inertia of the cross section of the cantilever. When we consider the plane stress problem, we set and . For the plane strain problem, we set and . Then, in both cases, the resulting exact stress components are as follows:We apply the exact displacement as essential boundary conditions at the boundary , and the exact traction at . The remaining boundaries at are free of traction.

Note that the exact displacement field is a cubic polynomial. By using basis functions of degree 3 and due to the simple geometry mappings, the exact solution u is in . Although the exact solution does not have peaks or singularities, we refine randomly chosen subdomains, in order to demonstrate the performance of the IETI-DP method in settings with hanging knots. In Fig. 14(b), the thick lines indicate the subdomain decomposition, while the thin lines schematically indicate how fine the discretizations of the respective subdomains are, and whether the setting is fully matching or not. The numerical tests confirm that, for all considered discretizations, the calculated numerical solution is exact (up to the accuracy to which we solve (36) for ).

The table presented in Fig. 14(d) shows that the condition numbers behave as expected. For cases which are not fully matching, the columns labeled and indicate the number of knot spans in one direction on the coarse and fine discretizations, respectively. The other columns are labeled as in the previous example. In the settings (C), (D), and (E) with hanging knots, it can be observed that the condition number of grows slowlier than in the unpreconditioned case. However, on coarse discretizations, both the absolute value of the condition number and the number of (P)CG iterations are larger in the preconditioned case than in the unpreconditioned case (cf. the discussion in the beginning of Section 5.3). The preconditioner performs much better in all cases with hanging knots.

Poisson problem on Yeti’s footprint

In our third example, we solve the Poisson problem (type (I) in Section 2.1 with ) on the physical domain resembling the footprint of a Yeti. The domain is shown in Fig. 15(a) and and consists of 21 subdomains. We set Dirichlet boundary conditions at the big toe, and Neumann boundary conditions everywhere else. The boundary conditions and the right hand side f are determined by the exact solutionwhere . This solution u is constructed in such a way that it has a peak at . We set and (see Fig. 15(a)).

Fig. 15

Yeti’s footprint with adaptive refinement.

Since the discussion of a posteriori error estimation is not in the scope of this paper, and because we know the exact solution u, we apply adaptive refinement based on the exact error. For each subdomain , we calculate the local error in the energy norm . Then we mark all subdomains for refinement, for whichholds. In each such refinement step, we apply uniform h-refinement on the marked subdomains. The global mesh after 5 such refinement steps is shown in Fig. 15(b). In Fig. 15(c), we compare the error obtained by global uniform refinement and by the described adaptive refinement. As expected, for a given number of DOF a more accurate solution can be achieved by adaptive, local h-refinement, as compared to global refinement in a fully matching setting.

In Fig. 15(d), the condition numbers and the (P)CG iteration numbers for the fully matching setting are presented. The column labeled “DOF” indicates the global number of DOF. In Fig. 15(e), the condition numbers and the (P)CG iterations in the adaptive refinement are shown.

Conclusion

We have proposed the IETI method which combines the ideas and advantages of IGA and the FETI method. We preserve the exact geometry representation from the coarsest discretization level, thereby eliminating the need for data transformation, and also eliminating consequent approximation errors in the geometry. At the same time, we apply the techniques from FETI methods to couple NURBS patches and to solve the presented model problems on complicated computational domains which consist of many NURBS patches and may contain holes.

We have discussed the coupling of interfaces with hanging knots, thereby introducing options for local refinement. These can be applied using only NURBS basis functions with tensor-product structure, without the need for involved local NURBS-refinement techniques. Numerical examples demonstrated the performance of the IETI-DP method and of the proposed preconditioners, both in fully matching settings, and in settings with hanging knots.

In this paper, we have assumed that the subdomain interfaces are either fully matching or that the discretization on one side is a refinement of the other. The coupling across fully matching interfaces can be extended to three-dimensional problems in a straightforward manner. Also, preconditioners from FETI-DP methods can be applied to a three-dimensional IETI method. The solver design and the choice of primal DOF (including edge averages which have been proven necessary in three dimensions [29,52]) are more involved, in particular in the presence of hanging vertices. The treatment of more general interfaces, including interfaces that are not necessarily geometrically conforming (i.e., that may have small gaps or overlaps) is open for future work. Another interesting issue for future research on the IETI method is the incorporation of fast iterative subdomain solvers, such as geometric multigrid solvers and solvers exploiting the tensor product structure, e.g., wavelet solvers.