Average surface and volume derived quantities in Cylindrical and Spherical

festim/__init__.py

@@ -74,14 +74,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

@@ -30,7 +30,7 @@ def __init__(self, field, surface) -> None:
 
     @property
     def allowed_meshes(self):
-        return ["cartesian"]
+        return ["cartesian", "spherical"]
 
     @property
     def export_unit(self):
@@ -51,3 +51,59 @@ 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
+
+    @property
+    def allowed_meshes(self):
+        return ["cylindrical"]
+
+    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
+
+
+AverageSurfaceSpherical = AverageSurface

festim/exports/derived_quantities/average_volume.py

@@ -19,6 +19,7 @@
             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
@@ -50,3 +51,94 @@
         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
+
+    @property
+    def allowed_meshes(self):
+        return ["cylindrical"]
+
+    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
+
+    @property
+    def allowed_meshes(self):
+        return ["spherical"]
+
+    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/simulation/test_initialise.py

@@ -117,7 +117,6 @@ def test_TXTExport_times_added_to_milestones(tmpdir):
         F.SurfaceFlux(field="solute", surface=1),
         F.TotalVolume(field="solute", volume=1),
         F.TotalSurface(field="solute", surface=1),
-        F.AverageSurface(field="solute", surface=1),
         F.AverageVolume(field="solute", volume=1),
         F.HydrogenFlux(surface=1),
         F.ThermalFlux(surface=1),
@@ -151,6 +150,62 @@ def test_cartesian_and_surface_flux_warning(quantity, sys):
         my_model.initialise()
 
 
+@pytest.mark.parametrize("sys", ["cartesian", "spherical"])
+def test_cylindrical_quantities_warning(sys):
+    """
+    Creates a Simulation object and checks that, if either a cartesian
+    or spherical meshes are given with a AverageSurfaceCylindrical, a warning is raised.
+
+    Args:
+        sys (str): type of the coordinate system
+    """
+    # build
+    my_model = F.Simulation()
+    my_model.mesh = F.MeshFromVertices([1, 2, 3], type=sys)
+    my_model.materials = F.Material(id=1, D_0=1, E_D=0)
+    my_model.T = F.Temperature(100)
+    my_model.dt = F.Stepsize(initial_value=3)
+    my_model.settings = F.Settings(
+        absolute_tolerance=1e-10, relative_tolerance=1e-10, final_time=4
+    )
+
+    derived_quantities = F.DerivedQuantities(
+        [F.AverageSurfaceCylindrical(field="solute", surface=1)]
+    )
+    my_model.exports = [derived_quantities]
+
+    # test
+    with pytest.warns(UserWarning, match=f"may not work as intended for {sys} meshes"):
+        my_model.initialise()
+
+
+def test_spherical_quantities_warning():
+    """
+    Creates a Simulation object and checks that, if a cylindrical
+    mesh is given with a AverageSurfaceSpherical, a warning is raised.
+    """
+    # build
+    my_model = F.Simulation()
+    my_model.mesh = F.MeshFromVertices([1, 2, 3], type="cylindrical")
+    my_model.materials = F.Material(id=1, D_0=1, E_D=0)
+    my_model.T = F.Temperature(100)
+    my_model.dt = F.Stepsize(initial_value=3)
+    my_model.settings = F.Settings(
+        absolute_tolerance=1e-10, relative_tolerance=1e-10, final_time=4
+    )
+
+    derived_quantities = F.DerivedQuantities(
+        [F.AverageSurfaceSpherical(field="solute", surface=1)]
+    )
+    my_model.exports = [derived_quantities]
+
+    # test
+    with pytest.warns(
+        UserWarning, match=f"may not work as intended for cylindrical meshes"
+    ):
+        my_model.initialise()
+
+
 @pytest.mark.parametrize(
     "value",
     [