LifeV_SOTA.tex 52.4 KB
 grandper committed Jul 01, 2010 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 %============================ % LifeV - State of the Art % July 2010 % % Author: Gwenol Grandperrin % % Contributors: % Paolo Crosetto % Simone Deparis % Gilles Fourestey % Cristiano Malossi % Samuel Quinodoz % Ricardo Ruiz % Paolo Tricerri %============================ \documentclass[11pt]{article} \usepackage{amssymb} \usepackage{amsmath} \usepackage{float} \usepackage{pifont} \usepackage{graphicx} \usepackage{subfigure} \usepackage{algorithme} \usepackage{longtable} %Pour aligner verticalement les tableaux \usepackage{array} \newcolumntype{M}[1]{>{\centering}m{#1}} \newcommand{\tick}{\ding{51}} \newcommand{\norm}[1]{\left\| #1 \right\|} \begin{document} \title{LifeV - State of the Art} \author{Gwenol Grandperrin, Paolo Crosetto, Simone Deparis,\\ Gilles Fourestey, Cristiano Malossi, Samuel Quinodoz, \\ Ricardo Ruiz, Paolo Tricerri} \date{\today} \maketitle  malossi committed Jul 20, 2010 43 \tableofcontents  grandper committed Jul 01, 2010 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83  %============================ % Introduction %============================ \section{Introduction} %LifeV is a finite element (FE) library. Nowadays super computers had become very efficient. However solving huge systems still requires research and development of efficient tools. One can think for example of the simulation of the entire cardiovascular system. To face such kind of challenges, the CMCS at EPFL in Switzerland, the MOX at Politechnico di Milano in Italy and REO group at INRIA in France have developed a C++ 3D finite element code called LifeV (prononced "life-five"). In 2006, CMCS started to parallelize the library. In this document we describe the actual level of both knowledge and development of LifeV. The content is actually made of contributions of the developers themselves. In Section~\ref{sec:lifev} we will describe the core of the library, i.e., the solvers, the finite elements, the boundary conditions and the tools available. Then in Section~\ref{sec:applications}, we present the current applications done with LiveV. Finally in the Section~\ref{sec:otherProjects}, one can find other projects where the LiveV library was one of the key aspect. %============================ % LifeV, a Finite Elements library %============================ \section{LifeV, a C++ finite element library} \label{sec:lifev} As mentioned in the introduction, LifeV is a finite element library which has been implemented in C++. The development started in 2002 and the code contains today around 50.000 lines. The code is free and available on http://www.lifev.org under LGPL license. In 2010, a new release of the code has been proposed where the whole code has been ported to a parallel version to achieve high performance computing on multiprocessors architectures. The largest number of degrees of freedom computed in a simulation with the new version is nine millions. All the efforts are now made to prepare LifeV for the architectures of tomorrow's supercomputers. % Note: the nine millions of DOF has been done with an ethiersteinman problem by Gilles Fourestey In Section~\ref{sec:solvers}, we describe first the physical solvers implemented in LifeV. Then we present the finite elements, the boundary conditions and the input/output management in Sections~\ref{sec:BC}, \ref{sec:FE} and \ref{sec:inOutput} respectively. % Solvers available in LifeV \subsection{Physical solvers available in LifeV} \label{sec:solvers} We refer to a physical solver class as a class that is able to build the discrete system associated to the approximation of a PDE. In this section, we go through the existing physical solvers in LifeV. Hence the reader will have a clear idea of the possible equations which can be solved. %\subsubsection{Harmonic Extension solver} \subsubsection{HarmonicExtensionSolver class} This solver deals with the Laplace equation: \begin{eqnarray} \Delta\eta & = & 0\\ +\text{b.c.}\nonumber \end{eqnarray} However more features have been implemented to use it as a tool in the Arbitrary Eulerian Lagrangian (ALE) framework (more informations can be found in~\cite{nobile} or~\cite{quarteroni}). The harmonic extension technique allows for example to move the points of a mesh in the latter framework.  grandper committed Jul 02, 2010 84   grandper committed Jul 01, 2010 85 The following table gives an overview of the main methods available in the HarmonicExtensionSolver class:  grandper committed Jul 02, 2010 86   grandper committed Jul 01, 2010 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 %\begin{table}[H] \begin{center} \begin{longtable}{p{3cm}p{7cm}} \hline \texttt{setUp(dataFile)} & Sets up the system solver with a GetPot datafile\footnote{GetPot is a parsing configuration files library available at http://getpot.sourceforge.net}\\ \texttt{initialize(disp)} & Initialization of the problem\\ \texttt{buildSystem()} & Builds the system\\ \texttt{updateSystem()} & Stores the displacement value in the old displacement variable\\ \texttt{iterate(BCHandler)} \\ & Updates the convective term, does the boundary conditions treatment and solves the linear system.\\ \texttt{dispOld()} & Returns the old value of the displacement\\ \texttt{dispDiff()} & Displacement difference (between the old and the new displacement)\\ \texttt{disp()} & Returns the displacement\\ \texttt{setDisplacement(disp)} \\ & Sets the value of the displacement with a vector\\ \texttt{setDispOld(disp)} & Sets the value of the old displacement\\ \texttt{getMap()} & Getter for the map\\ \texttt{comm()} & Returns the MPI communicator\\ \texttt{isLeader()} & Returns true on the Rank 0 process\\ \texttt{resetPrec(bool)} & Resets the preconditioner of the linear solver\\ \texttt{rescaleMatrix(dt)} & Multiplies the matrix by the scalar given as a parameter\\ \texttt{applyBoundaryConditions(rhs,BCHandler)} \\ & Applies the boundary conditions\\ \texttt{computeMatrix()} & Computes the matrix of the system\\ \texttt{updateDispDiff()} & Recomputes the value of the displacement difference\\ \hline \end{longtable} \end{center} %\caption{Main methods of the HarmonicExtensionSolver class} %\label{tab:harmonicExtensionMethods} %\end{table} %\subsubsection{Oseen solver} \subsubsection{Oseen class} We consider a flow of an incompressible and viscous fluid decribed by its velocity $\mathbf{u}$ and its pressure $p$. We assume that the reynolds number is low. This solver is designed to get the solution of the so-called Oseen problem: \begin{eqnarray} \alpha\mathbf{u}+\boldsymbol{\beta}\cdot\nabla\mathbf{u}-\nu\Delta\mathbf{u} + \nabla p & = & f\\ \nabla\cdot\mathbf{u} & = & 0\\ +\text{b.c.}\nonumber \end{eqnarray} where $\nu=\frac{\mu}{\rho}$ is the kinematic viscosity of the fluid, $\mu$ the viscosity, $\rho$ the density, and $\alpha$ and $\beta$ are two parameters. The main methods of the class are given in the following table: %\begin{table}[H] \begin{center} \begin{longtable}{p{3cm}p{8cm}} \hline \texttt{iterate(BCHandler)} \\ & Updates the convective term, does the boundary conditions treatment and solves the linear system.\\ \texttt{updateSystem(alpha,beta,rhs,matrix)} \\ & Updates the values of $\alpha$, $\boldsymbol{\beta}$, the rhs and the matrix \\ \texttt{solution()} & Returns the local solution vector\\ \texttt{residual()} & Returns the local residual vector\\ \texttt{reduceSolution(u,p)} \\ & Reduce the local solution into the global vectors for $\mathbf{u}$ and p\\ \texttt{setSourceTerm(s)} & Sets the source term functor\\ \texttt{sourceTerm()} & Returns the source term functor\\ \texttt{post\_proc()} & Returns the Post Processing stuff.\\ \texttt{area(flag)} & Computes the area on a given part of the boundary\\ \texttt{flux(flag,vector)} & Computes the flux on a given part of the boundary\\ \texttt{pressure(flag,vector)} \\ & Computes the pressure on a given part of the boundary\\ \texttt{LagrangeMultiplier(flag,BCHandler\_raw\_type)} \\ & Gets the Lagrange multiplier related to a flux imposed on a given part of the boundary\\ \texttt{flux(flag)} & Computes the flux on a given part of the boundary\\ \texttt{pressure(flag)} & Computes the pressure on a given part of the boundary.\\ \texttt{LagrangeMultiplier(flag)} \\ & Gets the Lagrange multiplier related to a flux imposed on a given part of the boundary\\ \texttt{postProcess(bool)} & Postprocessing\\ \hline \end{longtable} \end{center} %\caption{Main methods of the Oseen class} %\label{tab:oseenMethods} %\end{table} %\subsubsection{Navier-Stokes solver} \subsubsection{NavierStokesSolver class} The classical Navier-Stokes equations for an incompressible viscous flow reads: \begin{eqnarray} \frac{d}{dt}\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{u}-\nu\Delta\mathbf{u} + \nabla p & = & f \label{eqn:NS1}\\ \nabla\cdot\mathbf{u} & = & 0\label{eqn:NS2}\\ +\text{b.c.}\nonumber \end{eqnarray} where $\mathbf{u}$ is the fluid velocity, $p$ the pressure, $\nu=\frac{\mu}{\rho}$ the kinematic viscosity of the fluid, $\mu$ the viscosity, $\rho$ the density and $f$ the external forces. %The full implicit scheme is not yet implemented. The table below gives an overview of the main methods of the class. %\begin{table}[H] \begin{center} \begin{longtable}{p{3cm}p{8cm}} \hline \texttt{timeAdvance(source,time)} \\ & Updates the right hand side for time advancing\\ \texttt{iterate(time)} & Updates convective term, boundary conditions treatment and solve the linearized Navier-Stokes system\\ \texttt{eval(fx0,gx0,x0,status)} \\ & Sets the components of \texttt{fx0} relative to the degrees of freedom of $\mathbf{u}$ equal to the corresponding values of $\mathbf{u}$\\ \texttt{initialize(x0,t0,dt)} \\ & Initializes the values of the solution at the initial time step\\ \texttt{initialize(string)} \\ & Initializes the values of the solution at the initial time step with an external file\\ \texttt{initializeNull()} & Initializes the values of the solution at the initial time step with a previously computed solution\\ \texttt{initializeStokes(source,t0)} \\ & Initializes with steady Stokes solution (without convection)\\ \texttt{setUp(dataFile)} & Sets up the system solver with a GetPot datafile\\ \texttt{buildSystem()} & Builds the matrix of the linear system\\ \texttt{updateSystem()} & Updates the matrix of the linear system\\ \texttt{u()} & Returns the solution for $\mathbf{u}$\\ \texttt{p()} & Returns the solution for $p$\\ \texttt{setFluidBC()} & Returns true if the boundary conditions are set\\ \texttt{setFluidBC(BCHandler)} \\ & Sets the fluid BCHandler\\ \texttt{bcHandler()} & Returns the bcHandler used for the fluid boundary conditions\\ \texttt{setSourceTerm(s)} & Sets the source term functor\\ \texttt{sourceTerm()} & Gets the source term functor\\ \texttt{postProcess()} & Postprocessing\\ \texttt{dim\_u()} & Returns the number of degrees of freedom number for $\mathbf{u}$\\ \texttt{dim\_p()} & Returns the number of degrees of freedom number for $p$\\ \hline \end{longtable} \end{center} %\caption{Main methods of the NavierStokesSolver class} %\label{tab:navierStokesMethods} %\end{table} %\subsubsection{ALE Free Surface solver} % Gwenol Grandperrin, Samuel Quinodoz %To do %\subsubsection{Saint Venant-Kirchhof solver} \subsubsection{VenantKirchhofSolver class} \label{sec:venantKirchhofSolver} % Paolo Tricerri The Saint Venant-Krichhof model for elastic materials reads \begin{eqnarray} \rho\frac{\partial^2\boldsymbol{d}_s}{\partial t^2} - \nabla\cdot\sigma_s^o(\boldsymbol{d}_s) & = & \rho_0\mathbf{f} \label{eqn:VenantKirchhof}\\ +\text{b.c.}\nonumber \end{eqnarray} where $\sigma_s^o$ is the first Piola-Kirchhoff stress tensor. \sigma_s^o=\lambda tr(\epsilon)+2\mu_s\epsilon where $\lambda$ is the young modulus, $\mu_s$ the poisson coefficiant and \epsilon=\frac{\nabla\boldsymbol{d}_s+(\nabla\boldsymbol{d}_s)^T}{2}. The main methods of the class are the followings: %\begin{table}[H] \begin{center} \begin{longtable}{p{3cm}p{8cm}} \hline \texttt{updateSystem()} & Updates the right hand side for time advancing\\ \texttt{buildSystem(matrix)} \\ & Builds the linear system\\ \texttt{buildSystem()} & Builds the linear system\\ \texttt{iterate(sol)} & Solve the non-linear system\\ \texttt{iterate(BCHandler)} \\ & Updates the convective term, does the boundary conditions treatment and solves the linear system.\\ \texttt{iterateLin(BCHandler)} & \\ \texttt{getMap()} & \\ \texttt{getDisplayer()} & \\ \texttt{getMassStiff()} & \\ \texttt{dFESpace()} & \\ \texttt{BChandler()} & Returns the BCHandler\\ \texttt{residual()} & \\ \texttt{setBC(BCHandler)} & \\ \texttt{setSourceTerm(source)} & \\ \texttt{resetPrec(bool)} & \\ \texttt{setDisp(displacement)} & \\ \texttt{setRecur(recu)} & \\ \texttt{updateJacobian(sol,iter)} & \\ \texttt{solveJac(step,res,linear\_rel\_tol)} \\ & Solves the tangent problem for newton iterations\\ \texttt{solveJac(step,res,linear\_rel\_tol,BCHandler)} \\ & Solves the tangent problem with custom BC \\ \texttt{solveJacobian(Real)} & \\ \texttt{solveJacobian(Real,BCHandler)} & \\ \texttt{evalResidual(res,sol,iter)} \\ & Evaluates residual for newton iterations \\ \texttt{evalConstraintTensor()} & \\ \texttt{sourceTerm()} & \\ \texttt{disp()} & Returns the local displacement vector\\ \texttt{vel()} & Returns the local velocity vector\\ \texttt{rhsWithoutBC()} & Returns the local rhs vector without boundary conditions\\ \texttt{setUp(dataFile)} & Sets up the system solver with a GetPot datafile\\ \texttt{initialize(d0,w0)} & Initializes displacement and velocity with functions \\ \texttt{initialize(d0,w0)} & Initializes displacement and velocity with vectors \\ \texttt{initializeVel(w0)} & \\ \texttt{reduceSolution(disp,vel)} \\ & Reduces the local solution into the global vectors for $\mathbf{u}$ and $\boldsymbol{\eta}$\\ \texttt{comm()} & Returns the communicator\\ \texttt{rescaleMatrices()} & \\ \texttt{rescaleFactor()} & \\ \texttt{updateVel()} & \\ \texttt{offset()} & Returns the offset\\ \texttt{thickness()} & Returns the thickness\\ \texttt{density()} & Returns the density $\rho$\\ \texttt{rho()} & Returns $\rho$ \\ \texttt{young()} & Returns the young modulus $\lambda$ (stiffness of the material)\\ \texttt{poisson()} & Returns the poisson coefficiant $\mu_s$ (material compressibility)\\ \hline \end{longtable} \end{center} %\caption{Main methods of the VenantKirchhofSolver class} %\label{tab:venantKirchhofMethods} %\end{table} %\subsubsection{Monodomain solver} \subsubsection{MonodomainSolver class} % Ricardo Ruiz We consider the following equations: \begin{eqnarray} \frac{\partial v}{\partial t} - \nabla\cdot(M\nabla v) & = & f(v,w), \label{eqn:MonoDomain1}\\ \frac{\partial w}{\partial t} & = & g(v,w), \label{eqn:MonoDomain2}\\ + \text{Neumann b.c.}\nonumber \end{eqnarray} where $v$ is the electrical potential, $w$ is the gating variable, $M$ a given tensor and $f$, $g$ two given functional. This solver has been developped to simulate electricity propagation through the heart. In particular, this solver takes into account the orientation of the muscles fibers.  grandper committed Jul 02, 2010 308 The name of this solver has been chosen to contrast with the Bidomain solver where the equations are applied to \emph{two} superposed domains.  grandper committed Jul 01, 2010 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417  The main methods of the class are the followings: %\begin{table}[H] \begin{center} \begin{longtable}{p{3cm}p{8cm}} \hline \texttt{PDEiterate(BCHandler)} \\ & Updates sources, bc treatment and solve the monodomain system\\ \texttt{setUp(dataFile)} & Sets up the system solver with a GetPot datafile\\ \texttt{buildSystem()} & Builds time independent parts of PDE system\\ \texttt{updatePDESystem(alpha,sourceVec)} \\ & Updates time dependent parts of PDE system\\ \texttt{updatePDESystem(sourceVec)} \\ & Updates time dependent parts of PDE system\\ \texttt{initialize(source)} \\ & Initializes the system\\ \texttt{initialize(function)} \\ & Initializes the system\\ \texttt{initialize(vector)} \\ & Initializes the system\\ \texttt{solution\_u()} & Returns the local solution vector\\ \texttt{fiber\_vector()} & Returns the local fiber vector\\ \texttt{residual()} & Returns the local residual vector\\ \texttt{potentialFESpace()} \\ & Returns u finite element space\\ \texttt{setBC(BCHandler)} & Sets Monodomain boundary conditions\\ \texttt{postProcess(bool)} & Postprocessing\\ \texttt{resetPrec()} & Activates the reset of the preconditioner (which occurs during the solveSystem method)\\ \texttt{getRepeatedEpetraMap()} \\ & Returns the local repeated map\\ \texttt{getRepeatedEpetraMapVec()}\\ & Returns the local repeated map vector\\ \texttt{getMap()} & Returns the local map\\ \texttt{recomputeMatrix(bool)} \\ & Asks the solver to recompute the matrix\\ \texttt{matrMass()} & Returns the local mass matrix\\ \hline \end{longtable} \end{center} %\caption{Main methods of the MonodomainSolver class} %\label{tab:monodomainMethods} %\end{table} %\subsubsection{Bidomain solver} \subsubsection{BidomainSolver class} % Ricardo Ruiz We consider two superposed domains: an intra cellular domain $\Omega_i$ and an extra cellular domain $\Omega_e$. The equations of the problem read: \begin{eqnarray} \frac{\partial v}{\partial t} - \nabla\cdot(M_i\nabla u_i) & = & f(v,w),\\ \frac{\partial v}{\partial t} - \nabla\cdot(M_e\nabla u_e) & = & f(v,w),\\ \frac{\partial w}{\partial t} & = & g(v,w),\\ + \text{Neumann b.c.}\nonumber \end{eqnarray} where $u_i$ and $u_e$ are the potential in the intra cellular and extra cellular domain respectively, $v=u_i-u_e$, $M_i$ and $M_e$ are given tensors for the intra cellular and extra cellular domain respectively, and $f$, $g$ are two given functional. The main methods of the class are the followings: %\begin{table}[H] \begin{center} \begin{longtable}{p{3cm}p{8cm}} \hline \texttt{PDEiterate(BCHandler)} \\ & Updates sources, bc treatment and solve the bidomain system\\ \texttt{setUp(dataFile)} & Sets up the system solver with a GetPot datafile\\ \texttt{buildSystem()} & Builds time independent parts of PDE system\\ \texttt{updatePDESystem(alpha,sourceVec)} \\ & Updates time dependent parts of PDE system\\ \texttt{updatePDESystem(sourceVec)} \\ & Updates time dependent parts of PDE system\\ \texttt{initialize(source,source)} \\ & Initializes the system\\ \texttt{initialize(function,function)} \\ & Initializes the system\\ \texttt{initialize(vector,vector)} \\ & Initializes the system\\ \texttt{solution\_u()} & Returns the local solution vector for $v$\\ \texttt{solution\_ue()} & Returns the local solution vector for $u_e$\\ \texttt{solution\_uiue()} & Returns the global solution containing $v$ and $u_e$\\ \texttt{fiber\_vector()} & Returns the local fiber vector\\ \texttt{residual()} & Returns the local residual vector\\ \texttt{potentialFESpace()} \\ & Returns finite element space for $\mathbf{u}$\\ \texttt{setBC(BCHandler)} & Sets Bidomain boundary conditions\\ \texttt{resetPrec()} & Activates the reset of the preconditioner (which occurs during the updatePDESystem method)\\ \texttt{getRepeatedEpetraMap()} \\ & Returns the local repeated map\\ \texttt{getRepeatedEpetraMapVec()} \\ & Returns the local repeated map vector\\ \texttt{getMap()} & Returns the local map\\ \texttt{recomputeMatrix(bool)} \\ & Activates the rebuild of the matrix (which occurs during the solveSystem method)\\ \texttt{matrMass()} & Returns the local mass matrix\\ \texttt{bdf\_uiue()} & Returns a BdfT object\\ \hline \end{longtable} \end{center} %\caption{Main methods of the BidomainSolver class} %\label{tab:bidomainMethods} %\end{table} %\subsubsection{FSI solver} \subsubsection{FSISolver class} % Paolo Crosetto The FSI solver is able to solve the problem composed by the Navier-Stokes equations (\ref{eqn:NS1})-(\ref{eqn:NS2}) for the fluid and the St. Venant-Kirchhoff equations for the material in the Lagrangian formulation~(\ref{eqn:VenantKirchhof}). The coupling between the fluid and the structure is provided by three conditions: \begin{itemize} \item the continuity of the velocity at the fluid-structure interface \mathbf{u}\circ\mathcal{A}_t=\frac{d\mathbf{d}_s}{dt} \text{ on }\Gamma; where  simone committed Nov 30, 2010 418 419  \begin{array}{rrcl}  grandper committed Jul 01, 2010 420 421 \mathcal{A}_t: & \Omega_0 & \rightarrow & \Omega_t\nonumber\\ & \mathbf{x}_0 & \mapsto & \mathcal{A}_t(\mathbf{x}_0)=\mathbf{x}_0+\mathbf{d}_f(\mathbf{x}_0)  simone committed Nov 30, 2010 422 423 \end{array}  grandper committed Jul 01, 2010 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 is the ALE mapping from the reference domain $\Omega_0$ to the current domain $\Omega_t$. \item the continuity of stresses at the fluid-structure interface \sigma_s^0\mathbf{n}=J_s\sigma_fF_s^{-T}\mathbf{n} \text{ on }\Gamma, where $\mathbf{n}$ is the outward normal to the solid (or fluid) reference domain, $F_s=I+\nabla\mathbf{d}_s$ is the solid deformation gradient, $J_s$ its determinant and $\sigma_f=\mu\frac{\nabla\mathbf{u}+(\nabla\mathbf{u})^T}{2}-Ip$ is the Cauchy stress tensor; \item the geometric adherence \mathbf{d}_f = \mathbf{d}_s \text{ on }\Gamma. \end{itemize} This solver can treat the equations differently depending on different derivation of the class \texttt{FSIOperator}. Indeed, one can choose a monolitical\footnote{The monolitical approach is available only within the Mathcard project.} approach (where the the structure and the fluid part are solve at once in a big system) or a segregated approach (where the solution for the fluid and the structure are computed one after the other iteratively). The latter approach solves the solid and the fluid problem on different group of MPI processes. Moreover the convective term in the ALE Naviers-Stokes equation can treated explicitly or implicitly. The geometry can also be treated explicitly (we use the position of the mesh at the previous time step) or explicitly (we use the position of the mesh at the current time step). \begin{table}[H] \begin{center} \begin{tabular}{l|c|c|c|c|} & Geo. explicit & Geo. implicit & Conv. explicit & Conv. Implicit\\ \hline Monolitic & \tick & \tick & \tick & \tick \\ Segregated & \tick & & \tick & \\ \end{tabular} \end{center} \caption{Schemes available} \label{tab:FSISchemeAvailable} \end{table}% The main methods of the class are the followings: %\begin{table}[H] \begin{center} \begin{longtable}{p{3cm}p{8cm}} \hline \texttt{timeStep()} & Returns the time step\\ \texttt{timeEnd()} & Returns the final time\\ \texttt{timeInitial()} & Returns the initial time\\ \texttt{FSIOper()} & Returns the FSI operator\\ \texttt{displacement()} & Returns the displacement, which will be on the solid interface map \\ \texttt{isFluid()} & Returns true if the processor solves the fluid part of the problem\\ \texttt{isSolid()} & Returns true if the processor solves the solid part of the problem\\ \texttt{initSol(solInit)} & Sets the initial solution\\ \texttt{setSourceTerms(fluidSource,solidSource)} \\ & Sets the fluid and solid source term\\ \texttt{setFSIOperator(op)} \\ & Set the FSI operator to be used to couple the fluid and the structure\\ \texttt{setDataFromGetPot(dataFile)} & \\ \texttt{setup()} & \\ \texttt{initialize(u0,p0,d0,w0,df0)} & \\ \texttt{initialize(velName,pressName,velwName,depName,velSName,Tstart)} & \\ \texttt{initialize(u0,v0)} & \\ \texttt{iterate(time)} & \\ \texttt{showMe()} & Prints useful informations about the solver\\ \texttt{isMonolithic()} & Returns true if the monolithic approach is used to solve the problem\\ \hline \end{longtable} \end{center} %\caption{Main methods of the FSISolver class} %\label{tab:FSIMethods} %\end{table} Moreover some methods relative to the boundary conditions are implemented. They should be used carefully according to the type of solving method (monolithic (M), segregated with fixe point iterations (SFP) or segregated with Newton iterations (N)) and the type of boundary conditions as presented in the table below: %\begin{table}[H] \begin{center} \begin{longtable}{l|c|c|c|} Required methods & M & SFP & SN\\ \hline \texttt{setHarmonicExtensionBC(BCHandler)} & \tick & \tick & \tick \\ \texttt{setFluidBC(BCHandler)} & \tick & \tick & \tick \\ \texttt{setLineFluidBC(BCHandler)} & & & \tick \\ %\texttt{setInvLinFluidBC(BCHandler)} & & & \\ \texttt{setSolidBC(BCHandler)} & \tick & \tick & \tick \\ \texttt{setLinSolidBC(BCHandler)} & & & \tick \\ %\texttt{setInvLinSolidBC(BCHandler)} & & & \\ \texttt{setFluxBC(BCHandler)} & \tick & &\\ \texttt{setRobinBC(BCHandler)} & \tick & & \end{longtable} \end{center} %\caption{Main methods of the FSISolver class} %\label{tab:FSIMethods} %\end{table} % Boundary conditions \subsection{Boundary conditions} \label{sec:BC} No partial differential equation (PDE) is complete without suitable boundary conditions. In this section we describe the possibilities proposed in the library. LifeV manages the boundary conditions with two classes: \begin{itemize} \item The BCHandler class which stores the boundary conditions (flag of the boundary, type of the boundary condition to be applied and the value given through a function). \item The BCManage class which takes as argument the matrix $A$ and the right hand side $\mathbf{b}$ and the BCHandler. Then the boundary conditions are applied with respect to the BCHandler to modify $A$ and $\mathbf{b}$. \end{itemize} The implemented boundary conditions are the following: \begin{itemize} \item Dirichlet (type = essential) boundary conditions. Can be imposed by components. This conditions read \mathbf{u} = \mathbf{g} \text{ on }\Gamma, where $\mathbf{u}$ is the unknown, $\mathbf{g}$ a given function and $\Gamma$ the boundary of the domain. This condition is imposed by modifying the matrix $A=(a_{ij})$ and the right hand side (rhs) $\mathbf{b}=(b_j)$ in the following way. Suppose that we want to impose $u_i=c$. We set $b_i=c$ and $a_{ij}=\delta_{ij}$ where $\delta_{ij}=1$ if $i=j$ and $\delta_{ij}=0$ if $i\neq j$. \item Neumann (type = natural) boundary conditions. This condition derives naturally from the integration by parts in the weak formulation. It is imposed by modifying the right hand side only.  perego committed Dec 20, 2010 522 \item Robin boundary conditions.  grandper committed Jul 01, 2010 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 \item Dirichlet normal boundary condition. The value of the unknown in the normal direction can be imposed (i.e. $\mathbf{u}\cdot\mathbf{n}=g$, where $g$ is a known function). \item Directional boundary condition. Similar to the Dirichlet normal boundary condition but imposes the value in given direction. \end{itemize} % Finite Element (FE) \subsection{Finite Elements (FE)} \label{sec:FE} At the moment, the finite elements available in LifeV are only 3D tetrahedra elements, namely \begin{itemize} \item $\mathbb{P}0$ \item $\mathbb{P}1$ \item $\mathbb{P}1$Bubble \item $\mathbb{P}2$. \end{itemize} The degrees of freedom and their numbering is presented in Figure~\ref{fig:FE}. \begin{figure}[H] \begin{center}  grandper committed Jul 02, 2010 540 541 542 543 \subfigure[$\mathbb{P}0$]{\includegraphics[width=2cm]{images/fe/p03d.pdf}} \subfigure[$\mathbb{P}1$]{\includegraphics[width=2cm]{images/fe/p03d.pdf}} \subfigure[$\mathbb{P}1$Bubble]{\includegraphics[width=2cm]{images/fe/P1Bubble3D.pdf}} \subfigure[$\mathbb{P}2$]{\includegraphics[width=2cm]{images/fe/p23d.pdf}}  grandper committed Jul 01, 2010 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 \end{center} \caption{Overview of the 3D finite elements available} \label{fig:FE} \end{figure} Other elements and 2D elements are probably retrievable from LifeV-serial. Each finite element is an instance of the \texttt{RefFE} class. %\begin{table}[H] %\begin{center} %\begin{tabular}{m{2cm}|m{2cm}m{2cm}m{2cm}m{2cm}} % & \begin{center} 0D \end{center} & \begin{center} 1D \end{center} & \begin{center} 2D \end{center} & \begin{center} 3D \end{center} \\ %\hline %$\mathbb{P}$0 & \includegraphics[width=2cm]{images/FE/P00D.pdf} & \includegraphics[width=2cm]{images/FE/NA.pdf} & \includegraphics[width=2cm]{images/FE/P02D.pdf} & \includegraphics[width=2cm]{images/FE/P03D.pdf} \\ %$\mathbb{P}$1 & \includegraphics[width=2cm]{images/FE/NA.pdf} & \includegraphics[width=2cm]{images/FE/P11D.pdf} & \includegraphics[width=2cm]{images/FE/P12D.pdf} & \includegraphics[width=2cm]{images/FE/P13D.pdf} \\ %$\mathbb{P}$1Bubble & \includegraphics[width=2cm]{images/FE/NA.pdf} & \includegraphics[width=2cm]{images/FE/NA.pdf} & \includegraphics[width=2cm]{images/FE/NA.pdf} & \includegraphics[width=2cm]{images/FE/P1Bubble3D.pdf} \\ %$\mathbb{P}$2 & \includegraphics[width=2cm]{images/FE/NA.pdf} & \includegraphics[width=2cm]{images/FE/P21D.pdf} & \includegraphics[width=2cm]{images/FE/P22D.pdf} & \includegraphics[width=2cm]{images/FE/P23D.pdf} \\ %$\mathbb{P}$2tilde & \includegraphics[width=2cm]{images/FE/NA.pdf} & \includegraphics[width=2cm]{images/FE/NA.pdf} & \includegraphics[width=2cm]{images/FE/NA.pdf} & \includegraphics[width=2cm]{images/FE/P2Tilde3D.pdf} \\  grandper committed Jul 02, 2010 560 %RT0 & \includegraphics[width=2cm]{images/FE/NA.pdf} & \includegraphics[width=2cm]{images/FE/NA.pdf} & \includegraphics[width=2cm]{images/FE/NA.pdf} & \includegraphics[width=2cm]{images/FE/RT0Tetra3D.pdf}  grandper committed Jul 01, 2010 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 %\end{tabular} %\end{center} %\caption{Overview of the finite elements available} %\end{table} %The $\mathbb{P}$2tilde are Navier-Stokes $\mathbb{P}$2 basis oriented to the mass lumping. %\begin{table}[H] %\begin{center} %\begin{tabular}{m{1cm}|m{2cm}m{2cm}} % & \begin{center} 2D \end{center} & \begin{center} 3D \end{center} \\ %\hline %$\mathbb{Q}$0 & \includegraphics[width=2cm]{images/FE/Q02D.pdf} & \includegraphics[width=2cm]{images/FE/Q03D.pdf} \\ %$\mathbb{Q}$1 & \includegraphics[width=2cm]{images/FE/Q12D.pdf} & \includegraphics[width=2cm]{images/FE/Q13D.pdf} \\ %$\mathbb{Q}$2 & \includegraphics[width=2cm]{images/FE/Q22D.pdf} & \includegraphics[width=2cm]{images/FE/NA.pdf} \\ %RT0 & \includegraphics[width=2cm]{images/FE/NA} & \includegraphics[width=2cm]{images/FE/RT0Hexa3D.pdf} %\end{tabular} %\end{center} %\caption{Overview of the finite elements available} %\end{table} % Input/output in LifeV \subsection{Input/output in LifeV} \label{sec:inOutput} Scientific computing involves a lot of data: Physical parameters, methods or solvers settings, meshes of the domain(s), output of the solution, $\ldots$ In this section we go through the different kind of input and output data used in LifeV. \subsubsection{Input files} The most common way to pass data to LifeV are data files containing parameters and values stored. The informations are then retrieved in the program using GetPot\footnote{More informations about the format and GetPot are available on http://getpot.sourceforge.net}. \subsubsection{Mesh formats} LifeV contains parsers to load mesh from the INRIA mesh (\texttt{.mesh}) format. The following formats are probably retrievable from LifeV-serial: \begin{itemize} \item Mpp mesh (\texttt{.m++}); \item Gmsh mesh (\texttt{.msh}); \item Netgen mesh (\texttt{.vol}). \end{itemize} \subsubsection{3D Structured mesh generator} Instead of loading a mesh from a file, LifeV is able to generate a structured mesh with the function \texttt{structuredMesh3D}. This generator is useful to obtain quickly a mesh of a cuboid. One can also use it to obtain several refinements of a cuboid in order to evaluate the scalability of parallel algorithms or to test convergence rates on test problems. \subsubsection{Output files} When the computations are over, the user can collect the solution into output files. In the older versions of LifeV, the output were in the \texttt{vtk} format. The main disadvantage of this format is that each simulation generate hundreds of files! In the new version, all the outputs are stored in one unique file in the \texttt{hdf5} format. \subsubsection{Postprocessing} Since the outputs are written in the \texttt{vtk} or in the \texttt{hdf5} format, many possibilities are available to use the results. However most of the LifeV's users use Paraview\footnote{More details are available on http://www.paraview.org}. % Algorithms available in LifeV \subsection{Algorithms available in LifeV} LifeV does not only provide solvers but also tools and algorithms which helps to develop applications. \subsubsection{GeneralizedAitken class} \label{sec:aitken} % Simone Deparis This class is an implementation of the Aitken method which permits to compute automatically a relaxation parameter to accelerate the convergence of the fixed-point iterations. We give here a rough idea of the method. More details can be found in~\cite{deparis}. Let $\mathbf{f}$ be a vector map from $\mathbb{R}^d$ to $\mathbb{R}^d$. We want to solve the fixed point problem $\mathbf{f}(\mathbf{x})=\mathbf{x}$, which can be rewritten as $\mathcal{R}(\mathbf{x})=\mathbf{f}(\mathbf{x})-\mathbf{x} = 0$. If $\mathbf{x}^k$ denotes the approximation of $\mathbf{x}$ at the step $k$ of the algorithm (with $\mathbf{x}^k\neq\mathbf{x}$), then $\mathcal{R}(\mathbf{x}^k)=\epsilon\neq0$. We want to calculate $\omega^k$ such that we can find a new approximation $\mathbf{x}^{k+1}=\mathbf{x}^k+\omega^k\mathcal{R}(\mathbf{x}^k)$ optimizing the convergence rate of the series~$\{\mathbf{x}^k\}$. More precisely, we have the following algorithm. \begin{algorithme}[H] \begin{enumerate} \item Initial guess of $\mathbf{x}^0$ and $\omega^0$. \item Set $\mathbf{x}^{k+1}=\mathbf{x}^k-\omega^k\mathcal{R}(\mathbf{x}^k)$ \item Compute $\mathcal{R}(\mathbf{x}^{k+1})$ \item \IfThenElse{$\norm{\mathcal{R}(\mathbf{x}^{k+1})}<\epsilon$}{ Exit }{ Compute $\omega^{k+1}=\frac{(\mathcal{R}(\mathbf{x}^{k+1})-\mathcal{R}(\mathbf{x}^k))\cdot(\mathbf{x}^{k+1}-\mathbf{x}^k)}{\norm{\mathcal{R}(\mathbf{x}^{k+1})-\mathcal{R}(\mathbf{x}^k)}^2}$\\ $k=k+1$\\ Return to $2$.  grandper committed Jul 02, 2010 629 }  grandper committed Jul 01, 2010 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 \end{enumerate} \caption{Aitken acceleration method} \label{algoAitkenAccelerationMethod} \end{algorithme} \subsubsection{RCM reordering} % Paolo Crosetto This algorithm does a reordering of the nodes in such a way that the bandwidth of the matrix of the problem is minimized. The method is described in~\cite{saad}. The implemented code is tested via the test \texttt{test\_reorder}. \subsubsection{Eigen values computation} % Paolo Crosetto Using the Anazazi package from Trilinos\footnote{more informations available on http://trilinos.sandia.gov}, LifeV is able to compute eigen values from a matrix. In particular the test \texttt{test\_eigensolver} in the Mathcard testsuite compute the eigen values of a preconditioner. % Meshes available in LifeV \subsection{Meshes available in LifeV} The LifeV release proposes lots of meshes like cylinders or cubes used in examples. A particular mesh available is a mesh of the carotid artery (see Figure~\ref{fig:carotidMesh}). \begin{figure}[H] \centering \includegraphics[width=6cm]{images/geometriesAndMeshes/carotid.png} \caption{Mesh for blood flow simulations in the carotid artery} \label{fig:carotidMesh} \end{figure} %============================ % Applications %============================ \section{Applications} \label{sec:applications}  grandper committed Jul 02, 2010 659 In this section, we give an overview of the problem solved with LifeV.  grandper committed Jul 01, 2010 660 661 662 663 664 665 666 667  % Convergence: Eithersteinman \subsection{Convergence: Eithiersteinman} % Simone Deparis, Gwenol Grandperrin Eithier and Steinman have proposed an exact fully 3D Navier-Stokes solutions in~\cite{ethiersteinman}. This solution is used as a benchmark in LifeV to check the validity of the code. % Multiscale models \subsection{Multiscale models}  malossi committed Jul 20, 2010 668 In the contexte of blood flows simulation, we want to be able to split the cardiovascular system into several pieces. Each of the pieces can be veins, artery or organs. Moreover to be able to compute such system we want to introduce a 1D model for the smallest vessels. So at the end we will make the computation on a general network of 1D or 3D problems where the edges will be the coupling and the nodes the problems. The details about the theory of this section can be found in~\cite{malossi}.  grandper committed Jul 01, 2010 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718  \subsubsection{Equations} Currently, only the Oseen equations are available for the models in the network. Concerning the boundary conditions between the models one can choose between the two following possibilities: \begin{enumerate} \item flow-stress connection (see Figure~\ref{fig:multiscaleBCa}); \item stress-stress connection (see Figure~\ref{fig:multiscaleBCb}). \end{enumerate} \begin{figure}[H] \centering \subfigure[ ]{ \includegraphics[width=4cm]{images/multiscale/bca.pdf} \label{fig:multiscaleBCa} } \subfigure[ ]{ \includegraphics[width=4cm]{images/multiscale/bcb.pdf} \label{fig:multiscaleBCb} } \caption{Explicit method} \label{fig:multiscaleBC} \end{figure} If the stress and the flow are available the two possibilities are equivalent. \subsubsection{Solving methods} \begin{itemize} \item Method 1 (Explicit):\\ This method only works with simple geometries (no bifurcations) and a very small time-step. We first solve the problem by imposing zero Neumann boundary condition at the coupling boundaries. At each time step, the information is propagated through the boundaries (see Figures~\ref{fig:multiscaleExplicitMethoda} to~\ref{fig:multiscaleExplicitMethodc}). \begin{figure}[H] \centering \subfigure[t=0]{ \includegraphics[width=7cm]{images/multiscale/explicitMethoda.pdf} \label{fig:multiscaleExplicitMethoda} } \subfigure[t=1]{ \includegraphics[width=7cm]{images/multiscale/explicitMethodb.pdf} \label{fig:multiscaleExplicitMethodb} } \subfigure[t=2]{ \includegraphics[width=7cm]{images/multiscale/explicitMethodc.pdf} \label{fig:multiscaleExplicitMethodc} } \caption{Explicit method} \label{fig:multiscaleExplicitMethod} \end{figure} \item Method 2 (implicit)\\ We use an implicit Aitken method (see Section~\ref{sec:aitken}) to solve the system. With this method, there is no limitation for $\Delta t$ due to the network. The evaluation of the residual make this method quite cheap. The method works fine with simple geometries but fails when bifurcations are introduced. \item Method 3 (implicit)\\  malossi committed Jul 20, 2010 719 We use a Newton method to solve the system. We need to compute the jacobian, which is more expensive than the residual of Method 2. Nevertheless, this approach works fine also in presence of complex networks.  grandper committed Jul 01, 2010 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762  \end{itemize} Generally the solution is expensive, so it is better to reduce the number of iterations. \subsubsection{Applications} \begin{enumerate} \item The first application is a network of five cylinders connected to create a continuous vessel. A sinusoidal flow rate is imposed at one extremity. We obtain a Wonersley profile. Moreover an exact solution is available for this problem. \begin{figure}[H] \begin{center} \includegraphics[width=8cm]{images/multiscale/application1.pdf} \caption{Application 1} \label{fig:multiscaleApplication1} \end{center} \end{figure} \item In this application we introduce a bifurcation in the flow. This test allow to check if the flow is split correctly (the top and bottom branch should receive $\frac{1}{4}$ of the flow and the middle one $\frac{1}{2}$). \begin{figure}[H] \begin{center} \includegraphics[width=6cm]{images/multiscale/application2.pdf} \caption{Application 2} \label{fig:multiscaleApplication2} \end{center} \end{figure} \item We introduce a bifurcation in one of the models. Then we replace the bifurcation model by tree cylinders. We can then compare the accuracy of this substitution. \begin{figure}[H] \centering \subfigure[ ]{ \includegraphics[width=7cm]{images/multiscale/application3a.pdf} } \subfigure[ ]{ \includegraphics[width=7cm]{images/multiscale/application3b.pdf} } \caption{Application 3} \label{fig:multiscaleApplication3} \end{figure} \end{enumerate} \subsubsection{Future works}  malossi committed Jul 20, 2010 763 Cristiano Malossi is currently working on the implementation of 1D models and their integration into the models network. He is also collaborating with Paolo Crosetto for the integration of FSI models in the network.  grandper committed Jul 01, 2010 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954  % Aorta with Fluid Structure Interaction (FSI) and delective BC \subsection{Aorta with Fluid Structure Interaction (FSI) and delective BC} % Paolo Crosetto The theory concerning the resolution of FSI problem with parallel algorithm can be found in~\cite{crosetto}. A prior test is available (i.e. \texttt{test\_fsi}) which solves a FSI problem using the Newton method or an Aitken method. The convective term is treated explicitly but the geometry can be treated explicitly or implicitly. A simulation of an aorta is done using Navier-Stokes (derived from the Oseen problem). It is available in \texttt{CMCS applications/application}. \subsubsection{Future works} A test of scalability for the simulation of the aorta with the FSI equations is under construction. % Vessels deformation \subsection{Vessels deformation} % Paolo Tricerri & Margarita An important aspect of blood flows simulation is the vessels deformation. Nonlinear elasticity solver is already implemented in LifeV with the Saint Venant-Kirchhof model (see Section~\ref{sec:venantKirchhofSolver}). The \texttt{test\_nonlinearelasticity} test simulates the deformation of a cylinder. For this test we impose an homogenous Dirichlet boundary condition in the direction of the principal axis of the cylinder on $\Gamma_1$, an homogenous Neumann boundary condition on $\Gamma_2$ and a non homogenous Neumann boundary condition on $\Gamma_3$. \begin{figure}[H] \begin{center} \includegraphics[width=4cm]{images/monolitic/cylindredomain.pdf} \caption{Boundary conditions} \label{fig:monoliticDomain} \end{center} \end{figure} \subsubsection{Future works} In the serial version of LifeV, others models exists and should be ported to LiveV: \begin{enumerate} \item The weakly exponential model; \item The neo Hookean model. \end{enumerate} Moreover to obtain realistic model for the arterial wall deformation one should take into account the effects of elastin and collagen. Two tests (also in Life), \texttt{test\_structure\_1mec} and \texttt{test\_structure\_2mec}, take into account one or both substance respectively. % Heart electrical propagation \subsection{Heart electrical propagation} \label{sec:heartElectricalPropagation} % Ricardo Ruiz A human heart is stimulated with electrical signals. In this application we try to simulate the electrical propagation. A test case is available in the tests suite (i.e. \texttt{test\_heart}). An interesting domain is the truncated ellipsoidal domain (see Figure~\ref{fig:truncatedEllipsoidalDomain}). Indeed the fibers description are known and the problem becomes non-linear for this geometry. \begin{figure}[H] \begin{center} \includegraphics[width=4cm]{images/monoAndBidomain/truncatedellipsoidaldomain.pdf} \caption{Truncated ellipsoidal domain} \label{fig:truncatedEllipsoidalDomain} \end{center} \end{figure} % 1D tree \subsection{1D tree} % Gilles Fourestey, Tiziano Passerini The simulation of the entire cardiovascular system is very ambitious. A way to perform simulations is to consider a simple model: Each arteria, each vein and each capillary can be replaced by a 1D vessel. Putting all these models together we end up with a tree of vessels. %============================ %============================ %============================ %============================ % Other projects %============================ \section{Other projects} \label{sec:otherProjects} LifeV is not limited to an isolate development. Indeed several projects are currently involving the library: \begin{itemize} \item Mathcard project (EPFL-MOX) \item Synergia \end{itemize} The following sections proposed an overview of the applications and current works related to the above projects. % Monolitic (Mathcard) \subsection{Monolitic (Mathcard)} A \emph{monolitic} solver is implemented in LifeV. By \emph{monolitic} we mean that the solid and the fluid equations are solved in a unique matrix. The \texttt{test\_monolitic} tests a problem in a cylinder where the boundary are given in Figure~\ref{fig:monoliticDomain}. We impose an entering flux on $\Gamma_1$, Robin boundary condition on $\Gamma_2$ and absorbing Neumann boundary condition on $\Gamma_3$. The latter condition use Riemann invariant and the entering wave is imposed to $0$. \begin{figure}[H] \begin{center} \includegraphics[width=4cm]{images/monolitic/cylindredomain.pdf} \caption{Boundary conditions} \label{fig:monoliticDomain} \end{center} \end{figure} \subsubsection{Applications} As an applications, a simulation of blood flow in the aorta has been made using the FSI equations (see Figure~\ref{fig:aortaFSISimulation}). This application is available in \texttt{CMCS applications/application}. \begin{figure}[H] \centering \includegraphics[height=10cm]{images/geometriesAndMeshes/crosettoAorta.png} \caption{Aorta FSI simulation} \label{fig:aortaFSISimulation} \end{figure} % Description of an electro-mechanic coupling \subsection{Description of an electro-mechanic coupling (Mathcard)} % Ricardo Ruiz This application solves a problem composed by the equations~(\ref{eqn:MonoDomain1})-(\ref{eqn:MonoDomain2}) and the equation -\nabla\cdot\sigma_{os}=0, where $\sigma_{os}$ is the Piola-Kirschhof tensor. The implementation is given in the \texttt{test\_elecmech} test. However it is not ready yet. % Physiological geometries/meshes (Mathcard) \subsection{Physiological geometries/meshes (Mathcard)} %1 - artery.vtp %2 - bypass.vtp %3 - aorta_angio_IRM_pelvis_bifurcations2_mc_c1.vtp %4&5 - preop-afb.vtp %6 - postop-afb.vtp %7- virtual_aorta_1.vtp %8 - palmaz-stent.vtp %9 - palmaz-stent-expanded.vtp %10 - heart.mesh In order to do practical calculation, LifeV regrouped several physiological geometries for the Mathcard project. Since one of the practical interest is the simulation of the cardiovascular system, meshes simulating the vessels are available such as the bypass (proposed in~\cite{ku}) configuration where the blood flow is divided (see Figure~\ref{fig:vesselsGeometries1}), or an artery bifurcation (see Figure~\ref{fig:vesselsGeometries2}), or the aorta (proposed by \cite{simbios} and illustrated in~Figure~\ref{fig:vesselsGeometries3}). \begin{figure}[H] \centering \subfigure[Bypass]{ \includegraphics[height=3.5cm]{images/geometriesAndMeshes/bypass.png} \label{fig:vesselsGeometries1} } \subfigure[Artery]{ \includegraphics[height=3.5cm]{images/geometriesAndMeshes/artery.png} \label{fig:vesselsGeometries2} } \subfigure[Aorta]{ \includegraphics[height=8cm]{images/geometriesAndMeshes/virtualAorta.png} \label{fig:vesselsGeometries3} } \subfigure[Aorta bifurcations]{ \includegraphics[height=8cm]{images/geometriesAndMeshes/aortaIRM.png} \label{fig:vesselsGeometries4} } \caption{Vessels geometries} \label{fig:vesselsGeometries} \end{figure} The collection of geometries contains also particular situation like aneurysm geometry (before and after surgical operation, see Figures~\ref{fig:vesselsGeometries1} and~\ref{fig:vesselsGeometries2}; the two geometries are proposed in~\cite{wilson1} and~\cite{wilson2}) and medical device geometries like stents (compressed and extended version are showed on Figures~\ref{fig:vesselsGeometries3} and~\ref{fig:vesselsGeometries4}\footnote{Peter H. Feenstra of the Cardiovascular Biomechanics Research Laboratory at Stanford University is acknowledged for providing the geometry of the Palmaz-Schatz-like stent~(2006)}). \begin{figure}[H] \centering %\subfigure[ ]{ %\includegraphics[width=4cm]{images/geometriesAndMeshes/preopAFB1.png} %\label{fig:practicalInterestGeometries} %} \subfigure[Aneurysm before operation]{ \includegraphics[height=5cm]{images/geometriesAndMeshes/preopAFB2.png} \label{fig:practicalInterestGeometries1} } \subfigure[Aneurysm after operation]{ \includegraphics[height=5cm]{images/geometriesAndMeshes/postopAFB.png} \label{fig:practicalInterestGeometries2} } \subfigure[Palmaz stent]{ \includegraphics[height=3.7cm]{images/geometriesAndMeshes/stent.png} \label{fig:practicalInterestGeometries3} } \subfigure[Palmaz extended stent]{ \includegraphics[height=3.7cm]{images/geometriesAndMeshes/stentExpanded.png} \label{fig:practicalInterestGeometries4} } \caption{Practical interest of geometries} \label{fig:practicalInterestGeometries} \end{figure} Of course as presented in Section~\ref{sec:heartElectricalPropagation}, simulations can also target phenomena on organs like the heart (see Figure~\ref{fig:heartMesh}). \begin{figure}[H] \centering \subfigure[Front view]{ \includegraphics[height=5cm]{images/geometriesAndMeshes/heartfront.png} \label{fig:heartMesh1} } \subfigure[Top view]{ \includegraphics[height=5cm]{images/geometriesAndMeshes/hearttop.png} \label{fig:heartMesh2} } \caption{Mesh for the simulation of the heart} \label{fig:heartMesh} \end{figure} % Levelset (Synergia/bioreactor) \subsection{Level set (Synergia/bioreactor)} % Samuel Quinodoz The management of a level set function has been introduced in LifeV for the Synergia project. The aim is to simulate a bioreactor. The levelset is used to follow the deformation of the fluid in a tank full of water and nutrients. Let us consider a partition $\Omega_1(t)$, $\Omega_2(t)$ of a domain $\Omega(t)$ (which can depend on the time $t$) and let $\Gamma_I$ denotes the interface between the two subsets. The level set method describes the interface as the zero level set of a function $\phi$:  grandper committed Jul 02, 2010 955   grandper committed Jul 01, 2010 956 957 958 959 960 961 962 963 964 965 966 967 968 \label{levelsetfct} \left\{ \begin{array}{ll} \phi(\mathbf{x},t) > 0 & \mathbf{x}\in\Omega_1(t)\\ \phi(\mathbf{x},t) < 0 & \mathbf{x}\in\Omega_2(t)\\ \phi(\mathbf{x},t) = 0 & \mathbf{x}\in\Gamma_I(t) \end{array}\right. so that \Gamma_I(t) = \{\mathbf{x}\in\Omega(t)\mid\phi(\mathbf{x},t)=0\}. Let $H$ be the Heavyside function defined by  grandper committed Jul 02, 2010 969   grandper committed Jul 01, 2010 970 971 972 973 974 975 976 977 \label{heavysidefct} H(z)= \left\{ \begin{array}{ll} 0 & z<0\\ 1 & z>0 \end{array}\right. , such that the density and the viscosity quantities in $\Omega$ can be rewritten as  grandper committed Jul 02, 2010 978   grandper committed Jul 01, 2010 979 980 981 982 983 984 985 986 \label{heavysidefct} \begin{array}{ll} \rho = \rho_1+(\rho_2-\rho_1)H(z)\\ \mu = \mu_1+(\mu_2-\mu_1)H(z) \end{array}. The final system that describes the evolution of the fluid is composed of the Navier-Stokes equations~(\ref{eqn:NS1}) and~(\ref{eqn:NS2}) as well as the equation~(\ref{eqn:levelset}).  grandper committed Jul 02, 2010 987   grandper committed Jul 01, 2010 988 989 990 991 992 993 994 995 996 997 998 999 \label{eqn:levelset} \begin{array}{ll} \frac{\partial\phi}{\partial t}+\mathbf{u}\cdot\nabla\phi = 0 & \text{in }\Omega\times[0,T]\\ \phi=\phi_0 & \text{in }\Omega\text{ at }t=0\\ \phi=\phi_{in} & \text{on }\partial\Sigma_{in}\times[0,T] \end{array} where $\Sigma_{in} = \{(\mathbf{x},t)\in\partial\Omega\times(0,T)\mid\mathbf{u}(\mathbf{x},t)\cdot\mathbf{n}<0\}$. The main advantages of this method is that the interface remains sharp and handle any change in the topology of $\Gamma_I$. However it has drawbacks such as the distorsion of the interface caused by the numerical errors. To overcome this problem, one should use reinitialization of the interface. \subsubsection{Level set solver}  grandper committed Jul 02, 2010 1000 The level set solver solves the equations~(\ref{eqn:levelset}). It manages also the reinitialization of the interface by using the geometric approach. Methods to compute the surface and the volumes of the domains and to export the interface are available (only for $\phi$ approximate with $\mathbb{P}1$ finite element).