Merge pull request #125 from KVSlab/dynamic-meshing

Add meshing for moving domains and MovingAtrium.py

docs/_toc.yml

@@ -15,6 +15,7 @@ parts:
   - caption: Using VaMPy
     chapters:
     - file: preprocess
+    - file: movingpreprocess
     - file: simulation
     - file: postprocess
     - file: quantities
@@ -26,6 +27,7 @@ parts:
     - file: tutorials
     - file: artery
     - file: atrium
+    - file: movingatrium
 
   - caption: Python API
     chapters:

docs/api.rst

@@ -15,6 +15,11 @@ Pre-processing scripts
    :undoc-members:
    :show-inheritance:
 
+.. automodule:: vampy.automatedPreprocessing.moving_common
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 .. automodule:: vampy.automatedPreprocessing.DisplayData
    :members:
    :undoc-members:
@@ -96,6 +101,11 @@ CFD simulation scripts
    :undoc-members:
    :show-inheritance:
 
+.. automodule:: vampy.simulation.MovingAtrium
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 .. automodule:: vampy.simulation.Womersley
    :members:
    :undoc-members:

docs/artery.md

@@ -72,9 +72,9 @@ The centerline inside the surface model, colored by it\'s
 maximum inscribed sphere diameter.
 ```
 
-Meshing based on the centerline diameter is performed by supplying the
-`-m diameter` flag, and optionally the `--coarsening-factor` flag, as above. To generate a variable density mesh based on
-the centerline diameter, run the following command:
+Meshing based on the centerline diameter is performed by supplying the `-m diameter` flag, and optionally
+the `--coarsening-factor` flag, as above. To generate a variable density mesh based on the centerline diameter, run the
+following command:
 
 ``` console
 $ vampy-mesh -m diameter -c 1.2 -i C0001/model.vtp

docs/atrium.md

@@ -10,9 +10,8 @@ that VaMPy is also applicable to other vascular domains, not only tubular struct
 
 ```{attention} 
 Because VaMPy relies on `vtkPolyData` as input, the `.vtk` model needs
-to be converted to `.vtp` (or `.stl`) format, which can quicly be done in ParaView
-by using the `Extract Surface` filter, and saving the data as
-`LA_Endo_5.vtp`.
+to be converted to `.vtp` (or `.stl`) format, which can easily be performed in ParaView
+by using the `Extract Surface` filter, and saving the data as `LA_Endo_5.vtp`.
 ```
 
 ## Meshing an atrium with appendage refinement

docs/index.md

@@ -3,10 +3,12 @@
 The [Vascular Modeling Pypeline](https://github.com/KVSlab/vampy) (VaMPy) is a collection of automated scripts to
 prepare, run, and analyze vascular morphologies. This includes pre-processing scripts for generating a volumetric mesh,
 defining physiological boundary conditions, and inserting probes for sampling velocity and pressure. For the
-computational fluid dynamics (CFD) simulation, we have included two [Oasis](https://github.com/mikaem/Oasis) problem
-files for simulating flow in the [internal carotid artery](https://en.wikipedia.org/wiki/Internal_carotid_artery) and
-the [left atrium](https://en.wikipedia.org/wiki/Atrium_(heart)). There are also a variety of post-processing scripts,
-which computes wall shear stress-based metrics, more advanced turbulence metrics, and a variety of geometric parameters.
+computational fluid dynamics (CFD) simulation, we have included
+three [Oasis](https://github.com/mikaem/Oasis)/[OasisMove](https://github.com/KVSlab/OasisMove) problem files for
+simulating flow in the [internal carotid artery](https://en.wikipedia.org/wiki/Internal_carotid_artery) and
+the [left atrium](https://en.wikipedia.org/wiki/Atrium_(heart)), including moving domain simulations in the atrium.
+There are also a variety of post-processing scripts, which computes wall shear stress-based metrics, more advanced
+turbulence metrics, and a variety of geometric parameters.
 
 ## Contents
 

docs/installation.md

@@ -2,8 +2,9 @@
 
 `VaMPy` is a pure Python package that combines the morphological processing tools of
 [morphMan](https://github.com/KVSlab/morphMan) and [vmtk](http://www.vmtk.org/), the computational fluid dynamics
-solver [Oasis](https://github.com/mikaem/Oasis), and the finite element computing
-platform [FEniCS](https://fenicsproject.org/) into a powerful and user-friendly computational fluid dynamics pipeline.
+solver [Oasis](https://github.com/mikaem/Oasis)/[OasisMove](https://github.com/KVSlab/OasisMove), and the finite element
+computing platform [FEniCS](https://fenicsproject.org/) into a powerful and user-friendly computational fluid dynamics
+pipeline.
 
 ## Dependencies
 
@@ -13,15 +14,14 @@ The dependencies of `VaMPy` are:
 * FEniCS >= 2018.1
 * morphMan >= 1.2
 * vmtk >= 1.4.0
-* Oasis >= 2018.1
+* Oasis >= 2018.1 OR OasisMove >= 0.2.0
 * paramiko >= 3.0
 * cppimport >= 22.8.2
 
 ## Installing with `conda`
 
-To install `VaMPy` and all its dependencies to a *local environment*, we recommend using `conda`.
-Instructions for installing `VaMPy`
-with `conda` can be found [here](install:conda).
+To install `VaMPy` and all its dependencies to a *local environment*, we recommend using `conda`. Instructions for
+installing `VaMPy` with `conda` can be found [here](install:conda).
 
 ## Installing with Docker
 

docs/movingatrium.md

@@ -0,0 +1,122 @@
+(tutorial:movingatrium)=
+
+# Moving domain simulation of the left atrium
+
+The third tutorial considers a moving domain CFD simulation in a simplified left atrium geometry. This hand-crafted
+geometry, named `model.vtp`, can be found in the `models/moving_atrium` directory and was designed by the `VaMPy`
+authors using the [MeshMixer](https://meshmixer.com/) software. This tutorial aims to showcase `VaMPy`'s versatility,
+highlighting its capabilities beyond rigid wall simulations, particularly in the context of atrial simulations.
+Additionally, the tutorial introduces the concept of blood residence time ($T_R$) – a metric that captures blood flow
+stasis by solving a scalar transport equation.
+
+## Moving domain meshing of the left atrium
+
+The morphology of the simplified left atrium is depicted on the left side of {numref}`simple_atrium`. It features four
+pulmonary veins and a pouch-like structure representing the left atrial appendage. The left atrial appendage is
+particularly prone to blood clot formation, making it a significant region of interest within the left atrium. However,
+this tutorial focuses on showcasing the moving domain capabilities of `VaMPy`/`OasisMove`, so mesh refinement will not
+be discussed.
+
+To perform moving domain meshing for this model, execute the command below:
+
+```console
+$ vampy-mesh --moving-mesh -i models/moving_atrium/model.vtp -m constant -el 1.3 -fli 2 -flo 1.5 -at 
+```
+
+In this command, the `-fli` and `-flo` flags determine the length of the flow extensions at the inlets and outlet,
+respectively, and the `-at` flag is used to notify the pipeline that an atrium model is being meshed.
+
+Upon successful execution, the meshing command will generate the standard output files. Additionally, it will produce
+the displacement matrix, stored in `model_points.np`, and the deformed surfaces with flow extensions in `model_extended`
+. The latter is illustrated in the middle figure of {numref}`simple_atrium`. The final volumetric mesh, inclusive of
+four boundary layers, is displayed on the right side of {numref}`simple_atrium`.
+
+```{figure} figures/simpleatrium.gif
+---
+name: simple_atrium
+---
+From left to right: the surface model used in this tutorial, the deformed model with flow extensions, and the resulting volumetric mesh based on the provided meshing parameters.
+```
+
+## The displacement field and `model_points.np`
+
+The `model_points.np` file introduces a new input for the CFD simulation, describing the displacement field of each
+surface point. Each point within `model_points.np` is accomponied by its discrete displacement over a full cycle. For
+moving domain simulations driven by image-based motion, `OasisMove` utilizes this displacement field to define the wall
+boundary conditions. The process within `OasisMove` is as follows:
+
+1. **Load displacement matrix**: `OasisMove` identifies and loads the `model_points.np` file, storing the matrix in
+   the `points` variable.
+2. **Spline interpolation**: Using `SciPy`'s `splrep` method, spline interpolation is executed for each point, ensuring
+   a continuous description of the deformation over one cycle (see the left-hand side of {numref}`spline`).
+3. **Mesh velocity computation**: The first-order derivative of the spline is computed by `OasisMove`. This derivative
+   represents the mesh velocity and is set as the wall boundary condition (see the right-hand side of {numref}`spline`).
+
+To exemplify, we may consider the deformation field of the demo model located in the `moving_atrium` directory. Focusing
+on the surface point at index 1, the displacement curve, its spline-interpolated counterpart, and the mesh velocity for
+the $x$-component are illustrated in {numref}`spline`.
+
+```{figure} figures/spline.png
+---
+name: spline
+---
+On the left: the discrete displacement profile (x) for the surface point at index 1, together with its spline representation (blue curve). On the right: the first-order derivative, denoting the mesh velocity for that specific point.
+```
+
+---
+
+## Quantifying blood residence time
+
+In the majority of patients with atrial fibrillation, there's a risk of blood clot formation in the left atrial
+appendage. To further investigate this phenomenon using CFD simulations, it's essential to assess the duration a blood
+particle remains within the left atrium and its appendage. This duration can be quantified by a scalar metric termed the
+blood residence time, $T_R$.
+
+The blood residence time, $T_R$, is determined throughout the computational domain by solving a forced passive scalar
+equation, as formulated by Rossini et al.{cite}`rossini2016clinical`:
+
+```{math}
+:label: transport
+\frac{\partial T_R}{\partial t} + \mathbf u \cdot \nabla T_R = 1,
+```
+
+At the onset of the simulation, $T_R$ is initialized to zero everywhere in the domain, including at the pulmonary veins
+as a boundary condition. Within `OasisMove`, Equation {eq}`transport` is simulteneously solved as the Navier-Stokes
+equations, and incorperates a SUPG stabilization term to counteract the inherent instability of a purely advective
+transport equation.
+
+---
+
+## Moving domain CFD simulation
+
+With the generated mesh, standard CFD input files, and the displacement field from `model_points.np`, we can start the
+moving domain CFD simulation. It's important to note that such simulations are exclusively supported
+by [`OasisMove`](https://github.com/KVSlab/OasisMove). The specifics of our CFD problem are outlined in
+the [`MovingAtrium.py`](https://github.com/KVSlab/VaMPy/blob/master/src/vampy/simulation/MovingAtrium.py) file, found in
+the `simulation` directory. The inlet boundary conditions are defined by the normalized flow rate values located in
+the  `PV_values` file, and are adjusted by the model's mean flow rate. Meanwhile, the outlet is set with a zero-pressure
+boundary condition, and the wall boundary conditions are determined using the displacement field from `model_points.np`.
+
+To execute a moving domain CFD simulation of the left atrium model over four cardiac cycles, navigate to
+the `src/vampy/simulation` directory and run the following command:
+
+``` console
+$ oasismove NSfracStepMove problem=MovingAtrium mesh_path=../../../models/moving_atrium/model.xml.gz number_of_cycles=4 T=1000 dt=1
+```
+
+Upon completion of the simulation, the results, along with the associated mesh, are stored in a concise HDF5 format.
+Additionally, the residence time $T_R$ and the velocity field are saved in the XDMF files `blood.xdmf`
+and `velocity.xdmf`, respectively, which includes the mesh deformation. {numref}`moving_results` presents the volumetric
+representation of the blood residence time and velocity field throughout the four cardiac cycles. It's noteworthy to
+observe the accumulating $T_R$ within the left atrial appendage over time.
+
+<figure>
+<video controls style="width: 100%; height: auto;">
+    <source src="_static/videos/moving_results.mp4" type="video/mp4">
+</video>
+    <figcaption>From left to right: the volumetric depiction of the blood residence time, and corresponding velocity field over four cardiac cycles.</figcaption>
+</figure>
+
+```{bibliography}
+:filter: docname in docnames
+```

docs/movingpreprocess.md

@@ -0,0 +1,88 @@
+# Pre-processing for moving domains
+
+## Meshing for moving domain simulations
+
+In vascular modeling, understanding the dynamics of blood flow often needs the inclusion of vascular deformation. To
+address this, `VaMPy` incorporates a dynamic mesh or moving domain meshing pipeline. This feature enables automated
+pre-processing for both rigid and moving domains. Activating the moving domain meshing is straightforward and can be
+done through specific command line arguments. However, there are certain input requirements to be aware of.
+
+For moving domain meshing in `VaMPy`, users are expected to provide an input surface model, such as `model.vtp`, as they
+would typically do. In addition to this, users need to supply a series of deformed surface models. It's crucial that
+these deformed models maintain the same number of surface nodes/points as the original model, ensuring a one-to-one
+mapping between each model. As an illustrative example using the surface file `model.vtp`, the deformed surface models
+should be stored in a directory named `model_moved`. Each of these deformed models should be named in the
+format `model_moved_XX.vtp`, where `XX` represents an incrementing index. This would give us the following file
+structure:
+
+```
+moving_atrium
+├── model_moved
+│   ├── model_moved_01.vtp
+│   ├── model_moved_02.vtp
+│   ├── |
+│   ├── |
+│   ├── |
+│   └── model_moved_20.vtp
+└── model.vtp
+```
+
+To perform moving domain meshing, we add the `--moving-mesh` (`-mm`) flag:
+
+``` console
+$ vampy-mesh -i models/moving_atrium/model.vtp --moving-mesh -at 
+```
+
+Upon successful meshing, two new items will appear compared to the rigid domain meshing: the `model_points.np` file and
+the `model_extended` folder. This can be visualized in the subsequent file structure:
+
+```
+moving_atrium
+├── model_extended
+│   ├── model_moved_01.vtp
+│   ├── model_moved_02.vtp
+│   ├── |
+│   └── model_moved_20.vtp
+├── model_moved
+│   ├── model_moved_01.vtp
+│   ├── model_moved_02.vtp
+│   ├── |
+│   └── model_moved_20.vtp
+├── model.xml.gz
+├── model.vtu
+├── model_points.np
+├── model_probe_point.json
+├── model_info.json
+└── model.vtp
+```
+
+The `model_extended` folder contains the same amount of surface models as those in `model_moved`. However, these models
+come with the cylindrical flow extensions. For visual validation, you can use software like ParaView to ensure their
+suitability for simulation. Meanwhile, the `model_points.np` file captures the displacement field by tracking the
+movement of each point on the deformed surfaces. This file is crucial for CFD simulations, used to prescribe the wall
+boundary condition.
+
+## Clamping boundaries
+
+In case the input's deformed surface models exhibit deformations at the inlet and outlet boundaries, it's essential to
+adjust the flow extensions accordingly. To address this, we introduced the `--clamp-boundaries` (`-cl`) command line
+argument. This ensures that the original inlets and outlet boundaries maintain their image/surface-based motion, while
+the displacement is reduced towards the flow extension ends, which remain stationary or "clamped" in space.
+
+Two displacement reduction profiles are available: linear and sinusoidal. By default, the linear profile is applied.
+However, if you prefer the sinusoidal reduction, you can modify `moving_common.py` by replacing `profile="linear"`
+with `profile="sine"`. A side-by-side comparison of these profiles is illustrated in {numref}`clamp`, with the linear
+profile on the left and the sinusoidal on the right.
+
+```{figure} figures/clamp.png
+---
+name: clamp
+---
+On the left: A linear reduction between the image-based boundary and the end of the flow extension. On the right: A 
+sinusiodal profile is used to gradually reduce the movement between the image-based boundary and the end of the flow extension.
+```
+
+If the `--clamp boundaries` argument isn't specified, the deformation at the flow extension ends will mirror the
+boundary deformations observed in the image/surface-based motion.
+
+

docs/references.bib

@@ -10,4 +10,15 @@ @article{roney2021constructing
   pages={233--250},
   year={2021},
   publisher={Springer}
+}
+
+@article{rossini2016clinical,
+  title={A clinical method for mapping and quantifying blood stasis in the left ventricle},
+  author={Rossini, Lorenzo and Martinez-Legazpi, Pablo and Vu, Vi and Fernandez-Friera, Leticia and Del Villar, Candelas P{\'e}rez and Rodriguez-Lopez, Sara and Benito, Yolanda and Borja, Maria-Guadalupe and Pastor-Escuredo, David and Yotti, Raquel and others},
+  journal={Journal of biomechanics},
+  volume={49},
+  number={11},
+  pages={2152--2161},
+  year={2016},
+  publisher={Elsevier}
 }

docs/simulation.md

@@ -3,8 +3,8 @@
 ## Simulations in `Oasis`
 
 Following pre-processing, the next step of using the Vascular Modeling Pypeline is performing the computational fluid
-dynamics (CFD) simulations with `oasis`. Assuming `oasis` has been installed, start by navigating to the `simulation`
-folder:
+dynamics (CFD) simulations with [`Oasis`](https://github.com/mikaem/Oasis). Assuming `oasis` has been installed, start
+by navigating to the `simulation` folder:
 
 ``` console
 $ cd src/vampy/simulation
@@ -20,12 +20,28 @@ $ oasis NSfracStep problem=Artery mesh_path=../../../models/artery/artery.xml.gz
 Running the simulations will create the result folder `results_artery` (specific to the `Artery.py` problem) located
 inside `src/vampy/simulation`, with the results and corresponding mesh saved compactly in HDF5 format.
 
+## Simulations in `OasisMove`
+
+In case you decide to use [`OasisMove`](https://github.com/KVSlab/OasisMove) for CFD simulations, the command for
+running a (rigid wall) problem file is very similar to `Oasis`:
+
+``` console
+$ oasismove NSfracStepMove problem=Artery dynamic_mesh=False mesh_path=../../../models/artery/artery.xml.gz save_solution_after_cycle=0
+```
+
+For moving domain simulations (`MovingAtrium.py`), set `dynamic_mesh=True` and continue as usual:
+
+``` console
+$ oasismove NSfracStepMove problem=MovingAtrium dynamic_mesh=True mesh_path=[PATH_TO_ATRIUM_MODEL] 
+```
+
 ## Running parallel computational fluid dynamics in `Oasis`
 
-Oasis runs with [MPI](https://mpi4py.readthedocs.io/en/stable/), and problem files may be executed using MPI commands.
-To use MPI commands you will need to use an MPI interpreter, which is provided with the command `mpirun`, which is part
-of the Python MPI package *mpi4py*. On some systems this command is called `mpiexec` and *mpi4py* seems to include both.
-Assuming *mpi4py* is installed, you may run the CFD simulation in parallel using the `mpirun` command as follows:
+`Oasis`/`OasisMove` runs with [MPI](https://mpi4py.readthedocs.io/en/stable/), and problem files may be executed using
+MPI commands. To use MPI commands you will need to use an MPI interpreter, which is provided with the command `mpirun`,
+which is part of the Python MPI package *mpi4py*. On some systems this command is called `mpiexec` and *mpi4py* seems to
+include both. Assuming *mpi4py* is installed, you may run the CFD simulation in parallel using the `mpirun` command as
+follows:
 
 ``` console
 $ mpirun -np 4 oasis NSfracStep problem=Artery mesh_path=../../../models/artery/artery.xml.gz save_solution_after_cycle=0