Merge pull request #939 from jhdark/cyl_and_sph_total_quantities

Total surface and volume derived quantities for cylindrical and spherical meshes

festim/__init__.py

@@ -84,7 +84,16 @@
 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_surface import (
+    TotalSurface,
+    TotalSurfaceCylindrical,
+    TotalSurfaceSpherical,
+)
+from .exports.derived_quantities.total_volume import (
+    TotalVolume,
+    TotalVolumeCylindrical,
+    TotalVolumeSpherical,
+)
 from .exports.derived_quantities.total_volume import TotalVolume
 from .exports.derived_quantities.average_surface import (
     AverageSurface,

festim/exports/derived_quantities/total_surface.py

@@ -1,5 +1,6 @@
 from festim import SurfaceQuantity
 import fenics as f
+import numpy as np
 
 
 class TotalSurface(SurfaceQuantity):
@@ -58,3 +59,166 @@ def title(self):
 
     def compute(self):
         return f.assemble(self.function * self.ds(self.surface))
+
+
+class TotalSurfaceCylindrical(TotalSurface):
+    """
+    Computes the total value of a field on a given surface
+    int(f 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
+        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")
+        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
+
+    .. note::
+        Units are in H/m in 1D, H in 2D domains for hydrogen concentration
+        and K m in 1D, K m2 in 2D domains for temperature
+    """
+
+    def __init__(self, field, surface, azimuth_range=(0, 2 * np.pi)) -> None:
+        super().__init__(field=field, surface=surface)
+        self.r = None
+        self.azimuth_range = azimuth_range
+
+    @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(" m1", " m")
+        else:
+            return f"H m{dim-2}".replace(" m0", "")
+
+    @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):
+
+        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
+
+        tot_surf = f.assemble(self.function * self.r * self.ds(self.surface))
+        tot_surf *= self.azimuth_range[1] - self.azimuth_range[0]
+
+        return tot_surf
+
+
+class TotalSurfaceSpherical(TotalSurface):
+    """
+    Computes the total value of a field on a given surface
+    int(f ds)
+    ds is the surface measure in spherical coordinates.
+    ds = r**2 sin(theta) dtheta dphi
+
+    Args:
+        field (str, int):  the field ("solute", 0, 1, "T", "retention")
+        surface (int): the surface id
+        azimuth_range (tuple, optional): Range of the azimuthal angle
+            (phi) needs to be between 0 and 2 pi. Defaults to (0, 2 * np.pi)
+        polar_range (tuple, optional): Range of the polar angle
+            (theta) needs to be between 0 and pi. Defaults to (0, np.pi).
+
+    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
+
+    .. note::
+        Units are in H for hydrogen concentration
+        and K in 1D, K m in 2D domains for temperature
+    """
+
+    def __init__(
+        self, field, surface, azimuth_range=(0, 2 * np.pi), polar_range=(0, np.pi)
+    ) -> None:
+        super().__init__(field=field, surface=surface)
+        self.r = None
+        self.azimuth_range = azimuth_range
+        self.polar_range = polar_range
+
+    @property
+    def export_unit(self):
+        if self.field == "T":
+            return f"K m2"
+        else:
+            return "H"
+
+    @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 2 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] < 0 or value[1] > np.pi:
+            raise ValueError("Polar range must be between 0 and pi")
+        self._polar_range = value
+
+    @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
+
+        tot_surf = f.assemble(self.function * self.r**2 * self.ds(self.surface))
+
+        tot_surf *= (self.azimuth_range[1] - self.azimuth_range[0]) * (
+            np.cos(self.polar_range[0]) - np.cos(self.polar_range[1])
+        )
+
+        return tot_surf

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,163 @@ 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 for a given volume
+    int(f dx)
+    dx is the volume measure in cylindrical coordinates.
+    dx = r dr dtheta dz
+
+    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
+        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
+
+    .. note::
+        Units are in H/m in 1D and H in 2D for hydrogen concentration
+        and K m2 in 1D, K m3 in 2D domains for temperature
+    """
+
+    def __init__(self, field, volume, azimuth_range=(0, 2 * np.pi)) -> None:
+        super().__init__(field=field, volume=volume)
+        self.r = None
+        self.azimuth_range = azimuth_range
+
+    @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+1}"
+        else:
+            return f"H m{dim-2}".replace(" m0", "")
+
+    @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 2 pi")
+        self._azimuth_range = value
+
+    @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
+
+        tot_vol = f.assemble(self.function * self.r * self.dx(self.volume))
+        tot_vol *= self.azimuth_range[1] - self.azimuth_range[0]
+
+        return tot_vol
+
+
+class TotalVolumeSpherical(TotalVolume):
+    """Computes the total value of a field for a given volume
+    int(f dx)
+    dx is the volume measure in cylindrical coordinates.
+    dx = r**2 sin(theta) dtheta dphi dr
+
+    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 2 pi. Defaults to (0, 2 * np.pi)
+        polar_range (tuple, optional): Range of the polar angle
+            (theta) needs to be between 0 and pi. Defaults to (0, np.pi).
+
+    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 cylinder
+
+    .. note::
+        Units are in H for hydrogen concentration and K m2 for temperature
+    """
+
+    def __init__(
+        self, field, volume, azimuth_range=(0, 2 * np.pi), polar_range=(0, np.pi)
+    ) -> None:
+        super().__init__(field=field, volume=volume)
+        self.r = None
+        self.azimuth_range = azimuth_range
+        self.polar_range = polar_range
+
+    @property
+    def export_unit(self):
+        if self.field == "T":
+            return f"K m3"
+        else:
+            return f"H"
+
+    @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 polar_range(self):
+        return self._polar_range
+
+    @polar_range.setter
+    def polar_range(self, value):
+        if value[0] < 0 or value[1] > np.pi:
+            raise ValueError("Polar range must be between 0 and pi")
+        self._polar_range = value
+
+    @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
+
+        tot_vol = f.assemble(self.function * self.r**2 * self.dx(self.volume))
+
+        tot_vol *= (self.azimuth_range[1] - self.azimuth_range[0]) * (
+            np.cos(self.polar_range[0]) - np.cos(self.polar_range[1])
+        )
+
+        return tot_vol