Average surface and volume derived quantities in Cylindrical and Spherical

festim/__init__.py

@@ -75,14 +75,22 @@
 )
 from .exports.derived_quantities.hydrogen_flux import HydrogenFlux
 from .exports.derived_quantities.thermal_flux import ThermalFlux
-from .exports.derived_quantities.average_volume import AverageVolume
+from .exports.derived_quantities.average_volume import (
+    AverageVolume,
+    AverageVolumeCylindrical,
+    AverageVolumeSpherical,
+)
 from .exports.derived_quantities.maximum_volume import MaximumVolume
 from .exports.derived_quantities.minimum_volume import MinimumVolume
 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.average_surface import AverageSurface
+from .exports.derived_quantities.average_surface import (
+    AverageSurface,
+    AverageSurfaceCylindrical,
+    AverageSurfaceSpherical,
+)
 from .exports.derived_quantities.point_value import PointValue
 from .exports.derived_quantities.adsorbed_hydrogen import AdsorbedHydrogen
 

festim/exports/derived_quantities/average_surface.py

@@ -1,5 +1,6 @@
 from festim import SurfaceQuantity
 import fenics as f
+import numpy as np
 
 
 class AverageSurface(SurfaceQuantity):
@@ -47,3 +48,101 @@ def compute(self):
         return f.assemble(self.function * self.ds(self.surface)) / f.assemble(
             1 * self.ds(self.surface)
         )
+
+
+class AverageSurfaceCylindrical(AverageSurface):
+    """
+    Computes the average value of a field on a given surface
+    int(f ds) / int (1 * ds)
+    ds is the surface measure in cylindrical coordinates.
+    ds = r dr dtheta
+
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+
+    Attributes:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+        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 field
+        r (ufl.indexed.Indexed): the radius of the cylinder
+
+    Notes:
+        Units are in H/m3 for hydrogen concentration and K for temperature
+    """
+
+    def __init__(self, field, surface) -> None:
+        super().__init__(field=field, surface=surface)
+        self.r = None
+
+    def compute(self):
+
+        if self.r is None:
+            mesh = (
+                self.function.function_space().mesh()
+            )  # get the mesh from the function
+            rthetaz = f.SpatialCoordinate(mesh)  # get the coordinates from the mesh
+            self.r = rthetaz[0]  # only care about r here
+
+        # dS_z = r dr dtheta , assuming axisymmetry dS_z = theta r dr
+        # dS_r = r dz dtheta , assuming axisymmetry dS_r = theta r dz
+        # in both cases the expression with self.dx is the same
+
+        avg_surf = f.assemble(
+            self.function * self.r * self.ds(self.surface)
+        ) / f.assemble(1 * self.r * self.ds(self.surface))
+
+        return avg_surf
+
+
+class AverageSurfaceSpherical(AverageSurface):
+    """
+    Computes the average value of a field in a given volume
+    int(f ds) / int (1 * ds)
+    ds is the surface measure in cylindrical coordinates.
+    ds = r^2 sin(theta) dtheta dphi
+
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+
+    Attributes:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+        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 field
+        r (ufl.indexed.Indexed): the radius of the sphere
+
+    Notes:
+        Units are in H/m3 for hydrogen concentration and K for temperature
+    """
+
+    def __init__(self, field, surface) -> None:
+        super().__init__(field=field, surface=surface)
+        self.r = None
+
+    def compute(self):
+
+        if self.r is None:
+            mesh = (
+                self.function.function_space().mesh()
+            )  # get the mesh from the function
+            rthetaz = f.SpatialCoordinate(mesh)  # get the coordinates from the mesh
+            self.r = rthetaz[0]  # only care about r here
+
+        # dV_z = r dr dtheta , assuming axisymmetry dV_z = theta r dr
+        # dV_r = r dz dtheta , assuming axisymmetry dV_r = theta r dz
+        # in both cases the expression with self.dx is the same
+
+        avg_surf = f.assemble(
+            self.function * self.r**2 * self.ds(self.surface)
+        ) / f.assemble(1 * self.r**2 * self.ds(self.surface))
+
+        return avg_surf

festim/exports/derived_quantities/average_volume.py

@@ -1,5 +1,6 @@
 from festim import VolumeQuantity
 import fenics as f
+import numpy as np
 
 
 class AverageVolume(VolumeQuantity):
@@ -19,6 +20,7 @@ class AverageVolume(VolumeQuantity):
             file
         function (dolfin.function.function.Function): the solution function of
             the field
+        r (ufl.indexed.Indexed): the radius of the cylinder
 
     .. note::
         Units are in H/m3 for hydrogen concentration and K for temperature
@@ -46,3 +48,86 @@ def compute(self):
         return f.assemble(self.function * self.dx(self.volume)) / f.assemble(
             1 * self.dx(self.volume)
         )
+
+
+class AverageVolumeCylindrical(AverageVolume):
+    """
+    Computes the average value of a field in a given volume
+    int(f dx) / int (1 * dx)
+    dx is the volume measure in cylindrical coordinates.
+    dx = r dr dz dtheta
+
+    Note: for particle fluxes J is given in H/s, for heat fluxes J is given in W
+
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+
+    Attributes:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+        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 field
+        r (ufl.indexed.Indexed): the radius of the sphere
+    """
+
+    def __init__(self, field, volume) -> None:
+        super().__init__(field=field, volume=volume)
+        self.r = None
+
+    def compute(self):
+
+        if self.r is None:
+            mesh = (
+                self.function.function_space().mesh()
+            )  # get the mesh from the function
+            rthetaz = f.SpatialCoordinate(mesh)  # get the coordinates from the mesh
+            self.r = rthetaz[0]  # only care about r here
+
+        avg_vol = f.assemble(
+            self.function * self.r * self.dx(self.volume)
+        ) / f.assemble(1 * self.r * self.dx(self.volume))
+
+        return avg_vol
+
+
+class AverageVolumeSpherical(AverageVolume):
+    """
+    Computes the average value of a field in a given volume
+    int(f dx) / int (1 * dx)
+    dx is the volume measure in cylindrical coordinates.
+    dx = rho dtheta dphi
+
+    Note: for particle fluxes J is given in H/s, for heat fluxes J is given in W
+
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        volume (int): the volume id
+        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 field
+    """
+
+    def __init__(self, field, volume) -> None:
+        super().__init__(field=field, volume=volume)
+        self.r = None
+
+    def compute(self):
+
+        if self.r is None:
+            mesh = (
+                self.function.function_space().mesh()
+            )  # get the mesh from the function
+            rthetaphi = f.SpatialCoordinate(mesh)  # get the coordinates from the mesh
+            self.r = rthetaphi[0]  # only care about r here
+
+        avg_vol = f.assemble(
+            self.function * self.r**2 * self.dx(self.volume)
+        ) / f.assemble(1 * self.r**2 * self.dx(self.volume))
+
+        return avg_vol

test/unit/test_exports/test_derived_quantities/test_average_surface.py

@@ -1,6 +1,14 @@
-from festim import AverageSurface
+from festim import (
+    AverageSurface,
+    AverageSurfaceCylindrical,
+    AverageSurfaceSpherical,
+    x,
+    y,
+)
 import fenics as f
 import pytest
+import numpy as np
+from sympy.printing import ccode
 
 
 @pytest.mark.parametrize("field, surface", [("solute", 1), ("T", 2)])
@@ -43,3 +51,137 @@ def test_h_average(self):
         )
         computed = self.my_average.compute()
         assert computed == expected
+
+
+@pytest.mark.parametrize("radius", [1, 4])
+@pytest.mark.parametrize("r0", [3, 5])
+@pytest.mark.parametrize("height", [2, 7])
+@pytest.mark.parametrize("c_top", [8, 9])
+@pytest.mark.parametrize("c_bottom", [10, 11])
+def test_compute_cylindrical(r0, radius, height, c_top, c_bottom):
+    """
+    Test that AverageSurfaceCylindrical computes the value correctly on a hollow cylinder
+
+    Args:
+        r0 (float): internal radius
+        radius (float): cylinder radius
+        height (float): cylinder height
+        c_top (float): concentration top
+        c_bottom (float): concentration bottom
+    """
+    # 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)
+
+    top_surface = f.CompiledSubDomain(
+        f"on_boundary && near(x[1], {z1}, tol)", tol=1e-14
+    )
+
+    surface_markers = f.MeshFunction(
+        "size_t", mesh_fenics, mesh_fenics.topology().dim() - 1
+    )
+    surface_markers.set_all(0)
+    ds = f.Measure("ds", domain=mesh_fenics, subdomain_data=surface_markers)
+    # Surface ids
+    top_id = 2
+    top_surface.mark(surface_markers, top_id)
+
+    my_export = AverageSurfaceCylindrical("solute", top_id)
+    V = f.FunctionSpace(mesh_fenics, "P", 1)
+    c_fun = lambda z: c_bottom + (c_top - c_bottom) / (z1 - z0) * z
+    expr = f.Expression(
+        ccode(c_fun(y)),
+        degree=1,
+    )
+    my_export.function = f.interpolate(expr, V)
+    my_export.ds = ds
+
+    expected_value = c_bottom + (c_top - c_bottom) / (z1 - z0) * height
+
+    computed_value = my_export.compute()
+
+    assert np.isclose(computed_value, expected_value)
+
+
+@pytest.mark.parametrize("radius", [2, 4])
+@pytest.mark.parametrize("r0", [3, 5])
+@pytest.mark.parametrize("c_left", [8, 9])
+@pytest.mark.parametrize("c_right", [10, 11])
+def test_compute_spherical(r0, radius, c_left, c_right):
+    """
+    Test that AverageSurfaceSpherical computes the average value correctly
+    on a hollow sphere
+
+    Args:
+        r0 (float): internal radius
+        radius (float): sphere  radius
+        c_left (float): concentration left
+        c_right (float): concentration right
+    """
+    # creating a mesh with FEniCS
+    r1 = r0 + radius
+    mesh_fenics = f.IntervalMesh(10, r0, r1)
+
+    # marking physical groups (volumes and surfaces)
+    right_surface = f.CompiledSubDomain(
+        f"on_boundary && near(x[0], {r1}, tol)", tol=1e-14
+    )
+    surface_markers = f.MeshFunction(
+        "size_t", mesh_fenics, mesh_fenics.topology().dim() - 1
+    )
+    surface_markers.set_all(0)
+    # Surface ids
+    right_id = 2
+    right_surface.mark(surface_markers, right_id)
+    ds = f.Measure("ds", domain=mesh_fenics, subdomain_data=surface_markers)
+
+    my_export = AverageSurfaceSpherical("solute", right_id)
+    V = f.FunctionSpace(mesh_fenics, "P", 1)
+    c_fun = lambda r: c_left + (c_right - c_left) / (r1 - r0) * r
+    expr = f.Expression(
+        ccode(c_fun(x)),
+        degree=1,
+    )
+    my_export.function = f.interpolate(expr, V)
+
+    my_export.ds = ds
+
+    expected_value = c_left + ((c_right - c_left) / (r1 - r0)) * r1
+
+    computed_value = my_export.compute()
+
+    assert np.isclose(computed_value, expected_value)
+
+
+def test_average_surface_cylindrical_title_no_units_solute():
+    """A simple test to check that the title is set correctly in
+    festim.AverageSurfaceCylindrical with a solute field without units"""
+
+    my_export = AverageSurfaceCylindrical("solute", 4)
+    assert my_export.title == "Average solute surface 4"
+
+
+def test_average_surface_cylindrical_title_no_units_temperature():
+    """A simple test to check that the title is set correctly in
+    festim.AverageSurfaceCylindrical with a T field without units"""
+
+    my_export = AverageSurfaceCylindrical("T", 5)
+    assert my_export.title == "Average T surface 5"
+
+
+def test_average_surface_spherical_title_no_units_solute():
+    """A simple test to check that the title is set correctly in
+    festim.AverageSurfaceSpherical with a solute field without units"""
+
+    my_export = AverageSurfaceSpherical("solute", 6)
+    assert my_export.title == "Average solute surface 6"
+
+
+def test_average_surface_spherical_title_no_units_temperature():
+    """A simple test to check that the title is set correctly in
+    festim.AverageSurfaceSpherical with a T field without units"""
+
+    my_export = AverageSurfaceSpherical("T", 9)
+    assert my_export.title == "Average T surface 9"