Merge pull request #629 from jhdark/surface_flux_export

Surface Quantities: Surface Flux

festim/__init__.py

@@ -37,6 +37,8 @@
 
 from .stepsize import Stepsize
 
+from .exports.surface_quantity import SurfaceQuantity
+from .exports.surface_flux import SurfaceFlux
 from .exports.vtx import VTXExport
 from .exports.xdmf import XDMFExport
 

festim/exports/surface_flux.py

@@ -0,0 +1,42 @@
+from dolfinx import fem
+import festim as F
+import ufl
+
+
+class SurfaceFlux(F.SurfaceQuantity):
+    """Exports surface flux at a given subdomain
+
+    Args:
+        field (festim.Species): species for which the surface flux is computed
+        surface (festim.SurfaceSubdomain1D): surface subdomain
+        filename (str, optional): name of the file to which the surface flux is exported
+
+    Attributes:
+        field (festim.Species): species for which the surface flux is computed
+        surface (festim.SurfaceSubdomain1D): surface subdomain
+        filename (str): name of the file to which the surface flux is exported
+    """
+
+    def __init__(
+        self,
+        field: F.Species,
+        surface: F.SurfaceSubdomain1D,
+        filename: str = None,
+    ) -> None:
+        super().__init__(field, surface, filename)
+
+    def compute(self, n, ds):
+        """Computes the value of the surface flux at the surface
+
+        Args:
+            n (ufl.geometry.FacetNormal): normal vector to the surface
+            ds (ufl.Measure): surface measure of the model
+        """
+        self.value = fem.assemble_scalar(
+            fem.form(
+                -self.D
+                * ufl.dot(ufl.grad(self.field.solution), n)
+                * ds(self.surface.id)
+            )
+        )
+        self.data.append(self.value)

festim/exports/surface_quantity.py

@@ -0,0 +1,82 @@
+import festim as F
+import csv
+import os
+
+
+class SurfaceQuantity:
+    """Export SurfaceQuantity
+
+    Args:
+        field (festim.Species): species for which the surface flux is computed
+        surface (festim.SurfaceSubdomain1D): surface subdomain
+        filename (str, optional): name of the file to which the surface flux is exported
+
+    Attributes:
+        field (festim.Species): species for which the surface flux is computed
+        surface (festim.SurfaceSubdomain1D): surface subdomain
+        filename (str): name of the file to which the surface flux is exported
+        t (list): list of time values
+        data (list): list of values of the surface quantity
+    """
+
+    def __init__(self, field, surface, filename: str = None) -> None:
+        self.field = field
+        self.surface = surface
+        self.filename = filename
+
+        self.t = []
+        self.data = []
+
+    @property
+    def filename(self):
+        return self._filename
+
+    @filename.setter
+    def filename(self, value):
+        if value is None:
+            self._filename = None
+        elif not isinstance(value, str):
+            raise TypeError("filename must be of type str")
+        elif not value.endswith(".csv") and not value.endswith(".txt"):
+            raise ValueError("filename must end with .csv or .txt")
+        self._filename = value
+
+    @property
+    def surface(self):
+        return self._surface
+
+    @surface.setter
+    def surface(self, value):
+        if not isinstance(value, (int, F.SurfaceSubdomain1D)) or isinstance(
+            value, bool
+        ):
+            raise TypeError("surface should be an int or F.SurfaceSubdomain1D")
+        self._surface = value
+
+    @property
+    def field(self):
+        return self._field
+
+    @field.setter
+    def field(self, value):
+        # check that field is festim.Species
+        if not isinstance(value, (F.Species, str)):
+            raise TypeError("field must be of type festim.Species")
+
+        self._field = value
+
+    def write(self, t):
+        """If the filename doesnt exist yet, create it and write the header,
+        then append the time and value to the file"""
+
+        if not os.path.isfile(self.filename):
+            title = "Flux surface {}: {}".format(self.surface.id, self.field.name)
+
+            if self.filename is not None:
+                with open(self.filename, mode="w", newline="") as file:
+                    writer = csv.writer(file)
+                    writer.writerow(["t(s)", f"{title}"])
+
+        with open(self.filename, mode="a", newline="") as file:
+            writer = csv.writer(file)
+            writer.writerow([t, self.value])