From fb31540a39fdedd02f7da4af9bb29ef6209c697a Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 25 Jul 2024 23:19:08 +0200 Subject: [PATCH 01/29] first effort --- .../derived_quantities/average_surface.py | 136 +++++++++++++++++ .../derived_quantities/average_volume.py | 141 ++++++++++++++++++ 2 files changed, 277 insertions(+) diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py index a6d994c02..c4e030efb 100644 --- a/festim/exports/derived_quantities/average_surface.py +++ b/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): @@ -43,3 +44,138 @@ 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 + 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 + azimuth_range (tuple, optional): Range of the azimuthal angle + (theta) needs to be between 0 and 2 pi. Defaults to (0, 2 * np.pi) + + Notes: + Units are in H/m3 for hydrogen concentration and K 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 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 + + 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)) + + avg_surf *= self.azimuth_range[1] - self.azimuth_range[0] + + return avg_surf + + +class AverageSurfaceSpherical(AverageSurface): + """ + 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 dtheta + + Args: + 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 + azimuth_range (tuple, optional): Range of the azimuthal angle + (theta) needs to be between 0 and 2 pi. Defaults to (0, 2 * np.pi) + polar_range (tuple, optional): Range of the polar angle + (phi) needs to be between - pi and pi. Defaults to (-np.pi, np.pi) + + Notes: + Units are in H/m3 for hydrogen concentration and K for temperature + """ + + def __init__( + self, field, surface, azimuth_range=(0, 2 * np.pi), polar_range=(-np.pi, np.pi) + ) -> None: + super().__init__(field=field, surface=surface) + self.r = None + self.azimuth_range = azimuth_range + self.polar_range = polar_range + + @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] < -np.pi or value[1] > np.pi: + raise ValueError("Polar range must be between - pi and pi") + self._polar_range = value + + 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_vol = f.assemble( + self.function * self.r**2 * self.dx(self.volume) + ) / f.assemble(1 * self.r**2 * self.dx(self.volume)) + + avg_vol *= (self.polar_range[1] - self.polar_range[0]) * ( + -np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0]) + ) + + return avg_vol diff --git a/festim/exports/derived_quantities/average_volume.py b/festim/exports/derived_quantities/average_volume.py index 9c4140161..7f7dd850b 100644 --- a/festim/exports/derived_quantities/average_volume.py +++ b/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): @@ -42,3 +43,143 @@ 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 + 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 + azimuth_range (tuple, optional): Range of the azimuthal angle + (theta) needs to be between 0 and 2 pi. Defaults to (0, 2 * np.pi). + """ + + def __init__(self, field, volume, z, azimuth_range=(0, 2 * np.pi)) -> None: + super().__init__(field=field, volume=volume) + self.r = None + self.z = z + self.azimuth_range = azimuth_range + + @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 + + 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 + + values = f.assemble(self.function * self.r * self.z * self.dx(self.volume)) * ( + self.azimuth_range[1] - self.azimuth_range[0] + ) + + # volume = f.assemble(1 * self.r * self.z * self.dx(self.volume)) * ( + # self.azimuth_range[1] - self.azimuth_range[0] + # ) + volume = f.assemble(1 * self.r * self.z * self.dx(self.volume)) * ( + self.azimuth_range[1] - self.azimuth_range[0] + ) + + avg_vol = values / volume + avg_vol *= self.azimuth_range[1] - self.azimuth_range[0] + + 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 + azimuth_range (tuple, optional): Range of the azimuthal angle + (theta) needs to be between 0 and pi. Defaults to (0, np.pi) + polar_range (tuple, optional): Range of the polar angle + (phi) needs to be between - pi and pi. Defaults to (-np.pi, np.pi) + """ + + def __init__( + self, field, volume, azimuth_range=(0, np.pi), polar_range=(-np.pi, np.pi) + ) -> None: + super().__init__(field=field, volume=volume) + self.r = None + self.azimuth_range = azimuth_range + self.polar_range = polar_range + + @property + def polar_range(self): + return self._polar_range + + @polar_range.setter + def polar_range(self, value): + if value[0] < -np.pi or value[1] > np.pi: + raise ValueError("Polar range must be between - pi and pi") + self._polar_range = value + + @property + def azimuth_range(self): + return self._azimuth_range + + @azimuth_range.setter + def azimuth_range(self, value): + if value[0] < 0 or value[1] > np.pi: + raise ValueError("Azimuthal range must be between 0 and pi") + self._azimuth_range = value + + 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 + + values = ( + f.assemble(self.function * self.r**2 * self.dx(self.volume)) + * (self.polar_range[1] - self.polar_range[0]) + * (-np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0])) + ) + + volume = ( + f.assemble(1 * self.r**2 * self.dx(self.volume)) + * (self.polar_range[1] - self.polar_range[0]) + * (-np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0])) + ) + + avg_vol = values / volume + + return avg_vol From 5fb6032038addfabecb77619ac49315bf91960f1 Mon Sep 17 00:00:00 2001 From: J Dark Date: Thu, 25 Jul 2024 23:19:27 +0200 Subject: [PATCH 02/29] add new classes --- festim/__init__.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/festim/__init__.py b/festim/__init__.py index bd29c6078..aec818a1f 100644 --- a/festim/__init__.py +++ b/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.derived_quantities import DerivedQuantities From cd7c930d52fd494aca657f8b6c166440241c4824 Mon Sep 17 00:00:00 2001 From: JDark Date: Fri, 26 Jul 2024 14:07:33 +0200 Subject: [PATCH 03/29] remove comments and testing --- festim/exports/derived_quantities/average_volume.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/festim/exports/derived_quantities/average_volume.py b/festim/exports/derived_quantities/average_volume.py index 7f7dd850b..4dc70e1ce 100644 --- a/festim/exports/derived_quantities/average_volume.py +++ b/festim/exports/derived_quantities/average_volume.py @@ -95,15 +95,11 @@ def compute(self): self.azimuth_range[1] - self.azimuth_range[0] ) - # volume = f.assemble(1 * self.r * self.z * self.dx(self.volume)) * ( - # self.azimuth_range[1] - self.azimuth_range[0] - # ) volume = f.assemble(1 * self.r * self.z * self.dx(self.volume)) * ( self.azimuth_range[1] - self.azimuth_range[0] ) avg_vol = values / volume - avg_vol *= self.azimuth_range[1] - self.azimuth_range[0] return avg_vol From e82da3a1ab3af726713546dced38eb3e82ecf67a Mon Sep 17 00:00:00 2001 From: JDark Date: Fri, 26 Jul 2024 14:15:26 +0200 Subject: [PATCH 04/29] fix --- .../derived_quantities/average_surface.py | 29 ++++++++++++------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py index c4e030efb..a16d9bb1e 100644 --- a/festim/exports/derived_quantities/average_surface.py +++ b/festim/exports/derived_quantities/average_surface.py @@ -96,11 +96,14 @@ def compute(self): # 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)) + values = f.assemble(self.function * self.r * self.ds(self.surface)) * ( + self.azimuth_range[1] - self.azimuth_range[0] + ) - avg_surf *= self.azimuth_range[1] - self.azimuth_range[0] + surface_area = f.assemble(1 * self.r * self.ds(self.surface)) * ( + self.azimuth_range[1] - self.azimuth_range[0] + ) + avg_surf = values / surface_area return avg_surf @@ -170,12 +173,18 @@ def compute(self): # dV_r = r dz dtheta , assuming axisymmetry dV_r = theta r dz # in both cases the expression with self.dx is the same - avg_vol = f.assemble( - self.function * self.r**2 * self.dx(self.volume) - ) / f.assemble(1 * self.r**2 * self.dx(self.volume)) + values = ( + f.assemble(self.function * self.r**2 * self.dx(self.volume)) + * (self.polar_range[1] - self.polar_range[0]) + * (-np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0])) + ) - avg_vol *= (self.polar_range[1] - self.polar_range[0]) * ( - -np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0]) + surface_area = ( + f.assemble(1 * self.r**2 * self.dx(self.volume)) + * (self.polar_range[1] - self.polar_range[0]) + * (-np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0])) ) - return avg_vol + avg_surf = values / surface_area + + return avg_surf From 740e8e61b347c35c8b3085d1a36686f01d97f660 Mon Sep 17 00:00:00 2001 From: JDark Date: Fri, 26 Jul 2024 14:18:13 +0200 Subject: [PATCH 05/29] surface not volume --- festim/exports/derived_quantities/average_surface.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py index a16d9bb1e..827df5dac 100644 --- a/festim/exports/derived_quantities/average_surface.py +++ b/festim/exports/derived_quantities/average_surface.py @@ -174,13 +174,13 @@ def compute(self): # in both cases the expression with self.dx is the same values = ( - f.assemble(self.function * self.r**2 * self.dx(self.volume)) + f.assemble(self.function * self.r**2 * self.ds(self.surface)) * (self.polar_range[1] - self.polar_range[0]) * (-np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0])) ) surface_area = ( - f.assemble(1 * self.r**2 * self.dx(self.volume)) + f.assemble(1 * self.r**2 * self.ds(self.surface)) * (self.polar_range[1] - self.polar_range[0]) * (-np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0])) ) From b0ce4fb7f57318de05281bd7db9b73de9fe3632d Mon Sep 17 00:00:00 2001 From: JDark Date: Fri, 26 Jul 2024 14:21:11 +0200 Subject: [PATCH 06/29] fix documentation --- festim/exports/derived_quantities/average_surface.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py index 827df5dac..602625f5c 100644 --- a/festim/exports/derived_quantities/average_surface.py +++ b/festim/exports/derived_quantities/average_surface.py @@ -111,9 +111,9 @@ def compute(self): class AverageSurfaceSpherical(AverageSurface): """ 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 dtheta + int(f dS) / int (1 * dS) + dS is the volume measure in cylindrical coordinates. + dS = r dr dtheta Args: field (str, int): the field ("solute", 0, 1, "T", "retention") From a12c51916c9b68402c0f75e0e39e06df9d23fc2a Mon Sep 17 00:00:00 2001 From: JDark Date: Fri, 26 Jul 2024 14:27:04 +0200 Subject: [PATCH 07/29] fix azimuth range spherical and doc strings --- .../derived_quantities/average_surface.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py index 602625f5c..e39f927ad 100644 --- a/festim/exports/derived_quantities/average_surface.py +++ b/festim/exports/derived_quantities/average_surface.py @@ -49,9 +49,9 @@ def compute(self): 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 + 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") @@ -111,9 +111,9 @@ def compute(self): class AverageSurfaceSpherical(AverageSurface): """ Computes the average value of a field in a given volume - int(f dS) / int (1 * dS) - dS is the volume measure in cylindrical coordinates. - dS = r dr dtheta + 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") @@ -124,16 +124,16 @@ class AverageSurfaceSpherical(AverageSurface): function (dolfin.function.function.Function): the solution function of the field azimuth_range (tuple, optional): Range of the azimuthal angle - (theta) needs to be between 0 and 2 pi. Defaults to (0, 2 * np.pi) + (phi) needs to be between 0 and 2 pi. Defaults to (0, np.pi) polar_range (tuple, optional): Range of the polar angle - (phi) needs to be between - pi and pi. Defaults to (-np.pi, np.pi) + (theta) needs to be between - pi and pi. Defaults to (-np.pi, np.pi) Notes: Units are in H/m3 for hydrogen concentration and K for temperature """ def __init__( - self, field, surface, azimuth_range=(0, 2 * np.pi), polar_range=(-np.pi, np.pi) + self, field, surface, azimuth_range=(0, np.pi), polar_range=(-np.pi, np.pi) ) -> None: super().__init__(field=field, surface=surface) self.r = None From 307b95bc39d07f8b44332b297c31964fdfd363f8 Mon Sep 17 00:00:00 2001 From: JDark Date: Fri, 26 Jul 2024 14:28:33 +0200 Subject: [PATCH 08/29] no more than pi in range --- festim/exports/derived_quantities/average_surface.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py index e39f927ad..c7f51f1ad 100644 --- a/festim/exports/derived_quantities/average_surface.py +++ b/festim/exports/derived_quantities/average_surface.py @@ -146,7 +146,7 @@ def azimuth_range(self): @azimuth_range.setter def azimuth_range(self, value): - if value[0] < 0 or value[1] > 2 * np.pi: + if value[0] < 0 or value[1] > np.pi: raise ValueError("Azimuthal range must be between 0 and pi") self._azimuth_range = value From a83849e29f3cb5d3f6f12c7bf656a7ea60df22cd Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 26 Jul 2024 18:18:28 +0200 Subject: [PATCH 09/29] ranges cancel out surface --- .../derived_quantities/average_surface.py | 72 +++---------------- 1 file changed, 8 insertions(+), 64 deletions(-) diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py index c7f51f1ad..4d47b7101 100644 --- a/festim/exports/derived_quantities/average_surface.py +++ b/festim/exports/derived_quantities/average_surface.py @@ -61,27 +61,14 @@ class AverageSurfaceCylindrical(AverageSurface): file function (dolfin.function.function.Function): the solution function of the field - azimuth_range (tuple, optional): Range of the azimuthal angle - (theta) needs to be between 0 and 2 pi. Defaults to (0, 2 * np.pi) Notes: Units are in H/m3 for hydrogen concentration and K for temperature """ - def __init__(self, field, surface, azimuth_range=(0, 2 * np.pi)) -> None: + def __init__(self, field, surface) -> None: super().__init__(field=field, surface=surface) self.r = None - self.azimuth_range = azimuth_range - - @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 def compute(self): @@ -96,14 +83,9 @@ def compute(self): # dS_r = r dz dtheta , assuming axisymmetry dS_r = theta r dz # in both cases the expression with self.dx is the same - values = f.assemble(self.function * self.r * self.ds(self.surface)) * ( - self.azimuth_range[1] - self.azimuth_range[0] - ) - - surface_area = f.assemble(1 * self.r * self.ds(self.surface)) * ( - self.azimuth_range[1] - self.azimuth_range[0] - ) - avg_surf = values / surface_area + avg_surf = f.assemble( + self.function * self.r * self.ds(self.surface) + ) / f.assemble(1 * self.r * self.ds(self.surface)) return avg_surf @@ -123,42 +105,14 @@ class AverageSurfaceSpherical(AverageSurface): file function (dolfin.function.function.Function): the solution function of the field - azimuth_range (tuple, optional): Range of the azimuthal angle - (phi) needs to be between 0 and 2 pi. Defaults to (0, np.pi) - polar_range (tuple, optional): Range of the polar angle - (theta) needs to be between - pi and pi. Defaults to (-np.pi, np.pi) Notes: Units are in H/m3 for hydrogen concentration and K for temperature """ - def __init__( - self, field, surface, azimuth_range=(0, np.pi), polar_range=(-np.pi, np.pi) - ) -> None: + def __init__(self, field, surface) -> None: super().__init__(field=field, surface=surface) self.r = None - self.azimuth_range = azimuth_range - self.polar_range = polar_range - - @property - def azimuth_range(self): - return self._azimuth_range - - @azimuth_range.setter - def azimuth_range(self, value): - if value[0] < 0 or value[1] > 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] < -np.pi or value[1] > np.pi: - raise ValueError("Polar range must be between - pi and pi") - self._polar_range = value def compute(self): @@ -173,18 +127,8 @@ def compute(self): # dV_r = r dz dtheta , assuming axisymmetry dV_r = theta r dz # in both cases the expression with self.dx is the same - values = ( - f.assemble(self.function * self.r**2 * self.ds(self.surface)) - * (self.polar_range[1] - self.polar_range[0]) - * (-np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0])) - ) - - surface_area = ( - f.assemble(1 * self.r**2 * self.ds(self.surface)) - * (self.polar_range[1] - self.polar_range[0]) - * (-np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0])) - ) - - avg_surf = values / surface_area + 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 From 4bc4a6a41d754570abddf19c35cfa0597f4b9d4a Mon Sep 17 00:00:00 2001 From: J Dark Date: Fri, 26 Jul 2024 18:24:09 +0200 Subject: [PATCH 10/29] remove angles --- .../derived_quantities/average_volume.py | 74 ++----------------- 1 file changed, 8 insertions(+), 66 deletions(-) diff --git a/festim/exports/derived_quantities/average_volume.py b/festim/exports/derived_quantities/average_volume.py index 4dc70e1ce..7c87e1908 100644 --- a/festim/exports/derived_quantities/average_volume.py +++ b/festim/exports/derived_quantities/average_volume.py @@ -62,25 +62,11 @@ class AverageVolumeCylindrical(AverageVolume): file function (dolfin.function.function.Function): the solution function of the field - azimuth_range (tuple, optional): Range of the azimuthal angle - (theta) needs to be between 0 and 2 pi. Defaults to (0, 2 * np.pi). """ - def __init__(self, field, volume, z, azimuth_range=(0, 2 * np.pi)) -> None: + def __init__(self, field, volume) -> None: super().__init__(field=field, volume=volume) self.r = None - self.z = z - self.azimuth_range = azimuth_range - - @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 def compute(self): @@ -91,15 +77,9 @@ def compute(self): rthetaz = f.SpatialCoordinate(mesh) # get the coordinates from the mesh self.r = rthetaz[0] # only care about r here - values = f.assemble(self.function * self.r * self.z * self.dx(self.volume)) * ( - self.azimuth_range[1] - self.azimuth_range[0] - ) - - volume = f.assemble(1 * self.r * self.z * self.dx(self.volume)) * ( - self.azimuth_range[1] - self.azimuth_range[0] - ) - - avg_vol = values / volume + avg_vol = f.assemble( + self.function * self.r * self.dx(self.volume) + ) / f.assemble(1 * self.r * self.dx(self.volume)) return avg_vol @@ -121,39 +101,11 @@ class AverageVolumeSpherical(AverageVolume): file function (dolfin.function.function.Function): the solution function of the field - azimuth_range (tuple, optional): Range of the azimuthal angle - (theta) needs to be between 0 and pi. Defaults to (0, np.pi) - polar_range (tuple, optional): Range of the polar angle - (phi) needs to be between - pi and pi. Defaults to (-np.pi, np.pi) """ - def __init__( - self, field, volume, azimuth_range=(0, np.pi), polar_range=(-np.pi, np.pi) - ) -> None: + def __init__(self, field, volume) -> None: super().__init__(field=field, volume=volume) self.r = None - self.azimuth_range = azimuth_range - self.polar_range = polar_range - - @property - def polar_range(self): - return self._polar_range - - @polar_range.setter - def polar_range(self, value): - if value[0] < -np.pi or value[1] > np.pi: - raise ValueError("Polar range must be between - pi and pi") - self._polar_range = value - - @property - def azimuth_range(self): - return self._azimuth_range - - @azimuth_range.setter - def azimuth_range(self, value): - if value[0] < 0 or value[1] > np.pi: - raise ValueError("Azimuthal range must be between 0 and pi") - self._azimuth_range = value def compute(self): @@ -164,18 +116,8 @@ def compute(self): rthetaphi = f.SpatialCoordinate(mesh) # get the coordinates from the mesh self.r = rthetaphi[0] # only care about r here - values = ( - f.assemble(self.function * self.r**2 * self.dx(self.volume)) - * (self.polar_range[1] - self.polar_range[0]) - * (-np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0])) - ) - - volume = ( - f.assemble(1 * self.r**2 * self.dx(self.volume)) - * (self.polar_range[1] - self.polar_range[0]) - * (-np.cos(self.azimuth_range[1]) + np.cos(self.azimuth_range[0])) - ) - - avg_vol = values / volume + 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 From bb1a023b3afd9a8f46557ac8a011061901995983 Mon Sep 17 00:00:00 2001 From: JDark Date: Mon, 5 Aug 2024 09:08:11 +0200 Subject: [PATCH 11/29] update doc strings --- festim/exports/derived_quantities/average_surface.py | 12 ++++++++++++ festim/exports/derived_quantities/average_volume.py | 6 ++++++ 2 files changed, 18 insertions(+) diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py index cc67955ec..f4a4b477f 100644 --- a/festim/exports/derived_quantities/average_surface.py +++ b/festim/exports/derived_quantities/average_surface.py @@ -60,11 +60,16 @@ class AverageSurfaceCylindrical(AverageSurface): 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 @@ -82,6 +87,8 @@ def compute(self): ) # 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 + print(type(self.r)) + quit() # dS_z = r dr dtheta , assuming axisymmetry dS_z = theta r dr # dS_r = r dz dtheta , assuming axisymmetry dS_r = theta r dz @@ -104,11 +111,16 @@ class AverageSurfaceSpherical(AverageSurface): 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 diff --git a/festim/exports/derived_quantities/average_volume.py b/festim/exports/derived_quantities/average_volume.py index d8272111f..4b67dee34 100644 --- a/festim/exports/derived_quantities/average_volume.py +++ b/festim/exports/derived_quantities/average_volume.py @@ -20,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 @@ -61,11 +62,16 @@ class AverageVolumeCylindrical(AverageVolume): 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: From d4fe1c2b920baeeeaeb996fc9eda057d94e9ca1d Mon Sep 17 00:00:00 2001 From: JDark Date: Mon, 5 Aug 2024 09:14:43 +0200 Subject: [PATCH 12/29] remove print and quit --- festim/exports/derived_quantities/average_surface.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py index f4a4b477f..cd40b0cff 100644 --- a/festim/exports/derived_quantities/average_surface.py +++ b/festim/exports/derived_quantities/average_surface.py @@ -87,8 +87,6 @@ def compute(self): ) # 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 - print(type(self.r)) - quit() # dS_z = r dr dtheta , assuming axisymmetry dS_z = theta r dr # dS_r = r dz dtheta , assuming axisymmetry dS_r = theta r dz From d904cb3d0ae047f2434f4b2a7351d1af2bfecd08 Mon Sep 17 00:00:00 2001 From: JDark Date: Mon, 5 Aug 2024 14:59:01 +0200 Subject: [PATCH 13/29] added tests --- .../test_average_surface.py | 144 +++++++++++++++++- 1 file changed, 143 insertions(+), 1 deletion(-) diff --git a/test/unit/test_exports/test_derived_quantities/test_average_surface.py b/test/unit/test_exports/test_derived_quantities/test_average_surface.py index 556dd7865..be1e84a06 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_surface.py +++ b/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_cylindrical_flux_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_h_flux = AverageSurfaceCylindrical("solute", 4) + assert my_h_flux.title == "Average solute surface 4" + + +def test_cylindrical_flux_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_heat_flux = AverageSurfaceCylindrical("T", 5) + assert my_heat_flux.title == "Average T surface 5" + + +def test_spherical_flux_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_h_flux = AverageSurfaceSpherical("solute", 6) + assert my_h_flux.title == "Average solute surface 6" + + +def test_spherical_flux_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_heat_flux = AverageSurfaceSpherical("T", 9) + assert my_heat_flux.title == "Average T surface 9" From b96300df5ffa83a6d2fdacdb51787b774f983884 Mon Sep 17 00:00:00 2001 From: JDark Date: Mon, 5 Aug 2024 15:41:51 +0200 Subject: [PATCH 14/29] rename tests --- .../test_derived_quantities/test_average_surface.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/unit/test_exports/test_derived_quantities/test_average_surface.py b/test/unit/test_exports/test_derived_quantities/test_average_surface.py index be1e84a06..ea35a9a86 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_surface.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_surface.py @@ -155,7 +155,7 @@ def test_compute_spherical(r0, radius, c_left, c_right): assert np.isclose(computed_value, expected_value) -def test_cylindrical_flux_title_no_units_solute(): +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""" @@ -163,7 +163,7 @@ def test_cylindrical_flux_title_no_units_solute(): assert my_h_flux.title == "Average solute surface 4" -def test_cylindrical_flux_title_no_units_temperature(): +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""" @@ -171,7 +171,7 @@ def test_cylindrical_flux_title_no_units_temperature(): assert my_heat_flux.title == "Average T surface 5" -def test_spherical_flux_title_no_units_solute(): +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""" @@ -179,7 +179,7 @@ def test_spherical_flux_title_no_units_solute(): assert my_h_flux.title == "Average solute surface 6" -def test_spherical_flux_title_no_units_temperature(): +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""" From 41b90f13db82c6347e30e20e97d0a2e88c214fbf Mon Sep 17 00:00:00 2001 From: JDark Date: Mon, 5 Aug 2024 15:45:09 +0200 Subject: [PATCH 15/29] change export name --- .../test_average_surface.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/unit/test_exports/test_derived_quantities/test_average_surface.py b/test/unit/test_exports/test_derived_quantities/test_average_surface.py index ea35a9a86..c0489e704 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_surface.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_surface.py @@ -159,29 +159,29 @@ 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_h_flux = AverageSurfaceCylindrical("solute", 4) - assert my_h_flux.title == "Average solute surface 4" + 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_heat_flux = AverageSurfaceCylindrical("T", 5) - assert my_heat_flux.title == "Average T surface 5" + 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_h_flux = AverageSurfaceSpherical("solute", 6) - assert my_h_flux.title == "Average solute surface 6" + 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_heat_flux = AverageSurfaceSpherical("T", 9) - assert my_heat_flux.title == "Average T surface 9" + my_export = AverageSurfaceSpherical("T", 9) + assert my_export.title == "Average T surface 9" From 5944b682f8ad9dc5237bb944b7d90c89b9ccc940 Mon Sep 17 00:00:00 2001 From: JDark Date: Mon, 5 Aug 2024 15:45:19 +0200 Subject: [PATCH 16/29] test average volume --- .../test_average_volume.py | 125 +++++++++++++++++- 1 file changed, 124 insertions(+), 1 deletion(-) diff --git a/test/unit/test_exports/test_derived_quantities/test_average_volume.py b/test/unit/test_exports/test_derived_quantities/test_average_volume.py index 172f9b5e6..5bb300448 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_volume.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_volume.py @@ -1,6 +1,9 @@ -from festim import AverageVolume +from festim import AverageVolume, AverageVolumeCylindrical, AverageVolumeSpherical, x, y import fenics as f import pytest +import numpy as np +from sympy.printing import ccode +import math @pytest.mark.parametrize("field,volume", [("solute", 1), ("T", 2)]) @@ -43,3 +46,123 @@ def test_h_average(self): ) computed = self.my_average.compute() assert computed == expected + + +@pytest.mark.parametrize("radius", [2, 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 AverageVolumeCylindrical 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), 50, 50) + + volume_markers = f.MeshFunction("size_t", mesh_fenics, mesh_fenics.topology().dim()) + volume_markers.set_all(1) + dx = f.Measure("dx", domain=mesh_fenics, subdomain_data=volume_markers) + + my_export = AverageVolumeCylindrical("solute", 1) + V = f.FunctionSpace(mesh_fenics, "P", 1) + c_fun = lambda z: c_bottom + (c_top - c_bottom) / (height) * z + expr = f.Expression( + ccode(c_fun(y)), + degree=1, + ) + my_export.function = f.interpolate(expr, V) + my_export.dx = dx + + expected_value = (c_bottom + c_top) / 2 + + 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 AverageVolumeSpherical 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) + volume_markers = f.MeshFunction("size_t", mesh_fenics, mesh_fenics.topology().dim()) + volume_markers.set_all(1) + dx = f.Measure("dx", domain=mesh_fenics, subdomain_data=volume_markers) + + my_export = AverageVolumeSpherical("solute", 1) + 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.dx = dx + + expected_value = c_left + (3 * (c_right - c_left)) / (4 * (r1**3 - r0**3)) * ( + r1 + r0 + ) * (r1**2 + r0**2) + + computed_value = my_export.compute() + + assert np.isclose(computed_value, expected_value) + + +def test_average_volume_cylindrical_title_no_units_solute(): + """A simple test to check that the title is set correctly in + festim.AverageVolumeCylindrical with a solute field without units""" + + my_export = AverageVolumeCylindrical("solute", 4) + assert my_export.title == "Average solute volume 4" + + +def test_average_volume_cylindrical_title_no_units_temperature(): + """A simple test to check that the title is set correctly in + festim.AverageVolumeCylindrical with a T field without units""" + + my_export = AverageVolumeCylindrical("T", 5) + assert my_export.title == "Average T volume 5" + + +def test_average_volume_spherical_title_no_units_solute(): + """A simple test to check that the title is set correctly in + festim.AverageVolumeSpherical with a solute field without units""" + + my_export = AverageVolumeSpherical("solute", 6) + assert my_export.title == "Average solute volume 6" + + +def test_average_volume_spherical_title_no_units_temperature(): + """A simple test to check that the title is set correctly in + festim.AverageVolumeSpherical with a T field without units""" + + my_export = AverageVolumeSpherical("T", 9) + assert my_export.title == "Average T volume 9" From e747d52fbdd0e04321f779c25d251215bf20acb6 Mon Sep 17 00:00:00 2001 From: JDark Date: Mon, 5 Aug 2024 15:52:49 +0200 Subject: [PATCH 17/29] dont need haevy mesh --- .../test_exports/test_derived_quantities/test_average_volume.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit/test_exports/test_derived_quantities/test_average_volume.py b/test/unit/test_exports/test_derived_quantities/test_average_volume.py index 5bb300448..2bc67dc7b 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_volume.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_volume.py @@ -69,7 +69,7 @@ def test_compute_cylindrical(r0, radius, height, c_top, c_bottom): z0 = 0 z1 = z0 + height - mesh_fenics = f.RectangleMesh(f.Point(r0, z0), f.Point(r1, z1), 50, 50) + mesh_fenics = f.RectangleMesh(f.Point(r0, z0), f.Point(r1, z1), 10, 10) volume_markers = f.MeshFunction("size_t", mesh_fenics, mesh_fenics.topology().dim()) volume_markers.set_all(1) From 636a322c465b33b331df451fb0763df2ee8b90fc Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 10 Jan 2025 11:10:23 -0500 Subject: [PATCH 18/29] show units now default --- .../test_derived_quantities/test_average_surface.py | 8 ++++---- .../test_derived_quantities/test_average_volume.py | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/test/unit/test_exports/test_derived_quantities/test_average_surface.py b/test/unit/test_exports/test_derived_quantities/test_average_surface.py index 12d1d55be..f26363eab 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_surface.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_surface.py @@ -163,7 +163,7 @@ def test_average_surface_cylindrical_title_no_units_solute(): festim.AverageSurfaceCylindrical with a solute field without units""" my_export = AverageSurfaceCylindrical("solute", 4) - assert my_export.title == "Average solute surface 4" + assert my_export.title == "Average solute surface 4 (H m-3)" def test_average_surface_cylindrical_title_no_units_temperature(): @@ -171,7 +171,7 @@ def test_average_surface_cylindrical_title_no_units_temperature(): festim.AverageSurfaceCylindrical with a T field without units""" my_export = AverageSurfaceCylindrical("T", 5) - assert my_export.title == "Average T surface 5" + assert my_export.title == "Average T surface 5 (K)" def test_average_surface_spherical_title_no_units_solute(): @@ -179,7 +179,7 @@ def test_average_surface_spherical_title_no_units_solute(): festim.AverageSurfaceSpherical with a solute field without units""" my_export = AverageSurfaceSpherical("solute", 6) - assert my_export.title == "Average solute surface 6" + assert my_export.title == "Average solute surface 6 (H m-3)" def test_average_surface_spherical_title_no_units_temperature(): @@ -187,4 +187,4 @@ def test_average_surface_spherical_title_no_units_temperature(): festim.AverageSurfaceSpherical with a T field without units""" my_export = AverageSurfaceSpherical("T", 9) - assert my_export.title == "Average T surface 9" + assert my_export.title == "Average T surface 9 (K)" diff --git a/test/unit/test_exports/test_derived_quantities/test_average_volume.py b/test/unit/test_exports/test_derived_quantities/test_average_volume.py index f4c59f5eb..2e3ed0b54 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_volume.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_volume.py @@ -144,7 +144,7 @@ def test_average_volume_cylindrical_title_no_units_solute(): festim.AverageVolumeCylindrical with a solute field without units""" my_export = AverageVolumeCylindrical("solute", 4) - assert my_export.title == "Average solute volume 4" + assert my_export.title == "Average solute volume 4 (H m-3)" def test_average_volume_cylindrical_title_no_units_temperature(): @@ -152,7 +152,7 @@ def test_average_volume_cylindrical_title_no_units_temperature(): festim.AverageVolumeCylindrical with a T field without units""" my_export = AverageVolumeCylindrical("T", 5) - assert my_export.title == "Average T volume 5" + assert my_export.title == "Average T volume 5 (K)" def test_average_volume_spherical_title_no_units_solute(): @@ -160,7 +160,7 @@ def test_average_volume_spherical_title_no_units_solute(): festim.AverageVolumeSpherical with a solute field without units""" my_export = AverageVolumeSpherical("solute", 6) - assert my_export.title == "Average solute volume 6" + assert my_export.title == "Average solute volume 6 (H m-3)" def test_average_volume_spherical_title_no_units_temperature(): @@ -168,4 +168,4 @@ def test_average_volume_spherical_title_no_units_temperature(): festim.AverageVolumeSpherical with a T field without units""" my_export = AverageVolumeSpherical("T", 9) - assert my_export.title == "Average T volume 9" + assert my_export.title == "Average T volume 9 (K)" From 87c43121ff274bcf18367351cacacfd7aa82a0fe Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 10 Jan 2025 13:13:05 -0500 Subject: [PATCH 19/29] simplify avg sph surf class --- .../derived_quantities/average_surface.py | 46 ++----------------- 1 file changed, 3 insertions(+), 43 deletions(-) diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py index 965585a0b..4148d353e 100644 --- a/festim/exports/derived_quantities/average_surface.py +++ b/festim/exports/derived_quantities/average_surface.py @@ -105,48 +105,8 @@ def compute(self): 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 + Computes the average on a spherical "surface" in 1D. + Behaves identically to `AverageSurface`. """ - 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 + pass From d4b6979ae0a8eb6983d15da8d0fcbc2f3a845fd5 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 10 Jan 2025 13:13:11 -0500 Subject: [PATCH 20/29] fix tests --- .../test_average_surface.py | 28 +++++++---------- .../test_average_volume.py | 30 +++++-------------- 2 files changed, 18 insertions(+), 40 deletions(-) diff --git a/test/unit/test_exports/test_derived_quantities/test_average_surface.py b/test/unit/test_exports/test_derived_quantities/test_average_surface.py index f26363eab..d0a8cf416 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_surface.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_surface.py @@ -59,9 +59,7 @@ def test_h_average(self): @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): +def test_compute_cylindrical(r0, radius, height): """ Test that AverageSurfaceCylindrical computes the value correctly on a hollow cylinder @@ -69,19 +67,16 @@ def test_compute_cylindrical(r0, radius, height, c_top, c_bottom): 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 + z0, z1 = 0, 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 ) @@ -93,15 +88,16 @@ def test_compute_cylindrical(r0, radius, height, c_top, c_bottom): 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 + c_fun = lambda r: 2 * r + expr = f.Expression( - ccode(c_fun(y)), + ccode(c_fun(x)), degree=1, ) my_export.function = f.interpolate(expr, V) my_export.ds = ds - expected_value = c_bottom + (c_top - c_bottom) / (z1 - z0) * height + expected_value = 4 / 3 * (r1**3 - r0**3) / (r1**2 - r0**2) computed_value = my_export.compute() @@ -110,9 +106,7 @@ def test_compute_cylindrical(r0, radius, height, c_top, c_bottom): @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): +def test_compute_spherical(r0, radius): """ Test that AverageSurfaceSpherical computes the average value correctly on a hollow sphere @@ -120,8 +114,6 @@ def test_compute_spherical(r0, radius, c_left, c_right): 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 @@ -142,7 +134,7 @@ def test_compute_spherical(r0, radius, c_left, c_right): 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 + c_fun = lambda r: 4 * r expr = f.Expression( ccode(c_fun(x)), degree=1, @@ -151,7 +143,7 @@ def test_compute_spherical(r0, radius, c_left, c_right): my_export.ds = ds - expected_value = c_left + ((c_right - c_left) / (r1 - r0)) * r1 + expected_value = 4 * r1 computed_value = my_export.compute() diff --git a/test/unit/test_exports/test_derived_quantities/test_average_volume.py b/test/unit/test_exports/test_derived_quantities/test_average_volume.py index 2e3ed0b54..84f437d04 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_volume.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_volume.py @@ -54,9 +54,7 @@ def test_h_average(self): @pytest.mark.parametrize("radius", [2, 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): +def test_compute_cylindrical(r0, radius, height): """ Test that AverageVolumeCylindrical computes the value correctly on a hollow cylinder @@ -64,13 +62,10 @@ def test_compute_cylindrical(r0, radius, height, c_top, c_bottom): 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 + z0, z1 = 0, height mesh_fenics = f.RectangleMesh(f.Point(r0, z0), f.Point(r1, z1), 10, 10) @@ -80,15 +75,12 @@ def test_compute_cylindrical(r0, radius, height, c_top, c_bottom): my_export = AverageVolumeCylindrical("solute", 1) V = f.FunctionSpace(mesh_fenics, "P", 1) - c_fun = lambda z: c_bottom + (c_top - c_bottom) / (height) * z - expr = f.Expression( - ccode(c_fun(y)), - degree=1, - ) + c_fun = lambda r: 3 * r + expr = f.Expression(ccode(c_fun(x)), degree=1) my_export.function = f.interpolate(expr, V) my_export.dx = dx - expected_value = (c_bottom + c_top) / 2 + expected_value = 2 * (r1**3 - r0**3) / (r1**2 - r0**2) computed_value = my_export.compute() @@ -97,9 +89,7 @@ def test_compute_cylindrical(r0, radius, height, c_top, c_bottom): @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): +def test_compute_spherical(r0, radius): """ Test that AverageVolumeSpherical computes the average value correctly on a hollow sphere @@ -107,8 +97,6 @@ def test_compute_spherical(r0, radius, c_left, c_right): 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 @@ -121,7 +109,7 @@ def test_compute_spherical(r0, radius, c_left, c_right): my_export = AverageVolumeSpherical("solute", 1) V = f.FunctionSpace(mesh_fenics, "P", 1) - c_fun = lambda r: c_left + (c_right - c_left) / (r1 - r0) * r + c_fun = lambda r: 4 * r expr = f.Expression( ccode(c_fun(x)), degree=1, @@ -130,9 +118,7 @@ def test_compute_spherical(r0, radius, c_left, c_right): my_export.dx = dx - expected_value = c_left + (3 * (c_right - c_left)) / (4 * (r1**3 - r0**3)) * ( - r1 + r0 - ) * (r1**2 + r0**2) + expected_value = 3 * (r1**4 - r0**4) / (r1**3 - r0**3) computed_value = my_export.compute() From 152a9c5619213fd331629d03c60ea7f260dc9520 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 10 Jan 2025 13:21:34 -0500 Subject: [PATCH 21/29] remove unused imports --- .../test_derived_quantities/test_average_surface.py | 8 +------- .../test_derived_quantities/test_average_volume.py | 3 +-- 2 files changed, 2 insertions(+), 9 deletions(-) diff --git a/test/unit/test_exports/test_derived_quantities/test_average_surface.py b/test/unit/test_exports/test_derived_quantities/test_average_surface.py index d0a8cf416..df660951c 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_surface.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_surface.py @@ -1,10 +1,4 @@ -from festim import ( - AverageSurface, - AverageSurfaceCylindrical, - AverageSurfaceSpherical, - x, - y, -) +from festim import AverageSurface, AverageSurfaceCylindrical, AverageSurfaceSpherical, x import fenics as f import pytest import numpy as np diff --git a/test/unit/test_exports/test_derived_quantities/test_average_volume.py b/test/unit/test_exports/test_derived_quantities/test_average_volume.py index 84f437d04..a638d2f70 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_volume.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_volume.py @@ -1,9 +1,8 @@ -from festim import AverageVolume, AverageVolumeCylindrical, AverageVolumeSpherical, x, y +from festim import AverageVolume, AverageVolumeCylindrical, AverageVolumeSpherical, x import fenics as f import pytest import numpy as np from sympy.printing import ccode -import math @pytest.mark.parametrize("field,volume", [("solute", 1), ("T", 2)]) From 1e5f411e23b2ff1805713f2f2fcf0c1840217aa0 Mon Sep 17 00:00:00 2001 From: James Dark <65899899+jhdark@users.noreply.github.com> Date: Fri, 10 Jan 2025 14:41:06 -0500 Subject: [PATCH 22/29] Update festim/exports/derived_quantities/average_surface.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Rémi Delaporte-Mathurin <40028739+RemDelaporteMathurin@users.noreply.github.com> --- festim/exports/derived_quantities/average_surface.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py index 4148d353e..edc646e6e 100644 --- a/festim/exports/derived_quantities/average_surface.py +++ b/festim/exports/derived_quantities/average_surface.py @@ -103,10 +103,4 @@ def compute(self): return avg_surf -class AverageSurfaceSpherical(AverageSurface): - """ - Computes the average on a spherical "surface" in 1D. - Behaves identically to `AverageSurface`. - """ - - pass +AverageSurfaceSpherical = AverageSurface From 296d744f875c9ae38c73450b94daf43a25faeea7 Mon Sep 17 00:00:00 2001 From: jhdark Date: Fri, 10 Jan 2025 17:57:47 -0500 Subject: [PATCH 23/29] test avg surface on a rhombus --- .../test_average_surface.py | 65 ++++++++++++++++++- 1 file changed, 64 insertions(+), 1 deletion(-) diff --git a/test/unit/test_exports/test_derived_quantities/test_average_surface.py b/test/unit/test_exports/test_derived_quantities/test_average_surface.py index df660951c..f9142cd65 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_surface.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_surface.py @@ -1,4 +1,10 @@ -from festim import AverageSurface, AverageSurfaceCylindrical, AverageSurfaceSpherical, x +from festim import ( + AverageSurface, + AverageSurfaceCylindrical, + AverageSurfaceSpherical, + x, + y, +) import fenics as f import pytest import numpy as np @@ -174,3 +180,60 @@ def test_average_surface_spherical_title_no_units_temperature(): my_export = AverageSurfaceSpherical("T", 9) assert my_export.title == "Average T surface 9 (K)" + + +@pytest.mark.parametrize("j", [4, 5, 3]) +def test_average_surface_rhombus(j): + + # creating a mesh with FEniCS + # Define the number of divisions in the x and y directions + nx, ny = 10, 10 + mesh = f.RectangleMesh(f.Point(0, 0), f.Point(j, j), nx, ny) + + # Access mesh coordinates + coordinates = mesh.coordinates() + + # Rotation matrix for 45 degrees + theta = np.pi / 4 # 45 degrees in radians + rotation_matrix = np.array( + [[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]] + ) + + # Apply rotation to each coordinate + for i in range(len(coordinates)): + x_cord, y_cord = coordinates[i] + # Rotate the point (x, y) using the rotation matrix + new_coords = np.dot(rotation_matrix, np.array([x_cord, y_cord])) + coordinates[i] = new_coords + + # Create a new mesh with the rotated coordinates + mesh_rotated = f.Mesh(mesh) + + surface_markers = f.MeshFunction( + "size_t", mesh_rotated, mesh_rotated.topology().dim() - 1 + ) + surface_markers.set_all(0) + + # find all facets along y = x and mark them + surf_id = 2 + for facet in f.facets(mesh): + midpoint = facet.midpoint() + if np.isclose(midpoint[0], midpoint[1], atol=1e-10): + surface_markers[facet.index()] = surf_id + + ds = f.Measure("ds", domain=mesh_rotated, subdomain_data=surface_markers) + my_export = AverageSurface("solute", surf_id) + V = f.FunctionSpace(mesh_rotated, "P", 1) + + c_fun = lambda x, y: 2 * x + y + expr = f.Expression( + ccode(c_fun(x, y)), + degree=1, + ) + my_export.function = f.interpolate(expr, V) + my_export.ds = ds + + expected_value = (3 * j) / (2 * (2**0.5)) + computed_value = my_export.compute() + + assert np.isclose(computed_value, expected_value) From 58bd814cdb653030c74d329547d9237a7c35432b Mon Sep 17 00:00:00 2001 From: jhdark Date: Tue, 14 Jan 2025 09:54:16 -0500 Subject: [PATCH 24/29] add allowed meshes property --- festim/exports/derived_quantities/average_surface.py | 5 ++++- festim/exports/derived_quantities/average_volume.py | 9 ++++++++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py index edc646e6e..7579db678 100644 --- a/festim/exports/derived_quantities/average_surface.py +++ b/festim/exports/derived_quantities/average_surface.py @@ -1,6 +1,5 @@ from festim import SurfaceQuantity import fenics as f -import numpy as np class AverageSurface(SurfaceQuantity): @@ -83,6 +82,10 @@ 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: diff --git a/festim/exports/derived_quantities/average_volume.py b/festim/exports/derived_quantities/average_volume.py index d2b9f14fd..6428c0e5a 100644 --- a/festim/exports/derived_quantities/average_volume.py +++ b/festim/exports/derived_quantities/average_volume.py @@ -1,6 +1,5 @@ from festim import VolumeQuantity import fenics as f -import numpy as np class AverageVolume(VolumeQuantity): @@ -82,6 +81,10 @@ 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: @@ -121,6 +124,10 @@ 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: From e5f0a6a55242f188732cb498d02ff26cbaad918f Mon Sep 17 00:00:00 2001 From: jhdark Date: Tue, 14 Jan 2025 10:09:35 -0500 Subject: [PATCH 25/29] accept spherical meshes for evg surface --- festim/exports/derived_quantities/average_surface.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/festim/exports/derived_quantities/average_surface.py b/festim/exports/derived_quantities/average_surface.py index 7579db678..c56e0d793 100644 --- a/festim/exports/derived_quantities/average_surface.py +++ b/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): From ce2f1485d04ed605fc4281f501a7d6e40616ef59 Mon Sep 17 00:00:00 2001 From: jhdark Date: Tue, 14 Jan 2025 10:47:36 -0500 Subject: [PATCH 26/29] remove uneccesary test --- .../test_average_surface.py | 57 ------------------- 1 file changed, 57 deletions(-) diff --git a/test/unit/test_exports/test_derived_quantities/test_average_surface.py b/test/unit/test_exports/test_derived_quantities/test_average_surface.py index f9142cd65..d0a8cf416 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_surface.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_surface.py @@ -180,60 +180,3 @@ def test_average_surface_spherical_title_no_units_temperature(): my_export = AverageSurfaceSpherical("T", 9) assert my_export.title == "Average T surface 9 (K)" - - -@pytest.mark.parametrize("j", [4, 5, 3]) -def test_average_surface_rhombus(j): - - # creating a mesh with FEniCS - # Define the number of divisions in the x and y directions - nx, ny = 10, 10 - mesh = f.RectangleMesh(f.Point(0, 0), f.Point(j, j), nx, ny) - - # Access mesh coordinates - coordinates = mesh.coordinates() - - # Rotation matrix for 45 degrees - theta = np.pi / 4 # 45 degrees in radians - rotation_matrix = np.array( - [[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]] - ) - - # Apply rotation to each coordinate - for i in range(len(coordinates)): - x_cord, y_cord = coordinates[i] - # Rotate the point (x, y) using the rotation matrix - new_coords = np.dot(rotation_matrix, np.array([x_cord, y_cord])) - coordinates[i] = new_coords - - # Create a new mesh with the rotated coordinates - mesh_rotated = f.Mesh(mesh) - - surface_markers = f.MeshFunction( - "size_t", mesh_rotated, mesh_rotated.topology().dim() - 1 - ) - surface_markers.set_all(0) - - # find all facets along y = x and mark them - surf_id = 2 - for facet in f.facets(mesh): - midpoint = facet.midpoint() - if np.isclose(midpoint[0], midpoint[1], atol=1e-10): - surface_markers[facet.index()] = surf_id - - ds = f.Measure("ds", domain=mesh_rotated, subdomain_data=surface_markers) - my_export = AverageSurface("solute", surf_id) - V = f.FunctionSpace(mesh_rotated, "P", 1) - - c_fun = lambda x, y: 2 * x + y - expr = f.Expression( - ccode(c_fun(x, y)), - degree=1, - ) - my_export.function = f.interpolate(expr, V) - my_export.ds = ds - - expected_value = (3 * j) / (2 * (2**0.5)) - computed_value = my_export.compute() - - assert np.isclose(computed_value, expected_value) From 47f853e076193ce864cae468aca864aa8b2a7c88 Mon Sep 17 00:00:00 2001 From: jhdark Date: Tue, 14 Jan 2025 10:47:45 -0500 Subject: [PATCH 27/29] test allowed meshes --- test/simulation/test_initialise.py | 57 +++++++++++++++++++++++++++++- 1 file changed, 56 insertions(+), 1 deletion(-) diff --git a/test/simulation/test_initialise.py b/test/simulation/test_initialise.py index fdcb3a60d..3418a14d9 100644 --- a/test/simulation/test_initialise.py +++ b/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", [ From 301cba09ac941256ce952e8150150a7810458e77 Mon Sep 17 00:00:00 2001 From: jhdark Date: Tue, 14 Jan 2025 11:07:34 -0500 Subject: [PATCH 28/29] additional tests for avg volume --- test/simulation/test_initialise.py | 44 ++++++++++++++++++++++++++---- 1 file changed, 39 insertions(+), 5 deletions(-) diff --git a/test/simulation/test_initialise.py b/test/simulation/test_initialise.py index 3418a14d9..afa6ac1e3 100644 --- a/test/simulation/test_initialise.py +++ b/test/simulation/test_initialise.py @@ -150,8 +150,15 @@ def test_cartesian_and_surface_flux_warning(quantity, sys): my_model.initialise() +@pytest.mark.parametrize( + "quantity", + [ + F.AverageSurfaceCylindrical(field="solute", surface=1), + F.AverageVolumeCylindrical(field="solute", volume=1), + ], +) @pytest.mark.parametrize("sys", ["cartesian", "spherical"]) -def test_cylindrical_quantities_warning(sys): +def test_cylindrical_quantities_warning(quantity, sys): """ Creates a Simulation object and checks that, if either a cartesian or spherical meshes are given with a AverageSurfaceCylindrical, a warning is raised. @@ -169,9 +176,7 @@ def test_cylindrical_quantities_warning(sys): absolute_tolerance=1e-10, relative_tolerance=1e-10, final_time=4 ) - derived_quantities = F.DerivedQuantities( - [F.AverageSurfaceCylindrical(field="solute", surface=1)] - ) + derived_quantities = F.DerivedQuantities([quantity]) my_model.exports = [derived_quantities] # test @@ -179,7 +184,7 @@ def test_cylindrical_quantities_warning(sys): my_model.initialise() -def test_spherical_quantities_warning(): +def test_avg_surf_spherical_quantities_warning(): """ Creates a Simulation object and checks that, if a cylindrical mesh is given with a AverageSurfaceSpherical, a warning is raised. @@ -206,6 +211,35 @@ def test_spherical_quantities_warning(): my_model.initialise() +@pytest.mark.parametrize("sys", ["cartesian", "cylindrical"]) +def test_avg_volume_spherical_quantities_warning(sys): + """ + Creates a Simulation object and checks that, if either a cylindrical or + cartesian meshes are given with a AverageVolumeSpherical, 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.AverageVolumeSpherical(field="solute", volume=1)] + ) + my_model.exports = [derived_quantities] + + # test + with pytest.warns(UserWarning, match=f"may not work as intended for {sys} meshes"): + my_model.initialise() + + @pytest.mark.parametrize( "value", [ From 10fc3b00d217720472ffe0b5bfa29fc866b8a6de Mon Sep 17 00:00:00 2001 From: jhdark Date: Tue, 14 Jan 2025 11:25:47 -0500 Subject: [PATCH 29/29] simplify tests --- test/simulation/test_initialise.py | 90 ------------------- .../test_average_surface.py | 18 ++++ .../test_average_volume.py | 18 ++++ 3 files changed, 36 insertions(+), 90 deletions(-) diff --git a/test/simulation/test_initialise.py b/test/simulation/test_initialise.py index afa6ac1e3..c5c09aaa8 100644 --- a/test/simulation/test_initialise.py +++ b/test/simulation/test_initialise.py @@ -150,96 +150,6 @@ def test_cartesian_and_surface_flux_warning(quantity, sys): my_model.initialise() -@pytest.mark.parametrize( - "quantity", - [ - F.AverageSurfaceCylindrical(field="solute", surface=1), - F.AverageVolumeCylindrical(field="solute", volume=1), - ], -) -@pytest.mark.parametrize("sys", ["cartesian", "spherical"]) -def test_cylindrical_quantities_warning(quantity, 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([quantity]) - 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_avg_surf_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("sys", ["cartesian", "cylindrical"]) -def test_avg_volume_spherical_quantities_warning(sys): - """ - Creates a Simulation object and checks that, if either a cylindrical or - cartesian meshes are given with a AverageVolumeSpherical, 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.AverageVolumeSpherical(field="solute", volume=1)] - ) - my_model.exports = [derived_quantities] - - # test - with pytest.warns(UserWarning, match=f"may not work as intended for {sys} meshes"): - my_model.initialise() - - @pytest.mark.parametrize( "value", [ diff --git a/test/unit/test_exports/test_derived_quantities/test_average_surface.py b/test/unit/test_exports/test_derived_quantities/test_average_surface.py index d0a8cf416..17f47ac8d 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_surface.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_surface.py @@ -180,3 +180,21 @@ def test_average_surface_spherical_title_no_units_temperature(): my_export = AverageSurfaceSpherical("T", 9) assert my_export.title == "Average T surface 9 (K)" + + +def test_avg_surf_cylindrical_allow_meshes(): + """A simple test to check cylindrical meshes are the only + meshes allowed when using AverageSurfaceCylindrical""" + + my_export = AverageSurfaceCylindrical("solute", 2) + + assert my_export.allowed_meshes == ["cylindrical"] + + +def test_avg_surf_spherical_allow_meshes(): + """A simple test to check spherical meshes are one of the + meshes allowed when using AverageSurfaceSpherical""" + + my_export = AverageSurfaceSpherical("T", 6) + + assert "spherical" in my_export.allowed_meshes diff --git a/test/unit/test_exports/test_derived_quantities/test_average_volume.py b/test/unit/test_exports/test_derived_quantities/test_average_volume.py index a638d2f70..fc4f83059 100644 --- a/test/unit/test_exports/test_derived_quantities/test_average_volume.py +++ b/test/unit/test_exports/test_derived_quantities/test_average_volume.py @@ -154,3 +154,21 @@ def test_average_volume_spherical_title_no_units_temperature(): my_export = AverageVolumeSpherical("T", 9) assert my_export.title == "Average T volume 9 (K)" + + +def test_avg_vol_cylindrical_allow_meshes(): + """A simple test to check cylindrical meshes are the only + meshes allowed when using AverageVolumeCylindrical""" + + my_export = AverageVolumeCylindrical("solute", 3) + + assert my_export.allowed_meshes == ["cylindrical"] + + +def test_avg_vol_spherical_allow_meshes(): + """A simple test to check cylindrical meshes are the only + meshes allowed when using AverageVolumeSpherical""" + + my_export = AverageVolumeSpherical("solute", 5) + + assert my_export.allowed_meshes == ["spherical"]