diff --git a/festim/__init__.py b/festim/__init__.py index 10ce0105c..422b6f7bf 100644 --- a/festim/__init__.py +++ b/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 diff --git a/festim/exports/derived_quantities/total_volume.py b/festim/exports/derived_quantities/total_volume.py index fd770dde1..fda231981 100644 --- a/festim/exports/derived_quantities/total_volume.py +++ b/festim/exports/derived_quantities/total_volume.py @@ -1,5 +1,6 @@ from festim import VolumeQuantity import fenics as f +import numpy as np class TotalVolume(VolumeQuantity): @@ -58,3 +59,128 @@ def title(self): 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)) + ) diff --git a/test/unit/test_exports/test_derived_quantities/test_total_volume.py b/test/unit/test_exports/test_derived_quantities/test_total_volume.py index 5fd3a8cb3..87eec1602 100644 --- a/test/unit/test_exports/test_derived_quantities/test_total_volume.py +++ b/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", [