Merge pull request #735 from JonathanGSDUFOUR/Add_Quant(Issue638)

DerivedQuantities for fenicsx

festim/__init__.py

@@ -55,6 +55,13 @@
 from .exports.surface_quantity import SurfaceQuantity
 from .exports.volume_quantity import VolumeQuantity
 from .exports.total_volume import TotalVolume
+from .exports.average_volume import AverageVolume
+from .exports.maximum_volume import MaximumVolume
+from .exports.minimum_volume import MinimumVolume
+from .exports.total_surface import TotalSurface
+from .exports.maximum_surface import MaximumSurface
+from .exports.minimum_surface import MinimumSurface
+from .exports.average_surface import AverageSurface
 from .exports.surface_flux import SurfaceFlux
 from .exports.vtx import VTXExport, VTXExportForTemperature
 from .exports.xdmf import XDMFExport

festim/exports/average_surface.py

@@ -0,0 +1,31 @@
+import festim as F
+from dolfinx import fem
+import numpy as np
+
+
+class AverageSurface(F.SurfaceQuantity):
+    """Computes the average value of a field on a given surface
+
+    Args:
+        field (festim.Species): species for which the average surface is computed
+        surface (festim.SurfaceSubdomain): surface subdomain
+        filename (str, optional): name of the file to which the average surface is exported
+
+    Attributes:
+        see `festim.SurfaceQuantity`
+    """
+
+    @property
+    def title(self):
+        return f"Average {self.field.name} surface {self.surface.id}"
+
+    def compute(self, ds):
+        """
+        Computes the average value of the field on the defined surface
+        subdomain, and appends it to the data list
+        """
+
+        self.value = fem.assemble_scalar(
+            fem.form(self.field.solution * ds(self.surface.id))
+        ) / fem.assemble_scalar(fem.form(1 * ds(self.surface.id)))
+        self.data.append(self.value)

festim/exports/average_volume.py

@@ -0,0 +1,30 @@
+import festim as F
+from dolfinx import fem
+import numpy as np
+
+
+class AverageVolume(F.VolumeQuantity):
+    """Computes the average value of a field in a given volume
+
+    Args:
+        field (festim.Species): species for which the average volume is computed
+        volume (festim.VolumeSubdomain): volume subdomain
+        filename (str, optional): name of the file to which the average volume is exported
+
+    Attributes:
+        see `festim.VolumeQuantity`
+    """
+
+    @property
+    def title(self):
+        return f"Average {self.field.name} volume {self.volume.id}"
+
+    def compute(self, dx):
+        """
+        Computes the average value of solution function within the defined volume
+        subdomain, and appends it to the data list
+        """
+        self.value = fem.assemble_scalar(
+            fem.form(self.field.solution * dx(self.volume.id))
+        ) / fem.assemble_scalar(fem.form(1 * dx(self.volume.id)))
+        self.data.append(self.value)

festim/exports/maximum_surface.py

@@ -0,0 +1,27 @@
+import festim as F
+import numpy as np
+
+
+class MaximumSurface(F.SurfaceQuantity):
+    """Computes the maximum value of a field on a given surface
+
+    Args:
+        field (festim.Species): species for which the maximum surface is computed
+        surface (festim.SurfaceSubdomain): surface subdomain
+        filename (str, optional): name of the file to which the maximum surface is exported
+
+    Attributes:
+        see `festim.SurfaceQuantity`
+    """
+
+    @property
+    def title(self):
+        return f"Maximum {self.field.name} surface {self.surface.id}"
+
+    def compute(self):
+        """
+        Computes the maximum value of the field on the defined surface
+        subdomain, and appends it to the data list
+        """
+        self.value = np.max(self.field.solution.x.array[self.surface.indices])
+        self.data.append(self.value)

festim/exports/maximum_volume.py

@@ -0,0 +1,28 @@
+import festim as F
+import numpy as np
+
+
+class MaximumVolume(F.VolumeQuantity):
+    """Computes the maximum value of a field in a given volume
+
+    Args:
+        field (festim.Species): species for which the maximum volume is computed
+        volume (festim.VolumeSubdomain): volume subdomain
+        filename (str, optional): name of the file to which the maximum volume is exported
+
+    Attributes:
+        see `festim.VolumeQuantity`
+
+    """
+
+    @property
+    def title(self):
+        return f"Maximum {self.field.name} volume {self.volume.id}"
+
+    def compute(self):
+        """
+        Computes the maximum value of solution function within the defined volume
+        subdomain, and appends it to the data list
+        """
+        self.value = np.max(self.field.solution.x.array[self.volume.entities])
+        self.data.append(self.value)

festim/exports/minimum_surface.py

@@ -0,0 +1,27 @@
+import festim as F
+import numpy as np
+
+
+class MinimumSurface(F.SurfaceQuantity):
+    """Computes the minimum value of a field on a given surface
+
+    Args:
+        field (festim.Species): species for which the minimum surface is computed
+        surface (festim.SurfaceSubdomain): surface subdomain
+        filename (str, optional): name of the file to which the minimum surface is exported
+
+    Attributes:
+        see `festim.SurfaceQuantity`
+    """
+
+    @property
+    def title(self):
+        return f"Minimum {self.field.name} surface {self.surface.id}"
+
+    def compute(self):
+        """
+        Computes the minimum value of the field on the defined surface
+        subdomain, and appends it to the data list
+        """
+        self.value = np.min(self.field.solution.x.array[self.surface.indices])
+        self.data.append(self.value)

festim/exports/minimum_volume.py

@@ -0,0 +1,27 @@
+import festim as F
+import numpy as np
+
+
+class MinimumVolume(F.VolumeQuantity):
+    """Computes the minmum value of a field in a given volume
+
+    Args:
+        field (festim.Species): species for which the minmum volume is computed
+        volume (festim.VolumeSubdomain): volume subdomain
+        filename (str, optional): name of the file to which the minmum volume is exported
+
+    Attributes:
+        see `festim.VolumeQuantity`
+    """
+
+    @property
+    def title(self):
+        return f"Minimum {self.field.name} volume {self.volume.id}"
+
+    def compute(self):
+        """
+        Computes the minimum value of solution function within the defined volume
+        subdomain, and appends it to the data list
+        """
+        self.value = np.min(self.field.solution.x.array[self.volume.entities])
+        self.data.append(self.value)

festim/exports/surface_flux.py

@@ -4,29 +4,23 @@
 
 
 class SurfaceFlux(F.SurfaceQuantity):
-    """Exports surface flux at a given subdomain
+    """Computes the flux of a field on a given surface
 
     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
+        see `festim.SurfaceQuantity`
     """
 
-    def __init__(
-        self,
-        field: F.Species,
-        surface: F.SurfaceSubdomain1D,
-        filename: str = None,
-    ) -> None:
-        super().__init__(field, surface, filename)
+    @property
+    def title(self):
+        return f"{self.field.name} flux surface {self.surface.id}"
 
     def compute(self, ds):
-        """Computes the value of the surface flux at the surface
+        """Computes the value of the flux at the surface
 
         Args:
             ds (ufl.Measure): surface measure of the model

festim/exports/surface_quantity.py

@@ -66,12 +66,10 @@ def field(self, 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"""
-        title = "Flux surface {}: {}".format(
-            self.surface.id, self.field.name
-        )  # TODO this should be an attribute of the quantity
+
         if self.filename is not None:
             if self._first_time_export:
-                header = ["t(s)", f"{title}"]
+                header = ["t(s)", f"{self.title}"]
                 with open(self.filename, mode="w+", newline="") as file:
                     writer = csv.writer(file)
                     writer.writerow(header)

festim/exports/total_surface.py

@@ -0,0 +1,33 @@
+import festim as F
+from dolfinx import fem
+import ufl
+
+
+class TotalSurface(F.SurfaceQuantity):
+    """Computes the total value of a field on a given surface
+
+    Args:
+        field (`festim.Species`): species for which the total volume is computed
+        surface (`festim.SurfaceSubdomain`): surface subdomain
+        filename (str, optional): name of the file to which the total volume is exported
+
+    Attributes:
+        see `festim.SurfaceQuantity`
+    """
+
+    @property
+    def title(self):
+        return f"Total {self.field.name} surface {self.surface.id}"
+
+    def compute(self, ds: ufl.Measure):
+        """
+        Computes the total value of the field on the defined surface
+        subdomain, and appends it to the data list
+
+        Args:
+            ds (ufl.Measure): surface measure of the model
+        """
+        self.value = fem.assemble_scalar(
+            fem.form(self.field.solution * ds(self.surface.id))
+        )
+        self.data.append(self.value)

festim/exports/total_volume.py

@@ -4,24 +4,20 @@
 
 
 class TotalVolume(F.VolumeQuantity):
-    """Export TotalVolume
+    """Computes the total value of a field in a given volume
 
     Args:
-        field (`festim.Species`): species for which the total volume is computed
-        volume (`festim.VolumeSubdomain`): volume subdomain
+        field (festim.Species): species for which the total volume is computed
+        volume (festim.VolumeSubdomain): volume subdomain
         filename (str, optional): name of the file to which the total volume is exported
 
     Attributes:
         see `festim.VolumeQuantity`
     """
 
-    def __init__(
-        self,
-        field: F.Species,
-        volume: F.VolumeSubdomain,
-        filename: str = None,
-    ) -> None:
-        super().__init__(field, volume, filename)
+    @property
+    def title(self):
+        return f"Total {self.field.name} volume {self.volume.id}"
 
     def compute(self, dx: ufl.Measure):
         """

festim/exports/volume_quantity.py

@@ -26,6 +26,7 @@ def __init__(self, field, volume, filename: str = None) -> None:
 
         self.t = []
         self.data = []
+        self._first_time_export = True
 
     @property
     def filename(self):
@@ -58,13 +59,13 @@ def write(self, t):
         then append the time and value to the file"""
 
         if not os.path.isfile(self.filename):
-            title = "Total volume {}: {}".format(self.volume.id, self.field.name)
-
             if self.filename is not None:
-                with open(self.filename, mode="w", newline="") as file:
+                if self._first_time_export:
+                    header = ["t(s)", f"{self.title}"]
+                    with open(self.filename, mode="w+", newline="") as file:
+                        writer = csv.writer(file)
+                        writer.writerow(header)
+                    self._first_time_export = False
+                with open(self.filename, mode="a", 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])
+                    writer.writerow([t, self.value])

festim/hydrogen_transport_problem.py

@@ -567,7 +567,6 @@ def create_source_values_fenics(self):
         for source in self.sources:
             # create value_fenics for all F.ParticleSource objects
             if isinstance(source, F.ParticleSource):
-
                 source.create_value_fenics(
                     mesh=self.mesh.mesh,
                     temperature=self.temperature_fenics,
@@ -579,7 +578,6 @@ def create_flux_values_fenics(self):
         for bc in self.boundary_conditions:
             # create value_fenics for all F.ParticleFluxBC objects
             if isinstance(bc, F.ParticleFluxBC):
-
                 bc.create_value_fenics(
                     mesh=self.mesh.mesh,
                     temperature=self.temperature_fenics,
@@ -789,15 +787,25 @@ def post_processing(self):
         for export in self.exports:
             # TODO if export type derived quantity
             if isinstance(export, F.SurfaceQuantity):
-                export.compute(self.ds)
+                if isinstance(
+                    export, (F.SurfaceFlux, F.TotalSurface, F.AverageSurface)
+                ):
+                    export.compute(
+                        self.ds,
+                    )
+                else:
+                    export.compute()
                 # update export data
                 export.t.append(float(self.t))
 
                 # if filename given write export data to file
                 if export.filename is not None:
                     export.write(t=float(self.t))
             elif isinstance(export, F.VolumeQuantity):
-                export.compute(self.dx)
+                if isinstance(export, (F.TotalVolume, F.AverageVolume)):
+                    export.compute(self.dx)
+                else:
+                    export.compute()
                 # update export data
                 export.t.append(float(self.t))
 

festim/subdomain/surface_subdomain_1d.py

@@ -36,7 +36,7 @@ def locate_boundary_facet_indices(self, mesh):
         """
         assert mesh.geometry.dim == 1, "This method is only for 1D meshes"
         fdim = 0
-        indices = dolfinx.mesh.locate_entities_boundary(
+        self.indices = dolfinx.mesh.locate_entities_boundary(
             mesh, fdim, lambda x: np.isclose(x[0], self.x)
         )
-        return indices
+        return self.indices