Coverage for physioblocks / library / blocks / cavity.py: 100%
43 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-09 16:40 +0100
1# SPDX-FileCopyrightText: Copyright INRIA
2#
3# SPDX-License-Identifier: LGPL-3.0-only
4#
5# Copyright INRIA
6#
7# This file is part of PhysioBlocks, a library mostly developed by the
8# [Ananke project-team](https://team.inria.fr/ananke) at INRIA.
9#
10# Authors:
11# - Colin Drieu
12# - Dominique Chapelle
13# - François Kimmig
14# - Philippe Moireau
15#
16# PhysioBlocks is free software: you can redistribute it and/or modify it under the
17# terms of the GNU Lesser General Public License as published by the Free Software
18# Foundation, version 3 of the License.
19#
20# PhysioBlocks is distributed in the hope that it will be useful, but WITHOUT ANY
21# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
22# PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
23#
24# You should have received a copy of the GNU Lesser General Public License along with
25# PhysioBlocks. If not, see <https://www.gnu.org/licenses/>.
27"""
28Describes cavity blocks
29"""
31from dataclasses import dataclass
32from typing import Any
34import numpy as np
36from physioblocks.computing import Block, Expression, Quantity
37from physioblocks.registers import register_type
38from physioblocks.simulation import Time
40# Local id for the volume saved quantity
41CAVITY_VOLUME_LOCAL_ID = "volume"
43# Local id for the pressure dof in the cavity
44CAVITY_PRESSURE_LOCAL_ID = "pressure"
46# Local id for the pressure dof in the cavity
47CAVITY_DISP_LOCAL_ID = "disp"
49# Constant for the cavity block type id
50CAVITY_BLOCK_TYPE_ID = "spherical_cavity_block"
53@register_type(CAVITY_BLOCK_TYPE_ID)
54@dataclass
55class SphericalCavityBlock(Block):
56 r"""
57 Spherical cavity block implementation
59 .. tikz:: Spherical Cavity Scheme
61 \filldraw [fill=gray!20] (0, 0) circle [x radius = 3.0+0.2, y radius = 3.0+0.2];
62 \draw [thick, dashed] (0, 0) circle [x radius = 3.0, y radius = 3.0];
63 \filldraw [fill=red!25] (0, 0) circle [x radius = 3.0-0.2, y radius = 3.0-0.2];
64 \filldraw [fill=gray!40] (0, 0) circle [x radius = 1.8 + 0.5, y radius = 1.8 + 0.5];
65 \draw [dashed, thick] (0, 0) circle [x radius = 1.8, y radius = 1.8];
66 \filldraw [fill=red!25] (0, 0) circle [x radius = 1.8 - 0.5, y radius = 1.8 - 0.5];
67 \draw [thick][<->] (0, 0) -- node[above]{$R_0$} (-1.8, 0);
68 \draw [thick][<->] (0, 1.8-0.5) -- node[anchor=west, yshift=4]{$d_0$} (0, 1.8+0.5);
69 \draw [thick][dashed][<-] (-3.0, 0) -- node[anchor=north]{$y$} (-1.8, 0);
70 \draw [thick][<->] (0, 3.0-0.2) -- node[anchor=south, yshift=5]{$d(y)$} (0, 3.0+0.2);
71 \draw (3.0+0.2, 0) to[short, -*, i=$Q$] (5, 0) node[above]{$P$};
73 **Node 1:**
75 :math:`Q = - \frac {dV(y)}{dt}`
77 **Discretised:**
79 :math:`Q^{n + \frac{1}{2}} = - \frac {V(y^{n + 1}) - V(y^{n})}{\Delta t^n}`
81 .. note::
83 Notice that no dynamics is given on the **DOF** in the block.
84 **Model components** of the cavity can access the pressure and give it a dynamics.
86 """ # noqa: E501
88 disp: Quantity[np.float64]
89 """Displacement y"""
91 radius: Quantity[np.float64]
92 """Initial Sphere radius R0"""
94 thickness: Quantity[np.float64]
95 """Initial thickness d0"""
97 time: Time
98 """time"""
100 def initialize(self) -> None:
101 """
102 Initialize block attributes from current quantity values
103 """
104 self.inv_radius = 1.0 / self.radius.current
105 self.sphere_volume = (4.0 / 3.0) * np.pi * np.pow(self.radius.current, 3)
106 self.sphere_surface = 4.0 * np.pi * np.pow(self.radius.current, 2)
107 self.thickness_radius_ratio = self.thickness.current * self.inv_radius
108 self.half_thickness_radius_ratio = 0.5 * self.thickness_radius_ratio
110 def fluid_volume_current(self) -> Any:
111 """
112 Compute the current fluid volume in the cavity.
114 :return: the fluid volume
115 :rtype: np.float64
116 """
117 disp_cur_ratio = 1.0 + self.disp.current * self.inv_radius
118 return self.sphere_volume * np.pow(
119 disp_cur_ratio
120 - np.pow(disp_cur_ratio, -2) * self.half_thickness_radius_ratio,
121 3,
122 )
124 def fluid_volume_new(self) -> Any:
125 """
126 Compute the new fluid volume in the cavity.
128 :return: the fluid volume
129 :rtype: np.float64
130 """
131 disp_new_ratio = 1.0 + self.disp.new * self.inv_radius
132 return self.sphere_volume * np.pow(
133 disp_new_ratio
134 - np.pow(disp_new_ratio, -2) * self.half_thickness_radius_ratio,
135 3,
136 )
138 def cavity_flux(self) -> Any:
139 """
140 Compute the ventricule flux at local node 1.
142 :return: the flux
143 :rtype: np.float64
144 """
146 cavity_flux = -(
147 (self.fluid_volume_new() - self.fluid_volume_current()) * self.time.inv_dt
148 )
150 return cavity_flux
152 def dcavity_flux_ddisp(self) -> Any:
153 """
154 Compute the partial derivative of the cavity flux for ``disp``
156 :return: the flux partial derivative
157 :rtype: np.float64
158 """
160 disp_new_ratio = 1.0 + self.disp.new * self.inv_radius
162 return -(
163 self.time.inv_dt
164 * self.sphere_surface
165 * np.pow(
166 disp_new_ratio
167 - np.pow(disp_new_ratio, -2) * self.half_thickness_radius_ratio,
168 2,
169 )
170 * (1.0 + np.pow(disp_new_ratio, -3) * self.thickness_radius_ratio)
171 )
174# Define the cavity block flux expression.
175_cavity_block_flux_expr = Expression(
176 1,
177 SphericalCavityBlock.cavity_flux,
178 {CAVITY_DISP_LOCAL_ID: SphericalCavityBlock.dcavity_flux_ddisp},
179)
181SphericalCavityBlock.declares_flux_expression(
182 1, CAVITY_PRESSURE_LOCAL_ID, _cavity_block_flux_expr
183)
186# Define the cavity block volume of fluid expression.
187_cavity_block_fluid_volume_expr = Expression(
188 1, SphericalCavityBlock.fluid_volume_current
189)
190SphericalCavityBlock.declares_saved_quantity_expression(
191 CAVITY_VOLUME_LOCAL_ID, _cavity_block_fluid_volume_expr
192)