Classical and all-floating FETI methods for the simulation of arterial tissues

Classical and all-floating FETI methods for the simulation of arterial tissues

Christoph M. Augustin   christoph.augustin@medunigraz.at; Corresponding author Institute of Computational Mathematics, Graz University of Technology, Steyrergasse 30, 8010 Graz, AustriaInstitute of Biomechanics, Center of Biomedical Engineering, Graz University of Technology, Kronesgasse 5, 8010 Graz, Austria Institute of Biophysics, Medical University of Graz, Harrachgasse 21, 8010 Graz, Austria Gerhard A. Holzapfel Institute of Computational Mathematics, Graz University of Technology, Steyrergasse 30, 8010 Graz, AustriaInstitute of Biomechanics, Center of Biomedical Engineering, Graz University of Technology, Kronesgasse 5, 8010 Graz, Austria Institute of Biophysics, Medical University of Graz, Harrachgasse 21, 8010 Graz, Austria Olaf Steinbach Institute of Computational Mathematics, Graz University of Technology, Steyrergasse 30, 8010 Graz, AustriaInstitute of Biomechanics, Center of Biomedical Engineering, Graz University of Technology, Kronesgasse 5, 8010 Graz, Austria Institute of Biophysics, Medical University of Graz, Harrachgasse 21, 8010 Graz, Austria
Abstract

High-resolution and anatomically realistic computer models of biological soft tissues play a significant role in the understanding of the function of cardiovascular components in health and disease. However, the computational effort to handle fine grids to resolve the geometries as well as sophisticated tissue models is very challenging. One possibility to derive a strongly scalable parallel solution algorithm is to consider finite element tearing and interconnecting (FETI) methods. In this study we propose and investigate the application of FETI methods to simulate the elastic behavior of biological soft tissues. As one particular example we choose the artery which is – as most other biological tissues – characterized by anisotropic and nonlinear material properties. We compare two specific approaches of FETI methods, classical and all-floating, and investigate the numerical behavior of different preconditioning techniques. In comparison to classical FETI, the all-floating approach has not only advantages concerning the implementation but in many cases also concerning the convergence of the global iterative solution method. This behavior is illustrated with numerical examples. We present results of linear elastic simulations to show convergence rates, as expected from the theory, and results from the more sophisticated nonlinear case where we apply a well-known anisotropic model to the realistic geometry of an artery. Although the FETI methods have a great applicability on artery simulations we will also discuss some limitations concerning the dependence on material parameters.

This is the peer reviewed version of the following article: Augustin, C. M., Holzapfel, G. A. and Steinbach, O. (2014), Classical and all-floating FETI methods for the simulation of arterial tissues. Int. J. Numer. Meth. Engng., 99: 290–312. doi: 10.1002/nme.4674, which has been published in final form at onlinelibrary.wiley.com/doi/10.1002/nme.4674/abstract. This article may be used for non-commercial purposes in accordance with Wiley Terms and Conditions for Self-Archiving.
Keywords: artery, biological soft tissues, all-floating FETI, parallel computing

1 Introduction

The modeling of hyperelastic materials is realized by using a strain–energy function . For a comprehensive overview and the mathematical theory on elastic deformations, see, e.g., [1, 2, 3, 4]. A well established model for arterial tissues was introduced by Holzapfel et al. [5, 6]. This model was further developed and enlarged to collagen fiber dispersion in [6, 7, 8]; see [9] for the modeling of residual stresses in arteries which play also an important role in tissue engineering. An adequate model for the myocardium can be found in [10]. The fine mesh structure to model cardiovascular organs normally results in a very large number of degrees of freedom. The combination with the high complexity of the underlying partial differential equations demands fast solution algorithms and, conforming to up–to–date computer hardware architectures, parallel methods. One possibility to achieve these specifications is to use domain decomposition (DD) methods which acquired a lot of attention in the last years and resulted in the development of several overlapping as well as non–overlapping DD methods, see [11]. They all work according to the same principle: the computational domain is subdivided into a set of (overlapping or non–overlapping) subdomains . DD algorithms now decompose the large global problem into a set of smaller local problems on the subdomains, with suitable transmission or interface conditions. This yields a natural parallelization of the underlying problem. In addition to well established standard DD methods, other examples for more advanced domain decomposition methods are hybrid methods [12], mortar methods [13, 14, 15] and tearing and interconnecting methods [16].

In this paper we focus on the finite element tearing and interconnecting (FETI) method where the strategy is to decompose the computational domain into a finite number of non–overlapping subdomains. Therein the corresponding local problems can be handled efficiently by direct solvers. The reduced global system, that is related to discrete Lagrange multipliers on the interface, is then solved with a parallel Krylov space method to deduce the desired dual solution. This is, in the case of elasticity, the boundary stress and subsequentely, in a postprocessing step, we compute the primal unknown, i.e. the displacements, locally. For the global Krylov space method, such as the conjugate gradient (CG) or the generalized minimal residual (GMRES) method, we need to have a suitable preconditioning technique. Here we consider a simple lumped preconditioner and an almost optimal Dirichlet preconditioner, as proposed by Farhat et al. [17].

A variant of the classical FETI method is the all-floating tearing and interconnecting approach (AF-FETI) where, in contrast to the classical approach, the Dirichlet boundary acts as a part of the interface. It was introduced independently for the boundary element method by Steinbach and Of [18, 19] and as the Total-FETI (TFETI) method for finite elements by Dostál et al. [20]. This approach shows advantages in the implementation and, due to mapping properties of the involved operators, improves the convergence of the global iterative method for the considered problems. This behavior is illustrated with numerical examples, which are – to the best of our knowledge – the first application of all-floating FETI method to nonlinear and anisotropic biological materials.

An essential part of FETI methods is solving the local subproblems. Challenges occur with so-called floating subdomains which have no contribution to the Dirichlet boundary. These cases correspond to local Neumann problems and the solutions are – in the case of elasticity – only unique up to the six rigid body modes. For classical FETI it can happen that the kernel of the local operator is non-trivial and its dimension is lower than six. The problem to identify these kernels reliably causes trouble. One possibility to overcome this trouble is a modification of the classical approach, the dual-primal FETI (FETI-DP) method, cf. Farhat et al. [21] and Klawonn and Widlund [22]. In this variant some specific primal degrees of freedom are fixed. This yields solvable systems for all subdomains. Choosing the primal degrees of freedom may be very sophisticated [23]. This approach was already applied to model arterial tissues using FETI-DP by Klawonn and Rheinbach [24, 25], Brands et al. [26], Balzani et al. [28, 29] and Brinkhues et al. [27]. Note that for all-floating FETI the identification of the kernel of the local operators is no problem at all, since we treat all subdomains as floating subdomains, and hence have a kernel equal to six for all local operators. Moreover the resulting local systems are typically better conditioned than those arising in the FETI-DP approach, see Brzobohatỳ et al. [30]. All-floating FETI was used to model myocardial tissue in the preliminary work [31].

Both the classical FETI method, as well as all-floating FETI, need the construction of a generalized inverse matrix. This may be achieved using direct solvers with a sparsity preserving stabilization, see, e.g. [30], or stabilized iterative methods. For a mathematical analysis of FETI methods including convergence proofs for the classical one-level FETI method, see, e.g., [22, 32, 33].

2 Modeling Arterial Tissues

The deformation of a body is described by a function with the reference configuration at time and the current configuration at time . With this we introduce the displacement field in the reference configuration and the displacement field in the current configuration,

(1)

and the deformation gradient as, see, e.g., [2],

(2)

Moreover, we denote by the Jacobian of and by the right Cauchy–Green tensor. For later use, to model the nearly incompressible behavior of biological soft tissues, we introduce the following split of the deformation gradient in a volumetric and an isochoric part, compare Flory [34], i.e.

(3)

Consequently, this multiplicative split can be applied to other tensors such as the right Cauchy–Green tensor. Thus

(4)

As a starting point for the modeling of biological soft tissues the stationary equilibrium equations in the current configuration are considered to find a displacement field according to

(5)

where is the Cauchy stress tensor and is the body force at time .

In addition, we incorporate boundary conditions to describe displacements or normal stresses on the boundary , which is decomposed into disjoint parts such that . Dirichlet boundary conditions on correspond to a given displacement field , while Neumann boundary conditions on are identified physically with a given surface traction , where denotes the exterior normal vector at time .

The equilibrium equations and the boundary conditions may also be formulated in terms of the reference configuration, i.e.

(6)
(7)
(8)

where is the second Piola–Kirchhoff tensor and is the body force at time . In order to formulate the boundary conditions we introduce a prescribed displacement field , the exterior normal vector and the surface traction in the reference configuration.

Considering the study of the properties of soft biological soft tissues we have to deal with a nonlinear relationship between stress and strain, with large deformations and an anisotropic material. Since linear elasticity models are not adequate for treating such a complex behavior, we take a look at the more general concept of nonlinear elasticity.

The nonlinear stress-strain response is modeled via a constitutive equation that links the stress to a derivative of a strain-energy function , representing the elastic stored energy per unit reference volume. Derived from the Clausius–Duhem inequality, see [35, 36], we formulate the constitutive equations as

(9)

We make use of the Rivlin–Ericksen representation theorem [37] and its extension to anisotropic materials, cf. [38], to find a representation of the strain-energy function in terms of the principal invariants of .

Arteries are vessels that transport blood from the heart to the organs. In vivo the artery is a prestretched material under an internal pressure load. Healthy arteries are highly deformable composite structures and show a nonlinear stress-strain response with a typical stiffening effect at higher pressures. Reasons for this are the embedded collagen fibers which lead to an anisotropic mechanical behavior of arterial walls. We denote by and the predominant collagen fiber directions in the reference configuration. An important observation is that arteries do not change their volume within the physiological range of deformation, hence they are treated as a nearly incompressible material, see, e.g., [5]. In this work we focus on the in vitro passive behavior of the healthy artery, see Fig. 1.

1Endothelial cellInternal elastic laminaSmooth muscle cellCollagen fibrilsElastic fibrilsElastic laminaExternal elastic laminaBundles of collagen fibrils
Helically arranged fiber- reinforced medial layers
Composite reinforced by collagen fibers arranged in helical structures
Figure 1: Diagrammatic model of the major components of a healthy elastic artery, from [5]. The intima, the innermost layer is negligible for the modeling of healthy arteries, it plays a very important role in the modeling of diseased arteries, though. The two predominant directions of the collagen fibers in the media and the adventitia are indicated with black curves.

To capture the nearly incompressibility condition we remember the decomposition (3), which yields an additive split of the strain-energy function into a so-called volumetric and an isochoric part, i.e.

(10)

This procedure leads to constitutive equations in which the stress tensors are also additively decomposed into a volumetric and an isochoric part, i.e., cf. [2],

(11)

Here, the scalar-valued hydrostatic pressure is defined as

(12)

To capture the specifics of this fiber-reinforced composite, Holzapfel and Weizsäcker [39] and Holzapfel et al. [5] proposed an additional split of the strain-energy function into an isotropic and an anisotropic part so that the complete energy function can be written as

(13)

Following the classical approach we describe the volume changing part by

(14)

where , comparable to the bulk modulus in linear elasticity, serves as a penalty parameter to enforce the incompressibility constraint.

To model the isotropic non-collagenous matrix material the classical neo-Hookean model is used [2]. Thus

(15)

where is a stress-like material parameter and is the first principal invariant of the isochoric part of the right Cauchy–Green tensor. In (13), is associated with the deformation of the collagen fibers. According to [5], this transversely isotropic response is described by

(16)
(17)

with the invariants , and the material parameters and , which are both assumed to be positive. It is worth to mention that for the anisotropic responses, (16) and (17) only contribute for the cases or , respectively. This condition is explained with the wavy structure of the collagen fibers, which are regarded as not being able to support compressive stresses. Thus, the fibers are assumed to be active in tension () and inactive in compression (). This assumption is not only based on physical reasons but it is also essential for reasons of stability, see Holzapfel et al. [40].

The material parameters can be fitted to an experimentally observed response of the biological soft tissue. Following [5] we use the material parameters summarized in Table 1.

 kPa  kPa (-)
Table 1: Material parameters used in the numerical experiments; parameters taken from Holzapfel et al. [5].

Similar models can also be used for the description of other biological materials, e.g., for the myocardium, cf. [10].

3 Finite Element Approximation

3.1 Variational formulation of nonlinear elasticity problems

In this section we consider the variational formulation of the equilibrium equations (5) and (6) with the corresponding Dirichlet and Neumann boundary conditions. In particular, using spatial coordinates, the boundary value problem (5) is formally equivalent to the variational equations

(18)

valid for a smooth enough tensor field and all smooth enough vector fields , which vanish on , see, e.g., [1, Theorem 2.4-1]. Additionally,

(19)

and is the nonlinear operator in the current configuration which is induced by the stress tensor representation (11), and by using the related duality pairing . For later use, we introduce the corresponding terms in the reference configuration as and . Note that (18) formally corresponds to a variational formulation in linear elasticity. However, the integral and the involved terms have to be evaluated in the current configuration which comprises the nonlinearity of the system. If the test function is interpreted as the spatial velocity gradient, then is the rate of deformation tensor so that has the physical interpretation of the rate of internal mechanical work.

In terms of the reference configuration, the boundary value problem (6), (8) is formally equivalent to the variational equations

(20)

valid for a smooth enough tensor field and all smooth enough vector fields with on , see, e.g., [1, Theorem 2.6-1]. In (20) we use the definition of the directional derivative of the Green–Lagrange strain tensor, i.e.

(21)

which is also known as the variation or the material time derivative of the Green–Lagrange strain tensor in the literature.

It is important to note that results on existence of solutions in nonlinear elasticity can be stated given a polyconvex strain-energy function , which holds true for the anisotropic model (13) discussed in Section 2. For more details we refer to the results of Ball [41, 42], see also [1, 43] and Balzani et al. [44].

3.2 Linearization and discretization

In the following we confine ourselves to the reference configuration . The formulations in the current configuration can be deduced in an analogous way.

For the solution of the nonlinear system (20) we apply Newton’s method to obtain the recursion

(22)

with the tangential term , the displacement field of the -th Newton step , the increment and a suitable initial guess .

For the computational domain we consider an admissible decomposition into tetrahedral shape regular finite elements of mesh size , i.e. , and we introduce a conformal finite element space , , of piecewise polynomial continuous basis functions . Then the Galerkin finite element discretization of the linearized variational formulation (22) results in a system of algebraic equations to find , on such that

(23)

holds for all , on . Note that the initial guess has to satisfy an approximate Dirichlet boundary condition on to fulfill condition (7), where denotes a suitable approximation of the given displacement . For the computation of the tangential term we need to evaluate

(24)

For a more detailed presentation how to compute the tangential term, in particular the forth-order elasticity tensor we refer to [46, 45].

Note that the convergence rate of the Newton method is dependent on the initial guess, on the parameters used in the model and on the inhomogeneous Dirichlet and Neumann boundary conditions which influence .

In a time-stepping scheme we use zero for the initial guess, and the result of the -th time step as initial solution for the next step. The initial guess may also be the solution of a modified nonlinear elasticity problem such as the solution of the same nonlinear model but with modified parameters, e.g., a reduced penalty parameter , or modified boundary conditions, e.g., a reduced pressure on the surface. The latter is equivalent to an incremental load stepping scheme with a parameter , so that

(25)

Klawonn and Rheinbach [24] used a load stepping scheme of this kind, for more information on load stepping and global Newton methods, see [48, 47]. The standard finite element method (FEM) now yields a linear system of equations which is equivalent to the discretized variational formulation (23). Finally, we have to solve

(26)

with the solution vector in the -th Newton step and the increment . The tangent stiffness matrix is calculated according to

(27)

and the terms of the right hand side are constructed by

(28)

The additive split of the stress tensors (11) and the introduction of the hydrostatic pressure (12) leads to the additional equation

(29)

which has to be satisfied in a weak sense. For this we use the idea of static condensation where this volumetric variable is eliminated element-wise, see, e.g., [46]. This may be achieved in using discontinuous basis functions; in this paper we will concentrate on piecewise constants. In the case of tetrahedral elements, this approach leads to elements. Here is the order of the basis functions for the displacement field. It is known that linear finite elements are very prone to volumetric locking. Hence, for nearly incompressible materials piecewise quadratic elements () are a better choice, see Simo [49]. The resulting element is also the preferred choice to model nearly incompressible arterial materials in [24][29]. For the numerical results in this work (Section 5) we use both linear ( element) and quadratic ( element) ansatz functions for the displacement field and compare the results.

Note that due to the symmetry of the stress tensor and the major and minor symmetry properties of the elasticity tensor the operator is self-adjoint. We can also show, using the positive definiteness of the elasticity tensor, see [4], and the polyconvexity of the strain-energy function (Section 3.1), that this operator is -elliptic and bounded, see [4, 45]. With these properties of the operator we can state that the linearized system (23),(24) admits a unique solution . Furthermore, the tangent stiffness matrix is symmetric and positive definite.

Simulations with large deformations and the hence required derivative of the Neumann boundary conditions (8) would yield an additional non-symmetric mass matrix on the left hand side of (26). To stay with an symmetric system we neglect this matrix but compensate it with a surface update of the geometry after each Newton step. Thus, our whole system is symmetric and we can use the CG method as an iterative solver. Nonetheless, the FETI methods described in Section 4 also work for non-symmetric systems by using the GMRES method.

4 Finite Element Tearing and Interconnecting

To solve the linearized equations (26) arising in the Newton method we apply the finite element tearing and interconnecting approach [16], see also [24, 50, 51], and references given therein. The derivation of the FETI system for nonlinear mechanics will be performed in the reference configuration. In an analogous way this is also valid for the formulation in the current configuration.

Figure 2: Decomposition of a domain into four subdomains , .

For a bounded domain we introduce a non-overlapping domain decomposition

(30)

see Fig. 2. The local interfaces are given by for all . The skeleton of the domain decomposition (30) is denoted as

(31)

We assume that the finite element mesh matches the domain decomposition (30), i.e., we can reorder the degrees of freedom to rewrite the linear system (26) as

(32)

where the increments , the stiffness matrices and the terms on the right hand side , are related to the local degrees of freedom within the subdomain . All terms with an index correspond to degrees of freedom on the coupling boundary , see (31), while denote simple reordering matrices taking boolean values.

4.1 Classical FETI method

Starting from (32), the tearing is now carried out by

(33)

where is related to degrees of freedom on the coupling boundary . As the unknowns are typically not continuous over the interfaces we have to ensure the continuity of the solution on the interface, i.e.

(34)

This is done by applying the interconnecting

(35)

where the matrices are constructed from such that (34) holds. By using discrete Lagrange multipliers to enforce the constraint (35) we finally have to solve the linear system

(36)

4.2 all-floating FETI method

The idea of this special FETI method, cf., e.g., Of and Steinbach [19], is to treat all subdomains as floating subdomains, i.e. domains with no Dirichlet boundary conditions. In addition to the standard procedure of ‘gluing’ the subregions along the auxiliary interfaces, the Lagrange multipliers are now also used for the implementation of the Dirichlet boundary conditions, see Fig. 3. This simplifies the implementation of the FETI procedure since it is possible to treat all subdomains in the same way. In addition, some tests (Section 5) show more efficiency than the classical FETI approach and the asymptotic behavior improves. This is due to the mapping properties of the Steklov–Poincaré operator, see [19, Remark 1]. The drawback is an increasing number of degrees of freedom and Lagrange multipliers. Compare also to Dostál et al. [20] for the related Total-FETI method.

(a) (b)
Figure 3: Fully redundant classical FETI (a) and all-floating FETI (b) formulation: , , denote the local subdomains, the black dots correspond to the subdomain vertices and the dashed lines correspond to the constraints (34). The gray strip indicates Dirichlet boundary conditions. Note that the number of constraints for the all-floating approach rises with the number of vertices on the Dirichlet boundary.

If all regions are treated as floating subdomains the conformance of the Dirichlet boundary conditions is not given; they have to be enhanced in the system of constraints using the slightly modified interconnecting

(37)

where is a block matrix of the kind and the vector is of the form such that , if and only if is the index of a Dirichlet node of the subdomain , while equals the Dirichlet values corresponding to the vertices , see also [19].

For three-dimensional elasticity problems all subdomain stiffness matrices have now the same and known defect, which equals the number of six rigid body motions and which also simplifies the calculation of the later needed generalized inverse matrices . For all-floating FETI we finally get the linearized system of equations

(38)

4.3 Solving the FETI system

To solve the linearized systems (36) and (38) we follow the standard approach of tearing and interconnecting methods. For convenience we outline the procedure by means of the classical FETI formulation (Section 4.1). However the modus operandi is analogous for the all-floating approach.

First, note that in the case of a floating subdomain , i.e. , the local matrices are not invertible. Hence, we introduce a generalized inverse to represent the local solutions as

(39)

Here, correspond to the rigid body motions of elasticity and are unknown constants. For floating subdomains we additionally require the solvability conditions

(40)

In the case of a non-floating subdomain, i.e. , we may set . Note that it may happen that the kernel is non-trivial and its dimension is lower than 6. This is the case if the set is either a vertex or an edge. For classical FETI methods this requires the implementation of an effective method to identify these kernels reliably. Note that this is a key advantage of the all-floating FETI approach because all subdomains are here treated as floating subdomains, and hence we know the kernel of each local operator . With these kernels the solution of the local problems to find the generalized inverse can be reduced to sparse systems which are typically better conditioned as the systems arising from the FETI-DP method, see Brzobohatỳ et al. [30]. In Section 4.2 we comment on an all-floating approach where also Dirichlet boundary conditions are incorporated by using discrete Lagrange multipliers.

In general, the Schur complement system of (36) is constructed to obtain

(41)

This can be expressed as

(42)

with

(43)

and is constructed using for and . For the solution of the linearized system (42) the projection

(44)

is introduced. It now remains to consider the projected system

(45)

This can be solved by using a parallel iterative method with suitable preconditioning of the form

(46)

with modified jump operators which are obtained by multiplicity scaling, see [24, 51]. Since the local subproblems all yield symmetric tangent stiffness matrices , cf. Section 3, the matrix is also symmetric. This enables us to use the CG method as the global solver for (45). Be aware that the initial approximate solution has to satisfy the compatibility condition . A possible choice is

(47)

In a post processing we finally recover the vector of constants

(48)

and subsequently the desired solution (39).

4.4 Preconditioning

Following Farhat et al. [17] we apply either the lumped preconditioner

(49)

or the optimal Dirichlet preconditioner

(50)

where

(51)

is the Schur complement of the local finite element matrix . Alternatively, one may also use scaled hypersingular boundary integral operator preconditioners, as proposed in [52]. For comparison we employ an identity preconditioner which is constructed by using the identity matrix for in eq. (46).

5 Numerical Results

In this section some representative numerical examples for the finite element tearing and interconnecting approach for linear and nonlinear elasticity problems are presented. First, the FETI implementation is tested within linear elasticity. Here we are able to compare the computed results to a given exact solution. This enables us to show the efficiency of our implementation and also the convergence rates, as predicted from the theory. We compare the different preconditioning techniques and present differences between the classical FETI and the all-floating FETI approach.

Subsequently, we apply the FETI method to nonlinear elasticity problems. Thereby, we focus on the anisotropic model, as described in Section 2, and use a realistic triangulations of the aorta and a common carotid artery. As in the linear elastic case, different preconditioning techniques for the all-floating and for the classical FETI method are compared. In Section 5.3, we analyze the biomechanical behavior of an aorta up to an internal pressure of  mmHg and plot stress and displacement evolutions as a function of the internal pressure. Finally, in Section 5.4, we analyze our computational framework with respect to strong scaling properties.

The calculations were performed by using the VSC2-cluster (http://vsc.ac.at/) in Vienna. This Linux cluster features 1314 compute nodes, each with two AMD Opteron Magny Cours 6132HE (8 Cores, 2.2 GHz) processors and 8 x 4 RAM. This yields the total number of available processing units. As local direct solver we use Pardiso [53, 54], included in Intel’s Math Kernel Library (MKL).

5.1 Linear elasticity

In this section of numerical benchmarks we consider a linear elastic problem with the academic example of a unit cube which is decomposed into a certain number of subcubes. Dirichlet boundary conditions are imposed all over the surface . The parameters used are Young’s modulus GPa and Poisson’s ratio . The calculated solution is compared to the fundamental solution of linear elastostatics

(52)

for all , is an arbitrary point outside of the domain , and is the Kronecker delta, see [55]. The different strategies of preconditioning are compared and also the all-floating and classical FETI approaches. As global iterative method we use the CG method with a relative error reduction of . Under consideration is a linear elasticity problem using linear tetrahedral elements ( elements) with a uniform refinement over five levels () given a cube with subdomains.

Hence, the number of degrees of freedom associated with the coarsest mesh is for the all-floating FETI approach and for the classical FETI approach. The difference of the numbers is due to the decoupling of the Dirichlet boundary . For the finest mesh we have (all-floating) and (classical) degrees of freedom. The number of Lagrange multipliers varies between for level and for level 5. Again we have a higher number of Lagrange multipliers for the all-floating approach due to the decoupling of the Dirichlet boundary conditions. The computations were performed on VSC2 using processing units.

First note in Table 2 that for all examined settings, the L2 error, i.e.

(53)

where is the approximate and the exact solution, and the estimated order of convergence

(54)

behaves as predicted from the theory, i.e. it is of second order. As expected the least iteration numbers were observed for the optimal Dirichlet preconditioner. Nonetheless, since no additional time is required to compute the lumped preconditioner, in contrast to the more sophisticated Dirichlet preconditioner, this type of preconditioning yields comparable computational times for each level of refinement. As a comparison we also list the results of a very simple preconditioning technique, using the identity matrix for in (46), where almost no reduction of the condition numbers can be noticed.

Moreover, we observe that all-floating FETI yields better condition numbers for all preconditioners, and hence better convergence rates of the global conjugate gradient method. Although the global iterative method converges in less iterations for this approach, we achieve lower computational time for the classical FETI method for the linear elastic case with elements. This is mainly due to the larger expenditure of time to set up the all-floating FETI system, the larger coarse matrix , cf. (44), and due to the higher amount of Lagrange multipliers.

all-floating
identity prec. lumped prec. Dirichlet prec. L2 error eoc
1 61 it. 53.6 20.9 s 27 it. 10.3 19.7 s 21 it. 7.6 19.5 s 1.42E-04 -
2 71 it. 70.0 19.6 s 38 it. 19.7 18.8 s 26 it. 10.4 18.4 s 3.71E-05 1.94
3 88 it. 108.8 21.7 s 45 it. 26.1 22.3 s 27 it. 9.7 22.3 s 9.40E-06 1.98
4 119 it. 216.8 28.8 s 62 it. 53.2 26.4 s 32 it. 13.1 26.6 s 2.37E-06 1.99
5 160 it. 432.7 116.6 s 91 it. 126.2 99.0 s 37 it. 16.8 105.9 s 5.96E-07 1.99
classical
identity prec. lumped prec. Dirichlet prec. L2 error eoc
1 80 it. 98.2 7.1 s 35 it. 14.1 5.9 s 29 it. 10.0 5.9 s 1.47E-04 -
2 105 it. 161.4 7.8 s 58 it. 41.9 6.1 s 37 it. 16.4 5.8 s 3.72E-05 1.98
3 140 it. 295.7 9.3 s 85 it. 105.9 7.9 s 46 it. 25.4 7.7 s 9.41E-06 1.98
4 188 it. 580.9 15.2 s 125 it. 252.1 13.1 s 54 it. 35.8 12.2 s 2.37E-06 1.99
5 251 it. 1150.3 103.4 s 179 it. 555.7 88.2 s 60 it. 46.3 83.6 s 5.96E-07 1.99
Table 2: Iteration numbers (it.), condition numbers and computational time (in s) for each preconditioning technique using elements; is the level of uniform refinement. For the L2 error the definition is given in (53), while for the estimated error of convergence eoc the definition is given in (54).

From level 4, with a maximum of local degrees of freedom, to level 5, with a maximum of local degrees of freedom, we observe an increase in the local assembling and factorization time from approximately seconds up to about seconds for all kinds of preconditioners. This is mainly due to the higher memory requirements of the direct solver. Note also that the factorization of the local stiffness matrices by the direct solver is unfeasible, if the number of local degrees of freedom gets too large. The reason for that are memory limitations on the VSC2 cluster. A possibility to overcome this problem is the use of fast local iterative solvers, e.g., the CG method with a multigrid or a BPX preconditioner. Summing it up seems that the simple lumped preconditioner and the classical FETI approach appear to be favorable for this academic example, with very structured subdomains and the boundary . The latter yields a large number of floating subdomains for all-floating FETI which are non-floating for the classical FETI approach, and hence a much larger coarse matrix for all-floating FETI. The inversion of this matrix is the most time consuming part for the levels that also results in the higher computational time for all-floating FETI in these cases.

Next, we consider a linear elastic problem by using tetrahedral elements and quadratic ansatz functions, i.e. elements for the same mesh and parameter properties as above. The number of degrees of freedom now varies between (level ) and (level ) and the number of Lagrange multipliers between and . Note that for all preconditioning types and for both the all-floating and the classical FETI method the L2 error compared to the fundamental solution behaves as predicted from the theory as we get a cubic convergence rate, see Table 3.

For all-floating FETI we have the very interesting case that the global CG iteration numbers remain almost constant for the lumped preconditioner, and it even seems to be a decay for the identity and the Dirichlet preconditioner, if we increase the local degrees of freedom, i.e. increase the refinement level .

For the classical FETI approach the iteration numbers stay almost constant for the Dirichlet preconditioner and increase marginally for the other two preconditioning techniques. Concerning the computational time we have an analogous result as in the previous case with linear ansatz functions: the classical approach with the lumped preconditioner seems to be the best choice for this particular example.

all-floating
identity prec. lumped prec. Dirichlet prec. L2 error eoc
1 149 it. 444.7 23.3 s 73 it. 73.7 22.0 s 47 it. 36.7 18.7 s 1.13E-05 -
2 129 it. 330.8 21.9 s 75 it. 74.3 20.8 s 43 it. 27.7 19.3 s 1.44E-06 2.97
3 114 it. 210.3 30.3 s 73 it. 68.8 27.3 s 36 it. 16.6 28.5 s 1.81E-07 2.99
4 105 it. 167.8 99.8 s 69 it. 65.2 93.4 s 33 it. 14.4 90.2 s 2.26E-08 3.00
classical
identity prec. lumped prec. Dirichlet prec. L2 error eoc
1 120 it. 405.0 7.5 s 65 it. 48.9 6.9 s 40 it. 21.0 6.5 s 1.17E-05 -
2 108 it. 302.6 7.5 s 69 it. 57.6 6.7 s 41 it. 20.6 7.5 s 1.46E-06 3.00
3 112 it. 253.4 12.6 s 91 it. 116.2 11.7 s 42 it. 21.0 12.3 s 1.82E-07 3.01
4 136 it. 273.1 76.3 s 128 it. 262.8 77.3 s 48 it. 27.7 79.1 s 2.26E-08 3.01
Table 3: Iteration numbers (it.), condition numbers and computational time (in s) for each preconditioning technique using elements; is the level of uniform refinement. For the L2 error the definition is given in (53), while for the estimated error of convergence eoc the definition is given in (54).
Figure 4: Mesh of an aorta seen from above showing the brachiocephalic artery, and the left common carotid and subclavian arteries. The fine mesh consists of tetrahedrons and vertices, while colors indicate the displacement field with an internal pressure of  mmHg. Additionally, the splits show the decomposition of the mesh into subdomains (left). Coarser mesh consisting of tetrahedrons and vertices used in Section 5.3 with selected vertices A–E (right); colors show the distribution of the stress magnitude according to (56) with an internal pressure of  mmHg. For both images red indicates high and blue low values.
Figure 5: Mesh of a segment of a common carotid artery from two different points of view. The mesh consists of tetrahedrons and vertices. Color indicates the distribution of the stress magnitude according to (56) due to an internal pressure of  mmHg, red indicates high and blue low values. Additionally, the splits show the decomposition of the mesh into subdomains.
Figure 4: Mesh of an aorta seen from above showing the brachiocephalic artery, and the left common carotid and subclavian arteries. The fine mesh consists of tetrahedrons and vertices, while colors indicate the displacement field with an internal pressure of  mmHg. Additionally, the splits show the decomposition of the mesh into subdomains (left). Coarser mesh consisting of tetrahedrons and vertices used in Section 5.3 with selected vertices A–E (right); colors show the distribution of the stress magnitude according to (56) with an internal pressure of  mmHg. For both images red indicates high and blue low values.
Figure 6: Distribution of the stress magnitude inside the aorta (left); values of high stress in red and of low stress in blue. To the right the fiber directions (black curves) and the two layers (adventitia in red and media in orange) of the carotid artery are shown.

5.2 Arterial model on a realistic mesh geometry

In this section we present examples to show the applicability of the FETI approaches for biomechanical applications, in particular the inflation of an artery segment. We consider the mesh of an aorta and the mesh of a common carotid artery, see Figs. 5 and 5. The geometries are from AneuriskWeb [56] and Gmsh [57]. The generation of the volume mesh was performed using VMTK and Gmsh [57].

The fiber directions, see Fig. 6 (right), were calculated using a method described by Bayer et al. [58] for the myocardium. To adapt this method for the artery we first solved the Laplace equation on the domain with homogeneous Dirichlet boundary conditions on the inner surface and inhomogeneous Dirichlet boundary conditions on the outer wall. The gradient of the solution is used to define the transmural direction in each element. As a second step we repeat this procedure using homogeneous Dirichlet boundary conditions on the inlet surface and inhomogeneous boundary conditions on the outlet surfaces which yields the longitudinal direction . The cross product of these two vectors eventually provides the circumferential direction . With a rotation we get the two desired fiber directions and in the media and the adventitia, respectively. Thus,

(55)

The value for the angle are for the media and for the adventitia, taken from [5].

To describe the anisotropic and nonlinear arterial tissue, we use the material model (1317), with the parameters given in Table 1 and is varied. Dirichlet boundary conditions (7) are imposed on the respective intersection areas. We perform an inflation simulation on the artery segment where the interior wall is exposed to a constant pressure . This is performed using Neumann boundary conditions (8). If not stated otherwise, we present the results of one load step applying a rather low pressure of  mmHg. This is necessary to have a converging Newton method. Nonetheless, the material model as used is anisotropic. To simulate a higher pressure, an appropriate load stepping scheme, see (25), has to be used. However, this does not affect the number of local iterations significantly. As already mentioned in Section 4 we use the CG method as global iterative solver. Experiments with a standard non-symmetric nonlinear elasticity system and the necessary GMRES method as an iterative solver showed similar results, as presented in the following with the symmetric system. However, the memory requirements of the GMRES solver are much higher.

The local generalized pseudo-inverse matrices are realized with a sparsity preserving regularization by fixing nodes, see, e.g., [30], and the direct solver package Pardiso. The global nonlinear finite element system is solved by a Newton scheme, where the FETI approach is used in each Newton step. For the considered examples the Newton scheme needed four to six iterations. Due to the non-uniformity of the subdomains the efficiency of a global preconditioner becomes more important. It may happen that the decomposition of a mesh results in subdomains that have only a few points on the Dirichlet boundary. This negatively affects the convergence of the CG method using classical FETI, but does not affect the global iterative method of the all-floating approach at all. This is a major advantage of all-floating FETI since here all subdomains are treated the same, and hence all subdomains are stabilized. This behavior is observed for almost all settings for preconditioners and the penalty parameter as well as for linear and quadratic ansatz functions, see Tables 47.

all-floating
identity preconditioner lumped preconditioner Dirichlet preconditioner
10 1052 it. 57.6 s 160 it. 31.0 s 56 it. 22.8 s
100 1879 it. 94.6 s 305 it. 29.5 s 85 it. 25.4 s
1000 4122 it. 177.1 s 681 it. 48.8 s 209 it. 31.8 s
classical
identity preconditioner lumped preconditioner Dirichlet preconditioner
10 2056 it. 98.7 s 305 it. 35.5 s 117 it. 27.2 s
100 3711 it. 149.8 s 540 it. 35.5 s 144 it. 28.4 s
1000 8245 it. 327.8 s 1190 it. 60.9 s 263 it. 32.9 s
Table 4: Iteration numbers (it.) per Newton step and computational time (in s) per Newton step for the all-floating and the classical FETI approach with linear ansatz functions comparing the three considered preconditioners. The penalty parameter was varied from to  kPa. Mesh: mesh of the aorta subdivided in subdomains, computed with cores.
type identity preconditioner lumped preconditioner Dirichlet preconditioner
all-floating 10000 it. - s 1084 it. 100.6 s 497 it. 85.5 s
classical 5130 it. 357 s 1794 it. 200.2 s 588 it. 97.7 s
Table 5: Iteration numbers (it.) per Newton step and computational time (in s) per Newton step for the all-floating and the classical FETI approach with linear ansatz functions comparing the three considered preconditioners. The penalty parameter was set to  kPa. Mesh: mesh of the carotid artery with two layers (adventitia and media) subdivided in subdomains, computed with cores.
all-floating
identity preconditioner lumped preconditioner Dirichlet preconditioner
10 940 it. 491.1 s 283 it. 209.5 s 71 it. 157.3 s
100 1519 it. 1186.4 s 523 it. 332.0 s 105 it. 178.1 s
1000 3371 it. 2584.5 s 1372 it. 746.0 s 206 it. 282.7 s
classical
identity preconditioner lumped preconditioner Dirichlet preconditioner
10 1319 it. 654.2 s 333 it. 225.2 s 113 it. 188.4 s
100 2362 it. 1140.6 s 664 it. 402.6 s 110 it. 177.5 s
1000 5563 it. 4168.3 s 1742 it. 943.1 s 204 it. 280.1 s
Table 6: Iteration numbers (it.) per Newton step and computational time (in s) per Newton step for the all-floating and the classical FETI approach with quadratic ansatz functions comparing the three considered preconditioners. The penalty parameter was varied from to  kPa. Mesh: mesh of the aorta subdivided in subdomains, computed with cores.
type identity preconditioner lumped preconditioner Dirichlet preconditioner
all-floating 10000 it. - s 2163 it. 1133.9 s 674 it. 994.6 s
classical 6006 it. 2672.6 s 4798 it. 2306.8 s