Overview / Introduction

Fluid and structure mesh are first created separately, and then used in svFSI to perform Fluid-Structure-Interaction simulations with Arbitrary Lagrangian–Eulerian (ALE) coordinates.

Cardiovascular applications that involve significant deformation are difficult to simulate and model. Some examples include vascular simulations where the wall dilates to more than 10% of the vessel diameter, cardiac simulations where the heart wall displaces as it contracts and expands, or coronary simulations where the coronary vessels displace with the motion of the heart.

One way to computationally model situations like these is to use a fluid-structure interaction (FSI) solver. In FSI, separate domains are defined for the fluid part and solid parts of the computational geometry. The respective equations governing fluid flow (typically Navier-Stokes) are solved in the fluid domain, while the equations governing solid mechanics are solved in the solid domain. The two domains then interact through their interface where the solution variables (displacements, velocities, pressures, stresses) are required to match. This interface acts as a coupled boundary condition for both domains, as the solution of one domain will affect the solution in the other, and vice-versa.

SimVascular utilizes Arbitrary Lagrangian-Eularian (ALE) coordinates to perform FSI. As the name implies, this method is a combination of the Eularian and Lagrangian descriptions of motion that is particularly well suited for FSI problems. In the Eularian description, the computational mesh stays fixed and the motion is characterized as velocities flowing past the grid nodes. The Eularian description is typical for most fluid dynamics problems without FSI. While in the Lagrangian description, the mesh nodes move exactly with the motion of the fluid or solid and thus are often characterized as moving domain problems. Lagrangian problems are also often described in terms of the reference configuration (the initial computational geometry before any motion) and the current configuration (the current state of the geometry and mesh nodes). ALE combines these descriptions into a formulation that is convenient for describing FSI problems, where the mesh motion does not exactly match with the fluid motion, but instead a mesh velocity term is added to the convective term in the Navier-Stokes equations. More information can be found in the literature.

SimVascular provides an ALE-FSI solver in svFSI that has many features that make it convenient for modeling cardiovascular flows. First, svFSI is able to directly read in mesh formats that SimVascular outputs and use them directly in simulations. Second, svFSI provides capabilities to apply inflow and outflow boundary conditions that are common in cardiovascular modeling. These boundary conditions include specifying the flow rate at an inlet face with either a flat or parabolic profile, Neumann pressure boundary conditions, and outlet resistance boundary conditions that relate the flow and pressure.

Creating the mesh for ALE simulations

To run an FSI simulation we need a mesh for both the structural domain and the fluid domain. These two meshes must have their interface nodes coincide exactly in order to satisfy the interficial conditions that result from conservation of mass and momentum. The coincident nodes of the fluid mesh are mapped onto the corresponding nodes on the structural mesh and the solution of velocity, displacement, pressure, etc. are treated as equal in the structural and fluid domains.

The fluid domain geometry for patient-specific anatomies are generated using the usual SimVascular modeling pipeline. The Modeling Guide describes this process in more detail: http://simvascular.github.io/docsModelGuide.html. To create the geometry for the structural domain, we will make use of the boundary layer meshing feature in the Meshing module. Normally, the boundary layer meshing feature creates a thin layer of elements on the exterior of your fluid mesh to improve the accuracy of computations near the wall. Conceptually, this is very similar to how a boundary layer forms in fluid mechanics. Velocity gradients are usually high near the walls, so a fine mesh resolution is usually required to resolve this accurately.

The usual case for boundary layer meshing involves extruding this thin layer of elements inwards, starting from the walls and going into the direction of the vessel centers. To make a wall mesh, we will instead use the boundary layer meshing feature to extrude elements outwards to effectively make a new mesh with a specified thickness that surrounds our fluid domain. This new mesh will form the geometry of our structural domain.

Before we use the boundary layer meshing to extrude outwards, it is extremely important that we thoroughly smooth our model. The outward extrusion of elements has the potential to cause some elements to extrude into each other and overlap, which will cause the extrusion to fail. This is especially true at bifurcations, where the extruded elements from the two daughter branches could easily run into each other near the junction if not properly smoothed. We will therefore make thorough use of the local smoothing tools in the Models module before meshing. Below is an example of a bifurcation that would normally occur after lofting segmentations together without smoothing. If we try to extrude a boundary layer mesh from this, the elements at the sharp corner will intersect with each other and cause it to fail.

Sharp bifurcation that will cause problems with wall mesh extrusion.

After creating your model from the Segmentations, you should see a list of all the surfaces in your model in the “Face List” section. The first thing we should do is remesh the entire surface to aid in the surface smoothing. Click on the “Global Ops” menu to the right of the “Face List”, then click on the “Remesh” option to open up the remeshing tool. Click “Estimate Size” to choose an appropriate size for your model, then hit “Remesh” beneath the selected size. This should remesh the surface with triangles that should all be approximately the same size.

Remeshing in the Global Ops tab.

Now that we have remeshed our entire outer surface, we can start to smooth the walls. We only want to smooth the walls and leave the caps alone, so we will make use of the smoothing tools in the “Local Ops” tab. Use Ctrl + click to select all the “wall” surfaces in your model from the Face List, then click “Local Ops” to reveal the list of local operations. We will make use of the “Decimation” and “Cstr. Smooth” tools to perform our smoothing. Decimation will smooth the outer surface of the model by reducing the amount of surface triangles, while Cstr. Smooth will attempt to fillet over any sharp angles in the model while retaining the overall volume. Using these tools in tandem will iteratively smooth the model.

Before starting, check to confirm that all of the wall surfaces are selected in the Face List and none of the caps are selected. Now click the “Decimation” to open the Decimation menu, then click “Decimate Local” to decimate using the default settings.

Local decimation.

Now, click “Cstr. Local” to open up the constrained smoothing menu, then click “Smooth Local” to perform smoothing with the default settings. You should notice your model is a little smoother at the junctions.

Local constrained smoothing.

Repeat this cycle of Decimation and Constrained Smoothing several times to achieve a smoothed model. We recommend at least 10 repetitions of decimation and constrained smoothing. Repeated use of the Decimation feature will sometimes cause the surface of your model to lose so many triangles that the surface loses many of its key geometrical features. It is thus recommended to do a global remeshing every 5 decimations to ensure adequate surface accuracy. Once you are done smoothing, your bifurcations should be sufficiently smoothed and look something like this.

Smooth bifurcation.

Now that we have a smoothed model, we can open the “Meshes” module. Right-click on “Meshes” in the SV Data Manager and click “Create Mesh”. A new window will pop up asking you which model you would like to mesh, which meshing package to use, and the mesh name. Select the model you just smoothed in the previous step, select TetGen for the meshing package, and choose a name. You may leave the name section blank and SV will use the same name as your model. Once you have made these selections, the “SV Meshing” window should open and look like the image below.

SV Meshing.

We must first select an appropriate Global Max Edge Size for our model. This edge size will determine the thickness of our wall mesh. We will remesh the volume of the fluid and solid domains later, so select the edge size to be double what your desired thickness is. We select this to be double to cut down the meshing time. If you are not sure how to select the thickness of your model, choosing a thickness that is 10% of the mean radius in your model is a reasonable assumption used throughout the literature. Below the edge size selection, you should see a box to select boundary layer meshing. First, click the checkbox next to “Boundary Layer Meshing” to turn it on. Below this, you should notice three spaces to select parameters of the boundary layer meshing. Since we are producing this boundary layer for the wall mesh, we can use the same settings for these. “Portion Edge Size” determine the overall thickness of our boundary layer mesh, as a fraction of the Global Max Edge Size selected above. Since we selected a Global Max Edge Size to be double our desired thickness, we want this parameter to be 0.5. Next, you will have to choose the “Number of Layers” in your boundary layer mesh. Increasing this number will increase the accuracy of your structural domain calculations but also increase the number of elements and thus increase your cost. A reasonable number for this parameter is 2. Last, we must select the “Layer Decreasing Ratio”, which is a parameter that allows subsequent layers to be a smaller size than the one before it. Since this parameter does not matter too much for creating a wall mesh, we can select this to be 1.0 to make it so all our layers are the same size.

Below these parameter are three checkboxes. The “Extrude Boundary Layer Inward from Wall” checkbox will extrude the boundary layer mesh inwards if selected. Since we want to extrude the boundary layer mesh outwards, we will uncheck this box. The next box is a setting for “Use Constant Boundary Layer Thickness”, which will attempt to make the entire boundary layer the same thickness if selected. We recommend leaving this box unchecked so the boundary layer mesh can adaptively change thickness in areas of tricky geometry. The last checkbox, “Convert Boundary Layer to New Region/Domain”, is very important. We want to check this so that we will have a way to separate the fluid and structural domains later. Once you have set all of these settings, you boundary layer meshing box should look like this:

Boundary layer meshing settings.

Now, you are ready to run the mesher. Click “Run Mesher” near the top right of the “SV Meshing” window to run the mesher. If the meshing was successful, you should see a window pop up to inform you of the statistics of your mesh. If the meshing was unsuccessful, it is likely that your model may need to be smoothed more to avoid intersecting elements. Once you are successful in producing this boundary layer mesh, right-click it from the SV Data Manager, and click “Export Mesh-Complete” and choose a location to send the mesh.

If the number of elements is larger or smaller than you would like, feel free to adjust the global max edge size accordingly, making sure to adjust the boundary layer mesh settings to ensure you maintain the same wall mesh thickness. The “right” number of elements will depend on your geometry, but generally getting somewhere between 100,000 and 500,00 total elements will have decent accuracy for a reasonable cost. If you have access to a large parallel computing cluster where you can run the simulation in parallel with around 100 processors, we recommend about 1 million elements.

Navigating to the location where you saved the mesh-compelte, you should see that SimVascular has produced two folders. One folder contains the mesh for the fluid part (i.e. domain-1) and the other folder contains the mesh for the structural part of your domain (i.e. domain-2). These two folders are the inputs that are needed in the next step, which is setting up the simulation. We need to do one extra step in this phase before launching the svFSI plugin. If you go into the mesh-complete folder for the fluid domain, you should notice that there is a file called “wallscombined.vtp”. We will need to copy this file into the “mesh-surfaces” folder. Now go into the “mesh-surfaces” folder, and you should notice that there are many “wall"faces. We need to move these out of this folder. Select all of the "wall” faces, EXCEPT for “wallscombined.vtp” and move them to the parent folder. The wallscombined file will contain all the information from the other walls, but combined into 1 file to make assigning boundary conditions easier. In the end, your mesh-surfaces folder should contain .vtp files for each of the caps and the wallscombined.vtp. An example with an aorta model is shown below:

mesh-surfaces folder for the fluid domain.

Now exit from this folder, and enter the mesh-complete folder for the structural domain (i.e. domain-2). Inside this folder, you should notice two wallscombined (region0 and region1). We will need to copy both of these into the mesh-surfaces folder. The region0 file is the .vtp which designates the interface between the solid and fluid domains, while region1 is the outer wall. Navigate into the mesh-surfaces folder, and like we did in the fluid domain, remove all wall files EXCEPT for the two wallscombined files that we just put in here. Once you are done, your mesh-surfaces folder for the structural domain should look like below, where there is a surface for each cap in the model and the two walls_combined.

mesh-surfaces folder for the solid domain.

Basic GUI plugin use

The plugin provides a GUI with which to create input files for your svFSI simulation using SimVascular. In this section, we will describe how to set up a fluid-structure interaction simulation of flow in a thick-walled blood vessel. For simplicity in presentation, we assume that the vessel has only one inlet and outlet. Before setting up a fluid-structure interaction job, one should have a mesh for the fluid domain and a mesh for the annulus (artery wall) domain as described above. We have provided you with an example meshes for the fluid domain of a cylinder and the solid domain of the cylinder.

Initiate a job by pressing “New Job.” Click the job in the SV Data Manager, and ensure that it appears in the plugin window. The “Domains” panel allows one to add meshes to the simulation for the different physical domains. Click “Add Mesh-Complete” and select the fluid mesh. You should see three faces appear, one for in the inlet, outlet and fluid-solid interface. Ensure that the button for “fluid” is selected.

The domains panel.

Click “Add Mesh-Complete” and select the solid mesh. You should see four faces appear, one for in the solid inlet and outlet outlet (through which we expect no flow), one for the fluid-solid interface, at which the fluid and solid domains will be coupled, and one for the outer annulus wall, which the the edge of the physical domain of our problem. Ensure that the button for “solid” is selected.

The next panel “Physics,” which allows one to select the appropriate differential equations for your problem. Select “FSI” for a coupled fluid-structure interaction simulation and click the right arrow. This will also add “Mesh Motion,” based on the ALE method, to the domains. In the list of added equations, click “FSI” to set fluid-structure parameters. Suppose for now that we are working in CGS units. Set fluid density to 1.06 g/cm$^{3}$, viscosity to 0.04 g/(cm s), solid density to 1.06 g/cm$^{3}$, Set the elastic modulus to a nominal value of $7e7$ dynes/cm$^{2}$ (See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4009743/). We expect the artery wall to be nearly incompressible, and set the Poisson ratio to 0.49. Finally, select a “nHK” for a neo-Hookean solid.

The physics panel.

The mesh update equation determines the mesh dynamics in the ALE method. The parameter in the physics tab for this equation is the Poisson’s ratio, which should be set in the interval (0,.5), not too close to the boundary. A good default is .3.

The next tab is “output.” Select a quantity such as “WSS” for wall sheer stress, and click the right arrow to add it to variables that are output. Generally, for each piece of physics, there is a separate list of output variables. In the ALE method, since mesh motion is considered a separate domain and it has its own set of outputs.

The next tab is “boundary conditions.” Select “add,” which brings up a boundary condition menu. From there, select the type of boundary condition and other parameters accordingly. Set the lumen inlet to be a constant value Dirichlet condition of -10 cm$^{3}$/s. Select “impose flux” and a parabolic profile. Set the lumen outlet to be a constant value Neumann condition with value zero, corresponding to zero applied traction at the outlet. Set the annulus inlet and outlet to be zero Dirichlet boundary conditions, corresponding to zero displacement, which will hold them in place through the simulation. Set the annulus outer wall to a zero Neumann condition, which specifies no applied traction on the outer boundary of the physical domain.

Setting a Dirichlet boundary condition.

To couple the fluid and solid at their interface, set a Dirichlet boundary condition on “interface_solid.” Under options, select “projection” and pick the face “interface_fluid.” Do not apply a separate boundary condition for “interface_fluid,” as the coupling will handle the updates. (Note that one must click Dirichlet or Neumann to get the projection box to appear. Such a condition is not added to the input file.)

Setting a projection boundary condition at the fluid/solid interface.

For the mesh equation, set zero Dirichlet conditions on the lumen inlet and outlet. Do not set boundary conditions for the other faces.

The advanced tab options are described here, but changes are not needed for basic simulations. The button “coupled” allows the user to specify a lumped parameter network for model inlets and outlets. This requires compiling a separate executable and is beyond the scope of this tutorial. Tolerance controls when the nonlinear Newton’s iteration terminates. Residual dB reduction terminates also terminates the nonlinear Newton’s iteration, in this case when the residual decreases sufficiently from iteration to iteration. Min iterations and max iterations determines the range of nonlinear iterations allowed. Backflow stabilization sets the value of the coefficient that reduces the so-called “open-boundary instability.” This is required to be in (0,1), but is best left at the default value of .3.

A number of linear solver parameters are available; the default should suffice for now. Finally, there is “Remesher.” If this is off, then the ALE solver will not re-mesh the domains, regardless of mesh distortion. If this value is set to “TetGen,” then the ALE solver will re-mesh automatically when it decides it is necessary. Set the max edge sizes to be the same as your previous mesh spacings for each domain. Use default values for the remaining parameters.

The “Simulation Parameters” panel allows specification of general parameters for the solver. First, set the desired time step and number of steps. The value of the time step depends on your desired resolution in time and the maximum stable time step for your problem, which is determined by the spatial resolution of the mesh, underlying physics and choice of spatial discretization. The option “start from a previous simulation” allows restarting from a previous simulation. This is useful if there was a problem, or one simply wants to run the simulation for longer without re-running the first set of time steps. Under “Save simulation results,” the value start determines the first time step that is output. Increment determines how often output files are saved. Set this to one to output every time step, or to a larger value to only save a subset of the steps.

The simulation panel.

Finally, on the “Run Simulation” panel, click “Create Input File” to create the necessary files to run your simulation. If you wish to run with more than one processor, increase the slider accordingly. Finally, click “Run” to run your simulation.

SVFSI will create a directory called n-procs, where n is the number of MPI tasks for the simulation. This directory will contain .vtu files that with values of all requested fields, as well as a log file called history.dat and averages of various quantities over time.

Example simulation results

In this section, we show some example simulation results using the provided cylinder model. We will simply inflate the cylinder using a ramp for the pressure with a resistance outlet. Because the velocities and pressures in the domain are usually initialized to zero, it can be a good idea to slowly ramp up the pressure to avoid any sudden jump conditions that could produce oscillations. We apply a resistance boundary condition at the outlet to ensure a physiologic amount of flow. Below shows the inlet pressure vs. time, showing that we will ramp up slowly from 0 pressure to 80 mmHg over the course of a second.

Ramp pressure.

For reference, this simulation was run on a parallel cluster using 48 processors. Overall the simulation took about 10 hours to complete 1000 timesteps. The figure below shows the resulting displacement distribution. The wall has a thickness of approximately 0.016 cm, and we selected a physiologic Young’s modulus of about 1 MPa. The figure below shows the resulting displacement distribution across the cylinder.

Cylinder displacements.

In this appendix section, we discuss how to build svFSI from source.

Building svFSI from source

The source code for svFSI includes a build system in cmake. To build from source, one needs compilers for c,c++, fortran, an MPI compiler and runtime, and the lapack and blas libraries. LAPACK and BLAS must be installed in a standard location, which cmake will find automatically. For example, on OSX to install using the package manager brew (see https://brew.sh/) one might run

brew install cmake gcc open-mpi lapack 

or, on ubuntu, to install using apt-get run sudo apt-get install to add the following packages:

cmake
cmake-curses-gui
cmake-gui
gcc 
gfortran
libopenmpi-dev
libblas-dev
liblapack-dev
git

We recommend letting cmake find the default compilers for easy of building.

To build the svFSI solver, download the source code from

https://github.com/simvascular/svFSI

to a directory of your choosing. Make a build directory and go there.

mkdir svFSI\_build && cd svFSI\_build 

Initiate the cmake terminal interface to generate makefiles.

ccmake ../svFSI 

This will automatically search for compilers. Follow instructions if necessary. Push “c” for configure repeatedly until cmake presents the option “g” for generate. Hit “g” to create makefiles and exit. Run make:

make 

This will place the solver binary, called “svFSI” in a directory called svFSI_build/svFSI-build/bin. This will also create a script svFSI_build/svFSI-build/mysvfsi. If non-defuault compilers were passed to cmake, then runtime errors can occur, especially relating to libraries. If this occurs, try using the script, which will set paths appropriately.

  • Optional Trilinos package

The svFSI solver contains its own linear solver. Optionally, the advanced user may link to the Trilinos package. To accomplish this, the user must build Trilinos using cmake, then provide a path and change a single flag. Before building svFSI, download the Trininos source from github.

git clone https://github.com/trilinos/Trilinos.git

We recommend compiling the following branch, with which we have tested Trilinos.

git checkout remotes/origin/trilinos-release-12-10-branch

Make a build directory and go there.

mkdir Trilinos\_build && cd Trilinos\_build 

Run the cmake GUI to build Trilinos.

ccmake -DCMAKE_BUILD_TYPE=RELEASE \
      -DCMAKE_INSTALL_PREFIX:PATH=~/Trilinos_build \
      -DTPL_ENABLE_MPI:BOOL=ON \
      -DTrilinos_ENABLE_Amesos:BOOL=ON \
      -DTrilinos_ENABLE_AztecOO:BOOL=ON \
      -DTrilinos_ENABLE_Epetra:BOOL=ON \
      -DTrilinos_ENABLE_Ifpack:BOOL=ON \
      -DTrilinos_ENABLE_ML:Bool=ON \
      -DTrilinos_ENABLE_MueLu:Bool=ON \
      -DTrilinos_ENABLE_TESTS:BOOL=OFF \
      ../Trilinos

Note that CMAKEINSTALLPREFIX should point to your build directory. Run make to build the Trilinos source

make 

Run make install to place header and library files

make install 

This will create a directories lib/cmake/Trilinos and inside them a file called TrilinosConfig.cmake

lib/cmake/Trilinos/TrilinosConfig.cmake

If this directory does not exist or is empty, then an error has occured.

In the svFSI cmake source, change the SVUSETRILINOS option to on. This can be found in svFSI/Code/CMake/SimVascularOptions.cmake.

option(SV_USE_TRILINOS "Use Trilinos Library with svFSI" ON)

Now, when compiling svFSI, add

ccmake -DCMAKE_PREFIX_PATH=~/Trilinos_build/lib/cmake/Trilinos  ../svFSI

Run make.

make 

If Trilinos is found you should see output that says “Found Trilinos!” and displays the associated variables. The output should look something like this, wherein the elipsis contain many libraries and packages.

Found Trilinos!  Here are the details: 
   Trilinos_DIR = ~/Trilinos_build/lib/cmake/Trilinos
   Trilinos_VERSION = 12.10.1
   Trilinos_PACKAGE_LIST = MueLu;
   ...
   Trilinos_INCLUDE_DIRS = /Users/alex/sfw/Trilinos_build/include
   Trilinos_TPL_LIST = DLlib;LAPACK;BLAS;MPI;Pthread
   Trilinos_TPL_INCLUDE_DIRS = 
   Trilinos_TPL_LIBRARIES = /usr/lib/libdl.dylib
   ...
   Trilinos_BUILD_SHARED_LIBS = FALSE
   Trilinos_CXX_COMPILER = /usr/local/bin/mpicxx
   Trilinos_C_COMPILER = /usr/local/bin/mpicc
   Trilinos_Fortran_COMPILER = /usr/local/bin/mpif90
   Trilinos_CXX_COMPILER_FLAGS =  -std=c++11 -O3 -DNDEBUG
   Trilinos_C_COMPILER_FLAGS =  -O3 -DNDEBUG
   Trilinos_Fortran_COMPILER_FLAGS =  -O3
   Trilinos_LINKER = /usr/bin/ld
   Trilinos_EXTRA_LD_FLAGS = 
   Trilinos_AR = /usr/bin/ar
End of Trilinos details

Setting WITH_TRILINOS to true

If these are filled in, then cmake has found Trilinos and it will be linked. The Trilinos linear solvers are then available for use. In the linear solver options there are now the following Trilinos options for preconditioners.

Trilinos-Diagonal
Trilinos-BlockJacobi 
Trilinos-ILU
Trilinos-ILUT
Trilinos-IC
Trilinos-ML

This will automatically set the Trilinos linear solvers to be called, rather than the svFSILS linear solvers. The three options are GMRES,BICG,CG. There is currently no option to call Trilinos without a preconditioner. If one wishes to use the simplest possible preconditioning with Trilinos, use Trilinos-Diagonal.

For example, set

LS type: GMRES
{
  Preconditioner: Trilinos-ILU
  Tolerance:           1D-3
  Krylov space dimension: 200
}
  • Plugin installation

To build the svFSI plugin, download the source code from

    https://github.com/SimVascular/SimVascular-Plugin-svFSI

to a directory of your choosing. (Note that this repo is currently private, and the plugin is included with SimVascular). The plugin source code should not be in your SimVascular directory tree. Make a build directory and go there.

    mkdir SimVascular-Plugin-svFSI_build && cd SimVascular-Plugin-svFSI_build 

Initiate the cmake terminal interface to generate makefiles.

    ccmake ../SimVascular-Plugin-svFSI

This will automatically search for compilers. Follow instructions if necessary. Push “c” for configure repeatedly until cmake presents the option “g” for generate. Hit “g” to create makefiles and exit. Run make:

    make 

This will build the plugin. If passed, there should be a directory

    SimVascular-Plugin-svFSI/build/lib/plugins

Add that directory to the environment variable

    SV_PLUGIN_PATH

for example

    export SV_PLUGIN_PATH=$SV_PLUGIN_PATH:~/SimVascular-Plugin-svFSI_build/lib/plugins

That directory should contain a library file, a “.dylib” on mac. It should follow the naming convention

    libNAME.extension 

For example,

    liborg_sv_gui_qt_svfsi.dylib

Strip the “lib” prefix and the file extension, then add it to the variable

    SV_CUSTOM_PLUGINS 

as in

    export SV_CUSTOM_PLUGINS=$SV_CUSTOM_PLUGINS:org_sv_gui_qt_svfsi

Finally, if running on linux also set

    LD_LIBRARY_PATH

For example,

    export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:~/SimVascular-Plugin-svFSI_build/lib/plugins

SimVascular will look for these variables in the current environment, and assuming that everything has passed, load the plugin without any additional information. It is a good idea to use the same compilers for anything you build, which cmake will do be default if also building SimVascular from source. These instructions are generic for SimVascular plugin, and work for plugins other that the svFSI plugin.\

In this appendix section, we discuss an alternative method for producing a wall mesh using the software Meshmixer. For users with highly complex geometrical models with a large number of bifurcations, this method provides a higher degree of flexibility and robustness, at the cost of a higher number of steps and higher complexity.

Creating the mesh for ALE simulations using Meshmixer

To run an FSI simulation we need a mesh for both the structural domain and the fluid domain. These two meshes must have their interficial nodes coincide exactly in order to satisfy the interficial conditions that result from conservation of mass and momentum. The coincident nodes of the fluid mesh are mapped onto the corresponding nodes on the structural mesh and the solution of velocity, displacement, pressure, etc. are treated as equal in the structural and fluid domains.

We will now go through the procedure to create such a mesh step by step.

  1. Create the geometry for the fluid domain
  2. Create the geometry for the structural domain
  3. Create the mesh for the structural domain
  4. Create the mesh for the fluid domain

We will use both SimVascular and Meshmixer, which is a freely available software from https://www.meshmixer.com. Meshmixer is a versatile tool that is useful for preparing meshes. We will cover a subset of features in Meshmixer that is convenient for preparing meshes for FSI simulations, but more detailed documentation can be found at https://www.mmmanual.com.

Creating the geometry for the fluid domain

For most cardiovascular modeling applications, the geometry of the fluid domain is generated by segmenting blood vessels out of medical image data. This process is described on the main SimVascular documentation: http://simvascular.github.io/docsModelGuide.html. We now want delete all the caps off the model. Once you have done that, export it by right-clicking the model from the left-hand menu and selecting “Export as Solid Model”. When SimVascular prompts you for a name and location for the exported model, make sure to add an .stl extension to make sure the exported model is in .stl format. We will perform the next step in Meshmixer, and Meshmixer only takes in .stl format surfaces.

Exporting the fluid geometry from SimVascular as an .stl for Meshmixer.

Creating the solid model for the structural domain

In this step, we will import the fluid geometry into Meshmixer and utilize Meshmixer’s “Extrude” feature to create the wall geometry on top of the fluid. The resulting structural geometry resembles a “shell” that surrounds the fluid domain. We will also be using many of the geometry altering features within Meshmixer to smooth the structural geometry. Extrusions in Meshmixer often create intersecting elements and rough patches of geometry that are not suitable for meshing. These must be cleaned up in Meshmixer before meshing in SimVascular.

  1. Import the .stl file for the fluid domain into Meshmixer

When you open Meshmixer, you should see a button called “Import” that appears near the middle of the screen with a plus sign (+). Click this and select the fluid domain .stl that you created in the last step. After importing, you should see your fluid domain model appear in the main window as shown below.

Opening Meshmixer and importing a fluid model.

At this point, we are going to want to smooth all bifurcations and sharp angles in the fluid model. If we do not do this, extruding out the elements in these sections will cause the extruded elements to run into each other and intersect, which will cause problems when attempting to mesh. Ideally, we would like to turn all these sharp corners into rounded filets, as if you rolled a ball of clay into the sharp corners.

We will cover three main ways to smooth a model in Meshmixer: sculpting, local remeshing, and local smoothing. These all work well for smoothing a model, and your usage of them will depend on your personal preference and what works best with the models you are working with. We will first cover each of these tools in detail to discuss how to use them, then we will outline a general procedure that uses these tools that can be applied to most bifurcations. But keep in mind that a different combination of these tools might work better for you and your models.

A couple notes to keep in mind before modifying the geometry of your fluid geometry. If you do not like a particular change that you just applied, you can undo it using Ctrl + Z. You can adjust the center of rotation of the camera by hovering your cursor over a part of your model and hitting C. This is useful when working on specific parts of the model, and the camera rotation center is way off. Also be sure to not modify the geometry around the caps as much as possible. The cap geometry must be maintained in order to get a proper wall mesh around them.

a. Sculpting in Meshmixer

Using the sculpting tool in Meshmixer.

Sculpting is the most dynamic of the smoothing tools we will cover in this tutorial, and provides immediate feedback on the reshaped geometry. It allows users to point at a particular region of their geometry, and remodels it according to the currently selected sculpting brush. This is usually very convenient to do some quick touch ups to a small region that you noticed looks off. To activate the sculpting tool, selecting the “Sculpt” button from the left-hand menu, which has a chisel icon. You should see the sculpting options menu appear on the left hand side of the screen, which in the left most figure above. At the top you should see a slider that goes between “Volume” and “Surface”, which you should switch to “Surface” for this tutorial. After selecting “Surface” and moving your cursor over your geometry, you should notice a circle appear where your cursor is, and the circle should hug the surface of your geometry. If you left-click, the sculpting tool will remodel your geometry inside the circle according to your selected brush and parameters.

Below that is a button for “Brushes”, which allows you to select the sculpting shape that you would like to apply. You are of course free to experiment with the other brushes, but we have found the “BubbleSmooth” brush to work the best for smoothing sharp junctions in cardiovascular models. Below in the Properties window, you should see sliders for “Strength” and “Size”, two important parameters. Strength will determine how aggressively the sculpting tool will remodel your geometry, with higher numbers meaning more aggressive. When using the tool for the first time, we recommend using a low value of Strength, like 10, to see how it works. The Size parameter sets the radius of the sculpting circle. For smoothing the area between bifurcations, we recommend setting the Size such that the diameter of the sculpting circule is roughly the size of the space of the bifurcation.

b. Local remeshing in Meshmixer

Using the remeshing tool in Meshmixer.

Local remeshing will reconfigure the triangle elements in a given region such that the triangles are all have roughly equal size, shape, and aspect ratio.These triangles will generally be higher quality than the initial discretization that comes out of SimVascular, especially at the bifurcations.To remesh a region, we first have to “Select” it. Click on the “Select” button from the left hand menu, which should change your cursor into a selection tool that you can left-click on your model to start selecting regions. Before we do, it is useful to change the select type and size to be convenient for your model. The “Sphere Brush” will allow you to select all the cells inside a sphere, while the “Unwrap Brush” will select cells in a flat circule that hugs your surface. The “Size” parameter changes the size of the sphere or unwrap brush and should be set to a value that is convenient for the region you are remeshing. Once you have selected a brush and size, left-click on the regions in your model you want to remesh, like all the elements around a bifurcation, for example. After all the elements have been selected, hit R on your keyboard to initiate remeshing. Alternatively, you can click the “Remesh” button that appears from the “Edit…” menu after you have selected a region.

For larger sized regions, it might a few seconds to a few minutes to remesh. After the initial remeshing is complete, the remeshing properties menu should appear on the left hand side of the screen that will allow you to select various properties. Two convenient properties to set are the remeshing type, and the density/edge size. The first drop down menu in the remeshing properties menu allows you to select the remeshing type, of which “Adaptive Density” and “Target Edge Length” are common selections. Setting the mode to “Adaptive Density” will cause the remeshing edge size to be selected adaptively based on local geometrical features. The amount of elements triangles packed into the remeshing can be adjusted by adjusting the “Density” slider, which higher densities resulting in a higher number of elements. More elements typically results in a more faithful representation of your geometry, at an added cost of requiring more memory. Setting the mode to “Target Edge Length” will allow users to adjust the number of elements based on the maximum Edge Length slider. Smaller edge lengths will result in more elements.

Once you get a remesh you are satisfied with, click the “Accept” button at the bottom of the remeshing properties. This will commit your remesh to the model and take you back to the selection tool with your current region still selected.

c. Local Smoothing in Meshmixer

Using the smoothing tool in Meshmixer.

Local smoothing will attempt to smooth over any rough bumps or sharp corners automatically. This is usually useful to perform after remeshing a region to make sure the resulting remesh is smooth. Smoothing is achieved via the “Select” tool, just like remeshing. After selecting the region, press Ctrl + F to smooth the region. Alternatively, you can click “Deform…” then “Smooth” from the Select menu.

After the initial smoothing is complete, the smoothing properties menu will appear. Most of the time, the default settings for the smooth will provide adequate smoothing, though you are of course free to experiment and explore.

  1. Extrude the fluid domain

After the fluid domain surface has been sufficiently smoothed at all sharp corners and bifurcations, we are now ready to extrude out the wall surface.

a. Select all (e.g. by choosing select and then hit Ctrl+a / Cmd+a) b. Choose Edit > Extrude

c. Choose your offset and the offset direction “normal”. This will extrude the wall in normal direction by your given offset. The offset will correspond to the thickness of the vessel wall. Typically, the thickness of the vessel wall is determined by the size of the vessels you are modeling. This information can often be found in the literature, but in the absense of such data, setting the thickness to 10\% of the mean vessel diameter is a reasonable first approximation. If you prefer to have a variable thickness along the vessel, simply restrict your selection to the region you want to extrude by a given thickness. If you like the offset, hit “Accept”.

d. If you do a plane cut of your model now it should look something like this:

In green you can see the surface of the fluid, and in purple the extruded surface.

  1. Clean up the extruded surface

Even though we cleaned up the surface prior to extrusion, the extrusion will undoubtedly produce more intersecting, small, and poor quality elements. Use the sculpting, remeshing, and smoothing tools like we did earlier to clean up the model at the sharp corners and bifurcations again. One note when using the remeshing and smoothing tools: since there will be an interior fluid mesh along with the exterior surface mesh, it is important to select both when using the “Select” tool. In this scenario, it is often best to use the “Sphere Brush” when using the “Select” tool to get all the interior elements as well.

We will now outline a general recipe for extruding a fluid geometry to produce a structural geometry using the tools we just described. This is not the only way to do this, and if your geometry contains many complex features you will probably have to iterate several times to get an adequate structural geometry that can be used for FSI simulations.

  • After loading in the .stl for the fluid geometry, smooth the geometry at all the bifurcations by:
    • Use the “Sphere Brush” in the Select tool with a large enough size to select the entire region around a bifurcation.
    • Press “R” to remesh the region, select an appropriate edge size that captures the geometry well, then Accept.
    • While the region is still selected, Smooth it using Ctrl + F.
    • Use the “Sculpt” tool to further clean up the bifurcation to make it as smooth as possible.
    • Repeat for all other bifurcations.
  • Extrude your model. Press Ctrl + A and extrude to use a uniform thickness or select individual parts of your model and extrude seperately for variable thickness.
  • Smooth the post-extruded bifurcations using the same procedure as before.
    • Use the “Sphere Brush” in the Select tool with a large enough size to select the entire region around a bifurcation.
    • Press “R” to remesh the region, select an appropriate edge size that captures the geometry well, then Accept.
    • While the region is still selected, Smooth it using Ctrl + F.
    • Use the “Sculpt” tool to further clean up the bifurcation to make it as smooth as possible.
    • Repeat for all other bifurcations.

Below is an example of how a bifurcation looks before and after applying the above smoothing procedure. This example did not require much manual sculpting. But some bifurcations that produce sharp corners might require additional sculpting to ensure a smooth surface.

Smoothing bifurcations.

The solid model for the structural domain is finished. Export it as an stl file.

Create the solid domain mesh

Similar to the way we create the mesh for the fluid domain in SimVascular (see: http://simvascular.github.io/docsMeshing.html) we can now create the mesh for the solid domain.

  1. Import the .stl of the structural domain geometry into SimVascular

Right-click on the “Models” button and choose “Import solid Model”, then select the .stl model you just exported from Meshmixer:

Importing .stl model into SimVascular.
  1. Extract faces, rename them if wanted

Extracting faces allows users to reference specific parts of the exterior of the model for applying boundary conditions. It is important to have separate faces that represent each of the inlet and outlet faces, the inner wall, and the outer wall.

To extract faces from your .stl, double-click on your newly imported model in the “Models” section of the left-hand menu to bring up the SV Modeling window. Then, click “Face Ops” and “Extract”. This should reveal a dialogue box asking for a “Sep. Angle”, which is a parameter that determines how aggressively to search for faces. SimVascular will check the direction of the outward facing normals in each element of the surface, and if it finds the difference in the angle between two normals to be greater than the threshold, it will assign them different faces. We recommend a value of 70 for the Sep. Angle., but this depends on the geometry of your model. After extracting faces, you should see a list of “noname” faces appear on the “Face List”. Ideally, you should have a face for each of the inlets and outlets and the inner wall and outer wall, for a total of (number of inlets and outlet) + 2. If you click a noname face from the list, it will highlight in yellow on your 3D model, and the surface area of the face will be displayed in the bottom left corner of the SimVascular Window.

Now is a good time to click each face on the list so you can confirm its location on the model, then rename it appropriately for your convenience. An example naming convention could be to call the inner wall “wall_inner”, the outer wall “wall_outer”, and all of the inlets and outlets as “perimeter_X”, where X is the name of the particular vessel.

This step can be difficult as Meshmixer often produces poor quality triangles that do not extract correctly in SimVascular. You will sometimes see many more faces in the Face List than you expected. If some or all faces are not correctly detected (number of faces is not correct or not correctly distinguished between faces) you can try to use a different separation angle and/or use the merge tool in SimVascular, or improve the affected areas in meshmixer, e.g. by smoothing, sculpting or remeshing them.

This situation most likely happens due to poor quality triangles being produced on the wall surfaces from either the extrusion or remeshing. This can be addressed with the surface treatment tools within SimVascular. Use Ctrl + click to select all the surfaces in your structural model EXCEPT for the perimeter surfaces that surround your caps. Then, under the “Face Ops” tab, click the “Combine” button to combine all of these faces into one. This will allow us to apply surface treatment on all of the walls and poor triangles. We now follow the usual procedure for smoothing surfaces. Now, click “Local Ops”, then “Decimation”, then click “Decimate Local” to decimate some triangles on that surface. Thenm click “Cstr. Smooth”, then “Smooth Local” to smooth the wall surfaces. This will hopefully eliminate all of the poor quality triangles. We now need to re-extract the inner and outer wall surfaces into separate faces. Go back to “Face Ops” then “Extract”. Select an extraction angle (70 degrees is usually good) and click “Extract Faces”. Hopefully now when you do this, you should see a much smaller list of faces. If the number of faces equals the number of inlets and outlets + 2, then you are in great shape.

You will want separate faces for each of the inlets and outlet perimeters, as well as a separate face for the inner wall and a separate face for the outer wall. This is so that we can apply boundary conditions on all of these faces. Clicking on each face in your Face List should highlight their locations in your structural model in yellow. Do this for each face to confirm their location and to confirm that all faces are accounted for.

Extracting faces in SimVascular.
  1. Create the structural mesh

Right-click on “Meshes” in the Data Manager Window and select “Create mesh”, then choose the structural domain model. If you double click on the newly created mesh item, the “mesh interactor” should open up. If the surface mesh of your structural model from Meshmixer was of a good spatial resolution, and quality, you may want to toggle off the option “Surface remeshing” (this will also speed up the meshing step, especially for larger models).

Creating volumetric solid mesh.
  1. Run the mesher

This will create your solid domain mesh. Export the mesh by selecting “Export Mesh-Complete”. This will save a folder in the location of your choice that has all the mesh information. Inside this folder you will find a file called “mesh-complete.mesh.vtu” that contain the full volumetric mesh information of your model. You should also notice a folder called “mesh-surfaces” that has separate .vtp files for each of the extracted faces. Both the mesh-complete.mesh.vtu and all the .vtp mesh-surfaces are viewable in the free software Paraview.

Create the fluid domain mesh

When making the mesh for the fluid domain, we need to ensure that the nodes at the interface between the structural and fluid domains coincide. We will thus use the surface mesh of the inner wall of the solid domain that we just created in the previous step. You will find it in the exported mesh in the subdirectory “mesh-surfaces” of your exported solid mesh. If you used the suggested naming convention from the previous step, this inner wall will be named “wall_inner”.

  1. Import the interfacial mesh of the solid domain by right-clicking “Models” in the SV Data Manager, then click “Import Solid Model”. Select the surface representing the interface between the structural and fluid domains (i.e. the “inner wall” of your structural domain mesh). Rename it as you wish.

  2. Extract the face of the solid domain interface. You should get a single face which shows the entire wall. Rename it to something convenient like “wall_lumen”.

  3. Fill holes with ID’s.

We now need to fill the inlet and outlet holes and name them so we can apply boundary conditions. Under the “SV Modeling” window, click the “Face Ops” button, then click the newly revealed “Fill Holes w. IDs”. This should place “noname” faces at each inlet and outlet of your fluid domain. Rename them to something convenient to remember, such as by the name of the vessel they are representing.

  1. Remesh the inlets and outlets

Filling the inlet and outlet holes produced geometry on those surfaces, but it does not have a mesh appropriate for computational modeling. We thus need to remesh them. Make sure you are under the “Face Ops” tab by clicking “Face Ops” on the right side of the screen. Clicking the blue “Remesh” button should reveal the Remeshing dialogue, which asks you to select an appropriate size of the elements. The element size you select should reflect the kind of accuracy and cost you would like to use in the FSI simulation. A good starting point is to use the same size that you used for the structural domain solid mesh, though they do not have to be equal. Remesh each of the inlet and outlet faces but DO NOT remesh the wall interface. This will ensure that the nodes still coincide with the structural domain nodes after this remeshing.

Remeshing inlet and outlets of fluid domain.
  1. Create the fluid domain mesh.

In the SimVascular DataManager right click “Meshes” and create a new mesh from the fluid model. This time, it is REQUIRED to toggle off the option “Surface mesh”. This is important, because it will keep the surface nodes exactly where they are. This way the nodes on the fluid interface will coincide with the nodes on structural interface. Choose your refinement to be consistent with the size of the elements you chose for the wall and caps, then run the mesher. Once it finishes, right-click your fluid mesh and click “Export Mesh-Complete”. You are done!

Creating fluid domain mesh.

Restarting a simulation after re-meshing

If a simulation is to be restarted after re-meshing has occurred, the new mesh must be added to the input file. SVFSI will create a directory called n-procs, where n is the number of MPI tasks for the simulation, and we will refer to this as the run directory.

  1. In the run directory, find the stFile_n.bin. Make a symbolic link:

        ln -s stFile_n.bin stFile_last.bin
    

    Note that if the previous job completed smoothly, this file may already exist and no link is required to be made.

  2. If remeshing has occurred, the run directory will contain a hidden directory called “.remesh_tmp.dir”, which contains new mesh files that should be named according according to their initial name and the timestep of the last remesh. For example, lumen_21.vtu.

  3. Change the input.mfs file to set a restarted simulation and for the new mesh names.

    1. Set Continue previous simulation: 1
    2. All new paths and names for mesh files should be updated, including paths. For example, change to 120-procs/.remesh_tmp.dir/lumen_21.vtu.
    3. Set Check IEN order: F.
    4. Update the final time if desired.
  4. Run the simulation again.

Simulations with prescribed wall motion

In this section, we discuss how to set up a simulation with prescribed wall motion. For example, one may wish to extract the motion of the walls of the ventricle from scans, then reconstruct the flow fields by solving the Navier Stokes equations. This is accomplished by computing a displacement field for the prescribed motion of the mesh. This must be done offline by the user for their specific problem. This is typically done for a small subset of the total times, and displacements between the specified times are prescribed by using a piecewise linear interpolant.

For the wall motion file, the file format is as follows. First, specify the dimension of the problem (three), the number of times at which to specify the displacements, and the number of vertices in the moving wall mesh. Then specify the times at which the displacements occur. Next, for each vertex, specify its index and then the prescribed displacements for each time. Note that in the case of multiple moving faces, these numbers may not start at one for any given face, as indexing is global. If a vertex is on an edge between two faces, it should have the same index and displacement fields, specified redundantly in both files.

    Dimension n_times m_vertices
    t_1
    t_2
    ...
    t_n
    vertex_1_index
    displacement_1_vertex_1
    displacement_2_vertex_1
    ...
    displacement_n_vertex_1
    vertex_2_index 
    displacement_1_vertex_2
    displacement_2_vertex_2
    ...
    displacement_n_vertex_2
    ... 

For example,

    3   21   14907
    0.000000
    3.800E-2
    7.600E-2
    1.140E-1
    1.520E-1
    1.900E-1
    2.280E-1
    2.660E-1
    3.040E-1
    3.419E-1
    3.800E-1
    4.180E-1
    4.560E-1
    4.940E-1
    5.320E-1
    5.699E-1
    6.080E-1
    6.460E-1
    6.839E-1
    7.220E-1
    7.600E-1
    1
    0.000000   0.000000   0.000000   
    2.800E-1   -8.99E-2   -1.00E-2   
    8.799E-1   -2.09E-1   -8.10E-1   
    1.339999   -1.59E-1   -1.43000   
    1.509999   7.000E-2   -1.59000   
    1.310000   3.100E-1   -1.52000   
    1.069999   5.100E-1   -1.41000   
    1.009999   6.000E-1   -1.28000   
    9.699E-1   3.600E-1   -1.10000   
    1.049999   -1.09E-1   -8.80E-1   
    1.169999   -6.69E-1   -7.60E-1   
    1.169999   -1.17999   -7.90E-1   
    1.129999   -1.31999   -9.69E-1   
    1.049999   -1.25999   -1.34000   
    1.000000   -9.39E-1   -1.64000   
    9.899E-1   -4.29E-1   -1.96000   
    1.000000   3.000E-2   -2.15000   
    9.299E-1   3.600E-1   -1.90000   
    6.599E-1   4.200E-1   -1.14000   
    1.999E-1   1.500E-1   -3.50E-1   
    0.000000   0.000000   0.000000   
    2
    0.000000   0.000000   0.000000   
    2.700E-1   -8.99E-2   0.000000   
    8.500E-1   -1.99E-1   -6.80E-1   
    1.280000   -1.59E-1   -1.16999   
    1.430000   5.000E-2   -1.28000   
    1.230000   2.599E-1   -1.19999   
    1.010000   4.499E-1   -1.09000   
    9.500E-1   5.099E-1   -9.69E-1   
    9.100E-1   2.400E-1   -8.19E-1   
    1.000000   -2.59E-1   -6.70E-1   
    1.120000   -8.29E-1   -5.69E-1   
    1.130000   -1.31999   -5.89E-1   
    1.090000   -1.44999   -7.50E-1   
    1.000000   -1.36000   -1.06999   
    9.400E-1   -1.00999   -1.31000   
    9.300E-1   -4.69E-1   -1.56999   
    9.500E-1   0.000000   -1.72999   
    8.900E-1   3.299E-1   -1.53999   
    6.200E-1   3.999E-1   -9.39E-1   
    1.900E-1   1.399E-1   -2.90E-1   
    0.000000   0.000000   0.000000  
    ... 

Note that in this example, we wish the mesh motion to be periodic in time, and thus the final displacement is zero.

To set up the input file, set the equation to be FSI, even though there is not necessarily a structure, to use the ALE meshes rather than static meshes. For the moving wall, add the motion file when specifying the wall boundary condition, and turn on the option “Impose on state variable integral.”

    Add BC: moving_wall {
      Type: Dirichlet 
      Time dependence: General
      Temporal and spatial values file path: wall_motion.dat
      Profile: Flat
      Zero out perimeter: 1
      Impose flux: 0
      #---------------- Add this line to the moving boundary face -----------
      Impose on state variable integral: 1
    }

It is recommended to include remeshing if the wall motion is such that the domain undergoes large changes. For example, set

    Remesher: Tetgen {
      Max edge size: lumen { val: 3.0 }
      Min dihedral angle: 10.0
      Max radius ratio: 1.1
      Remesh frequency: 100
      Frequency for copying data: 1  
    }

The max edge size should be consistent with the original mesh size.

Under the mesh equation, we similarly add the motion file.

    Add equation: mesh {
       Coupled: 1
       Min iterations: 1
       Max iterations: 8
       Tolerance: 1e-3
       Residual dB reduction: -20
       Poisson ratio: 0
       Output: Spatial {
          Displacement: t
       }

       #---------- Add the BC for the moving_wall to the mesh equation as well ------
       Add BC: moving_wall {
          Type: Dirichlet 
          Time dependence: General
          Temporal and spatial values file path: wall_motion.dat
          Profile: Flat
          Zero out perimeter: 1
          Impose flux: 0
          #---------------- Add this line to the moving boundary face -----------
          Impose on state variable integral: 1
       }