diff --git a/docs/_static/videos/moving_results.mp4 b/docs/_static/videos/moving_results.mp4 new file mode 100644 index 00000000..67b1108a Binary files /dev/null and b/docs/_static/videos/moving_results.mp4 differ diff --git a/docs/_toc.yml b/docs/_toc.yml index c2dad6fa..6cde3f0c 100644 --- a/docs/_toc.yml +++ b/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: diff --git a/docs/api.rst b/docs/api.rst index 360b8fd5..360d45b7 100644 --- a/docs/api.rst +++ b/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: diff --git a/docs/artery.md b/docs/artery.md index 5a2221b9..69cc1440 100644 --- a/docs/artery.md +++ b/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 diff --git a/docs/atrium.md b/docs/atrium.md index 45a756ab..0d4b337f 100644 --- a/docs/atrium.md +++ b/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 diff --git a/docs/figures/clamp.png b/docs/figures/clamp.png new file mode 100644 index 00000000..07351f38 Binary files /dev/null and b/docs/figures/clamp.png differ diff --git a/docs/figures/simpleatrium.gif b/docs/figures/simpleatrium.gif new file mode 100644 index 00000000..e5919437 Binary files /dev/null and b/docs/figures/simpleatrium.gif differ diff --git a/docs/figures/spline.png b/docs/figures/spline.png new file mode 100644 index 00000000..3851556e Binary files /dev/null and b/docs/figures/spline.png differ diff --git a/docs/index.md b/docs/index.md index 50f75f9c..50d3c7e5 100644 --- a/docs/index.md +++ b/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 diff --git a/docs/installation.md b/docs/installation.md index c185fd5e..9b353e78 100644 --- a/docs/installation.md +++ b/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 diff --git a/docs/movingatrium.md b/docs/movingatrium.md new file mode 100644 index 00000000..abd0487f --- /dev/null +++ b/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. + +
+ +
From left to right: the volumetric depiction of the blood residence time, and corresponding velocity field over four cardiac cycles.
+
+ +```{bibliography} +:filter: docname in docnames +``` \ No newline at end of file diff --git a/docs/movingpreprocess.md b/docs/movingpreprocess.md new file mode 100644 index 00000000..3920efd2 --- /dev/null +++ b/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. + + diff --git a/docs/references.bib b/docs/references.bib index 474c4a2b..270d2a49 100644 --- a/docs/references.bib +++ b/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} } \ No newline at end of file diff --git a/docs/simulation.md b/docs/simulation.md index ef0d979e..9a6d5a88 100644 --- a/docs/simulation.md +++ b/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 diff --git a/docs/tutorials.md b/docs/tutorials.md index b5c27169..1c5accf9 100644 --- a/docs/tutorials.md +++ b/docs/tutorials.md @@ -1,19 +1,24 @@ # Tutorials These tutorials are meant to guide the user through the basic steps of performing a computational fluid dynamic (CFD) -simulation in a vascular body. The tutorials are divided into two main problems where we will first consider a model of -the [internal carotid artery](https://en.wikipedia.org/wiki/Internal_carotid_artery) (ICA) +simulation in a vascular body. The tutorials are divided into three main problems where we will first consider a model +of the [internal carotid artery](https://en.wikipedia.org/wiki/Internal_carotid_artery) (ICA) with an aneurysm, followed by a model of the [left atrium](https://en.wikipedia.org/wiki/Atrium_(heart)). For the ICA, we will demonstrate different approaches for meshing and pre-processing of the model, followed by a CFD simulation using the Oasis software, before post-processing the results. For the left atrium model, we will investigate local refinement -in a user-defined region, also followed by CFD simulation and post-processing of the results. +in a user-defined region, also followed by CFD simulation and post-processing of the results. In the third tutorial, +we'll explore a moving domain simulation for a simplified left atrium model, and examine blood residence time by +simultaneously solving a scalar transport equation during the simulation. We assume that you have VaMPy installed, meaning you have a working version -of [morphMan](https://github.com/KVSlab/morphMan) and [Oasis](https://github.com/mikaem/Oasis) on your computer. It is -also beneficial, but not necessary, that you are familiar with the [ParaView](https://www.paraview.org/) -visualization software, which we will frequently be using to visualize the geometries and results. +of [morphMan](https://github.com/KVSlab/morphMan) and [Oasis](https://github.com/mikaem/Oasis) on your computer. For the +moving domain tutorial you will need to install [OasisMove](https://github.com/KVSlab/OasisMove), a modified version +of `Oasis` for moving domain simulations. It is also beneficial, but not necessary, that you are familiar with +the [ParaView](https://www.paraview.org/) visualization software, which we will frequently be using to visualize the +geometries and results. ## Link to tutorials - [Internal carotid artery](tutorial:artery) - [Left atrium](tutorial:atrium) +- [Moving domain simulation in the Left atrium](tutorial:movingatrium) diff --git a/environment.yml b/environment.yml index 4de07df6..0c97cb9b 100644 --- a/environment.yml +++ b/environment.yml @@ -11,5 +11,6 @@ dependencies: - python >=3.8 - git - pip: + - scipy - cppimport - git+https://github.com/mikaem/Oasis diff --git a/models/atrium/README.md b/models/atrium/README.md index 905af2da..3e2926e5 100644 --- a/models/atrium/README.md +++ b/models/atrium/README.md @@ -4,10 +4,10 @@ The provided left atrium model is collected from the Left Atrial Segmentation Ch The model meshing and probe points can be computed using automatedPreProcessing with the following command: ``` -python automatedPreProcessing/automatedPreProcessing.py -m diameter -i test/Case_test_atrium/atrium.vtp -c 1.3 +python automatedPreProcessing/automatedPreProcessing.py -m diameter -i models/atrium/atrium.vtp -c 1.3 ``` -Note: Aneurysm workflow is intended for artery models. Therefore, boundary conditions for other models will be wrong. Adapting to other vascular models, such as the atrium, is work in progress. +Note: Boundary conditions for other than vascular models are still work in progress. --- diff --git a/models/moving_atrium/README.md b/models/moving_atrium/README.md new file mode 100644 index 00000000..15725229 --- /dev/null +++ b/models/moving_atrium/README.md @@ -0,0 +1,22 @@ +# Left atrium model for moving domain simulations + +The provided left atrium model (`model.vtp`) is manually created in MeshMixer by [Henrik A. Kjeldsberg](https://github.com/hkjeldsberg/) , and can be meshed +and simulated using VaMPy given that `OasisMove` is installed. + +In addition, the folder `model_moved` contains several displaced surface models, where `model_moved_00.vtp` +corresponds to the input model, `model.vtp`. These surface models are used to create the displacement matrix describing +the motion of the geometry. + +To mesh the model with dynamic/moving meshing, run the following command: + +``` +vampy-mesh -i models/moving_atrium/model.vtp -at -dm -cl -m constant -el 1.3 -bl False -fli 1 -flo 1 +``` + +Then, to run a moving simulation, navigate to `/simulation` and run the command + +``` +oasismove problem=MovingAtrium mesh_path=[PATH_TO_VAMPY]/models/moving_atrium/model.xml.gz +``` + + diff --git a/models/moving_atrium/model.vtp b/models/moving_atrium/model.vtp new file mode 100644 index 00000000..1a9b4718 Binary files /dev/null and b/models/moving_atrium/model.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_00.vtp b/models/moving_atrium/model_moved/model_moved_00.vtp new file mode 100644 index 00000000..1a9b4718 Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_00.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_01.vtp b/models/moving_atrium/model_moved/model_moved_01.vtp new file mode 100644 index 00000000..6957579f Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_01.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_02.vtp b/models/moving_atrium/model_moved/model_moved_02.vtp new file mode 100644 index 00000000..e87c247c Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_02.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_03.vtp b/models/moving_atrium/model_moved/model_moved_03.vtp new file mode 100644 index 00000000..e1b5baca Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_03.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_04.vtp b/models/moving_atrium/model_moved/model_moved_04.vtp new file mode 100644 index 00000000..3a2dc33c Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_04.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_05.vtp b/models/moving_atrium/model_moved/model_moved_05.vtp new file mode 100644 index 00000000..62a68180 Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_05.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_06.vtp b/models/moving_atrium/model_moved/model_moved_06.vtp new file mode 100644 index 00000000..f8571658 Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_06.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_07.vtp b/models/moving_atrium/model_moved/model_moved_07.vtp new file mode 100644 index 00000000..f25de340 Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_07.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_08.vtp b/models/moving_atrium/model_moved/model_moved_08.vtp new file mode 100644 index 00000000..9be3e173 Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_08.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_09.vtp b/models/moving_atrium/model_moved/model_moved_09.vtp new file mode 100644 index 00000000..590c2dfb Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_09.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_10.vtp b/models/moving_atrium/model_moved/model_moved_10.vtp new file mode 100644 index 00000000..913b32eb Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_10.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_11.vtp b/models/moving_atrium/model_moved/model_moved_11.vtp new file mode 100644 index 00000000..06646460 Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_11.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_12.vtp b/models/moving_atrium/model_moved/model_moved_12.vtp new file mode 100644 index 00000000..a7f7b5fd Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_12.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_13.vtp b/models/moving_atrium/model_moved/model_moved_13.vtp new file mode 100644 index 00000000..72aebc14 Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_13.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_14.vtp b/models/moving_atrium/model_moved/model_moved_14.vtp new file mode 100644 index 00000000..98d3e50b Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_14.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_15.vtp b/models/moving_atrium/model_moved/model_moved_15.vtp new file mode 100644 index 00000000..2f582175 Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_15.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_16.vtp b/models/moving_atrium/model_moved/model_moved_16.vtp new file mode 100644 index 00000000..fb5eac3e Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_16.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_17.vtp b/models/moving_atrium/model_moved/model_moved_17.vtp new file mode 100644 index 00000000..e7fa3404 Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_17.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_18.vtp b/models/moving_atrium/model_moved/model_moved_18.vtp new file mode 100644 index 00000000..60560af7 Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_18.vtp differ diff --git a/models/moving_atrium/model_moved/model_moved_19.vtp b/models/moving_atrium/model_moved/model_moved_19.vtp new file mode 100644 index 00000000..9163a075 Binary files /dev/null and b/models/moving_atrium/model_moved/model_moved_19.vtp differ diff --git a/pyproject.toml b/pyproject.toml index 454f1760..17e0473f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,7 +11,8 @@ readme = "README.md" dependencies = [ 'numpy', 'matplotlib', - 'cppimport' + 'cppimport', + 'scipy' ] [project.scripts] diff --git a/requirements.txt b/requirements.txt index 3e0cb0bf..47eb2953 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,4 +2,5 @@ jupyter-book matplotlib numpy ghp-import -cppimport \ No newline at end of file +cppimport +scipy \ No newline at end of file diff --git a/src/vampy/__init__.py b/src/vampy/__init__.py index ad77fd28..7ba1283a 100644 --- a/src/vampy/__init__.py +++ b/src/vampy/__init__.py @@ -30,6 +30,11 @@ except ModuleNotFoundError: print("WARNING: Oasis is not installed, running CFD is not available") +try: + from .simulation import MovingAtrium +except ModuleNotFoundError: + print("WARNING: OasisMove is not installed, running moving domain simulations (MovingAtrium) is not available") + meta = metadata("vampy") __version__ = meta["Version"] @@ -58,4 +63,5 @@ "simulation_common", "Probe", "Womersley", + 'MovingAtrium' ] diff --git a/src/vampy/automatedPreprocessing/automated_preprocessing.py b/src/vampy/automatedPreprocessing/automated_preprocessing.py index 76f2f4ba..4d6d0ae1 100644 --- a/src/vampy/automatedPreprocessing/automated_preprocessing.py +++ b/src/vampy/automatedPreprocessing/automated_preprocessing.py @@ -8,12 +8,14 @@ get_vtk_point_locator, extract_single_line, vtk_merge_polydata, get_point_data_array, smooth_voronoi_diagram, \ create_new_surface, compute_centers, vmtk_smooth_surface, str2bool, vmtk_compute_voronoi_diagram, \ prepare_output_surface, vmtk_compute_geometric_features + # Local imports from vampy.automatedPreprocessing import ToolRepairSTL +from vampy.automatedPreprocessing.moving_common import get_point_map, project_displacement, save_displacement from vampy.automatedPreprocessing.preprocessing_common import read_polydata, get_centers_for_meshing, \ dist_sphere_diam, dist_sphere_curvature, dist_sphere_constant, get_regions_to_refine, add_flow_extension, \ write_mesh, mesh_alternative, generate_mesh, find_boundaries, compute_flow_rate, setup_model_network, \ - radiusArrayName, scale_surface, get_furtest_surface_point, check_if_closed_surface + radiusArrayName, scale_surface, get_furtest_surface_point, check_if_closed_surface, remesh_surface from vampy.automatedPreprocessing.simulate import run_simulation from vampy.automatedPreprocessing.visualize import visualize_model @@ -22,7 +24,7 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f meshing_method, refine_region, is_atrium, add_flow_extensions, visualize, config_path, coarsening_factor, inlet_flow_extension_length, outlet_flow_extension_length, edge_length, region_points, compress_mesh, add_boundary_layer, scale_factor, resampling_step, - flow_rate_factor): + flow_rate_factor, moving_mesh, clamp_boundaries): """ Automatically generate mesh of surface model in .vtu and .xml format, including prescribed flow rates at inlet and outlet based on flow network model. @@ -51,6 +53,8 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f scale_factor (float): Scale input model by this factor resampling_step (float): Float value determining the resampling step for centerline computations, in [m] flow_rate_factor (float): Flow rate factor + moving_mesh (bool): Computes projected movement for displaced surfaces located in [filename_model]_moved folder + clamp_boundaries (bool): Clamps inlet(s) and outlet(s) if true """ # Get paths case_name = input_model.rsplit(path.sep, 1)[-1].rsplit('.')[0] @@ -65,6 +69,7 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f file_name_distance_to_sphere_diam = base_path + "_distance_to_sphere_diam.vtp" file_name_distance_to_sphere_const = base_path + "_distance_to_sphere_const.vtp" file_name_distance_to_sphere_curv = base_path + "_distance_to_sphere_curv.vtp" + file_name_distance_to_sphere_initial = base_path + "_distance_to_sphere_initial.vtp" file_name_probe_points = base_path + "_probe_point.json" file_name_voronoi = base_path + "_voronoi.vtp" file_name_voronoi_smooth = base_path + "_voronoi_smooth.vtp" @@ -76,14 +81,21 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f file_name_surface_name = base_path + "_remeshed_surface.vtp" file_name_xml_mesh = base_path + ".xml" file_name_vtu_mesh = base_path + ".vtu" + file_name_remeshed = base_path + "_remeshed.vtp" region_centerlines = None + # Dynamic mesh files + file_name_displacement_points = base_path + "_points.np" + folder_moved_surfaces = base_path + "_moved" + folder_extended_surfaces = base_path + "_extended" + # Open the surface file. print("--- Load model file\n") surface = read_polydata(input_model) # Scale surface if scale_factor is not None: + print(f"--- Scaling model by a factor {scale_factor}\n") surface = scale_surface(surface, scale_factor) resampling_step *= scale_factor @@ -146,8 +158,10 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f if refine_region: regions = get_regions_to_refine(capped_surface, region_points, base_path) for i in range(len(regions) // 3): - print("--- Region to refine ({}): {:.3f} {:.3f} {:.3f}" - .format(i + 1, regions[3 * i], regions[3 * i + 1], regions[3 * i + 2])) + print( + f"--- Region to refine ({i + 1}): " + + f"{regions[3 * i]:.3f} {regions[3 * i + 1]:.3f} {regions[3 * i + 2]:.3f}\n" + ) centerline_region, _, _ = compute_centerlines(source, regions, file_name_refine_region_centerlines, capped_surface, resampling=resampling_step) @@ -221,10 +235,9 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f if num_outlets != num_outlets_after: write_polydata(surface, file_name_surface_smooth) - print(("ERROR: Automatic clipping failed. You have to open {} and " + - "manually clipp the branch which still is capped. " + - "Overwrite the current {} and restart the script.").format( - file_name_surface_smooth, file_name_surface_smooth)) + print(f"ERROR: Automatic clipping failed. You have to open {file_name_surface_smooth} and " + + "manually clipp the branch which still is capped. " + + f"Overwrite the current {file_name_surface_smooth} and restart the script.") sys.exit(0) surface = surface_uncapped @@ -239,7 +252,7 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f surface = read_polydata(file_name_surface_smooth) elif smoothing_method in ["laplace", "taubin"]: - print("--- Smooth surface: {} smoothing\n".format(smoothing_method.capitalize())) + print(f"--- Smooth surface: {smoothing_method.capitalize()} smoothing\n") if not path.isfile(file_name_surface_smooth): surface = vmtk_smooth_surface(surface, smoothing_method, iterations=smoothing_iterations, passband=0.1, relaxation=0.01) @@ -252,6 +265,20 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f elif smoothing_method == "no_smooth" or None: print("--- No smoothing of surface\n") + if edge_length is not None and moving_mesh: + if path.exists(file_name_remeshed): + remeshed = read_polydata(file_name_remeshed) + else: + print("\n--- Remeshing surface for moving mesh\n") + surface = dist_sphere_constant(surface, centerlines, region_center, misr_max, + file_name_distance_to_sphere_initial, edge_length) + + remeshed = remesh_surface(surface, edge_length, "edgelengtharray") + remeshed = vtk_clean_polydata(remeshed) + write_polydata(remeshed, file_name_remeshed) + else: + remeshed = surface + # Add flow extensions if add_flow_extensions: if not path.isfile(file_name_model_flow_ext): @@ -263,14 +290,14 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f inlet_flow_extension_length, outlet_flow_extension_length = \ outlet_flow_extension_length, inlet_flow_extension_length - # Add extensions to inlet (artery) - surface_extended = add_flow_extension(surface, centerlines, is_inlet=True, - extension_length=inlet_flow_extension_length, - extension_mode=extension) + # Add extensions to inlet (artery) / outlet (atrium) + surface_extended = add_flow_extension(remeshed, centerlines, is_inlet=True, + extension_length=inlet_flow_extension_length) - # Add extensions to outlets (artery) + # Add extensions to outlets (artery) / inlets (Atrium) surface_extended = add_flow_extension(surface_extended, centerlines, is_inlet=False, - extension_length=outlet_flow_extension_length) + extension_length=outlet_flow_extension_length, + extension_mode=extension) surface_extended = vmtk_smooth_surface(surface_extended, "laplace", iterations=200) write_polydata(surface_extended, file_name_model_flow_ext) @@ -279,13 +306,26 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f else: surface_extended = surface - # Capp surface with flow extensions - capped_surface = vmtk_cap_polydata(surface_extended) + # Create displacement input file for moving mesh simulations + if moving_mesh: + print("--- Computing mesh displacement and saving points to file") + # Get a point mapper + distance, point_map = get_point_map(remeshed, surface_extended) + + # Project displacement between surfaces + points = project_displacement(clamp_boundaries, distance, folder_extended_surfaces, folder_moved_surfaces, + point_map, surface, surface_extended, remeshed, scale_factor) + + # Save displacement to numpy array + save_displacement(file_name_displacement_points, points) # Get new centerlines with the flow extensions if add_flow_extensions: if not path.isfile(file_name_flow_centerlines): print("--- Compute the model centerlines with flow extension.\n") + # Capp surface with flow extensions + capped_surface = vmtk_cap_polydata(surface_extended) + # Compute the centerlines. if has_outlet: inlet, outlets = get_centers_for_meshing(surface_extended, is_atrium, base_path, @@ -382,7 +422,8 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f file_name_centerlines, file_name_refine_region_centerlines, file_name_region_centerlines, file_name_distance_to_sphere_diam, file_name_distance_to_sphere_const, file_name_distance_to_sphere_curv, file_name_voronoi, file_name_voronoi_smooth, file_name_voronoi_surface, file_name_surface_smooth, - file_name_model_flow_ext, file_name_clipped_model, file_name_flow_centerlines, file_name_surface_name + file_name_model_flow_ext, file_name_clipped_model, file_name_flow_centerlines, file_name_surface_name, + file_name_remeshed, file_name_distance_to_sphere_initial ] for file in files_to_remove: if path.exists(file): @@ -534,6 +575,17 @@ def read_command_line(input_path=None): type=float, help="Flow rate factor.") + parser.add_argument('-mm', '--moving-mesh', + action="store_true", + default=False, + help="If true, assumes a dynamic/moving mesh and will perform computation of projection " + + "between moved surfaces located in the '[filename_model]_moved' folder.") + + parser.add_argument('-cl', '--clamp-boundaries', + action="store_true", + default=False, + help="Clamps boundaries at inlet(s) and outlet(s) if true. Only used for moving mesh.") + # Parse path to get default values if required: args = parser.parse_args() @@ -569,7 +621,8 @@ def verbose_print(*args): visualize=args.visualize, region_points=args.region_points, compress_mesh=args.compress_mesh, outlet_flow_extension_length=args.outlet_flowextension, add_boundary_layer=args.add_boundary_layer, scale_factor=args.scale_factor, resampling_step=args.resampling_step, - flow_rate_factor=args.flow_rate_factor) + flow_rate_factor=args.flow_rate_factor, moving_mesh=args.moving_mesh, + clamp_boundaries=args.clamp_boundaries) def main_meshing(): diff --git a/src/vampy/automatedPreprocessing/moving_common.py b/src/vampy/automatedPreprocessing/moving_common.py new file mode 100644 index 00000000..3a9d1814 --- /dev/null +++ b/src/vampy/automatedPreprocessing/moving_common.py @@ -0,0 +1,245 @@ +from os import listdir, mkdir + +from morphman.common import * +from scipy.interpolate import interp1d +from vtk.numpy_interface import dataset_adapter as dsa + +from vampy.automatedPreprocessing.preprocessing_common import scale_surface + + +############################################################## +# A Collection of utility scripts for moving mesh generation # +############################################################## + +def get_point_map(remeshed, remeshed_extended, profile="linear"): + """ + Create a map between original surface model and model with + flow extensions; remeshed_extended, for clamping the inlet(s) + and outlet. A linear or sinusodial profile for movement between endpoints + and moving domain is prescribed. + Args: + remeshed (vtkPolyData): Input surface + remeshed_extended (vtkPolyData): Input surface with flow extensions + profile (str): 'linear' or 'sine' for respective profiles + Returns: + distances (ndarray): Array of linear or sinusodial profile between boundary and moving mesh + point_map (ndarray): Mapping between surface IDs and points along extension + """ + remeshed = dsa.WrapDataObject(remeshed) + remeshed_extended = dsa.WrapDataObject(remeshed_extended) + + # Get lengths + num_re = remeshed.Points.shape[0] + num_ext = remeshed_extended.Points.shape[0] - remeshed.Points.shape[0] + + # Get locators + inner_feature = vtk_compute_connectivity(vtk_extract_feature_edges(remeshed.VTKObject)) + outer_feature = vtk_compute_connectivity(vtk_extract_feature_edges(remeshed_extended.VTKObject)) + locator_remeshed = get_vtk_point_locator(remeshed.VTKObject) + + n_features = outer_feature.GetPointData().GetArray("RegionId").GetValue(outer_feature.GetNumberOfPoints() - 1) + inner_features = np.array( + [vtk_compute_threshold(inner_feature, "RegionId", i, i + 0.1, source=0) for i in range(n_features + 1)]) + outer_features = np.array( + [vtk_compute_threshold(outer_feature, "RegionId", i, i + 0.1, source=0) for i in range(n_features + 1)]) + + inner_regions = [dsa.WrapDataObject(feature) for feature in inner_features] + inner_locators = [get_vtk_point_locator(feature) for feature in inner_features] + inner_points = [feature.GetNumberOfPoints() for feature in inner_features] + + outer_regions = [dsa.WrapDataObject(feature) for feature in outer_features] + outer_locators = [get_vtk_point_locator(feature) for feature in outer_features] + outer_points = [feature.GetNumberOfPoints() for feature in outer_features] + boundary_map = {i: j for i, j in zip(np.argsort(inner_points), np.argsort(outer_points))} + + # Get distance and point map + distances = np.zeros(num_ext) + point_map = np.zeros(num_ext) + + for i in range(num_ext): + point = remeshed_extended.Points[num_re + i] + tmp_id = -1 + tmp_p_norm = 1e16 + # Some hacks to find the correct corresponding points + for region_id in range(len(outer_features)): + region_id_out = boundary_map[region_id] + id_i = inner_locators[region_id].FindClosestPoint(point) + id_o = outer_locators[region_id_out].FindClosestPoint(point) + + p_i = inner_features[region_id].GetPoint(id_i) + p_o = outer_features[region_id_out].GetPoint(id_o) + + # Compute distances from current point to boundaries + p_a = np.array(p_i) - np.array(point) + p_b = np.array(p_o) - np.array(point) + p_norm = np.linalg.norm(p_b) + np.linalg.norm(p_a) + if p_norm < tmp_p_norm: + tmp_id = region_id + tmp_p_norm = p_norm + + regionId = tmp_id + regionId_out = boundary_map[regionId] + inner_id = inner_locators[regionId].FindClosestPoint(point) + outer_id = outer_locators[regionId_out].FindClosestPoint(point) + + dist_to_boundary = get_distance(point, outer_regions[regionId_out].Points[outer_id]) + dist_between_boundaries = get_distance(inner_regions[regionId].Points[inner_id], + outer_regions[regionId_out].Points[outer_id]) + if profile == "linear": + distances[i] = (dist_to_boundary / dist_between_boundaries) + elif profile == "sine": + distances[i] = easing_cos(dist_to_boundary, dist_between_boundaries) + point_map[i] = locator_remeshed.FindClosestPoint(inner_regions[regionId].Points[inner_id]) + + # Let the points corresponding to the caps have distance 0 + point_map = point_map.astype(int) + + return distances, point_map + + +def easing_cos(dist_b, dist_bb): + """ + A simple easing function between 0 and 1 for creating a sinusoidal profile + Args: + dist_b (ndarray): Distance to boundary for current point + dist_bb (ndarray): Distance between boundaries + + Returns: + distances (ndarray): Sinusoidal profile based on distance between point and boundary + """ + return - 0.5 * (np.cos(np.pi * dist_b / dist_bb) - 1) + + +def move_surface_model(surface, original, remeshed, remeshed_extended, distance, point_map, file_path, i, points, + clamp_boundaries): + """ + Computes and applies the displacement between the original and the given surface, then projects this + displacement onto a remeshed surface. The displaced points of the remeshed surface are then saved + to a specified file path. + + Args: + surface (vtkPolyData): The source surface used to compute the displacement. + original (vtkPolyData): The reference surface for computing the displacement. + remeshed (vtkPolyData): A remeshed version of the original surface. + remeshed_extended (vtkPolyData): An extended version of the remeshed surface. + distance (float): Scalar factor for the displacement. + point_map (ndarray): Mapping of the points from remeshed to remeshed_extended. + file_path (str): Path to save the displaced points of the remeshed surface. + i (int): Index to specify where in the 'points' array the data should be saved. + points (ndarray): An array where the displaced points will be saved. + clamp_boundaries (bool): If True, the displacement will be clamped by the distance in the extensions. + + Returns: + None: The function writes results to a file and modifies the points array in place. + """ + surface = dsa.WrapDataObject(surface) + original = dsa.WrapDataObject(original) + remeshed = dsa.WrapDataObject(remeshed) + remeshed_extended = dsa.WrapDataObject(remeshed_extended) + + if "displacement" in original.PointData.keys(): + original.VTKObject.GetPointData().RemoveArray("displacement") + + if "displacement" in remeshed_extended.PointData.keys(): + remeshed_extended.VTKObject.GetPointData().RemoveArray("displacement") + + # Get displacement field + original.PointData.append(surface.Points - original.Points, "displacement") + + # Get surface projection + projector = vmtkscripts.vmtkSurfaceProjection() + projector.Surface = remeshed_extended.VTKObject + projector.ReferenceSurface = original.VTKObject + projector.Execute() + + # New surface + new_surface = projector.Surface + new_surface = dsa.WrapDataObject(new_surface) + + # Manipulate displacement in the extensions + displacement = new_surface.PointData["displacement"] + if clamp_boundaries: + displacement[remeshed.Points.shape[0]:] = distance * displacement[point_map] + else: + displacement[remeshed.Points.shape[0]:] = displacement[point_map] + + # Move the mesh points + new_surface.Points += displacement + write_polydata(new_surface.VTKObject, file_path) + points[:, :, i] = new_surface.Points.copy() + new_surface.Points -= displacement + + +def save_displacement(file_name_displacement_points, points, number_of_points=100): + """ + Resample displacement points and write them + to numpy data array + Args: + file_name_displacement_points (str): Path to numpy point path + points (ndarray): Numpy array consisting of displacement surface points + number_of_points (int): Number of points to interpolate displacement points over + """ + points[:, :, -1] = points[:, :, 0] + + time = np.linspace(0, 1, points.shape[2]) + move = np.zeros((points.shape[0], points.shape[1], number_of_points + 1)) + time_r = np.linspace(0, 1, number_of_points + 1) + # Use interp1d if smooth displacement + x = interp1d(time, points[:, 0, :], axis=1) + y = interp1d(time, points[:, 1, :], axis=1) + z = interp1d(time, points[:, 2, :], axis=1) + + move[:, 0, :] = x(time_r) + move[:, 1, :] = y(time_r) + move[:, 2, :] = z(time_r) + + points = move + points.dump(file_name_displacement_points) + + +def project_displacement(clamp_boundaries, distance, folder_extended_surfaces, folder_moved_surfaces, point_map, + surface, surface_extended, remeshed, scale_factor): + """ + Projects the displacement of moved surfaces located in a specified folder onto an extended surface. The projected + surfaces are then saved to a designated output directory. + Only processes files in "folder_moved_surfaces" that end with ".vtp" or ".stl". + + Args: + clamp_boundaries (bool): If True, the displacement will be clamped by the distance in the extensions. + distance (float): Scalar factor for the displacement. + folder_extended_surfaces (str): Path to the directory where extended surfaces will be saved. + folder_moved_surfaces (str): Path to the directory containing the moved surfaces to be projected. + point_map (ndarray): Mapping of the points from remeshed to surface_extended. + surface (vtkPolyData): The reference surface for computing the displacement. + surface_extended (vtkPolyData): An extended version of the reference surface. + remeshed (vtkPolyData): A remeshed version of the original surface. + scale_factor (float): Factor to scale the surfaces. + + Returns: + array: A 3D numpy array where the displaced points will be saved. The shape is + (num_points, 3, num_surfaces). + """ + # Add extents to all surfaces + extended_surfaces = sorted([f for f in listdir(folder_moved_surfaces) if f.endswith(".vtp") or f.endswith(".stl")]) + if not path.exists(folder_extended_surfaces): + mkdir(folder_extended_surfaces) + + n_surfaces = len(extended_surfaces) + + print("--- Projecting surfaces\n") + points = np.zeros((surface_extended.GetNumberOfPoints(), 3, n_surfaces)) + for i in range(n_surfaces): + model_path = path.join(folder_moved_surfaces, extended_surfaces[i]) + if i == 0: + points[:, :, i] = dsa.WrapDataObject(surface_extended).Points + continue + + tmp_surface = read_polydata(model_path) + if scale_factor is not None: + tmp_surface = scale_surface(tmp_surface, scale_factor) + + new_path = path.join(folder_extended_surfaces, model_path.split("/")[-1]) + if not path.exists(new_path): + move_surface_model(tmp_surface, surface, remeshed, surface_extended, distance, point_map, new_path, i, + points, clamp_boundaries) + return points diff --git a/src/vampy/automatedPreprocessing/preprocessing_common.py b/src/vampy/automatedPreprocessing/preprocessing_common.py index d7415ecf..32efc885 100644 --- a/src/vampy/automatedPreprocessing/preprocessing_common.py +++ b/src/vampy/automatedPreprocessing/preprocessing_common.py @@ -1,8 +1,9 @@ import gzip -from os import path import json +from os import path import numpy as np +import vtkmodules.numpy_interface.dataset_adapter as dsa from morphman import vtk_clean_polydata, vtk_triangulate_surface, get_parameters, write_parameters, read_polydata, \ vmtkscripts, vtk, write_polydata, vtkvmtk, get_curvilinear_coordinate, vtk_compute_threshold, get_vtk_array, \ get_distance, get_number_of_arrays, vmtk_smooth_surface, get_point_data_array, create_vtk_array, \ @@ -15,6 +16,7 @@ # Global array names distanceToSpheresArrayName = "DistanceToSpheres" radiusArrayName = 'MaximumInscribedSphereRadius' +cellEntityArrayName = "CellEntityIds" def get_regions_to_refine(surface, provided_points, dir_path): @@ -876,3 +878,44 @@ def check_if_closed_surface(surface): cells = vtk_extract_feature_edges(surface) return cells.GetNumberOfCells() == 0 + + +def remesh_surface(surface, edge_length, element_size_mode="edgelength", exclude=None): + """ + Remeshes a given surface based on the specified parameters. + + Args: + surface (vtkPolyData): The input surface to be remeshed. + edge_length (float): The target edge length for remeshing. + element_size_mode (str, optional): Determines the method for sizing elements during remeshing. + exclude (list, optional): A list of entity IDs to be excluded during remeshing. Defaults to None. + + Returns: + remeshed_surface (vtkPolyData): The remeshed surface. + """ + surface = dsa.WrapDataObject(surface) + if cellEntityArrayName not in surface.CellData.keys(): + surface.CellData.append(np.zeros(surface.VTKObject.GetNumberOfCells()) + 1, cellEntityArrayName) + + remeshing = vmtkscripts.vmtkSurfaceRemeshing() + remeshing.Surface = surface.VTKObject + remeshing.CellEntityIdsArrayName = cellEntityArrayName + remeshing.TargetEdgeLength = edge_length + remeshing.MaxEdgeLength = 1e6 + remeshing.MinEdgeLength = 0.0 + remeshing.TargetEdgeLengthFactor = 1.0 + remeshing.TriangleSplitFactor = 5.0 + remeshing.ElementSizeMode = element_size_mode + if element_size_mode == "edgelength": + remeshing.TargetEdgeLengthArrayName = "" + else: + remeshing.TargetEdgeLengthArrayName = "Size" # Variable size mesh + + if exclude is not None: + remeshing.ExcludeEntityIds = exclude + + remeshing.Execute() + + remeshed_surface = remeshing.Surface + + return remeshed_surface diff --git a/src/vampy/simulation/Artery.py b/src/vampy/simulation/Artery.py index 9033fffa..530b4886 100755 --- a/src/vampy/simulation/Artery.py +++ b/src/vampy/simulation/Artery.py @@ -1,6 +1,5 @@ import json import pickle -from os import makedirs from pprint import pprint from dolfin import set_log_level @@ -9,7 +8,8 @@ from vampy.simulation.Probe import Probes # type: ignore from vampy.simulation.Womersley import make_womersley_bcs, compute_boundary_geometry_acrn -from vampy.simulation.simulation_common import get_file_paths, store_u_mean, print_mesh_information +from vampy.simulation.simulation_common import get_file_paths, store_u_mean, print_mesh_information, \ + store_velocity_and_pressure_h5, dump_probes # FEniCS specific command to control the desired level of logging, here set to critical errors set_log_level(50) @@ -124,11 +124,12 @@ def create_bcs(t, NS_expressions, V, Q, area_ratio, area_inlet, mesh, mesh_path, ds = Measure("ds", domain=mesh, subdomain_data=boundary) for i, ind in enumerate(id_out): dsi = ds(ind) - area_out.append(assemble(Constant(1.0, name="one") * dsi)) + area_out.append(assemble(Constant(1.0) * dsi)) bc_p = [] if MPI.rank(MPI.comm_world) == 0: print("=== Initial pressure and area fraction ===") + for i, ID in enumerate(id_out): p_initial = area_out[i] / sum(area_out) outflow = Expression("p", p=p_initial, degree=pressure_degree) @@ -136,7 +137,7 @@ def create_bcs(t, NS_expressions, V, Q, area_ratio, area_inlet, mesh, mesh_path, bc_p.append(bc) NS_expressions[ID] = outflow if MPI.rank(MPI.comm_world) == 0: - print(("Boundary ID={:d}, pressure: {:0.6f}, area fraction: {:0.4f}".format(ID, p_initial, area_ratio[i]))) + print(f"Boundary ID={ID}, pressure: {p_initial:.5f}, area fraction: {area_ratio[i]:0.5f}") # No slip condition at wall wall = Constant(0.0) @@ -219,12 +220,11 @@ def temporal_hook(u_, p_, mesh, tstep, dump_probe_frequency, eval_dict, newfolde diam_inlet = np.sqrt(4 * area_inlet[0] / np.pi) Re = U_mean * diam_inlet / nu print("=" * 10, "Time step " + str(tstep), "=" * 10) - print("Sum of Q_out = {:0.4f}, Q_in = {:0.4f}, mean velocity (inlet): {:0.4f}, Reynolds number (inlet): {:0.4f}" - .format(sum(Q_outs), Q_in, U_mean, Re)) + print(f"Sum of Q_out = {sum(Q_outs):.4f}, Q_in = {Q_in:.4f}, mean velocity (inlet): {U_mean:.4f}, " + + f"Reynolds number (inlet): {Re:.4f}") for i, out_id in enumerate(id_out): - print(("For outlet with boundary ID={:d}: target flow rate: {:0.4f} mL/s, " + - "computed flow rate: {:0.4f} mL/s, pressure updated to: {:0.4f}") - .format(out_id, Q_ideals[i], Q_outs[i], NS_expressions[out_id].p)) + print("For outlet with boundary ID={out_id}: target flow rate: {Q_ideals[i]:.4f} mL/s, " + + f"computed flow rate: {Q_outs[i]:.4f} mL/s, pressure updated to: {NS_parameters[out_id].p:.4f}") print() # Sample velocity and pressure in points/probes @@ -235,59 +235,11 @@ def temporal_hook(u_, p_, mesh, tstep, dump_probe_frequency, eval_dict, newfolde # Store sampled velocity and pressure if tstep % dump_probe_frequency == 0: - # Save variables along the centerline for CFD simulation - # diagnostics and light-weight post-processing - filepath = path.join(newfolder, "Probes") - if MPI.rank(MPI.comm_world) == 0: - if not path.exists(filepath): - makedirs(filepath) - - arr_u_x = eval_dict["centerline_u_x_probes"].array() - arr_u_y = eval_dict["centerline_u_y_probes"].array() - arr_u_z = eval_dict["centerline_u_z_probes"].array() - arr_p = eval_dict["centerline_p_probes"].array() - - # Dump stats - if MPI.rank(MPI.comm_world) == 0: - arr_u_x.dump(path.join(filepath, "u_x_%s.probes" % str(tstep))) - arr_u_y.dump(path.join(filepath, "u_y_%s.probes" % str(tstep))) - arr_u_z.dump(path.join(filepath, "u_z_%s.probes" % str(tstep))) - arr_p.dump(path.join(filepath, "p_%s.probes" % str(tstep))) - - # Clear stats - MPI.barrier(MPI.comm_world) - eval_dict["centerline_u_x_probes"].clear() - eval_dict["centerline_u_y_probes"].clear() - eval_dict["centerline_u_z_probes"].clear() - eval_dict["centerline_p_probes"].clear() + dump_probes(eval_dict, newfolder, tstep) # Save velocity and pressure for post-processing if tstep % save_solution_frequency == 0 and tstep >= save_solution_at_tstep: - # Assign velocity components to vector solution - assign(U.sub(0), u_[0]) - assign(U.sub(1), u_[1]) - assign(U.sub(2), u_[2]) - - # Get save paths - files = NS_parameters['files'] - p_path = files['p'] - u_path = files['u'] - file_mode = "w" if not path.exists(p_path) else "a" - - # Save pressure - viz_p = HDF5File(MPI.comm_world, p_path, file_mode=file_mode) - viz_p.write(p_, "/pressure", tstep) - viz_p.close() - - # Save velocity - viz_u = HDF5File(MPI.comm_world, u_path, file_mode=file_mode) - viz_u.write(U, "/velocity", tstep) - viz_u.close() - - # Accumulate velocity - u_mean0.vector().axpy(1, u_[0].vector()) - u_mean1.vector().axpy(1, u_[1].vector()) - u_mean2.vector().axpy(1, u_[2].vector()) + store_velocity_and_pressure_h5(NS_parameters, U, p_, tstep, u_, u_mean0, u_mean1, u_mean2) # Oasis hook called after the simulation has finished diff --git a/src/vampy/simulation/Atrium.py b/src/vampy/simulation/Atrium.py index 17dc27c9..378b4ff3 100644 --- a/src/vampy/simulation/Atrium.py +++ b/src/vampy/simulation/Atrium.py @@ -1,6 +1,5 @@ import json import pickle -from os import makedirs from pprint import pprint from dolfin import set_log_level, MPI @@ -9,7 +8,8 @@ from vampy.simulation.Probe import Probes # type: ignore from vampy.simulation.Womersley import make_womersley_bcs, compute_boundary_geometry_acrn -from vampy.simulation.simulation_common import store_u_mean, get_file_paths, print_mesh_information +from vampy.simulation.simulation_common import store_u_mean, get_file_paths, print_mesh_information, \ + store_velocity_and_pressure_h5, dump_probes # FEniCS specific command to control the desired level of logging, here set to critical errors set_log_level(50) @@ -116,19 +116,18 @@ def create_bcs(NS_expressions, mesh, mesh_path, nu, t, V, Q, id_in, id_out, **NS Q_scaled = Q_values * tmp_area / area_total # Create Womersley boundary condition at inlet - inlet = make_womersley_bcs(t_values, Q_scaled, nu, tmp_center, tmp_radius, tmp_normal, - V.ufl_element()) - NS_expressions["inlet_{}".format(ID)] = inlet + inlet = make_womersley_bcs(t_values, Q_scaled, nu, tmp_center, tmp_radius, tmp_normal, V.ufl_element()) + NS_expressions[f"inlet_{ID}"] = inlet # Initial condition for ID in id_in: for i in [0, 1, 2]: - NS_expressions["inlet_{}".format(ID)][i].set_t(t) + NS_expressions[f"inlet_{ID}"][i].set_t(t) # Create inlet boundary conditions bc_inlets = {} for ID in id_in: - bc_inlet = [DirichletBC(V, NS_expressions["inlet_{}".format(ID)][i], boundary, ID) for i in range(3)] + bc_inlet = [DirichletBC(V, NS_expressions[f"inlet_{ID}"][i], boundary, ID) for i in range(3)] bc_inlets[ID] = bc_inlet # Set outlet boundary conditions, assuming one outlet (Mitral Valve) @@ -240,59 +239,11 @@ def temporal_hook(mesh, dt, t, save_solution_frequency, u_, NS_expressions, id_i # Store sampled velocity and pressure if tstep % dump_probe_frequency == 0: - # Save variables along the centerline for CFD simulation - # diagnostics and light-weight post-processing - filepath = path.join(newfolder, "Probes") - if MPI.rank(MPI.comm_world) == 0: - if not path.exists(filepath): - makedirs(filepath) - - arr_u_x = eval_dict["centerline_u_x_probes"].array() - arr_u_y = eval_dict["centerline_u_y_probes"].array() - arr_u_z = eval_dict["centerline_u_z_probes"].array() - arr_p = eval_dict["centerline_p_probes"].array() - - # Dump stats - if MPI.rank(MPI.comm_world) == 0: - arr_u_x.dump(path.join(filepath, "u_x_%s.probes" % str(tstep))) - arr_u_y.dump(path.join(filepath, "u_y_%s.probes" % str(tstep))) - arr_u_z.dump(path.join(filepath, "u_z_%s.probes" % str(tstep))) - arr_p.dump(path.join(filepath, "p_%s.probes" % str(tstep))) - - # Clear stats - MPI.barrier(MPI.comm_world) - eval_dict["centerline_u_x_probes"].clear() - eval_dict["centerline_u_y_probes"].clear() - eval_dict["centerline_u_z_probes"].clear() - eval_dict["centerline_p_probes"].clear() + dump_probes(eval_dict, newfolder, tstep) # Save velocity and pressure for post-processing if tstep % save_solution_frequency == 0 and tstep >= save_solution_at_tstep: - # Assign velocity components to vector solution - assign(U.sub(0), u_[0]) - assign(U.sub(1), u_[1]) - assign(U.sub(2), u_[2]) - - # Get save paths - files = NS_parameters['files'] - p_path = files['p'] - u_path = files['u'] - file_mode = "w" if not path.exists(p_path) else "a" - - # Save pressure - viz_p = HDF5File(MPI.comm_world, p_path, file_mode=file_mode) - viz_p.write(p_, "/pressure", tstep) - viz_p.close() - - # Save velocity - viz_u = HDF5File(MPI.comm_world, u_path, file_mode=file_mode) - viz_u.write(U, "/velocity", tstep) - viz_u.close() - - # Accumulate velocity - u_mean0.vector().axpy(1, u_[0].vector()) - u_mean1.vector().axpy(1, u_[1].vector()) - u_mean2.vector().axpy(1, u_[2].vector()) + store_velocity_and_pressure_h5(NS_parameters, U, p_, tstep, u_, u_mean0, u_mean1, u_mean2) # Oasis hook called after the simulation has finished diff --git a/src/vampy/simulation/MovingAtrium.py b/src/vampy/simulation/MovingAtrium.py new file mode 100644 index 00000000..baec67b5 --- /dev/null +++ b/src/vampy/simulation/MovingAtrium.py @@ -0,0 +1,433 @@ +import json +import pickle +from os import makedirs +from pprint import pprint +from dolfin import set_log_level, UserExpression + +from oasismove.problems.NSfracStep import * +from oasismove.problems.NSfracStep.MovingCommon import get_visualization_writers +from scipy.interpolate import splrep, splev +from scipy.spatial import KDTree + +from vampy.simulation.Probe import Probes # type: ignore +from vampy.simulation.Womersley import make_womersley_bcs, compute_boundary_geometry_acrn +from vampy.simulation.simulation_common import store_u_mean, get_file_paths, print_mesh_information, \ + store_velocity_and_pressure_h5, dump_probes + +# FEniCS specific command to control the desired level of logging, here set to critical errors +set_log_level(50) + + +def problem_parameters(commandline_kwargs, scalar_components, Schmidt, NS_parameters, **NS_namespace): + """ + Problem file for running CFD simulation in left atrial models consisting of arbitrary number of pulmonary veins (PV) + (normally 3 to 7), and one outlet being the mitral valve. A Womersley velocity profile is applied at the inlets, + where the total flow rate is split between the area ratio of the PVs. The mitral valve is considered open with a + constant pressure of p=0. Flow rate and flow split values for the inlet condition are computed from the + pre-processing script automatedPreProcessing.py. The simulation is run for two cycles (adjustable), but only the + results/solutions from the second cycle are stored to avoid non-physiological effects from the first cycle. + One cardiac cycle is set to 0.951 s from [1], and scaled by a factor of 1000, hence all parameters are in + [mm] or [ms]. + The script has been adjusted to handle moving domains. + + [1] Hoi, Yiemeng, et al. "Characterization of volumetric flow rate waveforms at the carotid bifurcations of older + adults." Physiological measurement 31.3 (2010): 291. + """ + + if "restart_folder" in commandline_kwargs.keys(): + restart_folder = commandline_kwargs["restart_folder"] + mesh_path = commandline_kwargs["mesh_path"] + f = open(path.join(restart_folder, 'params.dat'), 'rb') + NS_parameters.update(pickle.load(f)) + NS_parameters['restart_folder'] = restart_folder + NS_parameters['mesh_path'] = mesh_path + globals().update(NS_parameters) + else: + # Override some problem specific parameters + # Parameters are in mm and ms + cardiac_cycle = float(commandline_kwargs.get("cardiac_cycle", 1000)) + number_of_cycles = int(commandline_kwargs.get("number_of_cycles", 2)) + track_blood = True + + NS_parameters.update( + # Moving atrium parameters + dynamic_mesh=True, # Run moving mesh simulation + compute_velocity_and_pressure=True, # Only solve mesh equations (For volume computation) + # Blood residence time + track_blood=track_blood, + # Backflow parameters + backflow_beta=0.2, + backflow_facets=[], + # Fluid parameters + nu=3.3018868e-3, # Viscosity [nu: 0.0035 Pa-s / 1060 kg/m^3 = 3.3018868E-6 m^2/s == 3.3018868E-3 mm^2/ms] + # Geometry parameters + id_in=[], # Inlet boundary ID + id_out=[], # Outlet boundary IDs + # Simulation parameters + cardiac_cycle=cardiac_cycle, # Run simulation for 1 cardiac cycles [ms] + T=cardiac_cycle * number_of_cycles, # Total simulation length + dt=0.1, # # Time step size [ms] + # Frequencies to save data + dump_probe_frequency=500, # Dump frequency for sampling velocity & pressure at probes along the centerline + save_solution_frequency=10, # Save frequency for velocity and pressure field + save_solution_frequency_xdmf=10, # Save frequency for velocity and pressure field + save_solution_after_cycle=0, # Store solution after 1 cardiac cycle + save_volume_frequency=1e10, # Save frequency for storing volume + save_flow_metrics_frequency=100, # Frequency for storing flow metrics + # Oasis specific parameters + checkpoint=1000, # Overwrite solution in Checkpoint folder each checkpoint + print_intermediate_info=500, + folder="results_moving_atrium", + mesh_path=commandline_kwargs["mesh_path"], + # Solver parameters + max_iter=2, + velocity_degree=1, + pressure_degree=1, + use_krylov_solvers=True, + krylov_solvers={'monitor_convergence': False, + 'report': False, + 'relative_tolerance': 1e-8, + 'absolute_tolerance': 1e-8} + ) + if track_blood: + if MPI.rank(MPI.comm_world) == 0: + print("-- Computing blood residence time --") + scalar_components += ["blood"] + + if MPI.rank(MPI.comm_world) == 0: + print("=== Starting simulation for MovingAtrium.py ===") + print("Running with the following parameters:") + pprint(NS_parameters) + + +def scalar_source(scalar_components, **NS_namespace): + """Return a dictionary of scalar sources.""" + return dict((ci, Constant(1)) for ci in scalar_components) + + +def mesh(mesh_path, **NS_namespace): + # Read mesh and print mesh information + if mesh_path.endswith(".xml") or mesh_path.endswith(".xml.gz"): + mesh = Mesh(mesh_path) + elif mesh_path.endswith(".h5"): + mesh = Mesh() + with HDF5File(MPI.comm_world, mesh_path, "r") as f: + f.read(mesh, "mesh", False) + + print_mesh_information(mesh) + + return mesh + + +class MeshMotionMapping(UserExpression): + def __init__(self, points, cycle, **kwargs): + self.motion_mapping = {} + self.counter = -1 + self.points = points + self.num_points = self.points.shape[0] + self.num_samples = self.points.shape[-1] + self.time = np.linspace(0, cycle, self.num_samples) + self.tree = self.create_kd_tree() + super().__init__(**kwargs) + + def get_motion_mapping(self): + return self.motion_mapping + + def create_kd_tree(self): + tree = KDTree(self.points[:, :, 0]) + return tree + + def eval(self, _, x): + self.counter += 1 + _, index = self.tree.query(x) + # FIXME: Set spline parameter objectively + s = 1E-2 + x_ = splrep(self.time, self.points[index, 0, :], s=s, per=True) + y_ = splrep(self.time, self.points[index, 1, :], s=s, per=True) + z_ = splrep(self.time, self.points[index, 2, :], s=s, per=True) + + self.motion_mapping[self.counter] = [x_, y_, z_] + + +class WallMotion(UserExpression): + def __init__(self, t, motion, max_counter, direction, cycle, **kwargs): + self.t = t + self.motion = motion + self.max_counter = max_counter + self.counter = -1 + self.direction = direction + self.cycle = cycle + super().__init__(**kwargs) + + def eval(self, values, _): + self.counter += 1 + # No motion for inlet/outlet + values[:] = splev(self.t % self.cycle, self.motion[self.counter][self.direction], der=1) + if self.counter == self.max_counter: + self.counter = -1 + + +def create_bcs(NS_expressions, dynamic_mesh, x_, cardiac_cycle, backflow_facets, mesh, mesh_path, nu, t, + V, Q, id_in, id_out, track_blood, **NS_namespace): + if mesh_path.endswith(".xml") or mesh_path.endswith(".xml.gz"): + mesh_filename = ".xml" + elif mesh_path.endswith(".h5"): + mesh_filename = ".h5" + + rank = MPI.rank(MPI.comm_world) + coords = ['x', 'y', 'z'] + # Variables needed during the simulation + boundary = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1, mesh.domains()) + if mesh_path.endswith(".h5"): + with HDF5File(MPI.comm_world, mesh_path, "r") as f: + f.read(boundary, "boundary") + + # Get IDs for inlet(s) and outlet(s) + info_path = mesh_path.split(mesh_filename)[0] + "_info.json" + with open(info_path) as f: + info = json.load(f) + + id_in[:] = info['inlet_ids'] + id_out[:] = info['outlet_id'] + id_wall = min(id_in + id_out) - 1 + + # Apply backflow at the outlet (Mitral valve) + backflow_facets[:] = info['outlet_id'] + + # Find corresponding areas + ds_new = Measure("ds", domain=mesh, subdomain_data=boundary) + area_total = 0 + for ID in id_in: + area_total += assemble(Constant(1.0) * ds_new(ID)) + + # Load normalized time and flow rate values + if dynamic_mesh: + flow_rate_path = mesh_path.split(mesh_filename)[0] + "_flowrate_moving.txt" + else: + flow_rate_path = mesh_path.split(mesh_filename)[0] + "_flowrate_rigid.txt" + + Q_mean = info['mean_flow_rate'] + t_values, Q_ = np.loadtxt(flow_rate_path).T + Q_values = Q_ * Q_mean # Scale normalized flow rate by mean flow rate + t_values *= 1000 # Scale time in normalised flow wave form to [ms] + + for ID in id_in: + tmp_area, tmp_center, tmp_radius, tmp_normal = compute_boundary_geometry_acrn(mesh, ID, boundary) + Q_scaled = Q_values * tmp_area / area_total + + # Create Womersley boundary condition at inlets + inlet = make_womersley_bcs(t_values, Q_scaled, nu, tmp_center, tmp_radius, tmp_normal, V.ufl_element()) + NS_expressions[f"inlet_{ID}"] = inlet + + # Initial condition + for ID in id_in: + for i in [0, 1, 2]: + NS_expressions[f"inlet_{ID}"][i].set_t(t) + + # Create inlet boundary conditions + bc_inlets = {} + for ID in id_in: + bc_inlet = [DirichletBC(V, NS_expressions[f"inlet_{ID}"][i], boundary, ID) for i in range(3)] + bc_inlets[ID] = bc_inlet + + # Set outlet boundary conditions, assuming one outlet (Mitral Valve) + bc_p = [DirichletBC(Q, Constant(0), boundary, ID) for ID in id_out] + + # Set wall boundary conditions + if dynamic_mesh: + # Moving walls + if rank == 0: + print("Loading displacement points") + points = np.load(mesh_path.split(mesh_filename)[0] + "_points.np", allow_pickle=True) + if rank == 0: + print("Creating splines for displacement") + + # Define wall movement + motion_mapping = MeshMotionMapping(points, cardiac_cycle, element=V.ufl_element()) + bc_tmp = DirichletBC(V, motion_mapping, boundary, id_wall) + bc_tmp.apply(x_["u0"]) + x_["u0"].zero() + + wall_counter_max = motion_mapping.counter + wall_motion = motion_mapping.get_motion_mapping() + + # Remove explicitly from memory. + del motion_mapping + del points + + if rank == 0: + print("Creating wall boundary conditions") + + for i, coord in enumerate(coords): + wall_ = WallMotion(t, wall_motion, wall_counter_max, i, cardiac_cycle, element=V.ufl_element()) + NS_expressions["wall_%s" % coord] = wall_ + + bc_wall = [DirichletBC(V, NS_expressions["wall_%s" % coord], boundary, id_wall) for coord in coords] + else: + # No slip on walls + bc_wall = [DirichletBC(V, Constant(0.0), boundary, id_wall)] * len(coords) + + # Create lists with all velocity boundary conditions + bc_u0 = [bc_inlets[ID][0] for ID in id_in] + [bc_wall[0]] + bc_u1 = [bc_inlets[ID][1] for ID in id_in] + [bc_wall[1]] + bc_u2 = [bc_inlets[ID][2] for ID in id_in] + [bc_wall[2]] + + if track_blood: + bc_blood = [DirichletBC(V, Constant(0.0), boundary, ID) for ID in id_in] + + return dict(u0=bc_u0, u1=bc_u1, u2=bc_u2, p=bc_p, blood=bc_blood) + else: + return dict(u0=bc_u0, u1=bc_u1, u2=bc_u2, p=bc_p) + + +def pre_solve_hook(u_components, id_in, id_out, dynamic_mesh, V, Q, cardiac_cycle, dt, + save_solution_after_cycle, mesh_path, mesh, newfolder, velocity_degree, + restart_folder, **NS_namespace): + id_wall = min(id_in + id_out) - 1 + # Extract diameter at mitral valve + if mesh_path.endswith(".xml") or mesh_path.endswith(".xml.gz"): + mesh_filename = ".xml" + elif mesh_path.endswith(".h5"): + mesh_filename = ".h5" + + info_path = mesh_path.split(mesh_filename)[0] + "_info.json" + with open(info_path) as f: + info = json.load(f) + + outlet_area = info['outlet_area'] + D_mitral = np.sqrt(4 * outlet_area / np.pi) + + # Create point for evaluation + n = FacetNormal(mesh) + eval_dict = {} + rel_path = mesh_path.split(mesh_filename)[0] + "_probe_point.json" + with open(rel_path, 'r') as infile: + probe_points = np.array(json.load(infile)) + + # Store points file in checkpoint + if MPI.rank(MPI.comm_world) == 0: + probe_points.dump(path.join(newfolder, "Checkpoint", "points")) + + eval_dict["centerline_u_x_probes"] = Probes(probe_points.flatten(), V) + eval_dict["centerline_u_y_probes"] = Probes(probe_points.flatten(), V) + eval_dict["centerline_u_z_probes"] = Probes(probe_points.flatten(), V) + eval_dict["centerline_p_probes"] = Probes(probe_points.flatten(), Q) + + if restart_folder is None: + # Get files to store results + files = get_file_paths(newfolder) + NS_parameters.update(dict(files=files)) + else: + files = NS_namespace["files"] + + # Save mesh as HDF5 file for post-processing + with HDF5File(MPI.comm_world, files["mesh"], "w") as mesh_file: + mesh_file.write(mesh, "mesh") + + # Create Probes path + probes_folder = path.join(newfolder, "Probes") + if MPI.rank(MPI.comm_world) == 0: + if not path.exists(probes_folder): + makedirs(probes_folder) + + # Define xdmf writers + viz_U, viz_blood = get_visualization_writers(newfolder, ['velocity', 'blood']) + + # Extract dof map and coordinates + coordinates = mesh.coordinates() + if velocity_degree == 1: + dof_map = vertex_to_dof_map(V) + else: + dof_map = V.dofmap().entity_dofs(mesh, 0) + + # Create vector function for storing velocity + Vv = VectorFunctionSpace(mesh, "CG", velocity_degree) + U = Function(Vv) + u_mean = Function(Vv) + u_mean0 = Function(V) + u_mean1 = Function(V) + u_mean2 = Function(V) + + # Time step when solutions for post-processing should start being saved + save_solution_at_tstep = int(cardiac_cycle * save_solution_after_cycle / dt) + + # Set mesh equation boundary condition + boundary = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1, mesh.domains()) + if dynamic_mesh: + # Set wall motion BCS + bc_mesh = dict((ui, []) for ui in u_components) + + # Zero wall velocity at inlet and outlet + noslip = Constant(0.0) + bc_out = [DirichletBC(V, noslip, boundary, ID) for ID in id_out] + bcu_in = [DirichletBC(V, noslip, boundary, ID) for ID in id_in] + + # Add wall movement to wall + bcu_wall_x = DirichletBC(V, NS_expressions["wall_x"], boundary, id_wall) + bcu_wall_y = DirichletBC(V, NS_expressions["wall_y"], boundary, id_wall) + bcu_wall_z = DirichletBC(V, NS_expressions["wall_z"], boundary, id_wall) + + bc_mesh["u0"] = [bcu_wall_x] + bc_out + bcu_in + bc_mesh["u1"] = [bcu_wall_y] + bc_out + bcu_in + bc_mesh["u2"] = [bcu_wall_z] + bc_out + bcu_in + else: + bc_mesh = {} + + return dict(outlet_area=outlet_area, id_wall=id_wall, D_mitral=D_mitral, n=n, eval_dict=eval_dict, U=U, + u_mean=u_mean, u_mean0=u_mean0, u_mean1=u_mean1, u_mean2=u_mean2, boundary=boundary, bc_mesh=bc_mesh, + coordinates=coordinates, viz_U=viz_U, save_solution_at_tstep=save_solution_at_tstep, + dof_map=dof_map, probes_folder=probes_folder, viz_blood=viz_blood) + + +def update_boundary_conditions(t, dynamic_mesh, NS_expressions, id_in, **NS_namespace): + # Update inlet condition + for ID in id_in: + for i in [0, 1, 2]: + NS_expressions[f"inlet_{ID}"][i].set_t(t) + + if dynamic_mesh: + # Update wall motion BCs + for coord in ["x", "y", "z"]: + NS_expressions[f"wall_{coord}"].t = t + + +def temporal_hook(mesh, id_wall, id_out, cardiac_cycle, dt, t, save_solution_frequency, u_, id_in, tstep, newfolder, + eval_dict, dump_probe_frequency, p_, save_solution_at_tstep, nu, D_mitral, U, u_mean0, u_mean1, + u_mean2, save_flow_metrics_frequency, save_volume_frequency, save_solution_frequency_xdmf, viz_U, + boundary, outlet_area, probes_folder, q_, viz_blood, track_blood, **NS_namespace): + if tstep % save_volume_frequency == 0: + compute_volume(mesh, t, newfolder) + + if tstep % save_solution_frequency_xdmf == 0: + for i in range(3): + assign(U.sub(i), u_[i]) + + viz_U.write(U, t) + if track_blood: + viz_blood.write(q_['blood'], t) + + if tstep % save_flow_metrics_frequency == 0: + h = mesh.hmin() + compute_flow_quantities(u_, D_mitral, nu, mesh, t, tstep, dt, h, outlet_area, boundary, id_out, id_in, id_wall, + period=cardiac_cycle, newfolder=newfolder, dynamic_mesh=False, write_to_file=True) + + # Sample velocity and pressure in points/probes + eval_dict["centerline_u_x_probes"](u_[0]) + eval_dict["centerline_u_y_probes"](u_[1]) + eval_dict["centerline_u_z_probes"](u_[2]) + eval_dict["centerline_p_probes"](p_) + + # Store sampled velocity and pressure + if tstep % dump_probe_frequency == 0: + dump_probes(eval_dict, newfolder, tstep) + + # Save velocity and pressure for post-processing + if tstep % save_solution_frequency == 0 and tstep >= save_solution_at_tstep: + store_velocity_and_pressure_h5(NS_parameters, U, p_, tstep, u_, u_mean0, u_mean1, u_mean2) + + +# Oasis hook called after the simulation has finished +def theend_hook(u_mean, u_mean0, u_mean1, u_mean2, T, dt, save_solution_at_tstep, save_solution_frequency, + **NS_namespace): + store_u_mean(T, dt, save_solution_at_tstep, save_solution_frequency, u_mean, u_mean0, u_mean1, u_mean2, + NS_parameters) diff --git a/src/vampy/simulation/__init__.py b/src/vampy/simulation/__init__.py index fffa452c..14933f18 100644 --- a/src/vampy/simulation/__init__.py +++ b/src/vampy/simulation/__init__.py @@ -5,4 +5,9 @@ from . import Womersley from . import Probe except ModuleNotFoundError: - print("WARNING: Oasis is not installed, running CFD is not available") \ No newline at end of file + print("WARNING: Oasis is not installed, running CFD is not available") + +try: + from . import MovingAtrium +except ModuleNotFoundError: + print("WARNING: OasisMove is not installed, running moving domain CFD is not available") diff --git a/src/vampy/simulation/simulation_common.py b/src/vampy/simulation/simulation_common.py index 8242fdff..2817e868 100644 --- a/src/vampy/simulation/simulation_common.py +++ b/src/vampy/simulation/simulation_common.py @@ -64,17 +64,17 @@ def print_mesh_information(mesh): if MPI.rank(MPI.comm_world) == 0: print("=== Mesh information ===") - print("X range: {} to {} (delta: {:.4f})".format(min(xmin), max(xmax), max(xmax) - min(xmin))) - print("Y range: {} to {} (delta: {:.4f})".format(min(ymin), max(ymax), max(ymax) - min(ymin))) - print("Z range: {} to {} (delta: {:.4f})".format(min(zmin), max(zmax), max(zmax) - min(zmin))) - print("Number of cells: {}".format(sum(num_cells))) - print("Number of cells per processor: {}".format(int(np.mean(num_cells)))) - print("Number of edges: {}".format(sum(num_edges))) - print("Number of faces: {}".format(sum(num_faces))) - print("Number of facets: {}".format(sum(num_facets))) - print("Number of vertices: {}".format(sum(num_vertices))) - print("Volume: {:.4f}".format(volume)) - print("Number of cells per volume: {:.4f}".format(sum(num_cells) / volume)) + print(f"X range: {min(xmin)} to {max(xmax)} (delta: {max(xmax) - min(xmin):.4f})") + print(f"Y range: {min(ymin)} to {max(ymax)} (delta: {max(ymax) - min(ymin):.4f})") + print(f"Z range: {min(zmin)} to {max(zmax)} (delta: {max(zmax) - min(zmin):.4f})") + print(f"Number of cells: {sum(num_cells)}") + print(f"Number of cells per processor: {int(np.mean(num_cells))}") + print(f"Number of edges: {sum(num_edges)}") + print(f"Number of faces: {sum(num_faces)}") + print(f"Number of facets: {sum(num_facets)}") + print(f"Number of vertices: {sum(num_vertices)}") + print(f"Volume: {volume:.4f}") + print(f"Number of cells per volume: {sum(num_cells) / volume:.4f}") def store_u_mean(T, dt, save_solution_at_tstep, save_solution_frequency, u_mean, u_mean0, u_mean1, u_mean2, @@ -108,3 +108,85 @@ def store_u_mean(T, dt, save_solution_at_tstep, save_solution_frequency, u_mean, # Save u_mean with HDF5File(MPI.comm_world, u_mean_path, "w") as u_mean_file: u_mean_file.write(u_mean, "u_mean") + + +def store_velocity_and_pressure_h5(NS_parameters, U, p_, tstep, u_, u_mean0, u_mean1, u_mean2): + """ + Store the velocity and pressure values to an HDF5 file. + + Args: + NS_parameters (dict): A dictionary containing the parameters for Navier-Stokes equations. + U (Function): A vector function space to assign the velocity components. + p_ (Function): The pressure function. + tstep (int): The current time step. + u_ (List[Function]): A list containing the velocity components. + u_mean0 (Function): The accumulated x-component of the velocity. + u_mean1 (Function): The accumulated y-component of the velocity. + u_mean2 (Function): The accumulated z-component of the velocity. + + Returns: + None + """ + # Assign velocity components to vector solution + for i in range(3): + assign(U.sub(i), u_[i]) + + # Get save paths + p_path = NS_parameters['files']['p'] + u_path = NS_parameters['files']['u'] + file_mode = "w" if not path.exists(p_path) else "a" + + # Save pressure + with HDF5File(MPI.comm_world, p_path, file_mode=file_mode) as viz_p: + viz_p.write(p_, "/pressure", tstep) + + # Save velocity + with HDF5File(MPI.comm_world, u_path, file_mode=file_mode) as viz_u: + viz_u.write(U, "/velocity", tstep) + + # Accumulate velocity + u_mean0.vector().axpy(1, u_[0].vector()) + u_mean1.vector().axpy(1, u_[1].vector()) + u_mean2.vector().axpy(1, u_[2].vector()) + + +def dump_probes(eval_dict, newfolder, tstep): + """ + Save variables along the centerline for CFD simulation diagnostics and light-weight post-processing. + + Args: + eval_dict (dict): Dictionary with probe arrays for velocity components and pressure along the centerline. + newfolder (str): The path to the folder where the probe data will be saved. + tstep (float): The current time step. + + Returns: + None + """ + # Construct the file path for probes + filepath = path.join(newfolder, "Probes") + + # Ensure the directory exists on the master process + if MPI.rank(MPI.comm_world) == 0 and not path.exists(filepath): + makedirs(filepath) + + # Extract probe arrays for each variable + variables = ["centerline_u_x_probes", "centerline_u_y_probes", "centerline_u_z_probes", "centerline_p_probes"] + arrs = {var: eval_dict[var].array() for var in variables} + + # Dump stats on the master process + if MPI.rank(MPI.comm_world) == 0: + filenames = { + "centerline_u_x_probes": f"u_x_{tstep}.probes", + "centerline_u_y_probes": f"u_y_{tstep}.probes", + "centerline_u_z_probes": f"u_z_{tstep}.probes", + "centerline_p_probes": f"p_{tstep}.probes" + } + for var, arr in arrs.items(): + arr.dump(path.join(filepath, filenames[var])) + + # Ensure all processes have dumped data before clearing + MPI.barrier(MPI.comm_world) + + # Clear stats + for var in variables: + eval_dict[var].clear()