From 82dbd1ffb77bec8b5389254e188e49d7892f48af Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 16:15:07 -0400 Subject: [PATCH 01/13] make sure post processing functions have a name --- src/festim/hydrogen_transport_problem.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 14d9c35b4..ee37837d2 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -646,6 +646,7 @@ def assign_functions_to_species(self): spe.sub_function_space = self.function_space.sub(idx) spe.sub_function = self.u.sub(idx) # TODO add this to discontinuous class spe.post_processing_solution = self.u.sub(idx).collapse() + spe.post_processing_solution.name = spe.name spe.collapsed_function_space, spe.map_sub_to_main_solution = ( self.function_space.sub(idx).collapse() ) From 14dcf9580d127503c90bb2ed2348c4829b403334 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 17:22:02 -0400 Subject: [PATCH 02/13] added test --- test/test_vtx.py | 81 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/test/test_vtx.py b/test/test_vtx.py index edbb2e531..335a9bda9 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -179,3 +179,84 @@ def test_filename_temp_raises_error_when_wrong_type(): """Test that the filename attribute for VTXTemperature export raises an error if the extension is not .bp""" with pytest.raises(TypeError): F.VTXTemperatureExport(1) + + +@pytest.mark.parametrize( + "expression", + [ + lambda x: x[0] + x[1] * 2, + lambda T: T + 1, + lambda c_A, c_B: c_A + c_B, + lambda c_A, c_B, x: c_A * c_B + x[0], + lambda c_A, T, x: c_A * T + x[0], + ], +) +def test_custom_field(tmp_path, expression): + """ + Test custom field export functionality. + This test checks that a custom field can be created with various types of + expressions. + """ + + my_model = F.HydrogenTransportProblem() + + mat = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + + vol = F.VolumeSubdomain(id=1, material=mat) + + top = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1)) + bottom = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[1], 0)) + left = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0)) + right = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[0], 1)) + + my_model.subdomains = [vol, top, bottom, left, right] + + dolfinx_mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) + my_model.mesh = F.Mesh(dolfinx_mesh) + + A = F.Species("A") + B = F.Species("B") + C = F.Species("C") + D = F.Species("D") + + my_model.species = [A, B, C, D] + + my_model.boundary_conditions = ( + [ + F.FixedConcentrationBC(species=A, subdomain=top, value=1), + F.FixedConcentrationBC(species=B, subdomain=left, value=1), + ] + + [ + F.FixedConcentrationBC(species=C, subdomain=surf, value=0) + for surf in [top, bottom, left, right] + ] + + [ + F.FixedConcentrationBC(species=D, subdomain=surf, value=0) + for surf in [top, bottom, left, right] + ] + ) + + my_model.reactions = [ + F.Reaction( + reactant=[A, B], product=[C], k_0=1, E_k=0, p_0=0, E_p=0, volume=vol + ), + F.Reaction(reactant=[C], product=[D], k_0=0.1, E_k=0, p_0=0, E_p=0, volume=vol), + ] + + my_model.temperature = 300 + + my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) + + custom_field = F.CustomField( + filename=tmp_path / "custom_field.bp", + expression=expression, + species_dependent_value={"c_A": A, "c_B": B}, + ) + + my_model.exports = [ + custom_field, + ] + + my_model.initialise() + + my_model.run() From 9f672f0925eb539a9ba783cc53f3ee4c9004e457 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 17:22:14 -0400 Subject: [PATCH 03/13] initial implementation --- src/festim/__init__.py | 7 ++- src/festim/exports/__init__.py | 3 +- src/festim/exports/vtx.py | 55 ++++++++++++++++++++++++ src/festim/helpers.py | 2 + src/festim/hydrogen_transport_problem.py | 21 +++++++++ 5 files changed, 86 insertions(+), 2 deletions(-) diff --git a/src/festim/__init__.py b/src/festim/__init__.py index 823d519d9..11db36aec 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -35,7 +35,12 @@ from .exports.total_surface import TotalSurface from .exports.total_volume import TotalVolume from .exports.volume_quantity import VolumeQuantity -from .exports.vtx import ExportBaseClass, VTXSpeciesExport, VTXTemperatureExport +from .exports.vtx import ( + ExportBaseClass, + VTXSpeciesExport, + VTXTemperatureExport, + CustomField, +) from .exports.xdmf import XDMFExport from .heat_transfer_problem import HeatTransferProblem from .helpers import ( diff --git a/src/festim/exports/__init__.py b/src/festim/exports/__init__.py index 31a287244..57c0e822f 100644 --- a/src/festim/exports/__init__.py +++ b/src/festim/exports/__init__.py @@ -10,12 +10,13 @@ from .total_surface import TotalSurface from .total_volume import TotalVolume from .volume_quantity import VolumeQuantity -from .vtx import ExportBaseClass, VTXSpeciesExport, VTXTemperatureExport +from .vtx import ExportBaseClass, VTXSpeciesExport, VTXTemperatureExport, CustomField from .xdmf import XDMFExport __all__ = [ "AverageSurface", "AverageVolume", + "CustomField", "ExportBaseClass", "MaximumSurface", "MaximumVolume", diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 7159360fd..cf2005d66 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -2,9 +2,11 @@ from pathlib import Path from dolfinx import fem, io +import ufl from festim.species import Species from festim.subdomain.volume_subdomain import VolumeSubdomain +import inspect class ExportBaseClass: @@ -172,3 +174,56 @@ def get_functions(self) -> list[fem.Function]: field.subdomain_to_post_processing_solution[self._subdomain] ) return outfiles + + +class CustomField(ExportBaseClass): + function: fem.Function + writer: io.VTXWriter + dolfinx_expression: fem.Expression + + def __init__( + self, + filename, + expression, + species_dependent_value=None, + times=None, + subdomain: VolumeSubdomain = None, + checkpoint=False, + ): + super().__init__( + filename=filename, + times=times, + ext=".bp", + ) + self.expression = expression + self.species_dependent_value = species_dependent_value or {} + self.checkpoint = checkpoint + self.subdomain = subdomain + + def set_dolfinx_expression( + self, + temperature: fem.Constant | fem.Function, + species: list, + time: fem.Constant, + ): + # check + arguments = inspect.signature(self.expression).parameters + kwargs = {} + if "t" in arguments: + kwargs["t"] = time + if "x" in arguments: + x = ufl.SpatialCoordinate(time._ufl_domain) + kwargs["x"] = x + if "T" in arguments: + kwargs["T"] = temperature + # check if there are other arguments and if they are in species_dependent_value + for arg in arguments: + if arg in self.species_dependent_value: + spe = self.species_dependent_value[arg] + kwargs[arg] = spe.post_processing_solution + # TODO handle discontinuous case and use self.subdomain + + self.dolfinx_expression = fem.Expression( + self.expression(**kwargs), + self.function.function_space.element.interpolation_points, + ) diff --git a/src/festim/helpers.py b/src/festim/helpers.py index 0a1ba5290..236d31159 100644 --- a/src/festim/helpers.py +++ b/src/festim/helpers.py @@ -34,6 +34,7 @@ def as_fenics_constant( ) +# TODO change this to accept species dependent values def as_mapped_function( value: Callable, function_space: fem.FunctionSpace | None = None, @@ -68,6 +69,7 @@ def as_mapped_function( return value(**kwargs) +# TODO change this to accept species dependent values def as_fenics_interp_expr_and_function( value: Callable, function_space: dolfinx.fem.function.FunctionSpace, diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index ee37837d2..f29b1f38f 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -449,6 +449,22 @@ def initialise_exports(self): else: adios4dolfinx.write_mesh(export.filename, mesh=self.mesh.mesh) + elif isinstance(export, exports.CustomField): + export.function = fem.Function(self.V_CG_1) + export.set_dolfinx_expression( + temperature=self.temperature_fenics, + species=self.species, + time=self.t, + ) + + export.writer = dolfinx.io.VTXWriter( + comm=export.function.function_space.mesh.comm, + filename=export.filename, + output=export.function, + engine="BP5", + ) + continue + elif isinstance(export, exports.SurfaceQuantity | exports.VolumeQuantity): # raise not implemented error if the derived quantity don't match the # type of mesh eg. SurfaceFlux is used with cylindrical mesh @@ -629,6 +645,7 @@ def define_function_spaces(self, element_degree: int = 1): ) self.V_DG_0 = fem.functionspace(self.mesh.mesh, element_DG0) self.V_DG_1 = fem.functionspace(self.mesh.mesh, element_DG1) + self.V_CG_1 = fem.functionspace(self.mesh.mesh, ("CG", 1)) self.u = fem.Function(self.function_space) self.u_n = fem.Function(self.function_space) @@ -982,6 +999,10 @@ def post_processing(self): self._get_temperature_field_as_function() ) export.writer.write(float(self.t)) + elif isinstance(export, exports.CustomField): + # update internal function + export.function.interpolate(export.dolfinx_expression) + export.writer.write(float(self.t)) # TODO if export type derived quantity if isinstance(export, exports.SurfaceQuantity): From 78e491a383882fb4f794b2c16de56da647b50c4c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 17:28:42 -0400 Subject: [PATCH 04/13] backwards compatibility --- src/festim/exports/vtx.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index cf2005d66..03d9d87c6 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -5,6 +5,7 @@ import ufl from festim.species import Species +from festim.helpers import get_interpolation_points from festim.subdomain.volume_subdomain import VolumeSubdomain import inspect @@ -225,5 +226,5 @@ def set_dolfinx_expression( self.dolfinx_expression = fem.Expression( self.expression(**kwargs), - self.function.function_space.element.interpolation_points, + get_interpolation_points(self.function.function_space.element), ) From d30f07cbef22bfd60b47ee52163fd39a269a4e0a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 23 Apr 2026 21:00:24 -0400 Subject: [PATCH 05/13] support for mixed domain --- src/festim/exports/vtx.py | 22 +++++-- src/festim/hydrogen_transport_problem.py | 35 +++++++++-- test/test_vtx.py | 80 ++++++++++++++++++++++++ 3 files changed, 127 insertions(+), 10 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 03d9d87c6..1cbce66e4 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -216,14 +216,28 @@ def set_dolfinx_expression( x = ufl.SpatialCoordinate(time._ufl_domain) kwargs["x"] = x if "T" in arguments: - kwargs["T"] = temperature + if ( + isinstance(temperature, fem.Function) + and self.subdomain.sub_T is not None + ): + kwargs["T"] = ( + self.subdomain.sub_T + ) # NOTE I'm not sure that sub_T is updated at every time step + else: + kwargs["T"] = temperature # check if there are other arguments and if they are in species_dependent_value for arg in arguments: if arg in self.species_dependent_value: spe = self.species_dependent_value[arg] - kwargs[arg] = spe.post_processing_solution - # TODO handle discontinuous case and use self.subdomain - + if spe.subdomain_to_post_processing_solution: + kwargs[arg] = spe.subdomain_to_post_processing_solution[ + self.subdomain + ] + else: + kwargs[arg] = spe.post_processing_solution + assert kwargs[arg] is not None, ( + f"Argument {arg} not found in species_dependent_value" + ) self.dolfinx_expression = fem.Expression( self.expression(**kwargs), get_interpolation_points(self.function.function_space.element), diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index f29b1f38f..4a529315a 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1117,6 +1117,7 @@ def __init__( self.interfaces = interfaces or [] self.surface_to_volume = surface_to_volume or {} self.subdomain_to_species = {} # maps subdomain to species defined in it + self.subdomain_to_V_CG1 = {} @property def method_interface(self): @@ -1357,6 +1358,10 @@ def define_function_spaces( u = dolfinx.fem.Function(V) u_n = dolfinx.fem.Function(V) + self.subdomain_to_V_CG1[subdomain] = dolfinx.fem.functionspace( + subdomain.submesh, ("CG", 1) + ) + # store attributes in the subdomain object subdomain.u = u subdomain.u_n = u_n @@ -1684,7 +1689,22 @@ def initialise_exports(self): export.filename, mesh=functions[0].function_space.mesh, ) + elif isinstance(export, exports.CustomField): + # need to find an appropriate function space on the right submesh + V = self.subdomain_to_V_CG1[export.subdomain] + export.function = fem.Function(V) + export.set_dolfinx_expression( + temperature=self.temperature_fenics, # need to pass the right temperature + species=self.species, + time=self.t, + ) + export.writer = dolfinx.io.VTXWriter( + comm=export.function.function_space.mesh.comm, + filename=export.filename, + output=export.function, + engine="BP5", + ) # compute diffusivity function for surface fluxes # for the discontinuous case, we don't use D_global as in # HydrogenTransportProblem @@ -1726,11 +1746,11 @@ def post_processing(self): continue # handle VTX exports if isinstance(export, exports.ExportBaseClass): - if not isinstance(export, exports.VTXSpeciesExport): - raise NotImplementedError( - f"Export type {type(export)} not implemented" - ) - if isinstance(export, exports.VTXSpeciesExport): + if isinstance(export, exports.CustomField): + # update internal function + export.function.interpolate(export.dolfinx_expression) + export.writer.write(float(self.t)) + elif isinstance(export, exports.VTXSpeciesExport): if export._checkpoint: for species in export.field: post_processing_solution = ( @@ -1746,7 +1766,10 @@ def post_processing(self): ) else: export.writer.write(float(self.t)) - + else: + raise NotImplementedError( + f"Export type {type(export)} not implemented" + ) # handle derived quantities if isinstance(export, exports.SurfaceQuantity): if isinstance( diff --git a/test/test_vtx.py b/test/test_vtx.py index 335a9bda9..9daf49118 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -260,3 +260,83 @@ def test_custom_field(tmp_path, expression): my_model.initialise() my_model.run() + + +@pytest.mark.parametrize( + "expression", + [ + lambda x: x[0] + x[1] * 2, + lambda T: T + 1, + lambda c_A, c_B: c_A + c_B, + lambda c_A, T: c_A * T, + # FIXME: the following don't work (easily) because of https://github.com/FEniCS/dolfinx/issues/3207 + # lambda c_A, c_B, x: c_A * c_B + x[0], + # lambda c_A, T, x: c_A * T + x[0], + # lambda T, x: T + x[0], + ], +) +def test_custom_field_discontinuous(tmp_path, expression): + """ + Test custom field export functionality. + This test checks that a custom field can be created with various types of + expressions. + """ + + my_model = F.HydrogenTransportProblemDiscontinuous() + + mat = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + + vol = F.VolumeSubdomain(id=1, material=mat) + + top = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1)) + bottom = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[1], 0)) + left = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0)) + right = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[0], 1)) + + my_model.subdomains = [vol, top, bottom, left, right] + + my_model.surface_to_volume = {top: vol, bottom: vol, left: vol, right: vol} + + dolfinx_mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) + my_model.mesh = F.Mesh(dolfinx_mesh) + + A = F.Species("A", subdomains=my_model.volume_subdomains) + B = F.Species("B", subdomains=my_model.volume_subdomains) + C = F.Species("C", subdomains=my_model.volume_subdomains) + D = F.Species("D", subdomains=my_model.volume_subdomains) + + my_model.species = [A, B, C, D] + + my_model.boundary_conditions = ( + [ + F.FixedConcentrationBC(species=A, subdomain=top, value=1), + F.FixedConcentrationBC(species=B, subdomain=left, value=1), + ] + + [ + F.FixedConcentrationBC(species=C, subdomain=surf, value=0) + for surf in [top, bottom, left, right] + ] + + [ + F.FixedConcentrationBC(species=D, subdomain=surf, value=0) + for surf in [top, bottom, left, right] + ] + ) + + my_model.temperature = lambda x: 300 + 100 * x[0] + + my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) + + custom_field = F.CustomField( + filename=tmp_path / "custom_field.bp", + expression=expression, + species_dependent_value={"c_A": A, "c_B": B}, + subdomain=vol, + ) + + my_model.exports = [ + custom_field, + ] + + my_model.initialise() + + my_model.run() From 00193d3ece222b6a4ac852c0a4ec001d06c9ad21 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 10:52:15 -0400 Subject: [PATCH 06/13] docstrings --- src/festim/exports/vtx.py | 45 +++++++++++++++++++++++++++++++++------ 1 file changed, 38 insertions(+), 7 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 1cbce66e4..ef927b6de 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -1,13 +1,15 @@ +import inspect import warnings +from collections.abc import Callable from pathlib import Path +from typing import Union from dolfinx import fem, io import ufl -from festim.species import Species from festim.helpers import get_interpolation_points +from festim.species import Species from festim.subdomain.volume_subdomain import VolumeSubdomain -import inspect class ExportBaseClass: @@ -178,18 +180,47 @@ def get_functions(self) -> list[fem.Function]: class CustomField(ExportBaseClass): + """Export a custom field to a VTX file + + Args: + filename: The name of the output file + expression: A function evaluating the custom field. Positional + arguments of the function can be "t" (time), "x" (spatial coordinate), + "T" (temperature), or any key from the `species_dependent_value` dictionary. + species_dependent_value: A dictionary mapping argument names + in `expression` to Species objects. Defaults to None. + times: if provided, the field will be exported at these timesteps. Otherwise + exports at all timesteps. Defaults to None. + subdomain: The volume subdomain on which the custom + field is evaluated. Defaults to None. + checkpoint: If True, the export will be a checkpoint file using + adios4dolfinx and won't be readable by ParaView. Default is False. + + Attributes: + filename: The name of the output file + expression: A function evaluating the custom field. + species_dependent_value: A dictionary mapping argument names to Species objects. + subdomain: The volume subdomain on which the custom field is evaluated. + checkpoint: If True, the export will be a checkpoint file. + times: if provided, the field will be exported at these timesteps. Otherwise + exports at all timesteps. + function: the function containing the custom field values + writer: The VTXWriter object used to write the file + dolfinx_expression: the dolfinx expression used to evaluate the function + """ + function: fem.Function writer: io.VTXWriter dolfinx_expression: fem.Expression def __init__( self, - filename, - expression, - species_dependent_value=None, - times=None, + filename: Union[str, Path], + expression: Callable, + species_dependent_value: Union[dict[str, Species], None] = None, + times: Union[list[float], list[int], None] = None, subdomain: VolumeSubdomain = None, - checkpoint=False, + checkpoint: bool = False, ): super().__init__( filename=filename, From 598103380e86ef38c019789993e9e3d585228d4d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 10:53:09 -0400 Subject: [PATCH 07/13] ruff --- src/festim/exports/vtx.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index ef927b6de..9b74dc95e 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -4,8 +4,8 @@ from pathlib import Path from typing import Union -from dolfinx import fem, io import ufl +from dolfinx import fem, io from festim.helpers import get_interpolation_points from festim.species import Species From dc5c21f8e984e3059b133519ba028775bc2d1132 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 10:55:45 -0400 Subject: [PATCH 08/13] removed unused argument --- src/festim/exports/vtx.py | 14 +++++++++++++- src/festim/hydrogen_transport_problem.py | 2 -- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 9b74dc95e..a3baa29b7 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -212,6 +212,10 @@ class CustomField(ExportBaseClass): function: fem.Function writer: io.VTXWriter dolfinx_expression: fem.Expression + expression: Callable + species_dependent_value: dict[str, Species] + subdomain: VolumeSubdomain + checkpoint: bool def __init__( self, @@ -235,9 +239,17 @@ def __init__( def set_dolfinx_expression( self, temperature: fem.Constant | fem.Function, - species: list, time: fem.Constant, ): + """ + Set the dolfinx expression used to evaluate the custom field. This is done by + evaluating the user-provided expression with the appropriate arguments and using + the result to create a dolfinx expression. + + Args: + temperature: The temperature field to use in the expression + time: The time to use in the expression + """ # check arguments = inspect.signature(self.expression).parameters kwargs = {} diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 4a529315a..b065d8ded 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -453,7 +453,6 @@ def initialise_exports(self): export.function = fem.Function(self.V_CG_1) export.set_dolfinx_expression( temperature=self.temperature_fenics, - species=self.species, time=self.t, ) @@ -1695,7 +1694,6 @@ def initialise_exports(self): export.function = fem.Function(V) export.set_dolfinx_expression( temperature=self.temperature_fenics, # need to pass the right temperature - species=self.species, time=self.t, ) From d660bcac6cc8ba53341bb6a34b3ff0cedeef13a0 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 11:06:14 -0400 Subject: [PATCH 09/13] added comments --- src/festim/exports/vtx.py | 27 ++++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index a3baa29b7..20ac72767 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -250,8 +250,16 @@ def set_dolfinx_expression( temperature: The temperature field to use in the expression time: The time to use in the expression """ - # check + # check if we are in a mixed domain/discontinuous case + mixed_domain = any( + spe.subdomain_to_post_processing_solution + for spe in self.species_dependent_value.values() + ) or (self.subdomain.sub_T if self.subdomain else None) + + # get the arguments of the user-provided expression arguments = inspect.signature(self.expression).parameters + + # create a dictionary mapping the arguments to the appropriate values kwargs = {} if "t" in arguments: kwargs["t"] = time @@ -259,20 +267,18 @@ def set_dolfinx_expression( x = ufl.SpatialCoordinate(time._ufl_domain) kwargs["x"] = x if "T" in arguments: - if ( - isinstance(temperature, fem.Function) - and self.subdomain.sub_T is not None - ): - kwargs["T"] = ( - self.subdomain.sub_T - ) # NOTE I'm not sure that sub_T is updated at every time step + if isinstance(temperature, fem.Function) and mixed_domain: + # fem.Function in mixed domain/discontinuous case, use sub_T + # NOTE I'm not sure that sub_T is updated at every time step + kwargs["T"] = self.subdomain.sub_T else: + # else use the provided temperature kwargs["T"] = temperature # check if there are other arguments and if they are in species_dependent_value for arg in arguments: if arg in self.species_dependent_value: spe = self.species_dependent_value[arg] - if spe.subdomain_to_post_processing_solution: + if mixed_domain: kwargs[arg] = spe.subdomain_to_post_processing_solution[ self.subdomain ] @@ -281,6 +287,9 @@ def set_dolfinx_expression( assert kwargs[arg] is not None, ( f"Argument {arg} not found in species_dependent_value" ) + + # evaluate the user-provided expression with the appropriate arguments and create a + # dolfinx.fem.Expression self.dolfinx_expression = fem.Expression( self.expression(**kwargs), get_interpolation_points(self.function.function_space.element), From 6f647488e671755b4c9f3cd501416ec9be456483 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 11:25:50 -0400 Subject: [PATCH 10/13] x is the same domain as self.function --- src/festim/exports/vtx.py | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 20ac72767..4b8fa1f22 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -264,7 +264,7 @@ def set_dolfinx_expression( if "t" in arguments: kwargs["t"] = time if "x" in arguments: - x = ufl.SpatialCoordinate(time._ufl_domain) + x = ufl.SpatialCoordinate(self.function.function_space.mesh) kwargs["x"] = x if "T" in arguments: if isinstance(temperature, fem.Function) and mixed_domain: @@ -288,9 +288,34 @@ def set_dolfinx_expression( f"Argument {arg} not found in species_dependent_value" ) + self.check_valid_inputs(kwargs, mixed_domain) + # evaluate the user-provided expression with the appropriate arguments and create a # dolfinx.fem.Expression self.dolfinx_expression = fem.Expression( self.expression(**kwargs), get_interpolation_points(self.function.function_space.element), ) + + def check_valid_inputs(self, kwargs: dict, mixed_domain: bool): + """ + Check if we are in the mixed domain/discontinuous case and if the user-provided + expression is valid in this case. + dolfinx.fem.Expression does not support co-dim 0 submeshes and time is defined + on the parent mesh, so we cannot have time-dependent custom fields in the mixed + domain/discontinuous case. + + When https://github.com/FEniCS/dolfinx/issues/3207 is resolved we should be + able to support this + """ + + # check the domain of all kwargs and check that they are the same + + if mixed_domain and "t" in kwargs: + raise NotImplementedError( + "Time-dependent custom fields are not implemented in the case of a " + "mixed domain/discontinuous case." + "dolfinx.fem.Expression does not support co-dim 0 submeshes and time is" + "defined on the parent mesh." + "See https://github.com/FEniCS/dolfinx/issues/3207 for more details." + ) From 3c8c002ba4abc7712fe3df4d72f46b030541beab Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 11:26:04 -0400 Subject: [PATCH 11/13] more tests now pass + check for error --- test/test_vtx.py | 71 +++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 67 insertions(+), 4 deletions(-) diff --git a/test/test_vtx.py b/test/test_vtx.py index 9daf49118..d6452ad1a 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -269,10 +269,9 @@ def test_custom_field(tmp_path, expression): lambda T: T + 1, lambda c_A, c_B: c_A + c_B, lambda c_A, T: c_A * T, - # FIXME: the following don't work (easily) because of https://github.com/FEniCS/dolfinx/issues/3207 - # lambda c_A, c_B, x: c_A * c_B + x[0], - # lambda c_A, T, x: c_A * T + x[0], - # lambda T, x: T + x[0], + lambda c_A, c_B, x: c_A * c_B + x[0], + lambda c_A, T, x: c_A * T + x[0], + lambda T, x: T + x[0], ], ) def test_custom_field_discontinuous(tmp_path, expression): @@ -340,3 +339,67 @@ def test_custom_field_discontinuous(tmp_path, expression): my_model.initialise() my_model.run() + + +@pytest.mark.parametrize( + "expression", + [lambda c_A, t: c_A * t, lambda t: t], +) +def test_custom_field_not_implemented_error(expression): + my_model = F.HydrogenTransportProblemDiscontinuous() + + mat = F.Material(D_0=1, E_D=0, K_S_0=1, E_K_S=0) + + vol = F.VolumeSubdomain(id=1, material=mat) + + top = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[1], 1)) + bottom = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[1], 0)) + left = F.SurfaceSubdomain(id=3, locator=lambda x: np.isclose(x[0], 0)) + right = F.SurfaceSubdomain(id=4, locator=lambda x: np.isclose(x[0], 1)) + + my_model.subdomains = [vol, top, bottom, left, right] + + my_model.surface_to_volume = {top: vol, bottom: vol, left: vol, right: vol} + + dolfinx_mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) + my_model.mesh = F.Mesh(dolfinx_mesh) + + A = F.Species("A", subdomains=my_model.volume_subdomains) + B = F.Species("B", subdomains=my_model.volume_subdomains) + C = F.Species("C", subdomains=my_model.volume_subdomains) + D = F.Species("D", subdomains=my_model.volume_subdomains) + + my_model.species = [A, B, C, D] + + my_model.boundary_conditions = ( + [ + F.FixedConcentrationBC(species=A, subdomain=top, value=1), + F.FixedConcentrationBC(species=B, subdomain=left, value=1), + ] + + [ + F.FixedConcentrationBC(species=C, subdomain=surf, value=0) + for surf in [top, bottom, left, right] + ] + + [ + F.FixedConcentrationBC(species=D, subdomain=surf, value=0) + for surf in [top, bottom, left, right] + ] + ) + + my_model.temperature = lambda x: 300 + 100 * x[0] + + my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) + + custom_field = F.CustomField( + filename="custom_field.bp", + expression=expression, + species_dependent_value={"c_A": A, "c_B": B}, + subdomain=vol, + ) + + my_model.exports = [ + custom_field, + ] + + with pytest.raises(NotImplementedError): + my_model.initialise() From a63c369a6066f9557f81b4042d36cba7a81eb607 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 24 Apr 2026 13:35:38 -0400 Subject: [PATCH 12/13] added not implemented error for implicit species --- src/festim/exports/vtx.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index 4b8fa1f22..d2572f3cb 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -8,7 +8,7 @@ from dolfinx import fem, io from festim.helpers import get_interpolation_points -from festim.species import Species +from festim.species import Species, ImplicitSpecies from festim.subdomain.volume_subdomain import VolumeSubdomain @@ -278,6 +278,11 @@ def set_dolfinx_expression( for arg in arguments: if arg in self.species_dependent_value: spe = self.species_dependent_value[arg] + if isinstance(spe, ImplicitSpecies): + raise NotImplementedError( + "Custom fields depending on implicit species are not" + "implemented yet." + ) if mixed_domain: kwargs[arg] = spe.subdomain_to_post_processing_solution[ self.subdomain From 78f648c10121b46afa1c0b35329696078c39a60c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Apr 2026 12:51:14 -0400 Subject: [PATCH 13/13] renamed export --- src/festim/__init__.py | 2 +- src/festim/exports/__init__.py | 9 +++++++-- src/festim/exports/vtx.py | 2 +- src/festim/hydrogen_transport_problem.py | 8 ++++---- test/test_vtx.py | 6 +++--- 5 files changed, 16 insertions(+), 11 deletions(-) diff --git a/src/festim/__init__.py b/src/festim/__init__.py index 84aec9c37..392c35621 100644 --- a/src/festim/__init__.py +++ b/src/festim/__init__.py @@ -39,7 +39,7 @@ ExportBaseClass, VTXSpeciesExport, VTXTemperatureExport, - CustomField, + CustomFieldExport, ) from .exports.xdmf import XDMFExport from .heat_transfer_problem import HeatTransferProblem diff --git a/src/festim/exports/__init__.py b/src/festim/exports/__init__.py index 57c0e822f..e18892474 100644 --- a/src/festim/exports/__init__.py +++ b/src/festim/exports/__init__.py @@ -10,13 +10,18 @@ from .total_surface import TotalSurface from .total_volume import TotalVolume from .volume_quantity import VolumeQuantity -from .vtx import ExportBaseClass, VTXSpeciesExport, VTXTemperatureExport, CustomField +from .vtx import ( + ExportBaseClass, + VTXSpeciesExport, + VTXTemperatureExport, + CustomFieldExport, +) from .xdmf import XDMFExport __all__ = [ "AverageSurface", "AverageVolume", - "CustomField", + "CustomFieldExport", "ExportBaseClass", "MaximumSurface", "MaximumVolume", diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py index d2572f3cb..30293869e 100644 --- a/src/festim/exports/vtx.py +++ b/src/festim/exports/vtx.py @@ -179,7 +179,7 @@ def get_functions(self) -> list[fem.Function]: return outfiles -class CustomField(ExportBaseClass): +class CustomFieldExport(ExportBaseClass): """Export a custom field to a VTX file Args: diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 2115d0762..2491da42c 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -448,7 +448,7 @@ def initialise_exports(self): else: adios4dolfinx.write_mesh(export.filename, mesh=self.mesh.mesh) - elif isinstance(export, exports.CustomField): + elif isinstance(export, exports.CustomFieldExport): export.function = fem.Function(self.V_CG_1) export.set_dolfinx_expression( temperature=self.temperature_fenics, @@ -997,7 +997,7 @@ def post_processing(self): self._get_temperature_field_as_function() ) export.writer.write(float(self.t)) - elif isinstance(export, exports.CustomField): + elif isinstance(export, exports.CustomFieldExport): # update internal function export.function.interpolate(export.dolfinx_expression) export.writer.write(float(self.t)) @@ -1707,7 +1707,7 @@ def initialise_exports(self): export.filename, mesh=functions[0].function_space.mesh, ) - elif isinstance(export, exports.CustomField): + elif isinstance(export, exports.CustomFieldExport): # need to find an appropriate function space on the right submesh V = self.subdomain_to_V_CG1[export.subdomain] export.function = fem.Function(V) @@ -1763,7 +1763,7 @@ def post_processing(self): continue # handle VTX exports if isinstance(export, exports.ExportBaseClass): - if isinstance(export, exports.CustomField): + if isinstance(export, exports.CustomFieldExport): # update internal function export.function.interpolate(export.dolfinx_expression) export.writer.write(float(self.t)) diff --git a/test/test_vtx.py b/test/test_vtx.py index 737bbdeac..da7923f3f 100644 --- a/test/test_vtx.py +++ b/test/test_vtx.py @@ -246,7 +246,7 @@ def test_custom_field(tmp_path, expression): my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) - custom_field = F.CustomField( + custom_field = F.CustomFieldExport( filename=tmp_path / "custom_field.bp", expression=expression, species_dependent_value={"c_A": A, "c_B": B}, @@ -324,7 +324,7 @@ def test_custom_field_discontinuous(tmp_path, expression): my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) - custom_field = F.CustomField( + custom_field = F.CustomFieldExport( filename=tmp_path / "custom_field.bp", expression=expression, species_dependent_value={"c_A": A, "c_B": B}, @@ -389,7 +389,7 @@ def test_custom_field_not_implemented_error(expression): my_model.settings = F.Settings(transient=False, atol=1e-9, rtol=1e-9) - custom_field = F.CustomField( + custom_field = F.CustomFieldExport( filename="custom_field.bp", expression=expression, species_dependent_value={"c_A": A, "c_B": B},