Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr

Computers and Mathematics with Applications xxx (xxxx) xxx

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr David Kamensky Department of Mechanical and Aerospace Engineering, University of California, San Diego, 9500 Gilman Drive, Mail Code 0411, La Jolla, CA 92093, USA

article

info

Article history: Available online xxxx Keywords: Isogeometric analysis Fluid–structure interaction Immersed boundary Open-source software FEniCS Heart valve

a b s t r a c t We recently developed the open-source library tIGAr, which extends the FEniCS finite element automation framework to isogeometric analysis. The present contribution demonstrates the utility of tIGAr in complex problems by applying it to immersogeometric fluid–structure interaction (FSI) analysis. This application is implemented as the new open-source library CouDALFISh (Coupling, via Dynamic Augmented Lagrangian, of Fluids with Immersed Shells, pronounced ‘‘cuttlefish’’), which uses the dynamic augmented Lagrangian (DAL) method to couple fluid and shell structure subproblems. The DAL method was introduced previously, over a series of papers largely focused on heart valve FSI, but an open-source implementation making extensive use of automation to compile numerical routines from high-level mathematical descriptions brings newfound transparency and reproducibility to these earlier developments on immersogeometric FSI analysis. The portions of CouDALFISh that do not use code generation also illustrate how a framework like FEniCS remains useful even when some functionality is outside the scope of its standard workflow. This paper summarizes the workings of CouDALFISh and documents a variety of benchmarks demonstrating its accuracy. Although the implementation emphasizes transparency and extensibility over performance, it is nonetheless demonstrated to be sufficient to simulate 3D FSI of an idealized aortic heart valve. Source code will be maintained at https://github.com/david-kamensky/CouDALFISh. © 2020 Elsevier Ltd. All rights reserved.

1. Introduction Beginning in 2014, a series of papers by Kamensky, Hsu, and other collaborators [1–3] introduced two concepts during an effort to develop a simulation framework for heart valve fluid–structure interaction (FSI): 1. Application of the term immersogeometric analysis (IMGA) to describe the use of immersed-boundary numerical methods (interpreted broadly) with exact design geometries. This extends isogeometric analysis (IGA) [4,5] to situations in which the partial differential equation (PDE) domain is difficult or impossible to parameterize in a conforming way.1

E-mail address: [email protected]. 1 Earlier work exists on combining immersed-boundary and isogeometric methods, e.g., [6], where the finite cell method [7] is used as the immersed-boundary approach. However, the term IMGA is first used to describe this trend in [1]. https://doi.org/10.1016/j.camwa.2020.01.023 0898-1221/© 2020 Elsevier Ltd. All rights reserved.

Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.

2

D. Kamensky / Computers and Mathematics with Applications xxx (xxxx) xxx

2. The dynamic augmented Lagrangian (DAL)2 method, which is an immersed-boundary method based on truncating the classical augmented Lagrangian iteration of Hestenes [10] and Powell [11] in time-dependent problems, to relax stability conditions on the spatial discretization of the Lagrange multiplier field. DAL-based IMGA of FSI has proven to be successful, and has been applied not only to heart valve FSI analysis [12–14], but also to optimization of other devices [15]. Extensions of the basic method are proposed in [9,16], and the convergence of one of the resulting DAL variants is analyzed a priori for an elliptic model problem in [8]. A great many other works have been published on immersed-boundary FSI analysis, but these are comprehensively reviewed elsewhere, e.g., [17–20], and we therefore focus on introducing the specific methods which the software described by this paper implements. Published work to-date on DAL-based IMGA has relied on closed-source prototype research codes to produce numerical results. If the source code were made openly-available, it would render results reproducible given the same input data, but this would provide little insight, as minimal consideration was given to software design, and the implementation details would be difficult to understand and adapt to new problems. (This becomes immediately clear during the training of new graduate students!) The idea that implementations can serve as a useful form of scholarly communication depends not only on reproducibility, but also transparency, which we define here as the ease with which other researchers – many of whom may be domain scientists rather than software specialists – can ascertain precisely what a given code does. Widespread transparency has long evaded the field of computational PDEs, since optimizations needed for sufficientlyfast solvers tend to obfuscate the function of software. However, a major breakthrough in this area has been the development of automatic code generation, where efficient numerical routines are compiled from high-level mathematical descriptions, rather than implemented and optimized by hand. This trend is exemplified by the FEniCS Project [21], which encompasses software to specify variational forms of PDEs in the high-level Unified Form Language (UFL) [22], compile these descriptions into numerical routines [23], and use these routines within the parallel finite element solver DOLFIN [24]. This software has been demonstrated to attain high levels of computational performance, even in problems with billions of degrees of freedom [25]. Another recent software project based on generating high-performance code from UFL is Firedrake [26–29]. Kamensky and Bazilevs [30] recently introduced the open-source library tIGAr [31], which allows users to perform IGA using FEniCS. Despite its recent release, tIGAr has already been successfully applied by a number of groups for research on diverse topics, including homogenized strain-gradient modeling of metamaterials [32], shock-capturing formulations based on novel entropy-like principles [33], residual-based shock-capturing in solids [34], divergence-conforming stabilized methods for incompressible flow [35], and machine-learning-based constitutive models for biological soft tissues [36]. In the present work, we couple a tIGAr-based thin shell discretization to a FEniCS-based fluid discretization, using the DAL method, to perform immersogeometric FSI analysis. This functionality is implemented within a library called CouDALFISh (pronounced ‘‘cuttlefish’’), which is an acronym for Coupling, via DAL, of Fluids and Immersed Shells. Special attention is paid to the modular design of the software solution, to permit straightforward interpretation, extension, and application to diverse fluid–thin structure interaction problems. The remainder of this paper is structured as follows. Section 2 briefly recapitulates the DAL method for FSI analysis and Section 3 reviews the main technical ideas behind the implementation of tIGAr. We then describe the design of CouDALFISh in Section 4. The accuracies of CouDALFISh and the modules used by it for the fluid and structure subproblems are verified by numerical results presented in Section 5. Lastly, the robustness of CouDALFISh is supported by its application to heart valve FSI analysis in Section 6. Conclusions and directions for future work are discussed in Section 7. The source code of CouDALFISh and the code and data for the applications demonstrated in this paper are openly available online [37]. 2. The DAL method for fluid–thin structure interaction This section describes the FSI temporal discretization used by CouDALFISh, which is based on the formulation of [9], with some minor modifications. We abuse notation somewhat in the interest of brevity; readers dissatisfied with any of the resulting ambiguities should consult the thoroughly documented source code and/or cited references. This discussion is primarily focused on the coupling between the fluid and shell structure, and we therefore introduce the following abbreviated notations for the subproblems. We write the time-discrete incompressible Navier–Stokes subproblem as: Find velocity un+1 ∈ V and pressure pn+1 ∈ Q such that for all v ∈ V and q ∈ Q, Fuf n un+1 , pn+1 , v, q = 0,

(

)

(1)

n

where the subscript u denotes the variational form’s dependence on the (known) previous time step’s velocity and V and Q are function spaces for the velocity and pressure. Likewise, for the shell structure, we denote its time-discrete subproblem by: Find displacement yn+1 ∈ Y such that for all z ∈ Y , Fyshn ,˙yn yn+1 , z = 0,

(

)

(2)

2 The method is introduced in [1] without a name; the first published use of ‘‘DAL’’ is [8], but the term was initially coined by J. A. Evans in discussions behind [9]. Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.

D. Kamensky / Computers and Mathematics with Applications xxx (xxxx) xxx

3

where Y is the displacement function space and yn and y˙ n are the known displacement and velocity from the previous time step. The formulas for F(f·) and F(sh ·,·) depend on the choice of time integrator and fluid and structure formulations; see [1] for examples of complete formulas. The DAL method approximately enforces the constraint3

∫

un+α ◦ φt n+α − y˙ n+α · δλ = 0 ∀δλ ∈ L,

)

( Γ0

(3)

where Γt is the midsurface of the shell structure at time t, φt maps Γ0 to Γt , n +α is some intermediate time level between n and n + 1 at which the constraint is collocated, y˙ n+α is the intermediate-level structure velocity, whose dependence on yn+1 , yn , and y˙ n varies with the details of the time integrator, and L is a function space on Γ0 against which the constraint is tested. As a starting point, consider the following augmented Lagrangian formulation: Find un+1 , pn+1 , yn+1 , and λn+1 ∈ L such that for all v, q, z, and δλ, F f + F sh +

∫ Γ0

∫

(v − z) · λn+1 un+α − y˙ n+α · δλ

(

−

∫Γ0 + Γ0

)

( ) β un+α − y˙ n+α · (v − z) = 0,

(4)

where β ≥ 0 is a penalty parameter, to augment constraint enforcement by the Lagrange multiplier λn+1 , and obvious function(al) arguments and compositions with φt n+α have been omitted for visual clarity. Remark 1. It is clear that β is a dimensional quantity. Guidelines for selecting it based on physical and discretization parameters are given in the cited references. In the limit of the spatial discretization length scale h → 0, analysis of parabolic model problems in [3,8] suggests that β should scale like C µ/h, where C is an O(1) dimensionless constant and µ is the dynamic viscosity of the fluid. However, the optimal choice, especially in the pre-asymptotic regime, remains an open question. To derive a DAL method similar to that of [9], we first regularize (4) by using a scalar Lagrange multiplier to enforce only the normal part of the constraint, and perturbing the constraint Lagrangian with a scalar parameter r ≥ 0: f

F +F

sh

∫ + Γ0

λn+1 (v − z) · nn+α

∫ (

(

n+α

u

−

∫Γ0 + Γ0

n+α

− y˙

)

·n

n+α

−

r λn+1

β

)

· δλ

( ) β un+α − y˙ n+α · (v − z) = 0,

(5)

where nn+α is some approximation of the normal vector field to the deformed shell structure midsurface Γt n+α . In this regularized problem, the tangential portion of the constraint is enforced entirely by the penalty term. If the space L = L2 (Γ0 ), one can also view the perturbed Lagrangian of the normal-direction constraint to be a penalty enforcement that gets stronger as r → 0 (cf. [38]). Note that, if the penalty scales inversely with the mesh refinement parameter h → 0, as suggested in Remark 1, this still enforces the original constraints in the converged limit, and, in that sense, does not represent a modification of the original problem (assuming that this scaling of the penalty is stable; see [3,8] for analysis using model problems). The defining feature of DAL (highlighted below in red) is to lag the Lagrange multiplier and update it based on the penalty term, removing it as an unknown from the implicitly-solved problem, thereby greatly simplifying the solution process: Find un+1 , pn+1 , and yn+1 such that for all v, q, and z, F f + F sh +

∫ Γ0

λn (v − z) · nn+α

∫ + Γ0

( ) β un+α − y˙ n+α · (v − z) = 0,

(6)

then find λn+1 such that for all δλ

( ( ))) ∫ ( ( ) r λn+1 λn+1 − λn + β un+α − y˙ n+α · nn+α − δλ = 0 . β Γ0

(7)

3 CouDALFISh imposes the constraint with a pulled-back fluid velocity on Γ , as written here, which differs from earlier work in [1] and other 0 previous references. We find the difference to be inconsequential in both benchmark problems and applications. Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.

4

D. Kamensky / Computers and Mathematics with Applications xxx (xxxx) xxx

For r > 0, the update (7) is implicit, but if L is piecewise-constant and one-point quadrature is used on its elements, a static condensation procedure provides the explicit update formula4

λn+1 =

1 1+r

) ) ( n ( λ + β un+α − y˙ n+α · nn+α .

(8)

In the particular implementation used by CouDALFISh, the implicit problem (6) is rendered ‘‘explicit in geometry’’, by using explicit predictors for the mapping φt n+α and the normal vector field nn+α . We find that this greatly increases the robustness of the nonlinear solution procedure. For time integration of the subproblems, CouDALFISh uses the second-order generalized-α method [39,40], to permit control over numerical dissipation of high-frequency solution components. The DAL algorithm is itself limited to first-order accuracy and previous error analysis has focused on the case of α = 1, i.e., backward Euler integration, for simplicity. However, the additional accuracy provided by the generalized-α method in subproblems improves the quality of flow features away from the immersed structure. We do not assume that the constructions of V or Q conform in any way to the immersed interface Γt , so the level of accuracy near the immersed surface will be dominated by spatial truncation error and there is no significant benefit to improving the temporal accuracy of DAL. Remark 2. A projection-based stabilization method, in which the stabilizing effect of the parameter r is applied only to the orthogonal complement to some coarse space for λn+1 is introduced in [16] and analyzed in the limit of r → ∞ for a scalar model problem in [8]. This projection-based method has various merits, as discussed at length in the cited works, but is not currently implemented in CouDALFISh. A weaker convergence analysis covering the un-stabilized case of r = 0 is given for a scalar model problem in [3, Section 3]. Remark 3. Although the methods are derived from unrelated lines of reasoning, there is a close relation between DAL and the staggered fluid–thin structure interaction scheme recently devised by Boilevin-Kayl et al. [41]. Boilevin-Kayl et al. also update the Lagrange multiplier explicitly [41, Algorithm 4.3, Step 3], and, up to minor differences, the scheme from the cited work can be interpreted as a clever choice of penalty parameter in DAL, to maintain stability while decoupling the implicit problem (6). If the structure is completely constrained by a Dirichlet boundary condition, the method of [41] is essentially equivalent to DAL (or an implicit variant of Goldstein et al.’s feedback method [42], as discussed at length in [3, Section 2.5.4]). However, in the context of simply applying Dirichlet boundary conditions to the fluid subproblem on arbitrary immersed boundaries, the intrinsic physical properties of these boundaries become irrelevant (if they are defined at all), and the argument of [41], based on structural inertia, no longer makes sense as a criterion for selecting the penalty. However, the analysis of [41] may provide guidance for choosing an optimal DAL penalty for soft, mobile immersed structures (cf. Remark 1). 3. Technical summary of tIGAr This section provides a brief overview of the tIGAr library, which allows for IGA using FEniCS. The library tIGAr is written in Python, and interfaces with the FEniCS project software through the latter’s Python bindings. The technical basis for tIGAr is the concept of extraction [43–45], namely, that smooth polynomial spline functions can be represented as linear combinations of less-smooth spline functions, including classical high-order finite element (FE) basis functions. Consider a matrix M with a column for each spline basis function, a row for each FE basis function, and entries selected such that each column provides the coefficients to express the corresponding spline basis function as a linear combination of FE basis functions. If the matrix AFE is the coordinate representation of a bilinear form in the finite element basis, then AIGA := MT AFE M is its representation in the spline basis. Similarly, for a linear form represented by the vector FFE , its spline basis coordinate representation is given by FIGA = MT FFE . FEniCS provides functionality to automate the assembly of AFE and FFE from mathematical descriptions of the corresponding forms. The additional library tIGAr provides functionality to apply (an arbitrary) M during solution procedures, alongside an abstract interface for generating M, subclasses of which may be implemented for a wide variety of spline types. Readers concerned with the additional technicalities of treating rational splines and spline-based geometrical mappings should refer to the original reference [30] and/or detailed documentation distributed with tIGAr. In summary, geometrical mappings and rational weights are incorporated into the definitions of variational forms, using user-defined extensions to UFL, in contrast to the more canonical approach of including them within shape function routines.5 4 Alternatively, one can consider the regularized constraint enforced ‘‘pointwise’’ on a quadrature rule, or leave L = L2 (Γ ) in the discrete problem 0 and confess to a quadrature crime, but this is somewhat at odds with the typical functional setting for the continuous problem (cf. [8]) in which the integrand of (7) is not well-defined at points. 5 This approach is made feasible by recent developments in UFL compilation technology, including the UFL Analyser and Compiler System (UFLACS), which became the default within FEniCS several years ago, and The Smart Form Compiler (TSFC) [27].

Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.

D. Kamensky / Computers and Mathematics with Applications xxx (xxxx) xxx

5

3.1. ShNAPr: Application to Kirchhoff–Love thin shells A nontrivial application of tIGAr documented in [30, Section 5.1] is to Kirchhoff–Love thin shell analysis, using the isogeometric formulations developed by Kiendl and co-workers [46–48], alongside independent parallel work by BuganzaTepole et al. [49]. The notorious complexity of shell formulations makes them an ideal showcase of the automatic code generation paradigm of the FEniCS Project; this is demonstrated in the classical finite element setting by the FEniCSShells library [50]. At the time of publication of [30], the Kirchhoff–Love shell demonstrations included with tIGAr (one of which is discussed at length in [30, Section 5.1]) already contained significant overlapping code with each other. As the present contribution involves yet more applications of this theory, we have consolidated these demonstrations into a separate library, called ShNAPr [51], which is a toolkit for developing Shell Nonlinear Analysis Programs using tIGAr. In particular, the redundant kinematic definitions between the St. Venant–Kirchhoff formulation of [46,47] and the hyperelastic formulations of [48,49] have been unified under the submodule ShNAPr.kinematics, while additional definitions for St. Venant–Kirchhoff and incompressible hyperelastic formulations have been placed in submodules ShNAPr.SVK and ShNAPr.hyperelastic. As many interesting applications of the Kirchhoff–Love shell theory (including the heart valve example shown in Section 6 of this paper) require robust and flexible (self-)contact of shell structures, ShNAPr also includes the submodule ShNAPr.contact, to support this need. The ShNAPr.contact submodule implements a meshfree discretization of the nonlocal (and arguably peridynamic [52]) contact formulation introduced in [53]. Benchmarking of ShNAPr is carried out in Section 5.2. 4. Design of CouDALFISh Like tIGAr, CouDALFISh is a Python library that interacts with FEniCS through the Python interface. The design of CouDALFISh assumes that the shell structure is discretized using ShNAPr (as summarized in Section 3.1), but leaves the discretization of the fluid subproblem flexible, assuming only that it is implemented using finite elements in FEniCS, accepting an arbitrary variational form specified in UFL at construction. The information regarding the shell formulation, fluid formulation, coupling, and current state of an FSI problem is encapsulated by the CouDALFISh library in the eponymous class, CouDALFISh.6 The main data members of a CouDALFISh instance are:

• A finite element mesh mesh_f for the fluid subproblem. • A tIGAr ExtractedSpline object named spline_sh encapsulating the geometry and function spaces for the shell structure.

• Generalized-α time integrator objects (as defined by the tIGAr submodule tIGAr.timeIntegration) for the fluid and structure, which permit access to the fluid and structure solution states at various time levels.

• An instance of ShellContactContext, from the ShNAPr submodule ShNAPr.contact, named contactContext_sh. While the purposes of including most of these data members are obvious, it is perhaps not intuitive that the data member contactContext_sh plays a critical role, regardless of whether there is contact in the shell structure subproblem. This is because the numerical integration rule used for the nonlocal contact method implemented by ShNAPr.contact is assumed by CouDALFISh to be the same rule used for the integrals over Γ0 in the DAL formulation. In particular, the numerical integration relies on the quadrature-based mass lumping trick described in [30, Section 5.3] to streamline the assembly procedure of the Lagrange multiplier and penalty forces acting on the shell structure. While quadrature rules based on a parameterization of an immersed surface are not optimal for integrating traces of functions on meshes into which the surfaces are immersed, the error resulting from such immersed surface quadrature tends to be relatively small in comparison to the effects of errors in volume quadrature [54]. Aside from its constructor, the main method of CouDALFISh that is intended for ‘‘public’’ use (although Python does not have formal public/private distinctions for class members) is the takeStep() method, which advances the state of the CouDALFISh by one time step, according to the DAL method outlined in Section 2. The loop over time steps is implemented in applications using CouDALFISh rather than within the library. This is to permit application-specific tasks within the time stepping loop, such as custom file output and/or extraction of quantities of interest. Most other methods of the class are mainly intended as ‘‘private’’ methods, to be called by takeStep() (although they could conceivably be used to implement alternative solution strategies). The most challenging part of the takeStep() method is of course the solution of the nonlinear algebraic problem resulting from spatial discretization of (6). For this, CouDALFISh applies the block iterative algorithm described in [3, Section 2.5.3]. The convergence behavior of this algorithm is analyzed for a model problem in [3, Section 4]. The fluid and shell subproblem tangent matrices needed by block iteration are obtained automatically using symbolic Gateaux differentiation functions in UFL, with the exception of the linearizations of the penalty terms, for which lumped-mass approximations are added directly to assembled matrices. Remark 4. CouDALFISh assumes that problems are posed in three spatial dimensions. This is an assumption inherited from its dependency ShNAPr, as shell theory is predicated on embedding 2D manifolds in 3D space. 2D problems may still be approximated, though, using a single layer of elements, as illustrated by the benchmark of Section 5.3.1. 6 To distinguish between the class and library, we typeset references to the class as inline code listings. Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.

6

D. Kamensky / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 1. Contours of velocity magnitude for the regularized lid-driven cavity problem, with magnitude-scaled arrow glyphs indicating direction.

5. Numerical testing The purpose of this section is to verify the correctness of CouDALFISh and the distinct fluid and shell structure libraries used to discretize its subproblems. As the benchmark problem definitions mostly come from published literature and the code and data for our computations are openly available, we focus on summarizing the main results of our testing and refer readers to these other resources for details. 5.1. Fluid benchmarks Although, as mentioned in Section 4, CouDALFISh is agnostic regarding the precise fluid formulation, computations presented in this paper use the small Python module VarMINT (Variational Multiscale Incompressible Navier–Stokes Toolkit) [55] to define the fluid formulation. This is a straightforward FEniCS implementation of the variational multiscale (VMS) formulation developed and widely applied by Bazilevs and collaborators [56–65].7 This section documents the results from several tests of the VarMINT module. Although the structure of VarMINT does not strictly require it, all results in this section use the same scalar finite element space for the pressure and each Cartesian component of the velocity. 5.1.1. Regularized lid-driven cavity We begin our numerical testing of VarMINT with steady Navier–Stokes flow. In particular, we use a variant of the regularized lid-driven cavity problem proposed in [68] and studied as well by [69]. We choose, a priori, the exact velocity solution given by [68], u(x) = 8x41 − 2x31 + x21

(

)(

4x32 − 2x2 e1 − 8 4x31 − 6x21 + 2x1

)

(

)(

x42 − x22 e2 ,

)

(9)

and the exact pressure solution p(x) = sin(π x1 ) sin(π x2 ),

(10)

which differs from that emerging from [68]’s complicated manual manipulations. The velocity field looks qualitatively like that of the classic lid-driven cavity problem, but without corner singularities. We automatically generate the corresponding source term from the strong form of the problem, using UFL, and apply the exact velocity solution as a weakly-enforced Dirichlet boundary condition [70] on a unit square domain. We consider unit length and velocity scales, and show results for Re = 100. A representative numerical solution is plotted in Fig. 1. Fig. 2 shows that, for Lagrange finite elements of varying polynomial degrees, the H 1 seminorm of the velocity error converges to zero at optimal rates with respect to h-refinement. 7 A more comprehensive open-source FEniCS-based solver using this formulation is IlliniFlow [66], which is described by Zhu and Yan in [67]. Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.

D. Kamensky / Computers and Mathematics with Applications xxx (xxxx) xxx

7

Fig. 2. Convergence of H 1 velocity error for the regularized lid-driven cavity problem, using elements of polynomial degree k.

Fig. 3. Convergence of H 1 velocity error at fixed time for the 2D Taylor–Green vortex problem, using elements of polynomial degree k.

5.1.2. 2D Taylor–Green vortex Moving now to a dynamic problem, we test VarMINT on the well-known 2D Taylor–Green vortex (as specified precisely in, e.g., [71, Section 9.10.1]), using the implicit midpoint rule for time discretization. Looking at the H 1 seminorm of velocity error at a fixed time, at Reynolds number 100 (based on a unit length scale, with the standard 2π × 2π domain of the cited reference), we see, in Fig. 3, that linear elements provide first-order convergence with respect to h-refinement and quadratic elements provide second-order convergence. In both cases we maintain a time step proportional to the spatial element size during refinement; this is an assumption behind the design of the stabilization parameters used in the formulation.8 See [56,74–78] for more discussion of stabilization parameter selection. Thus, having reached the limiting convergence rate of the time integrator and noting that refining faster in time than space is not appropriate for the stabilized formulation of VarMINT, we do not consider elements of cubic or higher degree in this problem. 5.2. Shell structure benchmarking We now perform basic verification of the shell structure formulation implemented by ShNAPr. The cursory testing considered here is augmented by further FSI tests in Section 5.3, which would also be unlikely to pass with a flawed implementation of the structure subproblem. 5.2.1. Scordelis–Lo roof We consider here the Scordelis–Lo roof example, as proposed originally in [79], and solved using isogeometric Kirchhoff–Love shells in [47, Section 6.2.1]. The problem consists of a cylindrical shell subject to gravity loading, and the quantity of interest is taken to be the vertical displacement at the midpoint of one of the straight boundaries. The full problem definition is readily available in the cited references and/or code example provided with ShNAPr. We compare solutions computed with ShNAPr using various numbers of elements and polynomial degrees to the reference value for the quantity of interest provided in [47, Section 6.2.1]. Fig. 4 demonstrates that discretizations of all polynomial degrees converge to this reference value with relatively-few elements. A representative displacement solution is plotted in Fig. 5. 8 This restriction can be circumvented by using modified stabilization [72] or dynamic subgrid models [73], although the standard definition is, in practice, robust over a large range of time step sizes. Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.

8

D. Kamensky / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 4. Convergence of vertical displacement at the midpoint of a straight edge in the Scordelis–Lo roof problem, for discretizations of varying polynomial degree p.

Fig. 5. Representative numerical solution of the Scordelis–Lo roof problem computed using ShNAPr, colored by the component of displacement parallel to gravity. Displacement is scaled by a factor of 20 for warping the initial geometry.

Remark 5. ShNAPr’s kinematic definitions assume geometric nonlinearity, but the Scordelis–Lo problem is defined as geometrically linear. To solve the geometrically linear problem with ShNAPr, we simply define the full nonlinear problem, but take only a single Newton iteration, starting from a zero initial guess. The effect of any geometric nonlinearity would be small, though, given the magnitude of the displacements relative to the cylinder’s radius of 25 and length of 50. 5.3. FSI benchmarks In this section, we test the ability of CouDALFISh to converge to the solutions of easy FSI benchmark problems. Section 5.3.1 looks at convergence more rigorously, measuring Sobolev norms of errors in approximations of a manufactured exact solution, while Section 5.3.2 looks at time histories of quantities of interest in a benchmark problem from the literature and compares results with a reference computation using boundary-fitted FSI techniques. 5.3.1. Manufactured solution This section demonstrates the convergence of CouDALFISh to the exact solution of a large-displacement fluid– membrane interaction problem derived in [8, Section 7]. Complete expressions for the fluid velocity, membrane displacement, and necessary source terms are evident from the source code of this demonstration, which is distributed with CouDALFISh, and a detailed procedure for deriving this solution is given in the cited reference. This manufactured solution includes both a jump in pressure across the immersed membrane and a jump in velocity gradient. These are the main difficult solution features that an effective immersed-boundary fluid–thin structure interaction method must be able to contend with in real problems. As noted in Remark 4, CouDALFISh assumes that problems are 3D, while this benchmark is 2D. We consider a 3D extrusion of the solution in the x3 direction and, as there is no variation of the exact solution in this direction, we only refine in the x1 and x2 directions to study convergence. A representative snapshot of a numerical solution is shown in Fig. 6. Note that the initial configuration of the membrane is flat; this example involves large deformations, in which the structure and its attendant low-regularity solution features pass through many elements of the background fluid mesh. The convergence of the L2 and H 1 norms of fluid velocity is shown in Fig. 7. The H 1 convergence rate is the best possible given that spatial regularity of the exact velocity solution (when considered as a function on the entire domain) is only H 3/2−ϵ for arbitrary ϵ > 0, due to the jump discontinuity in the velocity gradient. This is in agreement with the error analysis and previous numerical results of [8] (although some technical details of the numerical methods differ). Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.

D. Kamensky / Computers and Mathematics with Applications xxx (xxxx) xxx

9

Fig. 6. Numerical approximation of the FSI manufactured solution derived in [8, Section 7]. Pressure contours are shown on a slice and velocity is illustrated with magnitude-scaled arrow glyphs.

Fig. 7. Convergence of velocity error norms at fixed time for FSI manufactured solution derived in [8, Section 7], where the vertical offset C is chosen differently for each norm, to make both rates clearly visible in the same plot.

5.3.2. 2D valve We now consider a more qualitative benchmark, namely the 2D valve-like problem from [1, Section 4.7], variants of which have been used in earlier works, e.g., [80–82]. The problem consists of two cantilever beams extending from the walls of a channel and being driven by a periodic inflow boundary condition, at an approximate Reynolds number of 110. Complete problem data can be found in [1, Section 4.7.1]. As in [1], we extract histories of the displacement of a Lagrangian point at the tip of one of the beams, and compare with the same quantity of interest extracted from a highly-refined interface-tracking solution of the same problem, using fluid mesh deformation techniques described in [57,78,83–86],9 as discussed further in [1, Section 4.7.2]. A representative snapshot of a numerical solution computed using CouDALFISh is shown in Fig. 8. The time history of the x1 component of displacement for the tip of one of the beams is compared with the body-fitted reference computation in Fig. 9, for discretizations with 16, 24, and 32 fluid elements across the channel, showing clear convergence toward the reference data. The source of the reference results [1, Section 4.7.2] estimates that the beam tip displacements in the body-fitted computation are accurate to about 0.001 length units, which would be imperceptibly-small in Fig. 9. 6. Application to heart valve FSI Moving beyond mild benchmarks, we now apply CouDALFISh to simulate FSI in an aortic heart valve. The purpose of this is not to draw any scientific conclusions about valve function, but rather to demonstrate that CouDALFISh is capable of:

• Solving 3D problems of meaningful scale. • Remaining robust when there are large pressure jumps across the immersed structure. • Working in conjunction with structural self-contact.10 9 Code and data for the immersed computations are available at [37], but we are not authorized to distribute the source code for the body-fitted reference computation. Data from the reference computation is reused directly from [1]. 10 In this problem, contact happens to occur only between topologically-disconnected regions of the structure. However, that fact is not exploited by the ShNAPr.contact module, which treats all contact as self-contact and applies, without modification, to cases like [30, Figure 4]. Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.

10

D. Kamensky / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 8. A snapshot of an immersogeometric solution to the 2D valve benchmark, with a slice showing velocity magnitude.

Fig. 9. Time history of the x1 component of displacement for one of the leaflet tips. The body-fitted setup for the reference results is detailed in [1, Section 4.7.2].

Based on our experience in this application domain (reviewed in Section 1), these are the primary challenges from the perspective of numerical method development and implementation. While important for pursuing some scientific objectives, implementing new hyperelastic constitutive models is essentially trivialized by FEniCS’s UFL11 and designing and meshing a high-fidelity geometric model of the valve and aorta is time-consuming, but orthogonal to the development of CouDALFISh. Therefore, to avoid distracting from the objectives itemized above, and to produce a clear, compact, illustrative example that is easily distributed and understood, constitutive and geometrical modeling have been kept to a minimum. We summarize here the main modeling choices and numerical results, but refer readers to the heavily-annotated source code and input data for details. For the valve itself, we consider a B-spline model of a generic trileaflet valve, based on the geometry constructed from an industrial bioprosthetic heart valve in [1, Section 5.1], with some manual modifications to facilitate use with nonlocal contact. The simplest choice of a St. Venant–Kirchhoff material model is assumed, with parameters based on [1]. The section of the aorta surrounding the valve is defined geometrically as the union of several constructive solid geometry primitives, then meshed in a quasi-uniform fashion using CGAL [87], through the mshr module of FEniCS. While this entails some simplifications, it qualitatively captures the three-lobed sinus structure and allows the mesh to be defined directly in the code (thereby simplifying the distribution of this example), and is sufficient to accomplish the goals stated at the beginning of this section.12 The flow through the valve is driven by a pressure difference between the inflow and outflow (as designated based on the prevailing flow direction through the open valve, during systole). Because the domain is rigid, the absolute value of pressure is irrelevant, and the desired difference is applied as a traction at the inflow, in conjunction with a traction-free boundary condition at the outflow.13 The pressure applied at the inflow begins at 2 × 104 Ba, to push flow forward through 11 To be clear, we refer to the case where the appropriate model and parameter values are known. We do not claim that determining or calibrating an appropriate model is simplified by UFL. 12 CouDALFISh is not fundamentally limited to such simplified meshing. A FEniCS workflow that could be used with CouDALFISh is described in [88, Appendix A.3]. 13 We technically abuse terminology here, because the true traction differs from the applied boundary data where flow is entering the domain, due to backflow stabilization [89, Section 2.2.1] applied at both ends to maintain well-posedness.

Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.

D. Kamensky / Computers and Mathematics with Applications xxx (xxxx) xxx

11

Fig. 10. Several snapshots of the valve deformation, with fluid velocity magnitude plotted on a slice.

the valve (simulating systole), then switches to −105 Ba, to close the valve under physiological pressure (simulating diastole). To make this example friendlier to users with limited available hardware, the initial systolic phase is shortened to 0.1 s, roughly corresponding to exercise conditions (cf. [90]). A more complicated pressure waveform, interpolating tabulated data, would be straightforward to implement in FEniCS as a UserExpression, but would add unnecessary complexity, as it would not be relevant to the present goals. (In fact, the accelerated pulse and abrupt pressure change considered here is a much more strenuous test of the numerical methods, thereby better serving the current purposes.) Remark 6. A crucial implementation detail for this example is the configuration of the linear solver used to compute fluid increments within the block iteration. Three-dimensionality makes direct solvers impractically-slow, while the scaling of stabilization parameters needed to maintain robustness with large pressure jumps (as introduced in [1, Section 4.4.2], analyzed in [3, Appendix A], and found necessary for related methods in [41,91,92]) makes it essentially impossible to meet conventional convergence criteria for Krylov solvers using standard preconditioners. However, we find, empirically, that under-solving the fluid increments by taking a fixed number of GMRES [93] iterations – even with only marginal reduction in relative residual for the incremental linear system – is sufficient to converge the outer nonlinear iteration. Fig. 10 shows several snapshots of the coupled fluid and structure solutions, computed at a fluid resolution tractable for modern desktop workstations. Fig. 11 plots flow rate through the valve as a function of time, illustrating that the valve correctly cuts off flow during diastole (which is not a trivial result [1, Section 4.4.2]). The decaying oscillation in flow rate after closure is an expected feature of the solution, emanating from the physics of the problem, as discussed in [1, Section 5.4.4], although it is amplified artificially by the assumption of a rigid artery wall (cf. [2]), the sudden change from systole to diastole, and the absence of a resistance boundary condition at the outflow. We encourage readers to reproduce the computations themselves, using the publicly-available code and data, for more thorough inspection of the results. 7. Conclusion The implementation of an open, transparent, and extensible implementation of DAL-based IMGA is an important advance beyond earlier work on the method, which relied on closed and opaque research codes that required extensive modifications for new problems. The discretizations of the fluid and structure subproblems in CouDALFISh are generated automatically, using FEniCS and tIGAr, from high-level mathematical descriptions in UFL. The coupling between these Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.

12

D. Kamensky / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 11. Volumetric flow rate through the valve as a function of time. The artificially-severe step change from systolic to diastolic pressures at t = 0.1 s demonstrates the robustness of the method and implementation.

discretizations, however, leaves much room for further automation; at present it is implemented manually, which is amenable to neither interpretation at-a-glance nor easy modification of the coupling formulation. The implementation of meshfree nonlocal contact mechanics in ShNAPr is also manual. Further, the coupling and contact implementations introduce serial bottlenecks and collective communication, with only the fluid assembly and solution procedures leveraging the parallel capabilities of FEniCS. The difference in dimension between the fluid and shell structure subproblems implies that overall computational cost is overwhelmingly dominated by the fluid, so significant speedups are still attainable over fully serial computation, but a simple application of Amdahl’s law shows that the serial bottlenecks will dominate wall-clock time in a massively-parallel setting. We believe that the most productive path forward for addressing these limitations would be to devise and implement more flexible abstractions within the FEniCS framework for user-defined quadrature rules specified in physical space. The source code repositories for CouDALFISh [37], VarMINT [55], and ShNAPr [51] (along with their dependency, tIGAr [31]) will be updated at the cited URLs. References to specific code in this paper refer to the state of these repositories at the time of submission, which can be recovered using Git, if needed. Some numerical results shown here were computed at earlier states. We welcome any bug reports, suggestions, or contributions to these projects from the community. CRediT authorship contribution statement David Kamensky: Conceptualization, Methodology, Software, Validation, Investigation, Data curation, Writing - original draft, Writing - review & editing, Visualization, Supervision, Project administration. Acknowledgments This work was supported by start-up funds from the University of California, San Diego. FEniCS and tIGAr make extensive use of the library PETSc [94–96] (and its Python bindings [97]) for parallel linear algebra. Some examples from this paper construct spline geometry using igakit [98], which was originally developed for use with the open-source IGA library PetIGA [99,100]. References [1] D. Kamensky, M.-C. Hsu, D. Schillinger, J.A. Evans, A. Aggarwal, Y. Bazilevs, M.S. Sacks, T.J.R. Hughes, An immersogeometric variational framework for fluid–structure interaction: Application to bioprosthetic heart valves, Comput. Methods Appl. Mech. Engrg. 284 (2015) 1005–1053. [2] M.-C. Hsu, D. Kamensky, Y. Bazilevs, M.S. Sacks, T.J.R. Hughes, Fluid–structure interaction analysis of bioprosthetic heart valves: significance of arterial wall deformation, Comput. Mech. 54 (2014) 1055–1071. [3] D. Kamensky, M.-C. Hsu, Y. Yu, J.A. Evans, M.S. Sacks, T.J.R. Hughes, Immersogeometric cardiovascular fluid–structure interaction analysis with divergence-conforming B-splines, Comput. Methods Appl. Mech. Engrg. 314 (2017) 408–472, Special Issue on Biological Systems Dedicated to William S. Klug. [4] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (2005) 4135–4195. [5] J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs, Isogeometric Analysis: Toward Integration of CAD and FEA, Wiley, Chichester, 2009. [6] E. Rank, M. Ruess, S. Kollmannsberger, D. Schillinger, A. Düster, Geometric modeling, isogeometric analysis and the finite cell method, Comput. Methods Appl. Mech. Engrg. 249–252 (2012) 104–115. [7] J. Parvizian, A. Düster, E. Rank, Finite cell method, Comput. Mech. 41 (1) (2007) 121–133. [8] Y. Yu, D. Kamensky, M.-C. Hsu, X.Y. Lu, Y. Bazilevs, T.J.R. Hughes, Error estimates for projection-based dynamic augmented Lagrangian boundary condition enforcement, with application to fluid–structure interaction, Math. Models Methods Appl. Sci. 28 (12) (2018) 2457–2509.

Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.

D. Kamensky / Computers and Mathematics with Applications xxx (xxxx) xxx

13

[9] D. Kamensky, J.A. Evans, M.-C. Hsu, Stability and conservation properties of collocated constraints in immersogeometric fluid-thin structure interaction analysis, Commun. Comput. Phys. 18 (2015) 1147–1180. [10] M.R. Hestenes, Multiplier and gradient methods, J. Optim. Theory Appl. 4 (5) (1969) 303–320. [11] M.J.D. Powell, A method for nonlinear constraints in minimization problems, in: R. Fletcher (Ed.), Optimization, Academic Press, New York, 1969, pp. 283–298. [12] M.-C. Hsu, D. Kamensky, F. Xu, J. Kiendl, C. Wang, M.C.H. Wu, J. Mineroff, A. Reali, Y. Bazilevs, M.S. Sacks, Dynamic and fluid–structure interaction simulations of bioprosthetic heart valves using parametric design with T-splines and Fung-type material models, Comput. Mech. (2015) 1–15. [13] F. Xu, S. Morganti, R. Zakerzadeh, D. Kamensky, F. Auricchio, A. Reali, T.J.R. Hughes, M.S. Sacks, M.-C. Hsu, A framework for designing patientspecific bioprosthetic heart valves using immersogeometric fluid–structure interaction analysis, Int. J. Numer. Methods Biomed. Eng. 34 (4) (2018) e2938. [14] M.C.H. Wu, H.M. Muchowski, E.L. Johnson, M.R. Rajanna, M.-C. Hsu, Immersogeometric fluid–structure interaction modeling and simulation of transcatheter aortic valve replacement, Comput. Methods Appl. Mech. Engrg. 357 (2019) 112556. [15] M.C.H. Wu, D. Kamensky, C. Wang, A.J. Herrema, F. Xu, M.S. Pigazzini, A. Verma, A.L. Marsden, Y. Bazilevs, M.-C. Hsu, Optimizing fluid–structure interaction systems with immersogeometric analysis and surrogate modeling: Application to a hydraulic arresting gear, Comput. Methods Appl. Mech. Engrg. 316 (2017) 668–693, Special Issue on Isogeometric Analysis: Progress and Challenges. [16] D. Kamensky, J.A. Evans, M.-C. Hsu, Y. Bazilevs, Projection-based stabilization of interface lagrange multipliers in immersogeometric fluid– thin structure interaction analysis, with application to heart valve modeling, Comput. Math. Appl. 74 (9) (2017) 2068–2088, Advances in Mathematics of Finite Elements, honoring 90th birthday of Ivo Babuška. [17] C.S. Peskin, The immersed boundary method, Acta Numer. 11 (2002) 479–517. [18] R. Mittal, G. Iaccarino, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 239–261. [19] F. Sotiropoulos, X. Yang, Immersed boundary methods for simulating fluid–structure interaction, Prog. Aerosp. Sci. 65 (2014) 1–21. [20] B.E. Griffith, N.A. Patankar, Immersed methods for fluid–structure interaction, Annu. Rev. Fluid Mech. 52 (1) (2020) null. [21] A. Logg, K.-A. Mardal, G.N. Wells, et al., in: A. Logg, K.-A. Mardal, G.N. Wells (Eds.), Automated Solution of Differential Equations by the Finite Element Method, Springer, 2012. [22] M.S. Alnæs, A. Logg, K.B. Ølgaard, M.E. Rognes, G.N. Wells, Unified form language: A domain-specific language for weak formulations of partial differential equations, ACM Trans. Math. Software 40 (2) (2014) 9:1–9:37. [23] R.C. Kirby, A. Logg, A compiler for variational forms, ACM Trans. Math. Software 32 (3) (2006) 417–444. [24] A. Logg, G.N. Wells, DOLFIN: Automated finite element computing, ACM Trans. Math. Software 37 (2) (2010) 20:1–20:28. [25] C.N. Richardson, N. Sime, G.N. Wells, Scalable computation of thermomechanical turbomachinery problems, Finite Elem. Anal. Des. 155 (2019) 32–42. [26] G.-T. Bercea, A.T.T. McRae, D.A. Ham, L. Mitchell, F. Rathgeber, L. Nardi, F. Luporini, P.H.J. Kelly, A structure-exploiting numbering algorithm for finite elements on extruded meshes, and its performance evaluation in Firedrake, Geosci. Model Dev. 9 (10) (2016) 3803–3815. [27] M. Homolya, L. Mitchell, F. Luporini, D. Ham, TSFC: A structure-preserving form compiler, SIAM J. Sci. Comput. 40 (3) (2018) C401–C428. [28] M. Homolya, D. Ham, A parallel edge orientation algorithm for quadrilateral meshes, SIAM J. Sci. Comput. 38 (5) (2016) S48–S61. [29] A. McRae, G. Bercea, L. Mitchell, D. Ham, C. Cotter, Automated generation and symbolic manipulation of tensor product finite elements, SIAM J. Sci. Comput. 38 (5) (2016) S25–S47. [30] D. Kamensky, Y. Bazilevs, TIGAr: Automating isogeometric analysis with FEniCS, Comput. Methods Appl. Mech. Engrg. 344 (2019) 477–498. [31] https://github.com/david-kamensky/tIGAr. tIGAr source code. [32] H. Yang, B.E. Abali, D. Timofeev, W.H. Müller, Determination of metamaterial parameters by means of a homogenization approach based on asymptotic analysis, Contin. Mech. Thermodyn. (2019) http://dx.doi.org/10.1007/s00161-019-00837-4. [33] M.F.P. ten Eikelder, Y. Bazilevs, I. Akkerman, A theoretical framework for discontinuity capturing: Joining variational multiscale analysis and variation entropy theory, Comput. Methods Appl. Mech. Engrg. (2019) 112664. [34] Y. Bazilevs, D. Kamensky, G. Moutsanidis, S. Shende, Residual-based shock capturing in solids, Comput. Methods Appl. Mech. Engrg. 358 (2020) 112638. [35] J.A. Evans, D. Kamensky, Y. Bazilevs, Variational multiscale modeling with discretely divergence-free subscales, in: Society of Engineering Science Technical Meeting, 2019. [36] W. Zhang, T. Bui-Thanh, M.S. Sacks, A machine learning approach for soft tissue remodeling, in: Proceedings of FEniCS’19, 2019. [37] https://github.com/david-kamensky/CouDALFISh. CouDALFISh source code. [38] J.C. Simo, P. Wriggers, R.L. Taylor, A perturbed Lagrangian formulation for the finite element solution of contact problems, Comput. Methods Appl. Mech. Engrg. 50 (2) (1985) 163–180. [39] J. Chung, G.M. Hulbert, A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method, J. Appl. Mech. 60 (1993) 371–735. [40] K.E. Jansen, C.H. Whiting, G.M. Hulbert, A generalized-α method for integrating the filtered Navier-Stokes equations with a stabilized finite element method, Comput. Methods Appl. Mech. Engrg. 190 (2000) 305–319. [41] L. Boilevin-Kayl, M. Fernández, J.-F. Gerbeau, A loosely coupled scheme for fictitious domain approximations of fluid-structure interaction problems with immersed thin-walled structures, SIAM J. Sci. Comput. 41 (2) (2019) B351–B374. [42] D. Goldstein, R. Handler, L. Sirovich, Modeling a no-slip flow boundary with an external force field, J. Comput. Phys. 105 (2) (1993) 354–366. [43] M.J. Borden, M.A. Scott, J.A. Evans, T.J.R. Hughes, Isogeometric finite element data structures based on Bézier extraction of NURBS, Internat. J. Numer. Methods Engrg. 87 (2011) 15–47. [44] M.A. Scott, M.J. Borden, C.V. Verhoosel, T.W. Sederberg, T.J.R. Hughes, Isogeometric finite element data structures based on Bézier extraction of T-splines, Internat. J. Numer. Methods Engrg. 88 (2011) 126–156. [45] D. Schillinger, P.K. Ruthala, L.H. Nguyen, Lagrange extraction and projection for NURBS basis functions: A direct link between isogeometric and standard nodal finite element formulations, Internat. J. Numer. Methods Engrg. 108 (6) (2016) 515–534. [46] J. Kiendl, K.-U. Bletzinger, J. Linhard, R. Wüchner, Isogeometric shell analysis with Kirchhoff–Love elements, Comput. Methods Appl. Mech. Engrg. 198 (2009) 3902–3914. [47] J. Kiendl, Isogeometric Analysis and Shape Optimal Design of Shell Structures (PhD thesis), Lehrstuhl für Statik, Technische Universität München, 2011. [48] J. Kiendl, M.-C. Hsu, M.C.H. Wu, A. Reali, Isogeometric Kirchhoff–Love shell formulations for general hyperelastic materials, Comput. Methods Appl. Mech. Engrg. 291 (0) (2015) 280–303. [49] A. Buganza Tepole, H. Kabaria, K.-U. Bletzinger, E. Kuhl, Isogeometric Kirchhoff–Love shell formulations for biological membranes, Comput. Methods Appl. Mech. Engrg. (0) (2015). [50] J.S. Hale, M. Brunetti, S.P.A. Bordas, C. Maurini, Simple and extensible plate and shell finite element models through automatic code generation tools, Comput. Struct. 209 (2018) 163–181. [51] https://github.com/david-kamensky/ShNAPr. ShNAPr source code.

Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.

14

D. Kamensky / Computers and Mathematics with Applications xxx (xxxx) xxx

[52] D. Kamensky, M. Behzadinasab, J.T. Foster, Y. Bazilevs, Peridynamic modeling of frictional contact, J. Peridynamics Nonlocal Model. (2019). [53] D. Kamensky, F. Xu, C.-H. Lee, J. Yan, Y. Bazilevs, M.-C. Hsu, A contact formulation based on a volumetric potential: Application to isogeometric simulations of atrioventricular valves, Comput. Methods Appl. Mech. Engrg. 330 (2018) 522–546. [54] A. Stavrev, The Role of Higher-Order Geometry Approximation and Accurate Quadrature in NURBS Based Immersed Boundary Methods (Master’s thesis), Technische Universität München, Munich, Germany, 2012. [55] https://github.com/david-kamensky/VarMINT. VarMINT source code. [56] Y. Bazilevs, V.M. Calo, J.A. Cottrel, T.J.R. Hughes, A. Reali, G. Scovazzi, Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows, Comput. Methods Appl. Mech. Engrg. 197 (2007) 173–201. [57] Y. Bazilevs, V.M. Calo, T.J.R. Hughes, Y. Zhang, Isogeometric fluid–structure interaction: theory, algorithms, and computations, Comput. Mech. 43 (2008) 3–37. [58] I. Akkerman, Y. Bazilevs, V.M. Calo, T.J.R. Hughes, S. Hulshoff, The role of continuity in residual-based variational multiscale modeling of turbulence, Comput. Mech. 41 (2008) 371–378. [59] Y. Bazilevs, C. Michler, V.M. Calo, T.J.R. Hughes, Isogeometric variational multiscale modeling of wall-bounded turbulent flows with weakly enforced boundary conditions on unstretched meshes, Comput. Methods Appl. Mech. Engrg. 199 (2010) 780–790. [60] Y. Bazilevs, I. Akkerman, Large eddy simulation of turbulent Taylor–Couette flow using isogeometric analysis and the residual–based variational multiscale method, J. Comput. Phys. 229 (2010) 3402–3414. [61] K. Takizawa, Y. Bazilevs, T.E. Tezduyar, Space–time and ALE-VMS techniques for patient-specific cardiovascular fluid–structure interaction modeling, Arch. Comput. Methods Eng. 19 (2012) 171–225. [62] Y. Bazilevs, M.-C. Hsu, K. Takizawa, T.E. Tezduyar, ALE–VMS and ST–VMS methods for computer modeling of wind-turbine rotor aerodynamics and fluid–structure interaction, Math. Models Methods Appl. Sci. 22 (2012) 1230002. [63] Y. Bazilevs, M.-C. Hsu, M.A. Scott, Isogeometric fluid–structure interaction analysis with emphasis on non-matching discretizations, and with application to wind turbines, Comput. Methods Appl. Mech. Engrg. 249–252 (2012) 28–41. [64] M.-C. Hsu, I. Akkerman, Y. Bazilevs, Wind turbine aerodynamics using ALE–VMS: Validation and the role of weakly enforced boundary conditions, Comput. Mech. 50 (2012) 499–511. [65] A. Korobenko, M.-C. Hsu, I. Akkerman, Y. Bazilevs, Aerodynamic simulation of vertical-axis wind turbines, J. Appl. Mech. 81 (2014) 021011. [66] https://github.com/QimingZhu1992/IlliniFlow. IlliniFlow source code. [67] Q. Zhu, J. Yan, A moving-domain CFD solver in FEniCS with applications to tidal turbine simulations in turbulent flows, Comput. Math. Appl. (2019). [68] T.M. Shih, C.H. Tan, B.C. Hwang, Effects of grid staggering on numerical schemes, Internat. J. Numer. Methods Fluids 9 (2) (1989) 193–212. [69] T.M. van Opstal, J. Yan, C. Coley, J.A. Evans, T. Kvamsdal, Y. Bazilevs, Isogeometric divergence-conforming variational multiscale formulation of incompressible turbulent flows, Comput. Methods Appl. Mech. Engrg. 316 (2017) 859–879, Special Issue on Isogeometric Analysis: Progress and Challenges. [70] Y. Bazilevs, T.J.R. Hughes, Weak imposition of Dirichlet boundary conditions in fluid mechanics, Comput. & Fluids 36 (2007) 12–26. [71] J.A. Evans, Divergence-Free B-Spline Discretizations for Viscous Incompressible Flows (Ph.D. thesis), University of Texas at Austin, Austin, Texas, United States, 2011. [72] M.-C. Hsu, Y. Bazilevs, V.M. Calo, T.E. Tezduyar, T.J.R. Hughes, Improving stability of stabilized and multiscale formulations in flow simulations at small time steps, Comput. Methods Appl. Mech. Engrg. 199 (2010) 828–840. [73] R. Codina, J. Principe, O. Guasch, S. Badia, Time dependent subscales in the stabilized finite element approximation of incompressible flow problems, Comput. Methods Appl. Mech. Engrg. 196 (2007) 2413–2430. [74] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199–259. [75] T.E. Tezduyar, Stabilized finite element formulations for incompressible flow computations, Adv. Appl. Mech. 28 (1992) 1–44. [76] T.E. Tezduyar, Y. Osawa, Finite element stabilization parameters computed from element matrices and vectors, Comput. Methods Appl. Mech. Engrg. 190 (2000) 411–430. [77] M.-C. Hsu, I. Akkerman, Y. Bazilevs, Finite element simulation of wind turbine aerodynamics: Validation study using NREL phase VI experiment, Wind Energy (2014). [78] Y. Bazilevs, K. Takizawa, T.E. Tezduyar, Computational Fluid–Structure Interaction: Methods and Applications, Wiley, Chichester, 2013. [79] T. Belytschko, H. Stolarski, W.K. Liu, N. Carpenter, J.S.-J. Ong, Stress projection for membrane and shear locking in shell finite elements, Comput. Methods Appl. Mech. Engrg. 51 (1985) 221–258. [80] C. Hesch, A.J. Gil, A. Arranz Carreño, J. Bonet, On continuum immersed strategies for fluid-structure interaction, Comput. Methods Appl. Mech. Engrg. 247–248 (2012) 51–64. [81] A.J. Gil, A. Arranz Carreño, J. Bonet, O. Hassan, An enhanced immersed structural potential method for fluid–structure interaction, J. Comput. Phys. 250 (2013) 178–205. [82] T. Wick, Flapping and contact FSI computations with the fluid–solid interface-tracking/interface-capturing technique and mesh adaptivity, Comput. Mech. 53 (1) (2014) 29–43. [83] T. Tezduyar, S. Aliabadi, M. Behr, A. Johnson, S. Mittal, Parallel finite-element computation of 3D flows, Computer 26 (10) (1993) 27–36. [84] A.A. Johnson, T.E. Tezduyar, Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces, Comput. Methods Appl. Mech. Engrg. 119 (1994) 73–94. [85] K. Stein, T. Tezduyar, R. Benney, Mesh moving techniques for fluid–structure interactions with large displacements, J. Appl. Mech. 70 (2003) 58–63. [86] K. Stein, T.E. Tezduyar, R. Benney, Automatic mesh update with the solid-extension mesh moving technique, Comput. Methods Appl. Mech. Engrg. 193 (2004) 2019–2032. [87] The CGAL Project, CGAL User and Reference Manual, 4.13 ed., CGAL Editorial Board, 2018. [88] B.E. Abali, Computational Reality, Springer Singapore, 2017. [89] M. Esmaily-Moghadam, Y. Bazilevs, T.-Y. Hsia, I.E. Vignon-Clementel, A.L. Marsden, and Modeling of Congenital Hearts Alliance (MOCHA), A comparison of outlet boundary treatments for prevention of backflow divergence with relevance to blood flow simulations, Comput. Mech. 48 (2011) 277–291. [90] C.A. Taylor, T.J.R. Hughes, C.K. Zarins, Effect of excercise on hemodynamic conditions in the abdominal aorta, J. Vas. Surg. 29 (1999) 1077–1089. [91] H. Casquero, C. Bona-Casas, H. Gomez, NURBS-based numerical proxies for red blood cells and circulating tumor cells in microscale blood flow, Comput. Methods Appl. Mech. Engrg. 316 (2017) 646–667, Special Issue on Isogeometric Analysis: Progress and Challenges. [92] L. Boilevin-Kayl, M.A. Fernández, J.-F. Gerbeau, Numerical methods for immersed FSI with thin-walled structures, Comput. & Fluids 179 (2019) 744–763. [93] Y. Saad, M. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (1986) 856–869.

Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.

D. Kamensky / Computers and Mathematics with Applications xxx (xxxx) xxx

15

[94] S. Balay, S. Abhyankar, M.F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W.D. Gropp, D. Kaushik, M.G. Knepley, L.C. McInnes, K. Rupp, B.F. Smith, S. Zampini, H. Zhang, PETSc Web page, 2015, http://www.mcs.anl.gov/petsc. [95] S. Balay, S. Abhyankar, M.F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W.D. Gropp, D. Kaushik, M.G. Knepley, L.C. McInnes, K. Rupp, B.F. Smith, S. Zampini, H. Zhang, PETSc Users Manual, Technical Report ANL-95/11 - Revision 3.6, Argonne National Laboratory, 2015. [96] S. Balay, W.D. Gropp, L.C. McInnes, B.F. Smith, Efficient management of parallelism in object oriented numerical software libraries, in: E. Arge, A.M. Bruaset, H.P. Langtangen (Eds.), Modern Software Tools in Scientific Computing, Birkhäuser Press, 1997, pp. 163–202. [97] L.D. Dalcin, R.R. Paz, P.A. Kler, A. Cosimo, Parallel distributed computing using Python, Adv. Water Resour. 34 (9) (2011) 1124–1139, New Computational Methods and Software Tools. [98] igakit: Toolkit for IsoGeometric Analysis (IGA). https://bitbucket.org/dalcinl/igakit. [99] L. Dalcin, N. Collier, P. Vignal, A.M.A. Côrtes, V.M. Calo, PetIGA: A framework for high-performance isogeometric analysis, Comput. Methods Appl. Mech. Engrg. 308 (2016) 151–181. [100] A.F. Sarmiento, A.M.A. Côrtes, D.A. Garcia, L. Dalcin, N. Collier, V.M. Calo, PetIGA-MF: A multi-field high-performance toolbox for structure-preserving B-splines spaces, J. Comput. Sci. 18 (Supplement C) (2017) 117–131.

Please cite this article as: D. Kamensky, Open-source immersogeometric analysis of fluid–structure interaction using FEniCS and tIGAr, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.01.023.