Total Volume for Cylindrical and Spherical meshes

festim/__init__.py

@@ -81,7 +81,11 @@
 from .exports.derived_quantities.minimum_surface import MinimumSurface
 from .exports.derived_quantities.maximum_surface import MaximumSurface
 from .exports.derived_quantities.total_surface import TotalSurface
-from .exports.derived_quantities.total_volume import TotalVolume
+from .exports.derived_quantities.total_volume import (
+    TotalVolume,
+    TotalVolumeCylindrical,
+    TotalVolumeSpherical,
+)
 from .exports.derived_quantities.average_surface import AverageSurface
 from .exports.derived_quantities.point_value import PointValue
 from .exports.derived_quantities.adsorbed_hydrogen import AdsorbedHydrogen

festim/exports/derived_quantities/total_volume.py

         return f"H m{dim-3}".replace(" m0", "") 
 @property 
 def export_unit(self): 
     # obtain domain dimension 
     try: 
         dim = self.function.function_space().mesh().topology().dim() 
     except AttributeError: 
         dim = self.dx._domain._topological_dimension 
         # TODO we could simply do that all the time 
     # return unit depending on field and dimension of domain 
     if self.field == "T": 
         return f"K m{dim}".replace("1", "") 
     else: 
         return f"H m{dim-3}".replace(" m0", "") 
@@ -1,5 +1,6 @@
 from festim import VolumeQuantity
 import fenics as f
+import numpy as np
 
 
 class TotalVolume(VolumeQuantity):
@@ -58,3 +59,128 @@
 
     def compute(self):
         return f.assemble(self.function * self.dx(self.volume))
+
+
+class TotalVolumeCylindrical(TotalVolume):
+    """
+    Computes the total value of a field in a given volume
+    int(f r dx)
+
+    Args:
+        field (str, int): the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+        azimuth_range (tuple, optional): Range of the azimuthal angle
+            (theta) needs to be between 0 and 2 pi. Defaults to (0, 2 * np.pi).
+
+    Attributes:
+        field (str, int): the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+        export_unit (str): the unit of the derived quantity for exporting
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function of
+            the hydrogen solute field
+
+    .. note::
+        units are in H/m2 in 1D, H/m in 2D and H in 3D domains for hydrogen
+        concentration and K m in 1D, K m2 in 2D and K m3 in 3D domains for temperature
+
+    """
+
+    def __init__(self, field, volume, azimuth_range=(0, 2 * np.pi)):
+        super().__init__(field=field, volume=volume)
+        self.r = None
+        self.azimuth_range = azimuth_range
+
+    @property
+    def azimuth_range(self):
+        return self._azimuth_range
+
+    @azimuth_range.setter
+    def azimuth_range(self, value):
+        if value[0] < 0 or value[1] > 2 * np.pi:
+            raise ValueError("Azimuthal range must be between 0 and pi")
+        self._azimuth_range = value
+
+    @property
+    def allowed_meshes(self):
+        return ["cylindrical"]
+
+    def compute(self):
+        theta_0, theta_1 = self.azimuth_range
+        return (theta_1 - theta_0) * f.assemble(
+            self.function * self.r * self.dx(self.volume)
+        )
+
+
+class TotalVolumeSpherical(TotalVolume):
+    """
+    Computes the total value of a field in a given volume
+    int(f r dx)
+
+    Args:
+        field (str, int): the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+        azimuth_range (tuple, optional): Range of the azimuthal angle
+            (phi) needs to be between 0 and pi. Defaults to (0, np.pi).
+        polar_range (tuple, optional): Range of the polar angle
+            (theta) needs to be between - pi and pi. Defaults to (-np.pi, np.pi).
+
+    Attributes:
+        field (str, int): the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+        export_unit (str): the unit of the derived quantity for exporting
+        title (str): the title of the derived quantity
+        show_units (bool): show the units in the title in the derived quantities
+            file
+        function (dolfin.function.function.Function): the solution function of
+            the hydrogen solute field
+
+    .. note::
+        units are in H/m2 in 1D, H/m in 2D and H in 3D domains for hydrogen
+        concentration and K m in 1D, K m2 in 2D and K m3 in 3D domains for temperature
+
+    """
+
+    def __init__(
+        self, field, volume, azimuth_range=(0, np.pi), polar_range=(-np.pi, np.pi)
+    ):
+        super().__init__(field=field, volume=volume)
+        self.r = None
+        self.azimuth_range = azimuth_range
+        self.polar_range = polar_range
+
+    @property
+    def azimuth_range(self):
+        return self._azimuth_range
+
+    @azimuth_range.setter
+    def azimuth_range(self, value):
+        if value[0] < 0 or value[1] > np.pi:
+            raise ValueError("Azimuthal range must be between 0 and pi")
+        self._azimuth_range = value
+
+    @property
+    def polar_range(self):
+        return self._polar_range
+
+    @polar_range.setter
+    def polar_range(self, value):
+        if value[0] < -np.pi or value[1] > np.pi:
+            raise ValueError("Polar range must be between - pi and pi")
+        self._polar_range = value
+
+    @property
+    def allowed_meshes(self):
+        return ["spherical"]
+
+    def compute(self):
+        phi_0, phi_1 = self.azimuth_range
+        theta_0, theta_1 = self.polar_range
+
+        return (
+            (phi_1 - phi_0)
+            * (theta_1 - theta_0)
+            * f.assemble(self.function * self.r**2 * self.dx(self.volume))
+        )

test/unit/test_exports/test_derived_quantities/test_total_volume.py

@@ -1,8 +1,9 @@
-from festim import TotalVolume
+from festim import TotalVolume, TotalVolumeCylindrical, TotalVolumeSpherical
 import fenics as f
 import pytest
 from .tools import c_1D, c_2D, c_3D
 import pytest
+import numpy as np
 
 
 @pytest.mark.parametrize("field,volume", [("solute", 1), ("T", 2)])
@@ -46,6 +47,127 @@ def test_minimum(self):
         assert produced == expected
 
 
+@pytest.mark.parametrize("radius", [2, 3])
+@pytest.mark.parametrize("r0", [0, 2])
+@pytest.mark.parametrize("height", [2, 3])
+@pytest.mark.parametrize("azimuth_range", [(0, np.pi), (0, np.pi / 2)])
+def test_compute_cylindrical(r0, radius, height, azimuth_range):
+    """
+    Test that TotalVolumeCylindrical computes the volume correctly on a cylinder
+
+    Args:
+        r0 (float): internal radius
+        radius (float): cylinder radius
+        height (float): cylinder height
+        azimuth_range (tuple): range of the azimuthal angle
+    """
+    # creating a mesh with FEniCS
+    r1 = r0 + radius
+    z0 = 0
+    z1 = z0 + height
+    mesh_fenics = f.RectangleMesh(f.Point(r0, z0), f.Point(r1, z1), 10, 10)
+
+    # marking physical groups (volumes and surfaces)
+    volume_markers = f.MeshFunction("size_t", mesh_fenics, mesh_fenics.topology().dim())
+    volume_markers.set_all(1)
+
+    volume = 1
+    my_total = TotalVolumeCylindrical("solute", volume, azimuth_range)
+
+    dx = f.Measure("dx", domain=mesh_fenics, subdomain_data=volume_markers)
+
+    V = f.FunctionSpace(mesh_fenics, "P", 1)
+    r = f.interpolate(f.Expression("x[0]", degree=1), V)
+    c = 2 * r
+
+    my_total.dx = dx
+    my_total.function = c
+    my_total.r = f.Expression("x[0]", degree=1)
+
+    az0, az1 = azimuth_range
+
+    expected_value = f.assemble((az1 - az0) * c * r * dx(volume))
+    computed_value = my_total.compute()
+
+    assert np.isclose(expected_value, computed_value)
+
+
+@pytest.mark.parametrize(
+    "azimuth_range", [(-1, np.pi), (0, 3 * np.pi), (-1, 3 * np.pi)]
+)
+def test_azimuthal_range_cylindrical(azimuth_range):
+    """
+    Tests that an error is raised when the azimuthal range is out of bounds
+    """
+    with pytest.raises(ValueError):
+        TotalVolumeCylindrical("solute", 1, azimuth_range=azimuth_range)
+
+
+@pytest.mark.parametrize("radius", [2, 3])
+@pytest.mark.parametrize("r0", [0, 2])
+@pytest.mark.parametrize("azimuth_range", [(0, np.pi), (0, np.pi / 2)])
+@pytest.mark.parametrize("polar_range", [(-np.pi, np.pi), (np.pi / 2, -np.pi / 2)])
+def test_compute_spherical(r0, radius, azimuth_range, polar_range):
+    """
+    Test that TotalVolumeSpherical computes the volume correctly on a spherical mesh
+
+    Args:
+        r0 (float): internal radius
+        radius (float): sphere radius
+        azimuth_range (tuple): range of the azimuthal angle
+    """
+    # creating a mesh with FEniCS
+    r1 = r0 + radius
+    mesh_fenics = f.IntervalMesh(100, r0, r1)
+
+    # marking physical groups (volumes and surfaces)
+    volume_markers = f.MeshFunction("size_t", mesh_fenics, mesh_fenics.topology().dim())
+    volume_markers.set_all(1)
+
+    volume = 1
+    my_total = TotalVolumeSpherical("solute", volume, azimuth_range, polar_range)
+
+    dx = f.Measure("dx", domain=mesh_fenics, subdomain_data=volume_markers)
+
+    V = f.FunctionSpace(mesh_fenics, "P", 1)
+    r = f.interpolate(f.Expression("x[0]", degree=1), V)
+    c = 2 * r
+
+    my_total.dx = dx
+    my_total.function = c
+    my_total.r = f.Expression("x[0]", degree=1)
+
+    az0, az1 = azimuth_range
+    th0, th1 = polar_range
+
+    expected_value = f.assemble((az1 - az0) * (th1 - th0) * c * r**2 * dx(volume))
+    computed_value = my_total.compute()
+
+    assert np.isclose(expected_value, computed_value)
+
+
+@pytest.mark.parametrize(
+    "azimuth_range", [(-1, np.pi), (0, 3 * np.pi), (-1, 3 * np.pi)]
+)
+def test_azimuthal_range_spherical(azimuth_range):
+    """
+    Tests that an error is raised when the azimuthal range is out of bounds
+    """
+    with pytest.raises(ValueError):
+        TotalVolumeSpherical("solute", 1, azimuth_range=azimuth_range)
+
+
+@pytest.mark.parametrize(
+    "polar_range", [(0, 2 * np.pi), (-2 * np.pi, 0), (-2 * np.pi, 3 * np.pi)]
+)
+def test_polar_range_spherical(polar_range):
+    """
+    Tests that an error is raised when the polar range is out of bounds
+    """
+    with pytest.raises(ValueError):
+        TotalVolumeSpherical("solute", 1, polar_range=polar_range)
+
+
 @pytest.mark.parametrize(
     "function, field, expected_title",
     [