FEIBM_ANS
 All Files Pages
The step-feibm tutorial program
Table of contents
  1. Introduction
  2. The commented program
  1. Results
  2. The plain program

Introduction

We present the implementation of a solution scheme for fluid-structure interaction problems via the finite element software library deal.II. Specifically, we implement an immersed finite element method in which two independent discretizations are used for the fluid and immersed deformable body. In this type of formulation the support of the equations of motion of the fluid is extended to cover the union of the solid and fluid domains. The equations of motion over the extended solution domain govern the flow of a fluid under the action of a body force field. This body force field informs the fluid of the presence of the immersed solid. The velocity field of the immersed solid is the restriction over the immersed domain of the velocity field in the extended equations of motion. The focus of this paper is to show how the determination of the motion of the immersed domain is carried out in practice. We show that our implementation is automatically obtained from the choice of finite element spaces over the immersed solid and the extended fluid domains. We present a few results concerning the accuracy of the proposed method. For a careful discussion of the proposed approach see Heltai, L. and F. Costanzo (2012), "Variational implementation of immersed finite element methods," Computer Methods in Applied Mechanics and Engineering, 229–232, p. 110–127, which we will denoted by [HC2012].

Governing Equations

Fluid and Solid Domains. Referring to the figure below, $B_{t}$ represents the configuration of a regular body at time $t$. $B_{t}$ is a (possibly multiply connected) proper subset of a fixed control volume $\Omega$: $\Omega\setminus B_{t}$ is filled by a fluid and $B_{t}$ is an immersed body. Following common practice in continuum mechanics, we refer to $B_{t}$ as the current configuration of the immersed body. By contrast, we denote by $B$ the reference configuration of the immersed body. We denote the position of points in $B$ by $s$, whereas we denote the position at time $t$ of a point $P$ in $\Omega$ by $x_{P}(t)$. For the examples in the results section we take $B$ to coincide with the initial configuration $B_{0}$.

step-feibm.geometry.png

Motion of the Immersed Body. The motion of the solid is described by a function $\chi(s,t)$, which gives the position $x$ at time $t$ of a particle of the immersed domain with position $s$ in the reference configuration $B$, i.e., $x = \chi(s,t)$ for $s \in B$. The deformation gradient is the tensor $F = \nabla_{s} \chi(s,t)$, where $\nabla_{s}$ denotes the gradient operator relative to position in the reference configuration. We denote the determinant of $F(s,t)$ by $J(s,t)$ and we assume that $J(s,t) > 0$ for all $s \in B$ and $t$. From a practical viewpoint, as it is often done in solids mechanics, we describe the motion of the body through its displacement function. We denote the displacement by $w(s,t)$. The relation between $\chi(s,t)$ and $w$ is as follows:

\[ x = s + w(s,t). \]

From the above relation we have that $F(s,t) = I + \nabla_{s} w(s,t)$, where $I$ is the identity tensor.

Constitutive Equations. The fluid is assumed to be Newtonian so that its Cauchy stress is $\sigma_{f} = -p I + \eta [\nabla_{x} u + (\nabla_{x} u)^{T}]$, where $p$ is the Lagrange multiplier for the enforcement of incompressibility, $\eta$ is the dynamic viscosity, $\nabla_{x}$ is the gradient operator relative to position in the current configuration, and $u$ is the velocity field. For a Newtonian fluid, $p$ is also the pressure (mean normal stress). The immersed solid is assumed to be incompressible and viscoelastic, with the viscous response identical to that of the fluid. The elastic response is assumed to be admit a strain energy $W_{s}^{e}(F)$, which we assumed to be a convex function of the deformation gradient. Hence the Cauchy stress in the solid is given by $\sigma_{s} = - p I + \eta [\nabla_{x} u + (\nabla_{x} u)^{T}] + \sigma_{s}^{e}$, with $\sigma_{s}^{e} = J^{-1} P_{s}^{e} F^{-T}$, where $P_{s}^{e} = \partial W_{s}^{e}/\partial F$ is the first Piola-Kirchhoff stress tensor of the solid.

Velocity field and Displacement of the Immersed Body. The velocity field $u(x,t)$, with $x \in \Omega$, represents that velocity of the particle occupying the point $x$ at time $t$. As such, this field describes the velocity of the fluid for $x \in \Omega/B_{t}$ and the velocity of the solid for $x \in B_{t}$. Therefore, using the displacement function $w(s,t)$, for $s \in B$ and for al $t$, we have

\[ u(x,t)\big|_{x = s + w(s,t)} = \frac{\partial w(s,t)}{\partial t}. \]

Conditions at the Boundary of the Immersed Body. The boundary of the immersed body is viewed as a material surface. Therefore, the balance of linear momentum requires that the traction field be continuous across the boundary of the immersed body. In addition, we assume that there is no slip between the immersed body and the surrounding fluid.

Governing Equations: Strong Form. The motion of the system is governed by the following three equations, which, respectively, represent the balance of linear momentum, balance of mass accounting for incompressibility, and velocity compatibility:

\[ \nabla \cdot \sigma(x,t) + \rho b = \rho \biggl[\frac{\partial u(x,t)}{\partial t} + (\nabla u(x,t)) u(x,t) \biggr] \quad {\rm in}~\Omega,\quad \nabla \cdot u(x,t) = 0 \quad {\rm in}~\Omega, \quad u(x,t)\big|_{x = s + w(s,t)} = \frac{\partial w(s,t)}{\partial t} \quad {\rm in}~B, \]

where " $\nabla_{x} \cdot$" denotes the divergence operator (relative to position in the current configuration), $\rho$ is the density (here assumed to be a constant), $b$ is a (prescribed) body force field, and where $\sigma(x,t)$ is the Cauchy stress field in the entire domain $\Omega$, i.e., $\sigma(x,t) = \sigma_{f}(x,t)$ for $x \in \Omega/B_{t}$ and $\sigma(x,t) = \sigma_{s}(x,t)$ for $x \in B_{t}$.

As far as boundary conditions are concerned, we assume that a velocity and a traction distribution are prescribed on complementary subsets of the boundary of $\Omega$. Specifically, letting $\partial \Omega_{D} \cup \partial \Omega_{N} = \partial \Omega$, with $\partial \Omega_{D} \cap \partial \Omega_{N} = \emptyset$,

\[ u(x,t) = u_{g}(x,t)~{\rm for}~ x\in \partial \Omega_{D} \quad{\rm and}\quad \sigma(x,t) n(x,t) = \tau_{g}(x,t)~{\rm for}~x \in \partial \Omega_{N}, \]

where $u_{g}$ and $\tau_{g}$ are prescribed velocity and tractions distributions, and where $n$ denotes the outward unit normal to the $\partial \Omega$.

Weak Formulation

The primary unknowns of the problem are the velocity field $u(x,t)$ in $\Omega$, the displacement field $w(s,t)$ in $B$, and the Lagrange multiplier field $p$. For these we select appropriate function spaces $V$, $Y$, and $Q$ (for details see [HC2012]), respectively. With this in mind, the weak formulation from which the discrete formulation is derived is as follows:

\[ \int_{\Omega} \rho(\dot{u}(x,t) - b(x,t)] \cdot v(x,t) - \int_{\Omega} p(x,t) \, \nabla_{x} \cdot v(x,t) + \int_{\Omega} \eta [\nabla u(x,t) + (\nabla u(x,t))^{T}] \cdot \nabla v(x,t) - \int_{\partial \Omega_{N}} \tau_{g}(x,t) \cdot v(x,t) + \int_{B} P_{s}^{e}(s,t) F^{T}(s,t) \cdot [\nabla_{x} v(x,t)\big|_{x = s + w(s,t)}] = 0, \]

\[ \int_{\Omega} (\nabla_{x} \cdot u ) \; q = 0, \]

\[ \Phi_{B} \int_{B} \biggl(\frac{\partial w(s,t)}{\partial t} - u(x,t)\big|_{x = s + w(s,t)} \biggr) \cdot y = 0, \]

for all $v \in V_{0}$, $y \in Y$, and $q \in Q$. We observe that, since $F = I + \nabla_{s} w(s,t)$, $P_{s}^{e}$ is a function of the displacement gradient $\nabla_{s} w$. Finally, $\Phi_{B}$ is a constant needed to ensure that the dimensions of the equations are homogeneous with the dimensions of the rest of the equations.

The above weak formulation can be viewed as consisting of a typical formulation for the Navier-Stokes equations with some nonstandard terms. The latter are the last term on the left-hand of the first equation, which can be interpreted as a body force distribution that informs the fluid of the presence of the solid, and the very last equation, which is crucial for the tracking of the motion of the immersed body. What makes these terms unusual is the fact that, in the implementation of the discrete formulation, they require the evaluation of integrals over the triangulation of the domain $B$ of functions defined over the triangulation of the domain $\Omega$. In the next section we provide a discussion of how this is implemented in practice.

In [HC2012] we have discussed in detail the fact that the discrete formulation derived from the above weak formulation requires a careful treatment of the term

\[ \int_{B} P_{s}^{e}(s,t) F^{T}(s,t) \cdot [\nabla_{x} v(x,t)\big|_{x = s + w(s,t)}] \]

to yield a stable formulation. Specifically, it was shown that the above term should be viewed as the composition of two operators: (i) an "elastic" operator defined through the following term

\[ \int_{B} P_{s}^{e}(s,t) \cdot \nabla_{s} y, \]

and (ii) a "spread" operator whose definition stems from the last equation of the of the weak formulation (for details see, in particular, Remark 9 and Theorem 4 in [HC2012]). The code in the present example includes the determination of the elastic and spread operators as well as their composition.

Implementation

The implementation of the various terms that are common to the Navier-Stokes equations is done in a standard fashion and will not be discussed here. As alluded above, here we limit ourselves to the description of how to carry out the integration over $B$ of functions that are available through their finite element representation over the triangulation of $\Omega$.

Referring to the third equation in the weak formulation, let's consider the term

\[ \Phi_{B} \int_{B} u_{h}(x,t)\big|_{x = s + w_{h}(s,t)} \cdot y_{h}, \]

where $u_{h}$ denotes the finite element representation of $u$ given by interpolation functions supported over the triangulation $\Omega_{h}$ of $\Omega$, and where $y_{h}$ denotes shape functions supported over the triangulation $B_{h}$ of the domain $B$: the construction of the above term draws information from two independent triangulations. The above integral is computed by summing contributions from each cell $K$ of $B_{h}$. Each of these contributions is a sum over the $N_{Q}$ quadrature points. We illustrate this via the following figure:

step-feibm.integration.png

In our code, for each cell of the immersed body, we start by determining the position of the quadrature points of the element corresponding to the cell at hand. The position of the quadrature point is determined both relative to the reference unit element and relative to the global coordinate system adopted for the calculation, through the mappings:

\[ s_{K}: \hat{K} := [0,1]^{d} \mapsto K \in B_{h}, \quad{\rm and}\quad s + w_{h}: K \mapsto {\rm current~position~of~solid~cell}, \]

where $\hat{K}$ is the reference unit element and $d$ is the spatial dimension of the immersed solid. These maps allow us to determine the global coordinates of the quadrature points. These coordinates are then passed to a search algorithm that identifies the cells in $\Omega_{h}$ that contain the points in question. In turn, this identification allows us to evaluate the functions $v_h$. The overall operation is illustrated in the figure above where we show a cell of $B_{h}$ straddling four cells of $\Omega_{h}$ denoted fluid cells A–D. The quadrature points over the solid cell are denoted by filled circles. The contribution to the above integral due to the solid cell is then computed by summing the partial contributions corresponding to each of the fluid cells intersecting the solid cell in question. The implementation of an efficient search algorithm responsible for identifying the fluid cells intersecting an individual solid cell is the only technically challenging part of the procedure. We use the built in facilities of the deal.II library to perform this task. Once the fluid cells containing the quadrature points of a given solid cell are found, we determine the value of $v_{h}$ at the quadrature points using the interpolation infrastructure inherent in the finite element representation of fields defined over $\Omega_{h}$. The deal.II class we use for this implementation is the FEFieldFunction.

The commented program

Copyright (C) 2012 by Luca Heltai (1), Saswati Roy (2), and Francesco Costanzo (3)

(1) Scuola Internazionale Superiore di Studi Avanzati E-mail: luca..nosp@m.helt.nosp@m.ai@si.nosp@m.ssa..nosp@m.it (2) Center for Neural Engineering, The Pennsylvania State University E-Mail: sur16.nosp@m.4@ps.nosp@m.u.edu (3) Center for Neural Engineering, The Pennsylvania State University E-Mail: costa.nosp@m.nzo@.nosp@m.engr..nosp@m.psu..nosp@m.edu

This code was developed starting from the example step-33 of the deal.II FEM library.

This file is subject to QPL and may not be distributed without copyright and license information. Please refer to the webpage http://www.dealii.org/ -> License for the text and further information on this license.

Keywords: fluid-structure interaction, immersed method, finite elements, monolithic framework

Deal.II version: deal.II 7.2.pre

Include files

We include those elements of the deal.ii library whose functionality is needed for our purposes.

#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/point.h>
#include <deal.II/base/function.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/parsed_function.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/multithread_info.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/work_stream.h>
#include <deal.II/base/parallel.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/sparse_ilu.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/vector_view.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgp.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping_q_eulerian.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/fe_field_function.h>
#include <deal.II/numerics/data_out.h>

Elements of the C++ standard library

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <map>
#include <cmath>
#include <typeinfo>

Names imported into to the global namespace:

using namespace dealii;
using namespace ::Functions;
using namespace std;

This class collects all of the user-specified parameters, both pertaining to physics of the problem (e.g., shear modulus, dynamic viscosity), and to other numerical aspects of the simulation (e.g., names for the grid files, the specification of the boundary conditions). This class is derived from the ParameterHandler class in deal.II.

template <int dim>
class ProblemParameters :
public ParameterHandler
{
public:
ProblemParameters(int argc, char **argv);

Polynomial degree of the interpolation functions for the velocity of the fluid and the displacement of the solid. This parameters must be greater than one.

unsigned int degree;

Mass density of the fluid and of the immersed solid.

double rho;

Dynamic viscosity of the fluid and the immersed solid.

double eta;

Shear modulus of the neo-Hookean immersed solid.

double mu;

Time step.

double dt;

Final time.

double T;

Dimensional constant for the equation that sets the velocity of the solid provided by the time derivative of the displacement function equal to the velocity provided by the interpolation over the control volume.

double Phi_B;

Displacement of the immersed solid at the initial time.

ParsedFunction<dim> W_0;

Velocity over the control volume at the initial time.

ParsedFunction<dim> u_0;

Dirichlet boundary conditions for the control volume.

ParsedFunction<dim> u_g;

Body force field.

ParsedFunction<dim> force;

Mesh refinement level of the control volume.

unsigned int ref_f;

Mesh refinement level for the immersed domain.

unsigned int ref_s;

Maps of boundary value functions: 1st: a boundary indicator; 2nd: a boundary value function.

map<unsigned char, const Function<dim> *> boundary_map;

Maps of boundary value functions for homogeneous Dirichlet boundary values: 1st: a boundary indicator; 2nd: a zero boundary value function.

map<unsigned char, const Function<dim> *> zero_boundary_map;

Vector of flags for distinguishing between velocity and pressure degrees of freedom.

vector<bool> component_mask;

Map storing the boundary conditions: 1st: a boundary degree of freedom; 2nd: the value of field corresponding to the given degree of freedom.

map<unsigned int, double> boundary_values;

Flag to indicate whether or not the Newton iteration scheme must update the Jacobian at each iteration.

bool update_jacobian_continuously;

Flag to indicate whether or not the time integration scheme must be semi-implicit.

bool semi_implicit;

Flag to indicate how to deal with the non-uniqueness of the pressure field.

bool fix_pressure;

Flag to indicate whether homogeneous Dirichlet boundary conditions are applied.

bool all_DBC;

When set to true, an update of the system Jacobian is performed at the beginning of each time step.

bool update_jacobian_at_step_beginning;

Name of the mesh file for the solid domain.

string solid_mesh;

Name of the mesh file for the fluid domain.

string fluid_mesh;

Name of the output file.

string output_name;

The interval of timesteps between storage of output.

int output_interval;

Flag to indicating the use the spread operator.

bool use_spread;

List of available consitutive models for the elastic stress of the immersed solid:

INH_0: incompressible neo-Hookean with $P_{s}^{e} = \mu (F - F^{-T})$.

INH_1: incompressible neo-Hookean with $P_{s}^{e} = \mu F$.

CircumferentialFiberModel: $P_{s}^{e} = \mu F (e_{\theta} \otimes e_{\theta}) F^{-T}$.

enum MaterialModel {INH_0 = 1, INH_1, CircumferentialFiberModel};

Variable to identify the constitutive model for the immersed solid.

unsigned int material_model;

String to store the name of the finite element that will be used for the pressure field.

string fe_p_name;

Variable to store the center of the ring with circumferential fibers

Point<dim> ring_center;
};

Class constructor: the name of the input file is immersed_fem.prm. If the file does not exist at run time, it is created, and the simulation parameters are given default values.

template <int dim>
ProblemParameters<dim>::ProblemParameters(int argc, char **argv) :
W_0(dim),
u_0(dim+1),
u_g(dim+1),
force(dim+1),
component_mask(dim+1, true)
{

Declaration of parameters for the ParsedFunction objects in the class.

this->enter_subsection("W0");
ParsedFunction<dim>::declare_parameters(*this, dim);
this->leave_subsection();
this->enter_subsection("u0");
ParsedFunction<dim>::declare_parameters(*this, dim+1);
this->leave_subsection();
this->enter_subsection("ug");
ParsedFunction<dim>::declare_parameters(*this, dim+1);
this->declare_entry("Function expression", "if(y>.99, 1, 0); 0; 0",
Patterns::Anything());
this->leave_subsection();
this->enter_subsection("force");
ParsedFunction<dim>::declare_parameters(*this, dim+1);
this->leave_subsection();

Declaration of the class parameters and assignment of default values.

this->declare_entry (
"Velocity finite element degree",
"2",
Patterns::Integer(2,10)
);
this->declare_entry ("Fluid refinement", "3", Patterns::Integer());
this->declare_entry ("Solid refinement", "2", Patterns::Integer());
this->declare_entry ("Delta t", ".1", Patterns::Double());
this->declare_entry ("Final t", "1", Patterns::Double());
this->declare_entry ("Update J cont", "false", Patterns::Bool());
this->declare_entry (
"Force J update at step beginning",
"false",
Patterns::Bool()
);
this->declare_entry ("Density", "1", Patterns::Double());
this->declare_entry ("Viscosity", "1", Patterns::Double());
this->declare_entry ("Elastic modulus", "1", Patterns::Double());
this->declare_entry ("Phi_B", "1", Patterns::Double());
this->declare_entry ("Semi-implicit scheme", "true", Patterns::Bool());
this->declare_entry ("Fix one dof of p", "false", Patterns::Bool());
this->declare_entry (
"Solid mesh",
"meshes/solid_square.inp",
Patterns::Anything()
);
this->declare_entry (
"Fluid mesh",
"meshes/fluid_square.inp",
Patterns::Anything()
);
this->declare_entry ("Output base name", "out/square", Patterns::Anything());
this->declare_entry ("Dirichlet BC indicator", "0", Patterns::Integer(0,254));
this->declare_entry ("All Dirichlet BC", "true", Patterns::Bool());
this->declare_entry (
"Interval (of time-steps) between output",
"1",
Patterns::Integer()
);
this->declare_entry ("Use spread operator","true", Patterns::Bool());
this->declare_entry (
"Solid constitutive model",
"INH_0",
Patterns::Selection ("INH_0|INH_1|CircumferentialFiberModel"),
"Constitutive models available are: "
"INH_0: incompressible neo-Hookean with "
"P^{e} = mu (F - F^{-T}); "
"INH_1: incompressible Neo-Hookean with P^{e} = mu F; "
"CircumferentialFiberModel: incompressible with "
"P^{e} = mu F (e_{\\theta} \\otimes e_{\\theta}) F^{-T}; "
"this is suitable for annular solid comprising inextensible "
"circumferential fibers"
);
this->declare_entry (
"Finite element for pressure",
"FE_DGP",
Patterns::Selection("FE_DGP|FE_Q"),
"Select between FE_Q (Lagrange finite element space of "
"continuous, piecewise polynomials) or "
"FE_DGP(Discontinuous finite elements based on Legendre "
"polynomials) to approximate the pressure field"
);
this->enter_subsection (
"Equilibrium Solution of Ring with Circumferential Fibers"
);
this->declare_entry ("Inner radius of the ring", "0.25", Patterns::Double());
this->declare_entry ("Width of the ring", "0.0625", Patterns::Double());
this->declare_entry (
"Any edge length of the (square) control volume",
"1.",
Patterns::Double()
);
this->declare_entry (
"x-coordinate of the center of the ring",
"0.5",
Patterns::Double()
);
this->declare_entry (
"y-coordinate of the center of the ring",
"0.5",
Patterns::Double()
);
this->leave_subsection ();

Specification of the parmeter file. If no parameter file is specified in input, use the default one, else read each additional argument.

if(argc == 1)
this->read_input ("immersed_fem.prm");
else
for(int i=1; i<argc; ++i)
this->read_input(argv[i]);

Reading in the parameters.

this->enter_subsection ("W0");
W_0.parse_parameters (*this);
this->leave_subsection ();
this->enter_subsection ("u0");
u_0.parse_parameters (*this);
this->leave_subsection ();
this->enter_subsection ("ug");
u_g.parse_parameters (*this);
this->leave_subsection ();
this->enter_subsection ("force");
force.parse_parameters (*this);
this->leave_subsection ();
ref_f = this->get_integer ("Fluid refinement");
ref_s = this->get_integer ("Solid refinement");
dt = this->get_double ("Delta t");
T = this->get_double ("Final t");
update_jacobian_continuously = this->get_bool ("Update J cont");
update_jacobian_at_step_beginning = this->get_bool (
"Force J update at step beginning"
);
rho = this->get_double ("Density");
eta = this->get_double ("Viscosity");
mu = this->get_double ("Elastic modulus");
Phi_B = this->get_double ("Phi_B");
semi_implicit = this->get_bool ("Semi-implicit scheme");
fix_pressure = this->get_bool ("Fix one dof of p");
solid_mesh = this->get ("Solid mesh");
fluid_mesh = this->get ("Fluid mesh");
output_name = this->get ("Output base name");
unsigned char id = this->get_integer ("Dirichlet BC indicator");
all_DBC = this->get_bool ("All Dirichlet BC");
output_interval = this->get_integer (
"Interval (of time-steps) between output"
);
use_spread =this->get_bool ("Use spread operator");
if(this->get("Solid constitutive model") == string("INH_0"))
material_model = INH_0;
else if (this->get("Solid constitutive model") == string("INH_1"))
material_model = INH_1;
else if (this->get("Solid constitutive model")
== string("CircumferentialFiberModel"))
material_model = CircumferentialFiberModel;
else
cout
<< " No matching constitutive model found! Using INH_0."
<< endl;
component_mask[dim] = false;
static ZeroFunction<dim> zero (dim+1);
zero_boundary_map[id] = &zero;
boundary_map[id] = &u_g;
degree = this->get_integer ("Velocity finite element degree");
fe_p_name = this->get ("Finite element for pressure");
fe_p_name +="<dim>(" + Utilities::int_to_string(degree-1) + ")";
this->enter_subsection (
"Equilibrium Solution of Ring with Circumferential Fibers"
);
ring_center[0] = this->get_double ("x-coordinate of the center of the ring");
ring_center[1] = this->get_double ("y-coordinate of the center of the ring");
this->leave_subsection();

The following lines help keeping track of what prm file goes

with a specific output. Therefore, they are here for

convenience and not for any specific computational need.

ofstream paramfile((output_name+"_param.prm").c_str());
this->print_parameters(paramfile, this->Text);
}

This class provides the exact distribution of the Lagrange multiplier for enforcing incompressibility both in the fluid and in the immersed domains for a specific problem. The latter concerns the equilibrium of a circular cylinder immersed in an incompressible Newtonian fluid. The immersed domain is assumed to be an incompressible elastic material with an elastic response proportional to the stretch of elastic fibers wound in the circumferential direction (hoop). The constitutive equations of the cylinder correspond to the choice of "Solid constitutive model" as "CircumferentialFiberModel". Finally, we refer to this particular problem as the "Hello world" problem for immersed methods. We learned this expression from our colleague Boyce E. Griffith, currently at the Leon H. Charney Division of Cardiology, Department of Medicine, NYU School of Medicine, New York University.

template<int dim>
class ExactSolutionRingWithFibers :
public Function<dim>
{
public:

No default constructor is defined. Simulation objects must be initialized by assigning the simulation parameters, which are elements of objects of type ProblemParameters.

ExactSolutionRingWithFibers (ProblemParameters<dim> &par);
void vector_value (const Point <dim> &p,
Vector <double> &values) const;
void vector_value_list(const vector< Point<dim> > &points,
vector <Vector <double> > &values ) const;

Inner radius of the ring.

double R;

Width of the ring.

double w;

Edge length of the square control volume.

double l;

Center of the ring.

Point<dim> center;
private:
ProblemParameters<dim> &par;
};

Class constructor.

template<int dim>
ExactSolutionRingWithFibers<dim>::ExactSolutionRingWithFibers (
ProblemParameters<dim> &prm
)
:
Function<dim>(dim+1),
par(prm)
{
par.enter_subsection (
"Equilibrium Solution of Ring with Circumferential Fibers"
);
R = par.get_double ("Inner radius of the ring");
w = par.get_double ("Width of the ring");
l = par.get_double ("Any edge length of the (square) control volume");
center[0] = par.get_double ("x-coordinate of the center of the ring");
center[1] = par.get_double ("y-coordinate of the center of the ring");
par.leave_subsection ();
}

It provides the Lagrange multiplier (pressure) distribution in the control volume and in the ring at equilibrium.

template<int dim>
void
ExactSolutionRingWithFibers<dim>::vector_value
(
const Point <dim> &p,
Vector <double> &values
) const
{
double r = p.distance (center);
values = 0.0;
double p0 = -numbers::PI*par.mu/(2*l*l)*w*(2*R+w);
if (r >= (R+w))
values(dim) = p0;
else if (r >= R)
values(dim) = p0 + par.mu*log((R + w)/r);
else
values(dim) = p0 + par.mu*log(1.0 + w/R);
}

It provides the Lagrange multiplier (pressure) distribution in the control volume and in the ring at equilibrium.

template<int dim>
void
ExactSolutionRingWithFibers<dim>::vector_value_list
(
const vector< Point<dim> > &points,
vector <Vector <double> > &values
) const
{
Assert (points.size() == values.size(),
ExcDimensionMismatch(points.size(), values.size()));
for (unsigned int i = 0; i < values.size(); ++i)
vector_value (points[i], values[i]);
}

It defines simulations objects. The only method in the public interface is run(), which is invoked to carry out the simulation.

template <int dim>
class ImmersedFEM
{
public:

No default constructor is defined. Simulation objects must be initialized by assigning the simulation parameters, which are elements of objects of type ProblemParameters.

ImmersedFEM(ProblemParameters<dim> &par);
~ImmersedFEM();
void run ();
private:

The parameters of the problem.

ProblemParameters<dim> &par;

Vector of boundary indicators. The type of this vector matches the return type of the function Triangulation< dim, spacedim >::get_boundary_indicator().

vector<unsigned char> boundary_indicators;

Triangulation over the control volume (fluid domain). Following deal.II conventions, a triangulation pertains to a manifold of dimension dim embedded in a space of dimension spacedim. In this case, only a single dimensional parameter is specified so that the dimension of the manifold and of the containing space are the same.

Triangulation<dim> tria_f;

Triangulations of the immersed domain (solid domain). Following deal.II conventions, a triangulation pertains to a manifold of dimension dim embedded in a space of dimension spacedim. While in this case the two dimension parameters are set equal to each other, it is possible to formulate problems in which the immersed domain is a manifold of dimension lower than that of the containing space.

Triangulation<dim, dim> tria_s;

FESystem for the control volume. It consists of two fields: velocity (a vector field of dimension dim) and pressure (a scalar field). The meaning of the parameter dim is as for the Triangulation<dim> tria_f element of the class.

FESystem<dim> fe_f;

A variable to check whether the pressure field is approximated using the FE_DGP elements.

bool dgp_for_p;

This is the FESystem for the immersed domain.

FESystem<dim, dim> fe_s;

The dof_handler for the control volume.

DoFHandler<dim> dh_f;

The dof_handler for the immersed domain.

DoFHandler<dim, dim> dh_s;

The triangulation of for the immersed domain defines the reference configuration of the immersed domain. As the immersed domain moves through the fluid, it is important to be able to conveniently describe quantities defined over the immersed domain according to an Eulerian view. It is therefore convenient to define a MappingQEulerian object that will support such a description.

MappingQEulerian<dim, Vector<double>, dim> * mapping;

The quadrature object for the control volume.

QGauss<dim> quad_f;

The quadrature object for the immersed domain.

QTrapez<1> qtrapez;
QIterated<dim> quad_s;

Constraints matrix for the control volume.

ConstraintMatrix constraints_f;

Constraints matrix for the immersed domain.

ConstraintMatrix constraints_s;

Sparsity pattern.

BlockSparsityPattern sparsity;

Jacobian of the residual.

BlockSparseMatrix<double> JF;

Object of BlockSparseMatrix<double> type to be used in place of the real Jacobian when the real Jacobian is not to be modified.

BlockSparseMatrix<double> dummy_JF;

State of the system at current time step: velocity, pressure, and displacement of the immersed domain.

BlockVector<double> current_xi;

State of the system at previous time step: velocity, pressure, and displacement of the immersed domain.

BlockVector<double> previous_xi;

Approximation of the time derivative of the state of the system.

BlockVector<double> current_xit;

Current value of the residual.

BlockVector<double> current_res;

Newton iteration update.

BlockVector<double> newton_update;

Vector to compute the average pressure when the average pressure is set to zero.

Vector<double> pressure_average;

Vector to represent a uniform unit pressure.

Vector<double> unit_pressure;

Number of degrees of freedom for each component of the system.

unsigned int n_dofs_u, n_dofs_p, n_dofs_up, n_dofs_W, n_total_dofs;

A couple of vectors that can be used as temporary storage. They are defined as a private member of the class to avoid that the object is allocated and deallocated when used, so to gain in efficiency.

Vector<double> tmp_vec_n_total_dofs;
Vector<double> tmp_vec_n_dofs_up;

Matrix to be inverted when solving the problem.

SparseDirectUMFPACK JF_inv;

Scalar used for conditioning purposes.

double scaling;

Variable to keep track of the previous time.

double previous_time;

The first dof of the pressure field.

unsigned int constraining_dof;

A container to store the dofs corresponding to the pressure field.

set<unsigned int> pressure_dofs;

Storage for the elasticity operator of the immersed domain.

Vector <double> A_gamma;

Mass matrix of the immersed domain.

SparseMatrix<double> M_gamma3;

Inverse of M_gamma3.

SparseDirectUMFPACK M_gamma3_inv;

M_gamma3_inv * A_gamma.

Vector <double> M_gamma3_inv_A_gamma;

Area of the control volume.

double area;

File stream that is used to output a file containing information about the fluid flux, area and the centroid of the immersed domain over time.

ofstream global_info_file;

Function declarations

void create_triangulation_and_dofs ();
void apply_constraints (vector<double> &local_res,
FullMatrix<double> &local_jacobian,
const Vector<double> &local_up,
const vector<unsigned int> &dofs);
void compute_current_bc (const double time);
void apply_current_bc (
BlockVector<double> &vec,
const double time);
void assemble_sparsity (Mapping<dim, dim> &mapping);
void get_area_and_first_pressure_dof ();
void residual_and_or_Jacobian (
BlockVector<double> &residual,
BlockSparseMatrix<double> &Jacobian,
const BlockVector<double> &xit,
const BlockVector<double> &xi,
const double alpha,
const double t
);
void distribute_residual (
Vector<double> &residual,
const vector<double> &local_res,
const vector<unsigned int> &dofs_1,
const unsigned int offset_1
);
void distribute_jacobian (
SparseMatrix<double> &Jacobian,
const FullMatrix<double> &local_Jac,
const vector<unsigned int> &dofs_1,
const vector<unsigned int> &dofs_2,
const unsigned int offset_1,
const unsigned int offset_2
);
void distribute_constraint_on_pressure (
Vector<double> &residual,
const double average_pressure
);
void distribute_constraint_on_pressure (
SparseMatrix<double> &jacobian,
const vector<double> &pressure_coefficient,
const vector<unsigned int> &dofs,
const unsigned int offset
);
void localize (
Vector<double> &local_M_gamma3_inv_A_gamma,
const Vector<double> &M_gamma3_inv_A_gamma,
const vector<unsigned int> &dofs
);
void get_Agamma_values (
const FEValues<dim,dim> &fe_v_s,
const vector< unsigned int > &dofs,
const Vector<double> &xi,
Vector<double> &local_A_gamma
);
void get_Pe_F_and_DPeFT_dxi_values (
const FEValues<dim,dim> &fe_v_s,
const vector< unsigned int > &dofs,
const Vector<double> &xi,
const bool update_jacobian,
vector<Tensor<2,dim,double> > &Pe,
vector<Tensor<2,dim,double> > &F,
vector< vector<Tensor<2,dim,double> > > & DPe_dxi
);
void calculate_error () const;
unsigned int n_dofs() const {
return n_total_dofs;
};
void output_step (
const double t,
const BlockVector<double> &solution,
const unsigned int step_number,
const double h
);
template<class Type>
inline void set_to_zero (Type &v) const;
template<class Type>
inline void set_to_zero (Table<2,Type> &v) const;
template<class Type>
inline void set_to_zero(vector<Type> &v) const;
double norm(const vector<double> &v);
};

Constructor: Initializes the FEM system of the control volume; Initializes the FEM system of the immersed domain; Initializes, corresponding dof handlers, and the quadrature rule; It runs the create_triangulation_and_dofs function.

template <int dim>
ImmersedFEM<dim>::ImmersedFEM (ProblemParameters<dim> &par)
:
par (par),
fe_f (
FE_Q<dim>(par.degree),
dim,
*FETools::get_fe_from_name<dim>(par.fe_p_name),
par.degree-1
),
fe_s (FE_Q<dim, dim>(par.degree), dim),
dh_f (tria_f),
dh_s (tria_s),
quad_f (par.degree+2),
quad_s (qtrapez, 4*(par.degree+8))
{
if(par.degree <= 1)
cout
<< " WARNING: The chosen pair of finite element spaces is not stable."
<< endl
<< " The obtained results will be nonsense."
<< endl;
if( Utilities::match_at_string_start(par.fe_p_name, string("FE_DGP")))
dgp_for_p = true;
else dgp_for_p = false;
create_triangulation_and_dofs ();
global_info_file.open((par.output_name+"_global.gpl").c_str());
}

Distructor: deletion of pointers created with new and closing of the record keeping file.

template <int dim>
ImmersedFEM<dim>::~ImmersedFEM ()
{
delete mapping;
global_info_file.close();
}

Determination of the current value of time dependent boundary values.

template <int dim>
void
ImmersedFEM<dim>::compute_current_bc (const double t)
{
par.u_g.set_time(t);
VectorTools::interpolate_boundary_values (
StaticMappingQ1<dim>::mapping,
dh_f,
par.boundary_map,
par.boundary_values,
par.component_mask
);

Set to zero the value of the first dof associated to the pressure field.

if(par.fix_pressure == true) par.boundary_values[constraining_dof] = 0;
}

Application of time dependent boundary conditions.

template <int dim>
void
ImmersedFEM<dim>::apply_current_bc
(
BlockVector<double> &vec,
const double t
)
{
compute_current_bc(t);
map<unsigned int, double>::iterator it = par.boundary_values.begin(),
itend = par.boundary_values.end();
if(vec.size() != 0)
for(; it != itend; ++it)
vec.block(0)(it->first) = it->second;
else
for(; it != itend; ++it)
constraints_f.set_inhomogeneity(it->first, it->second);
}

Defines the triangulations for both the control volume and the immersed domain. It distributes degrees of freedom over said triangulations. Both grids are assumed to be available in UCD format. The naming convention is as follows: fluid_[dim]d.inp for the control volume and solid_[dim]d.inp for the immersed domain. This function also sets up the constraint matrices for the enforcement of Dirichlet boundary conditions. In addition, it sets up the framework for enforcing the initial conditions.

template <int dim>
void
ImmersedFEM<dim>::create_triangulation_and_dofs ()
{
if(par.material_model == ProblemParameters<dim>::CircumferentialFiberModel)
{

This is used only by the solution of the problem with the immersed domain consisting of a circular cylinder. We only implemented this in two dimensions.

Assert(dim == 2, ExcNotImplemented());
ExactSolutionRingWithFibers<dim> ring(par);

Construct the square domain for the control volume using the parameter file.

GridGenerator::hyper_cube (tria_f, 0., ring.l);

Construct the hyper shell using the parameter file.

GridGenerator::hyper_shell(tria_s, ring.center,
ring.R, ring.R+ring.w);
static const HyperShellBoundary<dim> shell_boundary(ring.center);
tria_s.set_boundary(0, shell_boundary);
}
else
{

As specified in the documentation for the "GridIn" class the triangulation corresponding to a grid needs to be empty at this time.

GridIn<dim> grid_in_f;
grid_in_f.attach_triangulation (tria_f);
{
ifstream file (par.fluid_mesh.c_str());
Assert (file, ExcFileNotOpen (par.fluid_mesh.c_str()));

A grid in ucd format is expected.

grid_in_f.read_ucd (file);
}
GridIn<dim, dim> grid_in_s;
grid_in_s.attach_triangulation (tria_s);
ifstream file (par.solid_mesh.c_str());
Assert (file, ExcFileNotOpen (par.solid_mesh.c_str()));

A grid in ucd format is expected.

grid_in_s.read_ucd (file);
}
cout
<< "Number of fluid refines = "
<< par.ref_f
<< endl;
tria_f.refine_global (par.ref_f);
cout
<< "Number of active fluid cells: "
<< tria_f.n_active_cells ()
<< endl;
cout
<< "Number of solid refines = "
<< par.ref_s
<< endl;
tria_s.refine_global (par.ref_s);
cout
<< "Number of active solid cells: "
<< tria_s.n_active_cells ()
<< endl;

Initialization of the boundary_indicators vector.

boundary_indicators = tria_f.get_boundary_indicators ();

Distribution of the degrees of freedom. Both for the solid and fluid domains, the dofs are renumbered first globally and then by component.

dh_f.distribute_dofs (fe_f);
DoFRenumbering::boost::Cuthill_McKee (dh_f);

Consistently with the fact that the various components of the system are stored in a block matrix, now renumber velocity and pressure component wise.

vector<unsigned int> block_component (dim+1,0);
block_component[dim] = 1;
DoFRenumbering::component_wise (dh_f, block_component);
vector<unsigned int> dofs_per_block (2);
DoFTools::count_dofs_per_block (dh_f, dofs_per_block, block_component);

Accounting of the number of degrees of freedom for the fluid domain on a block by block basis.

n_dofs_u = dofs_per_block[0];
n_dofs_p = dofs_per_block[1];
n_dofs_up = dh_f.n_dofs ();

Simply distribute dofs on the solid displacement.

dh_s.distribute_dofs (fe_s);
DoFRenumbering::boost::Cuthill_McKee (dh_s);

Determine the total number of dofs.

n_dofs_W = dh_s.n_dofs ();
n_total_dofs = n_dofs_up+n_dofs_W;
cout
<< "dim (V_h) = "
<< n_dofs_u
<< endl
<< "dim (Q_h) = "
<< n_dofs_p
<< endl
<< "dim (Z_h) = "
<< dh_s.n_dofs ()
<< endl
<< "Total: "
<< n_total_dofs
<< endl;
vector<unsigned int> all_dofs (2);
all_dofs[0] = n_dofs_up;
all_dofs[1] = n_dofs_W;

Re-initialization of the BlockVectors containing the values of the degrees of freedom and of the residual.

current_xi.reinit (all_dofs);
previous_xi.reinit (all_dofs);
current_xit.reinit (all_dofs);
current_res.reinit (all_dofs);
newton_update.reinit (all_dofs);

Re-initialization of the average and unit pressure vectors.

pressure_average.reinit (n_dofs_up);
unit_pressure.reinit (n_dofs_up);

Re-initialization of temporary vectors.

tmp_vec_n_total_dofs.reinit(n_total_dofs);
tmp_vec_n_dofs_up.reinit(n_dofs_up);

We now deal with contraint matrices.

{
constraints_f.clear ();
constraints_s.clear ();

Enforce hanging node constraints.

DoFTools::make_hanging_node_constraints (dh_f, constraints_f);
DoFTools::make_hanging_node_constraints (dh_s, constraints_s);

To solve the problem we first assemble the Jacobian of the residual using zero boundary values for the velocity. The specification of the actual boundary values is done later by the apply_current_bc function.

VectorTools::interpolate_boundary_values (
StaticMappingQ1<dim>::mapping,
dh_f,
par.zero_boundary_map,
constraints_f,
par.component_mask);
}

Determine the area (in 2D) of the control volume and find the first dof pertaining to the pressure.

get_area_and_first_pressure_dof ();
constraints_f.close ();
constraints_s.close ();

The following matrix plays no part in the formulation. It is defined here only to use the VectorTools::project function in initializing the vectors previous_xi.block(0) and unit_pressure.

ConstraintMatrix cc;
cc.close();

Construction of the initial conditions.

if(fe_f.has_support_points())
{
VectorTools::interpolate (dh_f, par.u_0, previous_xi.block(0));
VectorTools::interpolate (
dh_f,
ComponentSelectFunction<dim>(dim, 1., dim+1),
unit_pressure
);
}
else
{
VectorTools::project (dh_f, cc, quad_f, par.u_0, previous_xi.block(0));
VectorTools::project (
dh_f,
cc,
quad_f,
ComponentSelectFunction<dim>(dim, 1., dim+1),
unit_pressure
);
}
if(fe_s.has_support_points())
VectorTools::interpolate (dh_s, par.W_0, previous_xi.block(1));
else
VectorTools::project (dh_s, cc, quad_s, par.W_0, previous_xi.block(1));
mapping = new MappingQEulerian<dim, Vector<double>, dim> (par.degree,
previous_xi.block(1), dh_s);

We now deal with the sparsity patterns.

{
BlockCompressedSimpleSparsityPattern csp (2,2);
csp.block(0,0).reinit (n_dofs_up, n_dofs_up);
csp.block(0,1).reinit (n_dofs_up, n_dofs_W );
csp.block(1,0).reinit (n_dofs_W , n_dofs_up);
csp.block(1,1).reinit (n_dofs_W , n_dofs_W );

As stated in the documentation, now we must call the function csp.collect_sizes.() since have changed the size of the sub-objects of the object csp.

csp.collect_sizes();
Table< 2, DoFTools::Coupling > coupling(dim+1,dim+1);
for(unsigned int i=0; i<dim; ++i)
{

Velocity is coupled with pressure.

coupling(i,dim) = DoFTools::always;

Pressure is coupled with velocity.

coupling(dim,i) = DoFTools::always;
for(unsigned int j=0; j<dim; ++j)

The velocity components are coupled with themselves and each other.

coupling(i,j) = DoFTools::always;
}

The pressure is coupled with itself.

coupling(dim, dim) = DoFTools::always;

Find the first pressure dof. Then tell all the pressure dofs that they are related to the first pressure dof.

set<unsigned int>::iterator it = pressure_dofs.begin();
constraining_dof = *it;
for(++it; it != pressure_dofs.end(); ++it)
{
csp.block(0,0).add(constraining_dof, *it);
}
DoFTools::make_sparsity_pattern (dh_f,
coupling,
csp.block(0,0),
constraints_f,
true);
DoFTools::make_sparsity_pattern (dh_s, csp.block(1,1));
sparsity.copy_from (csp);
assemble_sparsity(*mapping);
}

Here is the Jacobian matrix.

JF.reinit(sparsity);

Boundary conditions at t = 0.

apply_current_bc(previous_xi, 0);

Resizing other containers concerning the elastic response of the immersed domain.

A_gamma.reinit(n_dofs_W);
M_gamma3_inv_A_gamma.reinit(n_dofs_W);

Creating the mass matrix for the solid domain and storing its inverse.

ConstantFunction<dim> phi_b_func (par.Phi_B, dim);
M_gamma3.reinit (sparsity.block(1,1));

Using the deal.II in-built functionality to create the mass matrix.

MatrixCreator::create_mass_matrix (dh_s, quad_s, M_gamma3, &phi_b_func);
M_gamma3_inv.initialize (M_gamma3);
}

Relatively standard way to determine the sparsity pattern of each block of the global Jacobian.

template <int dim>
void
ImmersedFEM<dim>::assemble_sparsity (Mapping<dim, dim> &immersed_mapping)
{
FEFieldFunction<dim, DoFHandler<dim>, Vector<double> > up_field (dh_f, tmp_vec_n_dofs_up);
vector< typename DoFHandler<dim>::active_cell_iterator > cells;
vector< vector< Point< dim > > > qpoints;
vector< vector< unsigned int> > maps;
vector< unsigned int > dofs_f(fe_f.dofs_per_cell);
vector< unsigned int > dofs_s(fe_s.dofs_per_cell);
typename DoFHandler<dim,dim>::active_cell_iterator
cell = dh_s.begin_active(),
endc = dh_s.end();
FEValues<dim,dim> fe_v(immersed_mapping, fe_s, quad_s,
update_quadrature_points);
CompressedSimpleSparsityPattern sp1(n_dofs_up, n_dofs_W);
CompressedSimpleSparsityPattern sp2(n_dofs_W , n_dofs_up);
for(; cell != endc; ++cell)
{
fe_v.reinit(cell);
cell->get_dof_indices(dofs_s);
up_field.compute_point_locations (fe_v.get_quadrature_points(),
cells, qpoints, maps);
for(unsigned int c=0; c<cells.size(); ++c)
{
cells[c]->get_dof_indices(dofs_f);
for(unsigned int i=0; i<dofs_f.size(); ++i)
for(unsigned int j=0; j<dofs_s.size(); ++j)
{
sp1.add(dofs_f[i],dofs_s[j]);
sp2.add(dofs_s[j],dofs_f[i]);
}
}
}
sparsity.block(0,1).copy_from(sp1);
sparsity.block(1,0).copy_from(sp2);
}

Determination of the volume (area in 2D) of the control volume and identification of the first dof associated with the pressure field.

template <int dim>
void
ImmersedFEM<dim>::get_area_and_first_pressure_dof ()
{
area = 0.0;
typename DoFHandler<dim,dim>::active_cell_iterator
cell = dh_f.begin_active (),
endc = dh_f.end ();
FEValues<dim,dim> fe_v (fe_f,
quad_f,
update_values |
update_JxW_values);
vector<unsigned int> dofs_f(fe_f.dofs_per_cell);

Calculate the area of the control volume.

for(; cell != endc; ++cell)
{
fe_v.reinit (cell);
cell->get_dof_indices (dofs_f);
for(unsigned int i=0; i < fe_f.dofs_per_cell; ++i)
{
unsigned int comp_i = fe_f.system_to_component_index(i).first;
if(comp_i == dim)
{
pressure_dofs.insert(dofs_f[i]);
if (dgp_for_p) break;
}
}
for(unsigned int q=0; q<quad_f.size(); ++q) area += fe_v.JxW(q);
}

Get the first dof pertaining to pressure.

constraining_dof = *(pressure_dofs.begin());
}

Assemblage of the various operators in the formulation along with their contribution to the system Jacobian.

template <int dim>
void
ImmersedFEM<dim>::residual_and_or_Jacobian
(
BlockVector<double> &residual,
BlockSparseMatrix<double> &jacobian,
const BlockVector<double> &xit,
const BlockVector<double> &xi,
const double alpha,
const double t
)
{

Determine whether or not the calculation of the Jacobian is needed.

bool update_jacobian = !jacobian.empty();

Reset the mapping to NULL.

if(mapping != NULL) delete mapping;

In a semi-implicit scheme, the position of the immersed body coincides with the position of the body at the previous time step.

if(par.semi_implicit == true)
{
if(fabs(previous_time - t) > 1e-12)
{
previous_time = t;
previous_xi = xi;
}
mapping = new MappingQEulerian<dim, Vector<double>, dim> (par.degree,
previous_xi.block(1),
dh_s);
}
else
mapping = new MappingQEulerian<dim, Vector<double>, dim> (par.degree,
xi.block(1),
dh_s);

In applying the boundary conditions, we set a scaling factor equal to the diameter of the smallest cell in the triangulation.

scaling = GridTools::minimal_cell_diameter(tria_f);

Initialization of the residual.

residual = 0;

If the Jacobian is needed, then it is initialized here.

if(update_jacobian)
{
jacobian.clear();
assemble_sparsity(*mapping);
jacobian.reinit(sparsity);
}

Evaluation of the current values of the external force and of the boundary conditions.

par.force.set_time(t);
compute_current_bc(t);

Computation of the maximum number of degrees of freedom one could have on a fluid-solid interaction cell. Rationale the coupling of the fluid and solid domains is computed by finding each of the fluid cells that interact with a given solid cell. In each interaction instance we will be dealing with a total number of degrees of freedom that is the sum of the dofs of the current solid cell and the dofs of the current fluid cell in the list of fluid cells interacting with the solid cell in question.

unsigned int n_local_dofs = fe_f.dofs_per_cell + fe_s.dofs_per_cell;

Storage for the local dofs in the fluid and in the solid.

vector< unsigned int > dofs_f(fe_f.dofs_per_cell);
vector< unsigned int > dofs_s(fe_s.dofs_per_cell);

FEValues for the fluid.

FEValues<dim> fe_f_v (fe_f,
quad_f,
update_values |
update_gradients |
update_JxW_values |
update_quadrature_points);

Number of quadrature points on fluid and solid cells.

const unsigned int nqpf = quad_f.size();
const unsigned int nqps = quad_s.size();

The local residual vector: the largest possible size of this vector is n_local_dofs.

vector<double> local_res(n_local_dofs);
vector<Vector<double> > local_force(nqpf, Vector<double>(dim+1));
FullMatrix<double> local_jacobian;
if(update_jacobian) local_jacobian.reinit(n_local_dofs, n_local_dofs);

Since we want to solve a system of equations of the form $f(\xi', \xi, t) = 0$, we need to manage the information in $\xi'$ as though it were independent of the information in $\xi$. We do so by defining a vector of local degrees of freedom that has a length equal to twice the total number of local degrees of freedom. This information is stored in the vector local_x.