An accurate, robust, and efficient finite element framework with applications to anisotropic, nearly and fully incompressible elasticity

Fiber-reinforced soft biological tissues are typically modeled as hyperelastic, anisotropic, and nearly incompressible materials. To enforce incompressibility a multiplicative split of the deformation gradient into a volumetric and an isochoric part is a very common approach. However, the finite element analysis of such problems often suffers from severe volumetric locking effects and numerical instabilities. In this paper, we present novel methods to overcome volumetric locking phenomena for using stabilized P1–P1 elements. We introduce different stabilization techniques and demonstrate the high robustness and computational efficiency of the chosen methods. In two benchmark problems from the literature as well as an advanced application to cardiac electromechanics, we compare the approach to standard linear elements and show the accuracy and versatility of the methods to simulate anisotropic, nearly and fully incompressible materials. We demonstrate the potential of this numerical framework to accelerate accurate simulations of biological tissues to the extent of enabling patient-specific parameterization studies, where numerous forward simulations are required.


Introduction
Computer models of biological tissues, e.g., the simulation of vessel mechanics or cardiac electro-mechanics (EM), aid in understanding the biomechanical function of the organs and show promise to be a powerful tool for predicting therapeutic responses. Advanced applications include simulations to assess passive filling properties and active response to pacing therapies [1], simulations of growth and remodeling processes occurring in the failing heart or arteries [2][3][4][5], as well as rupture risk assessment in arterial aneurysms [6,7]. Here, predictions of in-silico models are often based on the computation of local stresses, hence, an accurate computation of strain and stress is indispensable to achieve mechanistically sound predictions. To build confidence in simulation outcomes proper model calibration, validation, uncertainty quantification, and sensitivity analyses are required. Additionally, high computational efficiency and excellent strong scaling properties are crucial to perform simulations with highlyresolved, complex, or heterogeneous geometries within tractable time frames; to facilitate model personalization using a large number of forward simulations; and to simulate tissue behavior over a broad range of experimental protocols and extended observation periods.
In-silico models of cardiac tissue and vessel walls are typically based on the theory of hyperelasticity and properties of soft tissues include a nonlinear relationship between stress and strain with large deformations and a nearly-incompressible, anisotropic materials [8][9][10]. Commonly, the resulting non-linear formulations are approximately solved using a finite element (FE) approach [11][12][13][14]. However, volumetric locking phenomena, that are resulting in ill-conditioned global stiffness matrices are frequently encountered. Here and in the following, we use the definition of locking given in the seminal work of Babuška and Suri [15], meaning the dependance of convergence of a FE solution on a critical parameter -in our applications the bulk modulus κ -and the lack thereof in the limit. In fact, this is one of the classical problems of modeling nearly incompressible hyperelasticity [16,17]. Volumetric locking, often completely invalidating FE solutions, is in particular prevalent for fiber-reinforced soft biological tissues due to a high stiffness in the preferred fiber directions compared to a relatively soft matrix material [18] and is thus extensively studied in recent publications [19][20][21][22].
Typically, the modeling of (nearly) incompressible elastic materials involves a split of the deformation gradient into a volumetric and an isochoric part [23]. Here, volumetric locking phenomena are very common when using unstable approximation pairs such as Q1-P0 or P1-P0 elements, i.e., when linear shape functions are the choice to approximate the displacement field u and the hydrostatic pressure p is statically condensed from the system of equations on the element level. It is well known that in such cases solution algorithms are likely to show very low convergence rates, and that variables of interest such as stresses can be inaccurate [20].
To some degree volumetric locking problems in anisotropic hyperelasticity for these simple elements can be reduced by using augmented Lagrangian methods [24,25], formulations with an unsplit deformation gradient for the anisotropic contribution [26,27], and methods with simplified kinematics for the anisotropic contributions [22]. Another possibility to obtain more accurate results is the use of higher order polynomials to approximate the displacement [28][29][30][31]. However, the incompressibility constraint is still modeled by a penalty formulation, hence, volumetric locking may still be an issue and the modeling of fully incompressible materials is not possible. Additionally, already for quadratic ansatz functions the considerable larger amount of degrees of freedom increases computational cost significantly.
A more sophisticated approach -also allowing the modeling of fully incompressible materials -is the reformulation of the underlying equations into a saddle point problem by introducing the hydrostatic pressure p as an additional unknown to the system. Here, from mathematical theory, approximation pairs for u and p have to fulfill the Ladyzhenskaya-Babuŝka-Brezzi (LBB) or inf-sup conditions [32][33][34] to guarantee stability. A popular choice are quadratic ansatz functions for the displacement and linear ansatz functions for the pressure, i.e., the Taylor-Hood element [35,36]. Though stable, this element leads to a vast increase in degrees of freedom and the development of scalable solvers is very difficult. Consequently, this entails a high computational burden; especially for applied problems in the field of tissue mechanics with highly resolved geometries.
A computationally potentially more favorable choice are equal order pairs with a stabilization, widely used for linear and isotropic elasticity [37][38][39][40][41][42][43]. Yet, their extension to non-linear problems is challenging [22,44,45]. In the specific case of modeling biological tissues Hu-Washizu-based formulations are often used, e.g., [46][47][48][49]. However, especially for problems undergoing large strains, this mixed three field approach shows limited performance and robustness [22]. Further, variational multiscale (VMS) formulations have been used for nonlinear solid mechanics [39,50,51]. VMS is frequently used in computational fluid mechanics with the key advantage that these formulations offer both stabilization as well as benefits in turbulence modeling. The main idea behind VMS is to split FE solutions into resolvable and unresolvable parts, where the unresolved parts of the solution are modeled with the residual of the physics. On the other hand, VMS formulations are introduced at the discrete FE level, and thus one usually obtains modified variational formulations that include mesh-depended parameters. For solid mechanics, introducing VMS results in a modified FE formulation that can encode anisotropic material behavior very well [39]. However, this comes at the cost of a very intrusive change into the formulation as the consistent linearization introduces non-standard 6th order tensors which need to be implemented into existing code frameworks.
A very promising and efficient stabilization approach for nearly incompressible elasticity problems is a variant of the MINI element [52], originally established for computational fluid dynamics problems. This element is modified for the application of incompressible hyperelasticity and a bubble function is included in the ansatz space of displacements. To improve efficiency, the support of this bubble can be eliminated from the system of equation using static condensation. First uses of MINI elements have been reported [53,54] though still using a piecewise constant ansatz for the hydrostatic pressure. Even more efficient and notably simpler to implement is a pressure projection method originally introduced for the Stokes problem [55]. To the best of our knowledge the here proposed methods were not yet applied in this form to anisotropic and nearly incompressible materials.
A significant advantage of both stabilization techniques is that these do not rely on artificial stabilization parameters that may influence the numerical solution. We illustrate in different benchmarks that the same setting can be used for a large variety of tissue mechanics problems allowing for a one-for-all approach. By comparing to other methods previously reported in the literature, we show that our methods are suitable to compute accurate strain and stress fields and outperform existing contributions in terms of efficiency.
The paper is outlined as follows: in Section 2 we recall the mathematical background of modeling anisotropic, nearly incompressible elasticity and introduce the theoretical framework of our stabilization techniques. Subsequently, Section 3 documents three benchmarks problems to show the applicability of the stabilized P1-P1 elements in different scenarios. For each benchmark we give a detailed problem description and discuss results and computational efficiency by comparing to the literature and analyzing strong-scaling properties.
To show the usefulness of the presented methods to clinically relevant problems, we present a 3D EM model of the heart that is coupled to a 0D model of the circulatory system. This constitutes one of the most complete model of cardiac EM in the literature to date as all components, i.e., electrophysiology, cellular dynamics, active stress, passive tissue mechanics, pre-and afterload, are based on physiological, state-of-the-art models. Finally, Section 4 concludes the paper with a brief summary and all required equations to implement the methods in a software framework are given in Appendix A.

Almost incompressible nonlinear elasticity
Let Ω 0 ⊂ R 3 denote the reference configuration and let Ω t ⊂ R 3 denote the current configuration of the domain of interest. Assume that the boundary of Ω 0 is decomposed into ∂Ω 0 = Γ D,0 ∪ Γ N,0 with |Γ D,0 | > 0. Here, Γ D,0 describes the Dirichlet part of the boundary and Γ N,0 describes the Neumann part of the boundary, respectively. Further, let N be the unit outward normal on ∂Ω 0 . The nonlinear mapping φ : X ∈ Ω 0 → x ∈ Ω t , defined by φ := X + u(X, t), with displacement u, maps points in the reference configuration to points in the current configuration. Following standard notation, we introduce the deformation gradient F, the Jacobian J , and the left Cauchy-Green tensor C as F := Grad φ = I + Grad u, J := det(F), C := F ⊤ F.
Here, Grad(•) denotes the gradient with respect to the reference coordinates X ∈ Ω 0 . The displacement field u is sought as infimizer of the functional over all admissible fields u with u = g D on Γ D,0 , where, Ψ denotes the strain energy function; ρ 0 denotes the material density in reference configuration; f denotes a volumetric body force; g D denotes a given boundary displacement; and h denotes a given follower load [56,57] defined as with given external constant surface pressure p ext ∈ R + . Note that in general p ext is assumed constant over Γ N,0 but may vary over time. For ease of presentation it is assumed that ρ 0 is constant and f , and g D do not depend on u. For changes of the presented formulation to model transient simulations including inertial terms see [58]. Further, for a discussion of the existence of infimizers in the case of follower loads see the pioneering works of Ball [56],Ciarlet [57].
In this study, we consider nearly incompressible materials, meaning that J ≈ 1. A possibility to model this behavior was originally proposed by Flory [23] using a split of the deformation gradient F -hence referred to as Flory split -such that Here, F vol describes the volumetric change while F describes the isochoric change. By setting F vol := J 1 3 I and F := J − 1 3 F we get det(F) = 1 and det(F vol ) = J . Analogously, by setting C vol := J 2 3 I and C := J − 2 3 C, we have C = C vol C. Assuming a hyperelastic material, the Flory split also postulates an additive decomposition of the strain energy function The function U (J ) will be used in the form where κ > 0 denotes the bulk modulus. In the literature many different choices for the functions Θ(J ) are proposed, see e.g [59][60][61] for examples and related discussion. We want to emphasize that for all choices of Θ it holds that Θ(1) = 0 and Θ ′ (1) = 1. For studying also the limit case κ → ∞ we will consider a reformulation of Eq. (1) as perturbed Lagrangian-multiplier (PL) functional, see [62][63][64][65][66][67] as well as [68,Chapter 8] for details. Introducing the pressure p we seek stationary points (u, The above PL formulation is a special case assuming the form (4) for U (J ), see also [63] and [39,Section 3.3]. Variable p is not the physical hydrostatic pressure P hydro but corresponds to the applied hydrostatic pressure [69,Chapter 7,(7.4.44)]. The two pressures are linked by showing that in the incompressible limit the absolute values of both coincide. To guarantee well-definedness, we assume that with H 1 (Ω 0 ) being the standard Sobolev space of square integrable functions having a square integrable gradient, and Q = L 2 (Ω 0 ). For a more in-depth discussion we refer to [56,57]. To solve for stationary points of Eq. (5) we calculate the variations with respect to test functions v and q. This results in the following non-linear variational problem, find (u, p) ∈ V g D × Q such that R inc (u, p; q) = 0, for all (v, q) ∈ V 0 × Q. Here, where a isc (u; v) := Components of the second Piola-Kirchhoff stress tensor are computed as with Dev(•) being the deviatoric operator in the Lagrangian description, see Eq. (C.1). When modeling electrically active tissue, we consider an additive decomposition of the isochoric part of the stress tensor. The total stress tensor is now given by the additive decomposition with the active stress S a defined as in [70].
To simulate the effect of the circulatory system, these equations are coupled to a 0D lumped model as in [71]. The corresponding nonlinear variational problem reads as find (u, p) ∈ V g D × Q and p CAV ∈ R n CAV such that for all (v, q) ∈ V 0 × Q, and i = 1, . . . , n CAV , where p CAV = { p CAV,1 , p CAV,2 , . . . , p CAV,n CAV } denote the cavity pressures, and n CAV denotes the number of cavities in the model. Here, the variations read as assuming summation over double occurring indices in (13), where with Γ CAV,i,0 denoting the closed surface of the ith cavity in reference configuration. The expression for cavity volume V CAV,i follows from applying Nanson's formula to the definition of cavity volume in the current configuration For a more detailed account on the coupling of nonlinear elastic equations with 0D lumped parameter models, and the computation of V CS , the volume as predicted by the 0D model of the circulatory system for the intra-cavitary pressure p CAV,i , we refer to [71].

Consistent linearization
In the subsequent section we will use • k to indicate a quantity that depends on a previous Newton iterate. For the subsequent discretization we need the consistent linearization of (6)-(7) and (10)- (12) and we obtain the following linear saddle-point problem: for each (u k , p k ) ∈ V g D × Q, find increments (∆u, ∆ p) ∈ V 0 × Q such that where a k (∆u, ∆v) := using the following abbreviations and C isc given in (C.3). For the deviation of term (17) see Appendix A, other terms in (15)- (16) have been discussed previously, see [58]. Note, that by choosing U (J ) as 1 2 Θ 2 (J ) we automatically obtain -in the absence of follower loads -a symmetric system, similar to [66]. However, due to the PL formulation we can avoid issues when U ′′ (J ) = 0.
In the case of an attached circulatory system we obtain the following linearized system: find increments where e CAV,i ( p k CAV,i ; ∆ p CAV ) := and d CAV defined as in (A.3). The term (21) depends on the chosen model for the circulatory system and a detailed discussion is out of the scope of this work. For a detailed derivation of the explicit representation of the compliance matrix (21) stemming from the model used in Section 3 we refer to [71].

Finite element discretization
Here we provide a summary of the FE discretization used in the subsequent results. The framework builds upon methods previously introduced for isotropic, passive mechanics in [58]. In the following, we extend this approach to anisotropic tissues also allowing for complex EM simulations that are coupled to a 0D system of the circulatory system.
Let T h be a FE partitioning of Ω consisting of isoparametric tetrahedral and/or hexahedral elements. We assume standard regularity assumptions [57] and invertibility of the isoparametric mapping F K from the reference elementK to a physical element K ∈ T h . For tetrahedral elements this poses no additional restrictions, for hexahedral elements we refer to [72] for details. Let furtherP 1 andQ 1 denote the space of lowest order linear/trilinear FE functions on the reference tetrahedron/hexahedron. The discrete analogue to (15)-(16) reads as: with V and Q denoting suitable spaces of polynomials which will be specified in the subsequent sections.

Pressure projection stabilized equal order pair
The pressure projection stabilization was originally introduced for solving Stokes problems [55] and has also been applied in the context of linear elasticity [73,74]. Recently, we extended its use to isotropic, nonlinear elasticity [58]. A similar approach can be used for anisotropic materials, we setV :=Q :=P 1 /Q 1 for tetrahedral or hexahedral elements. To ensure stability, we have to modify the definition of the residuals in (7) and (11) tõ where the projection operator Π h is defined elementwise as In contrast to [58], the parameter µ * is no longer an arbitrary value but set to |K | 1/3 ; a choice that showed excellent results for all discussed anisotropic problems in Section 3 as well as isotropic benchmarks in [58]. We note, that the integral in (24) has to be understood as a sum over all the elements of the triangulation of the domain Ω 0 . For a more comprehensive overview and implementation details we refer to [58].

MINI element
One of the earliest strategies in constructing a stable FE pairing for discrete saddle-point problems arising from Stokes equations is the MINI element, dating back to the works of Brezzi et al. see for example [52,75]. Briefly, the strategy is to enrich the basis of lowest order elements by adding a higher degree polynomial with support restricted to the interior of the element. Thus, for the tetrahedral reference elementK we definê where (ξ, η, ζ ) ∈K see also [76]. For the hexahedral reference elementK = [−1, 1] 3 we definê for (ξ, η, ζ ) ∈K and (α, β) ∈ [1,8] denoting the indices of two ansatz functions for diagonal opposite nodes in K , see [58]. The choice of the two numbers α, β depends on the nodal ordering of the reference hexahedron. A visualization of the ansatz functions used in this work is given in Figs. 1 and 2.
Classical results [76] guarantee the stability of the MINI element for tetrahedral meshes in the almost incompressible linear elastic case. For hexahedral elements we were able to prove stability in the almost incompressible linear elasticity case provided an enrichment like (25) of the displacement ansatz space by two bubble functions see [58]. Due to the compact support of the bubble functions, static condensation can be applied to remove the interior degrees of freedom during assembly. Static condensation can be done by standard procedures [76] with the  exception of follower loads which is discussed in Appendix B. As a result, these degrees of freedom are not needed to be considered in the full global stiffness matrix assembly which is a key advantage of the MINI element.

Material models
Arterial and myocardial tissue as modeled in Section 3 is considered as a non-linear, hyperelastic, nearly incompressible, and anisotropic material with a layered organization of fibers. To model this behavior in our benchmark problems we used strain energy functions of the form (3). First, single Fung-type exponential models of the form with a > 0 kPa being a scaling parameter, a i are fiber directions, and Q is a function in terms of scalar strain components. The specific form of U (J ) and Q will be discussed later in the benchmark section. Second, separated, invariant-based Fung-type exponential models, also often referred to as Holzapfel-Ogden-type models. Here, we compared the standard formulation -using anisotropic splitting (AS) - with the formulation using an unsplit deformation gradient for the anisotropic contribution -without anisotropic split (WAS) - introduced in [77,78]. The specific form of the volumetric, U (J ), isotropic, Ψ iso , and anisotropic, Ψ aniso /Ψ aniso , contributions will be discussed later for each of the benchmark problems. Note that results on existence of solutions in nonlinear elasticity as well as stability of fiber-reinforced materials are highly dependent on polyconvex strain-energy functions. This is not guaranteed for above mentioned models, however, we use the same models and parameters for our benchmark problems as reported in previous works, i.e., the setup from Gültekin et al. [79] for the benchmark in Section 3.1 and the setup from Land et al. [80] for the benchmark in Section 3.2. For the third numerical example in Section 3.3 we use a state-of-the-art Holzapfel-Ogden type material which are most commonly used to describe passive tissue mechanics.
For more details regarding the existence of solutions and the stability of the material models we refer to [56,57,78,81,82].

Numerical examples
Biomechanical applications often require highly resolved meshes and thus efficient and massively parallel solution algorithms for the linearized system of equations become an important factor to deal with the resulting computational load. Extending our previous implementations for cardiac electro-mechanics (EM) [83] we implemented the stabilization techniques in the FE framework Cardiac Arrhythmia Research Package (CARPentry) [84,85], built upon extensions of the openCARP EP framework [86] (http://www.opencarp.org). We solve the stabilized saddlepoint problem (22)-(23) by using a GMRES method with a block preconditioner based on a smoothed aggregation algebraic multigrid (GAMG) approach which is incorporated in PETSc [87,88].
In all of the following benchmark problems our goal was to study the performance and accuracy of different FE discretizations, namely (i) Q1/P1-P0-AS: discretization with piecewise linear displacements and piecewise constant pressure using the strain energy function (27); (ii) Q1/P1-P0-WAS: discretization with piecewise linear displacement and piecewise constant pressure using the strain energy function (28); (iii) Projection: equal order discretization with piecewise linear displacements and pressure, stabilized as described in Section 2.3.1 using the strain energy function (27); (iv) MINI: discretization using MINI elements as described in Section 2.3.2 using the strain energy function (27).
Meshes for Sections 3.1 and 3.2 were created using Gmsh [89]. Tetrahedral meshes were created as tetrahedralization of the hexahedral mesh and we use ℓ to index the various uniform refinement levels of the meshes.

Extension, inflation and torsion of a simplified artery model
Simulation setup. First, we show the applicability of our proposed methods to a benchmark problem from Gültekin et al. [79] where a simplified artery model is represented by a thick-walled cylindrical tube. The dimensions of this idealized geometry with its centerline on the z-axis are as follows: height H = 10 mm, inner radius R 1 = 8 mm, and outer radius R 2 = 10 mm. Two symmetric families of fibers, f 0 and s 0 are immersed in the tissue, having an angle of 40 • with circumferential θ-axis, see Fig. 3A. As for loading, a monotonically increasing displacement up to 2 mm superimposed by a monotonically increasing torsion up to 60 • is applied on the top of the tube (marked blue in Fig. 3B). Additionally, a linearly increasing pressure (follower load) up to 500 mmHg is applied on the inside of the tube (marked red in Fig. 3B). Finally, the lower part of the tube is clamped at zero displacement. The material is described by the strain-energy function (27), Ψ AS , with and invariants and analogously with Ψ aniso (C, f 0 , s 0 ) for the WAS formulation, Ψ WAS , (28). Material parameters were taken from [79, Table 1], i.e., κ = 5000 kPa, µ = 10 kPa, k 1 = 500 kPa, and k 2 = 2.0. In case of the stabilized equal-order elements (Projection and MINI), we set 1/κ = 0 to render the material incompressible. To assess mesh convergence simulations were performed on seven discretization levels, see Table 1.

Results.
A comparison of the radial, σ rr , the circumferential, σ θ θ , and the axial, σ zz , components of the Cauchy stress tensor is shown in Fig. 4 for the finest discretization level ℓ = 7. We see that with the exception of the lowest order discretizations with anisotropic splitting (Q1/P1-P0-AS) the stress distribution is very similar and also matches results in Gültekin et al. [79]. The observation that simulations with Q1/P1-P0-AS are not accurate for this benchmark problem is further emphasized in Figs. 5 and 6. Here, Fig. 5 shows the displacements (u x , u y , u z ) and Fig. 6 shows the stress components (σ rr , σ θ θ , σ zz ) at the evaluation points A and B over the discretization levels. In agreement with [79] the lowest order discretization with anisotropic splitting converges to a lower value than the other discretization types hinting at possible locking phenomena. All other formulations perform similarly well. Here, discrepancies at the finest level ℓ = 7 are rather due to differences in the meshes for tetrahedral and hexahedral grids. Fig. 7 shows a distribution of the Jacobian det(F) on the finest level ℓ = 7. Unsurprisingly, with a mean value µ close to 1 the saddle-point formulations (Projection, MINI) satisfy incompressibility better than the penalty formulations (P1/Q1-P0-AS, P1/Q1-P0-WAS). While the AS formulation led to a small increase in volume (µ > 1), the WAS formulation resulted in a slightly reduced volume (µ < 1).
Numerical performance. Computational times for the simulation using different element types are given in Fig. 8; left, for the coarse problem (ℓ = 1) and right, for the finest grid (ℓ = 7). For all cases we used a relative error reduction of ϵ = 10 −8 for the GMRES linear solver and a relative error reduction of ϵ = 10 −6 for the residual of the Newton method. With a load stepping scheme the applied traction and displacement was increased to arrive at the final prescribed displacement and an inner pressure of 500 mmHg. Using 25 loading steps for the coarse problem required a total number of 100 linear solving steps. For the fine problem we chose 250 loading steps and required 1000 linear solving steps. Strong scaling was achieved for the coarse problem up to 16 cores on a desktop machine (AMD Ryzen Threadripper 2990X), see Fig. 8, left. Here, computational times were averaged over 5 runs Discussion. In accordance with Gültekin et al. [79] we have shown for this benchmark that the concept Q1/P1-P0-WAS is able to match quasi-incompressible responses compared to the gold standard of locking-free elements. For this extreme loading case standard Q1/P1-P0-AS elements cannot reproduce accurate stress distributions; not even on very fine grids, see Fig. 6(b).
Further, this benchmark highlights the high efficiency of our methods. Computational times of the benchmark as presented in [79] -on a 3.0 GHz CPU unit and the same hexahedral mesh -were 11 s for Q1-P0-WAS elements and 24 s, i.e., only 2 to 3 times slower, for locking-free Hex Projection elements, see Fig. 8, left. Similarly, for tetrahedral elements as well as different refinement levels, simulations are only to a relatively small factor slower when using locking-free elements (Fig. 8). Further, as illustrated in this figure, strong scaling properties were excellent for both: computations using the coarsest grid on a desktop machine, as well as computations using the finest grid on a supercomputer. Due to a higher number of linear iterations and higher matrix assembly times, simulations with hexahedral meshes were more expensive compared to simulations with tetrahedral grids. However, as also observed, e.g., by Chamberland et al. [90] hexahedral elements were slightly more accurate than their tetrahedral equivalent. In realworld applications the minor gain in accuracy has to be balanced with the higher cost of hex elements in terms of computing and mesh generation when geometrically detailed anatomical models are used.

Inflation and active contraction of a simplified ventricle
Simulation setup. To verify our EM setup we repeated the inflation and active contraction benchmark from Land et al. [80]. We generated the reference geometry of an idealized ventricle as a tetrahedral mesh of a truncated ellipsoid and prescribed a local orthonormal coordinate system with fiber, f 0 , sheet, s 0 , and sheet-normal, n 0 , directions according to this paper. We constructed three levels of refinement, see Table 2 for discretization details.
As material we used the transversely isotropic law by Guccione et al. [91], based on Eq. (26), with Θ(J ) := ln(J ) and Constitutive parameters were a = 2 kPa, b f = 8, b t = 2, and b fs = 4. In the benchmark paper the material is considered to be fully incompressible, hence, we chose 1/κ = 0 for the saddle-point formulation (Projection,  MINI). For the penalty formulation (P1-P0 elements) we chose κ = 1000 kPa which was the best trade-off between convergence of the solver for all three levels in Table 2, near incompressibility, and minimization of locking effects.  As the material law above does not allow for a WAS formulation we repeated the benchmark using a separated, invariant-based Fung-type exponential model as in Eq. (27). In particular, we chose a Holzapfel-Ogden material [92] of the form with invariants such that contributions of compressed fibers are excluded, and the interaction-invariant Analogously, we used the constitutive equation above with Ψ aniso (C, f 0 , s 0 , n 0 ) for the WAS formulation. Material parameters were taken from [93], a = 0.809 kPa, b = 7.474, a f = 1.911 kPa, b f = 22.063, a n = 0.227 kPa, b n = 34.802, a fs = 0.547 kPa, and b fs = 5.691, fitted to human myocardial experiments in [94].
Results. For the transversely isotropic law (26), we compared our results to selected reference solutions from the benchmark paper [80], namely, the result from IBM with the Cardioid framework [95] using P2-P1 elements and the result from Simula with FEniCS [96] using two-dimensional P2-P1 elements. First, the final location of the apex is measured and, second, circumferential, longitudinal, and radial strains at the endocardium, epicardium, and midwall are calculated on points along apex-to-base lines, see [80] for more details. Results show that the apex location Fig. 9(a) and strains Fig. 10 are very similar for the finest level (ℓ = 3) for all chosen element types. For the level ℓ = 2 the strain solution using simple P1-P0 elements is not converged showing differences to the benchmark solutions especially in boundary regions at the apex (p1) and the base (p10), see Fig. 10(a). We repeated simulations as above measuring the final apex location and calculating strains along apex-to-base lines using the orthotropic law (29). We compared results using P1-P0-WAS, P1-P0-AS, projection-stabilized, and MINI elements in Fig. 9(b) and Fig. 11. Apex locations are very similar for all element types, however, strains are different, especially in boundary regions close to the apex and the base. We can see in Fig. 11 that even for the Fig. 10. Ellipsoid benchmark, Guccione material: longitudinal (LONG), circumferential (CIRC), and radial (TRANS) strains at endocardium, epicardium, and midwall. Index of points increases from the apex to the base. Own results (in black) are compared to benchmark results (in gray) presented in [80]. Discussion. For the transversely-isotropic Guccione material model strains with the presented projection-stabilized and MINI elements match results using higher order P2-P1 elements even on coarser grids. Here, also linear tetrahedral elements seem to be accurate given a fine enough discretization; this behavior was also observed in [80].
In contrast, for the orthotropic Holzapfel-Ogden material, we see in Fig. 11 that P1-P0 elements cannot always accurately reproduce strains even on the finest grid. On the other hand we can assume that strains are accurate and almost converged for projection-stabilized and MINI elements as results for levels 2 and 3 are very similar. Interestingly, for this benchmark, we see no difference in accuracy between the standard P1-P0 and the P1-P0-WAS formulation. Most likely the reason for this is that parameters fitted to human myocardial data [93] are not as stiff in fiber direction compared to the more extreme artificial benchmark case in Section 3.1. Overall, both approaches with simple linear elements fail to match results from gold-standard elements, especially in boundary regions.

3D-0D closed-loop model of the heart and circulation
Simulation setup. Finally, we show the applicability of our method to an advanced model of computational cardiac EM. Here, a 3D model of bi-ventricular EM is coupled to the physiologically comprehensive 0D CircAdapt model representing atrial mechanics and closed-loop circulation, see Fig. 12. In the present paper, the myocardium of the ventricles was modeled as a nonlinear hyperelastic, (nearly) incompressible and orthotropic material as in Eq. (27). In particular, for this application, we chose the model proposed by Gültekin et al. [97] Θ(J ) := ln(J ), with modified unimodular fourth-invariants to support dispersion of fibers , i ∈ f, s and standard invariants I 1 := tr(C), I 8fs = f 0 · C s 0 .
Analogously, we used the constitutive equation above with Ψ aniso (C, f 0 , s 0 ) for the WAS formulation.
A reaction-eikonal model [85] was used to generate electrical activation sequences. Cellular dynamics were described by the Grandi-Pasqualini-Bers model [98] of human ventricular electrophysiology coupled to the Land-Niederer model [99] of human active stress generation to account for length and velocity dependence of active stress generation, see also Augustin et al. [83] for more details on this strong coupling as well as Regazzoni and Quarteroni [100] for implementation details on the velocity-dependent active stress model. The active stress tensor is computed according to [101] as where S a is the scalar valued active stress generated by the cardiac myocytes and κ f is the same dispersion parameter as above.
For the time integration of Cauchy's equation of motion we used a variant of the generalized-α integrator [102] with spectral radius ρ ∞ = 0 and damping parameters β mass = 0.1 ms −1 , β stiff = 0.1 ms.
Bi-ventricular finite element models. The bi-ventricular geometry was created according to [71] with an average spatial resolution of 1.3 mm for the LV and 1.2 mm for the RV. The resulting mesh used for simulations consisted of 557 316 tetrahedral elements and 111 234 nodes. Fiber and sheet directions were computed by a rule-based method [103] with fiber angles changing linearly from −60 • at the epicardium to +60 • at the endocardium [104]. The valve openings were closed by solid lids as described in [71] to generate the closed surfaces as needed for Eq. (14).
Boundary conditions. Boundary conditions on the epicardium were modeled using spatially varying normal Robin boundary conditions [105] to simulate in-vivo constraints imposed by the pericardium. The basal cut plane was constrained by omni-directional spring type boundary conditions. The 3D ventricular PDE model was coupled to Red colors indicate oxygenated and blue colors de-oxygenated blood. For more information on the 3D-0D coupling see [71]. the 0D ODE model CircAdapt [106] representing cardiovascular system dynamics according to Augustin et al. [71]. An ex-vivo setting without pericardial boundary conditions is used to calibrate parameters replicating ex-vivo passive inflation experiments by Klotz et al. [107].
Parameterization. Passive material parameters a = 0.4, b = 6.55, a f = 3.05, b f = 29.05, a s = 1.25, b s = 36.65, a fs = 0.15, and b fs = 6.28 were taken from [97]. Dispersion parameters have been identified previously by mechanical experiments on passive cardiac tissue by Sommer et al. [94] and are set to κ f = 0.08 and κ s = 0.09.
To eliminate potential differences in stress/strain due to parameterization we fitted parameters to achieve similar pressure-volume (PV) loops and a similar end-diastolic PV relationship (EDPVR) for all element types. First, initial passive material parameters above were fitted to the empiric description of the EDPVR by Klotz et al. [107]. For each element type we used a backward displacement algorithm and boundary conditions replicating experiments in [107] according to Marx et al. [108], see Fig. 13. This fitting resulted in multiplicative scaling factors of 0.4529 (P1-P0 elements) and 0.9582 (locking-free elements) for the stress-like material parameters (a, a f , a s , a fs ); and in multiplicative scaling factors of 1.0322 (P1-P0 elements) and 0.7981 (locking-free elements) for the dimensionless parameters (b, b f , b s , b fs ). The list of fitted passive parameters is given in Table 3.
Active stress parameters were fitted using locking-free elements to reach a target peak pressure of 105 mmHg in the LV. Using the same active stress parameters, the simulations with P1-P0 elements resulted in slightly higher peak pressure, 109.2 mmHg, see Fig. 15 (a). The end-diastolic state remained almost unchanged while the ejection fraction increased slightly. We attribute this to a higher contractility of the elements when the tissue is not modeled as a fully incompressible continuum. To achieve similar PV loops for P1-P0 elements, simulations were repeated with reduced active tension T • ref . See Table 3 for a summary of all active stress parameters. Simulations parameters of the circulatory system were set as in Augustin et al. [71] with a cycle length of 0.585 s. Fig. 13 (a) show that the passive parameterization -that was performed individually for all element types -allowed to reach the predicted stress-free volume and the given end-diastolic Fig. 13. 3D-0D model of the heart: Passive model calibration and unloading following Marx et al. [108]. (a) An ex-vivo setup was used to calibrate the passive material parameters to fit the EDPVR of Klotz et al. [107]; (b) Calibrated as well as default material parameters, see Table 3, were used in an in-vivo setup to generate prestress using a backward-displacement scheme. Shown are the final loading curves. Note that loading curves for P1-P0-WAS and P1-P0-AS coincide and are thus only given once. The same applies for projection-stabilized and MINI elements. The marker "A" denotes the time-point for stress plots in Figs. 16 and 18.  Table 3. (a) Converged left-ventricular PV loops. While the choice of κ has a high influence on the outcome when using P1-P0 elements it is negligible when using the proposed stabilized elements. Note that PV loops for P1-P0-WAS and P1-P0-AS coincide and are thus only given once. The same applies for projection-stabilized and MINI elements. (c) Change in tissue volume over the last two beats. As expected, higher values of κ lead to less myocardial volume change. The choice of 1/κ = 0 for projection-stabilized and MINI elements renders the material incompressible and the myocardial volume does not deviate from its initial value. volume almost perfectly while reproducing the shape of the Klotz EDPVR curve. Boundary conditions for this experiment correspond to the ex-vivo setting described above replicating passive inflation experiments [107]. With fitted material parameters we repeat the backward displacement algorithm to get a stress-free reference geometry with in-vivo boundary conditions to model the constraints imposed by the pericardium. Loading this stress-free configuration to end-diastolic pressure results in loading curves as shown in Fig. 13 (b). When using fitted parameters, these loading curves are almost identical for all element types. Comparing to results obtained using default parameters from Gültekin et al. [97], we get quite similar loading curves using stabilized elements. This is expected as the values of the fitted parameters, see Table 3, are close to the default values obtained from triaxial shear tests [97]. Further, we can see that the choice of κ has little influence on the loading for stabilized elements: using κ = 650 kPa or 1/κ = 0 results in coinciding curves. On the other hand, for P1-P0 elements the choice of κ has a clear impact: in line with expectations, higher values of κ result in a stiffer behavior of the tissue. Looking at stresses after the initial loading phase, marked by the letter "A" in Fig. 13 (b), we also see that for stabilized elements fitted and default parameters produce similar results and also the influence of the choice of κ is negligible, see Figs. 16 and 18. For P1-P0 elements, however, the stress fields are generally less smooth compared to stabilized elements and the dependency on κ is clearly visible.

Results. First, simulation results in
The pre-stressed configuration at marker "A" matches the geometry obtained from imaging and serves as the starting point for EM heart beat experiments as described in the following. To get to a converged solution of the closed-loop 3D-0D system we simulated 30 heart beats for each FE setting: 18 initialization beats with 1 Newton step corresponding to a semi-implicit (linearly-implicit) discretization method [109]; 10 beats with 2 Newton steps which is required to get a correct update due to the velocity dependence of the active-stress model; and two final beats with a fully converged Newton with a relative error reduction of the residual of ϵ = 10 −6 . First, using default passive parameters we again see a clear influence of the bulk modulus κ when using P1-P0 elements while for stabilized elements PV loops almost coincide for different values of κ, see Fig. 14 (a). Shown is element-wise passive stress which is not smoothed at the (a) loaded state using fitted parameters given in Table 3; (b) loaded state using default parameters given in Table 3 and a bulk modulus of κ = 650 kPa for all element types; (c) loaded state using default parameters given in Table 3 and a bulk modulus of κ = 5000 kPa for P1-P0 elements and 1/κ = 0 for projection-stabilized and MINI elements.
Looking at myocardial mass in Fig. 14 (b) it shows that for projection-stabilized and MINI elements the mass remains unchanged at the initial value of 133.56 cm 3 all the time when the tissue is modeled to be fully incompressible (1/κ = 0). Using a smaller value of κ = 650 kPa the tissue behaves nearly incompressible and especially during the ejection phase the myocardial mass decreases slightly by around 1.2 %. This is about the same volume change that we see for P1-P0-AS elements using a penalty formulation with the same κ = 650 kPa. Compared to this, P1-P0-WAS elements result in a minimally larger volume change of around 1.4 %.
Further, see Fig. 15 (a), for a comparison of the final three PV loops using fitted passive but the same active stress parameters which resulted in higher pressures for P1-P0 elements. In Fig. 15 (b)-(d) we show traces for the refitted active stress parameters as described above. Here, the final three PV loops in Fig. 15 (b) coincide. Hence, the solution is converged and there is also no difference between the simulation with two Newton steps and the fully converged Newton method. This also holds true for pressures in the ventricles and adjacent arteries, see Fig. 15 (c), in-and outflow traces Fig. 15 (d), and strains/stresses.
In contrast to pressure, volume, and flow traces, stresses show a very different pattern when comparing lockingfree to simple P1-P0 elements, see Fig. 17. In this plot we show element-wise, total first principal stress at three Fig. 17. 3D-0D model of the heart: Snapshots of first principal stress values at end-diastole (first row), peak systole (second row), and at mitral valve opening (third row), see Fig. 15 "B", "C", "D" for a visualization of the time-points. Shown is element-wise total stress which is not smoothed. Gray outlines indicate the end-diastolic configuration. Compared are P1-P0-AS elements (first column), P1-P0-WAS elements (second column) stabilized P1-P1 elements (third column), and MINI elements (fourth column). time points marked by B, C, and D in the PV loop in Fig. 15 (b) which represent B, the most expanded (enddiastole), C, the highest total stress (peak systole), and D, the most contracted (beginning of filling phase) states of the ventricles. Especially at end-diastole Fig. 17 (a) and the beginning of the filling phase Fig. 17 (c) where passive stress dominates over active stress we see a distinct checkerboard pattern for P1-P0 elements while solutions for projection-stabilized and MINI elements are smooth. Also in violin plots showing the stress distribution over the whole tissue domain we see a clear difference for these time points, see Fig. 19 (a) and (c). On the other hand, the stress distribution is very similar for all element types at peak-systole where active stress dominates, see Fig. 17 (b) and Fig. 19 (b).
Numerical performance. Computational times for the simulation using different element types are given in Table 4. Simulations were performed on 128 cores of Archer2 and we distinguish between solver-time, t s,• , the accumulated time of the linear solver (GMRES), and assembly-time, t a,• , the accumulated time of matrix and vector assembling of the linearized system (18)- (20).
In total, for a full simulation with loading, 18 initialization beats with 1 Newton step, 10 initialization beats with 2 Newton steps, and 2 final beats with a fully converging Newton method the computational costs were around 3.5 h for P1-P0 elements, 13 h for projection stabilized elements, and 17.5 h for MINI elements, see Table 4 for exact values. Here, in addition to GMRES solver and assembly times, also input-output times, the solution of the Fig. 18. 3D-0D model of the heart: Violin plots of the first principal stress distribution after loading (a) loaded state using fitted parameters given in Table 3; (b) loaded state using default parameters given in Table 3 and a bulk modulus of κ = 650 kPa for all element types; (c) loaded state using default parameters given in Table 3 and a bulk modulus of κ = 5000 kPa for P1-P0 elements and 1/κ = 0 for projection-stabilized and MINI elements. Table 3 Summary of electrical and mechanical material parameters.   Discussion. In this benchmark, we show one of the most complete model of cardiac EM that is currently available: (i) Cardiac electrophysiology was modeled by a reaction-Eikonal model which predicts potential fields with high fidelity even on coarser grids [85]. (ii) Cellular dynamics were modeled by the physiological Grandi-Pasqualini-Bers model [98] which is coupled to the Land-Niederer model [99], representing human ventricular electrophysiology and active stress generation. This allows for strong coupling, i.e., to account for length and velocity effects on the cytosolic calcium transient, using an approach as described in Augustin et al. [83]. (iii) Passive tissue mechanics was modeled by the recent Holzapfel-Ogden type model [97] and active stress according to [101] using a recent approach by Regazzoni and Quarteroni [100] to avoid oscillations. Note that both, passive and active stress computation, account for fiber dispersion, hence, allowing to model the active tension generated by dispersed fibers. (iv) Spatially varying Robin boundary conditions were included to model the effect of the pericardium [105]. (v) The 3D PDE Table 4 Summary of computational times on 128 cores of Archer2 for the different finite element types. Given are degrees of freedom (DOF) for the finest resolution bi-ventricular geometry as well as solver, assembly, and total computational times for one beat using one Newton iteration (t s,1 , t a,1 ) two Newton iterations (t s,2 , t a,2 ) and a fully converged Newton solution (t s,c , t a,c ). T b,1 , T b,2 , and T b,c correspond to the total simulation time per heart beat for one Newton iteration, two Newton iteration, and fully converged Newton scenarios. Timings refer to a single heart beat lasting 0.585 s at a time step size of 1 ms. In addition, the times required for the loading phase, T ld , and the total simulation times including 0D solution, IO, and postprocessing are presented. Note that for P1-P0 elements the hydrostatic pressure p is statically condensed on the element level and thus not considered for counting DOFs. model was coupled to the physiologically comprehensive 0D closed-loop model CircAdapt of the cardiovascular system. This allows to replicate physiological behaviors under experimental standard protocols altering loading conditions and contractility [71]. (vi) Simulations were performed using locking-free finite elements -as presented in this paper -to accurately compute stress distributions. All computational models of cardiac EM presented in the literature so far, e.g., [28,71,95,[110][111][112][113][114][115][116], are missing one or mostly more of the above points. While the importance of model components (i)-(v) was discussed extensively in references above we could show in this paper that also (vi) is necessary to compute accurate stress fields in the tissue. Depending on the application this could be a very critical modeling component, e.g., for the accurate prediction of rupture risks or for the estimation of growth and remodeling based on stress. Note that also for this benchmark we see no substantial difference between simulation outcomes using standard P1-P0-AS and the P1-P0-WAS formulation; both approaches fail to match stress results from gold-standard elements.
While stress fields differ vastly the PV loops computed with locking-free and simple P1-P0 elements are very similar; at least if the bulk modulus κ is not too large, see Fig. 14 (a), or if active and passive tissue parameters are fitted independently for each element type, see Fig. 15 (b). For the fitting of passive parameters to the Klotz curve, which is state-of-the-art in recently published EM models [14,28,71] it shows that using P1-P0 elements correspond to a softer material. Here, the fitting compensates volumetric locking effects to a certain degree. On the other hand, the reference peak tension parameter (T • ref ) had to be slightly reduced to reach the same target value as projection stabilized and MINI elements. As PV loops for default parameters coincide in Fig. 14 (a), we attribute the higher active stress generation in this case to the softer passive material and hence, a reduced resistance of the ground matrix material to compression in fiber direction. However, when using default parameters or fitted parameters for active and passive tissue properties, PV traces for the different element types are almost identical. This shows that simple P1-P0 elements can predict most simulation outputs as good as gold standard formulations and are thus adequate for cardiac EM simulations when stresses are not a primary quantity of interest. In this case, the computational efficiency of P1-P0 elements might trump the numerical accuracy of locking-free elements.
High resolution EM models require efficient numerical solvers to limit the computational cost that results from a high number of degrees of freedom to capture anatomical details as well as a high number of time steps. Strong scaling characteristics of our EM framework were reported in detail previously [83,117]. In the present work, we showed that strong scaling is preserved when using locking free elements and coupling to a 0D model of the circulatory system.
Using the advanced approach presented in this paper, the time needed for the passive filling of the bi-ventricular model of the heart is very low. For the projection stabilized element loading times are less than half a minute on 128 cores@2.25 GHz of Archer2 (Table 4) for the simulation with 444 936 degrees of freedom. Even on 8 cores the passive filling could be achieved within 5 minutes (Fig. 20). In comparison, a recent work [21] reports compute times for a similar passive inflation scenario using locking-free elements that were around 162 min (203 214 degrees of freedom, 8 cores@2.8 GHz). Fast loading times are crucial for the estimation of the stress-free reference configuration using fixed-point iterations [118,119].
Computational cost for one heart beat -using grids with a comparable number of elements and nodes but, in general, P1-P0 elements -range from 1.8 to 24 h in previous studies [95,111,112,114,115]. In this work, we could show that even with locking-free elements one heart beat can be simulated within 27 min on 128 cores of Archer2, see Table 4, and within 11 min using 1024 cores, see Fig. 20. Still, P1-P0 elements are computationally less expensive with one heart beat in around 7 min using 128 cores and 3.2 min using 1024 cores. Using a coarser mesh as in [71] fast computational times are also possible on desktop machines with 2.5 min for one heart beat on 32 cores using P1-P0 elements and 7.15 min using locking-free elements. This computational efficiency is of paramount importance for future parameterization studies where numerous forward simulations have to be carried out to personalize models to patient data.
Limitations. First, we set an arbitrary number of 30 heart beats which was more than enough to reach a limit cycle in all experiments. However, an automatic stopping criterion could be used that stops the simulation after reaching the limit cycle. Simulation times could be further reduced by accelerating the convergence to a limit cycle using data-driven 0D emulators [120] or by tuning the 0D CircAdapt model to predict PV traces from the 3D-0D model more effectively, hence, reducing the number of beats needed to a converged solution.
Second, the parameterization of the model has room for improvement. Passive parameters were fitted to the empiric Klotz curve using end-diastolic volume and pressure and active stress was fitted using a target peak pressure value. Other model components such as the reaction-Eikonal model and CircAdapt were not parameterized to fit experimental data, default values from the literature were used. The personalization of the complete model to patientspecific data is not within the scope of this contribution, however, the computational efficiency of the model is of crucial importance for parameter identification studies that often require a large number of forward simulations.
Third, no independent validation of the model is performed. This could be done by comparing displacements or strains predicted by the model to observations from cine MRI or 3D tagged MRI data, see, e.g., [14,114,121]. However, in this work, we focused on showing advantages of locking-free elements for applied simulations using an advanced setup which is necessary to replicate physiological behavior. A rigorous, independent validation as described above will be the focus of future studies.
Finally, locking-free formulations as presented in this paper require the solution of a block system, which in turn necessitates suitable preconditioning for computational efficiency. This is not a trivial task, however, preconditioners used for simulations in this paper are publicly available through the open-source software framework PETSc.

Conclusion
In this study, we introduced stabilization techniques that accelerate simulations of elastic materials, in particular, we apply the methods to model nearly and fully incompressible fiber-reinforced solids such as arterial wall or myocardial tissue. A MINI element formulation and a simple and computationally efficient technique based on a local pressure projection were presented. Both methods were applied for the first time for simulations of anisotropic materials and showed to be an excellent choice when the use of higher order or Taylor-Hood elements is not desired. This is the case, e.g., for detailed, high-resolution problem domains that results in a high number of degrees of freedom. We showed that both approaches are very versatile and can be applied to stationary and transient problems as well as hexahedral and tetrahedral grids without modifications. It is worth noting that all required implementations are purely on the element level, thus, facilitating an inclusion in existing FE codes. Furthermore, solvers and preconditioners used to solve the linearized block system of equations are available through the open-source software package PETSc [87].
We showed the robustness and accuracy of the chosen approaches in two benchmark problems from the literature: first, a thick-walled cylindrical tube representing arterial tissue and second, an ellipsoid representing LV myocardial tissue. Additionally, in a third application of the stabilization approaches, we presented a complex 3D-0D model of the ventricles. This constitutes a highly advanced computational EM model of the heart where all components are captured by physiological, state-of-the-art models. We could show that for the first time accurate and physiologically detailed cardiovascular simulations are feasible within a clinically tractable time frame.
Computational efficiency of the methods is unprecedented in the literature and the framework shows excellent strong scaling on desktop and HPC architectures. The high versatility of the one-fits-all approach allows the simulation of nearly and fully incompressible fiber-reinforced materials in many different scenarios. Overall, this offers the possibility to perform accurate simulations of biological tissues in clinically tractable time frames, also enabling parameterization studies where numerous forward simulations have to be carried out to personalize models to patient data.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Linearization
We will give a short summary of the linearization of a cavity volume V CAV defined by x · n ds x .
Using Nanson's formula and x = X + u we can rewrite this as

Appendix B. Static condensation for inhomogeneous Neumann boundary condition
While homogeneous Neumann boundary conditions do not alter the process of static condensation, the procedure needs to be adapted to for inhomogeneous ones. First, looking at the definition of the nonlinear residual R vol in (6) we see that this can be split as where R vol,Ω 0 holds all the terms coming from integration over the domain Ω 0 and R vol,Γ N ,0 holds all the terms coming from integration over the Neumann surfaces. Next, note that the bubble functionsψ B for tetrahedral elements as well as their hexahedral counterpartsψ B,1 ,ψ B,2 have compact support in the FE interior. However, the respective gradients do not vanish on the FE boundary. Consider an arbitrary FE K ∈ T h with K ∩ Γ N ,0 ̸ = ∅. The gradient of a bubble function occurs in bilinear-form a k N in (17), for the argument ∆u ∈ V h . This yields a non-zero contribution to the element-stiffness-matrix, whereas there is no contribution from R vol,Γ N ,0 to the total element residual vector. Using the decomposition of local degrees of freedom into exterior, E and interior, I it follows that the local block system can be written in the following form ⎛ The interior degrees of freedom can be statically condensed. On element level this leads to the static condensed system II R vol,I ,R Γ ,vol := R Γ ,E − K Γ ,EI K −1 II R vol,I , R inc := R inc,E − C I K −1 II R inc,I . The individual element matrices/vectorsÃ,Ã Γ ,R, andR Γ can be assembled into a global stiffness matrix through loops over volume elements and surface elements respectively. In the case of an attached circulatory system a static condensation can be performed in an analogous way ⎛ ⎜ ⎜ ⎝ where we assumed n CAV = 1 for brevity of presentation. The generalization to n CAV > 1 is straightforward. Static condensation of all interior degrees of freedom leads to

Appendix C. Tensor calculus
We use the following results from tensor calculus, for more details we refer to, e.g., [68,122].