From 031620af46e3218c9fff057eb8313c92af4484b5 Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Tue, 2 Jun 2026 09:46:47 -0600
Subject: [PATCH 01/14] new interface residual export

---
 src/festim/exports/vtx.py | 30 +++++++++++++++++++++++++++++-
 1 file changed, 29 insertions(+), 1 deletion(-)

diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py
index 26e67f1f0..26a68f386 100644
--- a/src/festim/exports/vtx.py
+++ b/src/festim/exports/vtx.py
@@ -345,6 +345,35 @@ def check_valid_inputs(self, kwargs: dict):
             )
 
 
+class VTXInterfaceResidualExport(ExportBaseClass):
+    """Export the concentration jump (residual) at an interface to a VTX file.
+
+    The residual is defined as ``c_0 - c_1``, where ``c_0`` and ``c_1`` are the
+    concentrations of the species on each side of the interface, interpolated
+    onto a submesh built from the interface facets.
+
+    Args:
+        field: The species whose interface residual is exported.
+        filename: The name of the output file.
+        interface: The interface between the two subdomains.
+        times: if provided, the field will be exported at these timesteps.
+            Otherwise exports at all timesteps. Defaults to None.
+
+    Attributes:
+        field: The species to export.
+        interface: The interface between the two subdomains.
+        function: The residual function on the interface submesh (set during
+            initialisation).
+        writer: The VTXWriter object used to write the file (set during
+            initialisation).
+    """
+
+    def __init__(self, field, filename, interface, times=None):
+        super().__init__(filename, ".bp", times)
+        self.field = field
+        self.interface = interface
+
+
 class ReactionRateExport(CustomFieldExport):
     """Export a reaction rate to a VTX file
 
@@ -370,7 +399,6 @@ def __init__(
         subdomain: VolumeSubdomain | None = None,
         checkpoint: bool = False,
     ):
-
         reactant_names = [reactant.name for reactant in reaction.reactant]
         if isinstance(reaction.product, list):
             product_names = [product.name for product in reaction.product]

From ecb87a59585210e84988b4e93637f402b624a6df Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Tue, 2 Jun 2026 09:46:57 -0600
Subject: [PATCH 02/14] add to inits

---
 src/festim/__init__.py         | 1 +
 src/festim/exports/__init__.py | 2 ++
 2 files changed, 3 insertions(+)

diff --git a/src/festim/__init__.py b/src/festim/__init__.py
index b2095f861..4a7b5244a 100644
--- a/src/festim/__init__.py
+++ b/src/festim/__init__.py
@@ -42,6 +42,7 @@
     CustomFieldExport,
     ExportBaseClass,
     ReactionRateExport,
+    VTXInterfaceResidualExport,
     VTXSpeciesExport,
     VTXTemperatureExport,
 )
diff --git a/src/festim/exports/__init__.py b/src/festim/exports/__init__.py
index ad3988734..597ba7c09 100644
--- a/src/festim/exports/__init__.py
+++ b/src/festim/exports/__init__.py
@@ -16,6 +16,7 @@
     CustomFieldExport,
     ExportBaseClass,
     ReactionRateExport,
+    VTXInterfaceResidualExport,
     VTXSpeciesExport,
     VTXTemperatureExport,
 )
@@ -38,6 +39,7 @@
     "SurfaceQuantity",
     "TotalSurface",
     "TotalVolume",
+    "VTXInterfaceResidualExport",
     "VTXSpeciesExport",
     "VTXTemperatureExport",
     "VolumeQuantity",

From 54d355b2ef9212fc8047dfb3d9d19f7dece12edd Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Tue, 2 Jun 2026 09:54:28 -0600
Subject: [PATCH 03/14] export initialisation and writing

---
 src/festim/hydrogen_transport_problem.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index 1e6682ccb..886a82b47 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -1738,6 +1738,9 @@ def initialise_exports(self):
                     self.temperature_fenics,
                     engine="BP5",
                 )
+            elif isinstance(export, exports.VTXInterfaceResidualExport):
+                export.initialise(self.temperature_fenics)
+
             elif isinstance(export, exports.CustomFieldExport):
                 # need to find an appropriate function space on the right submesh
                 V = self.subdomain_to_V_CG1[export.subdomain]
@@ -1853,6 +1856,8 @@ def post_processing(self):
                         export.writer.write(float(self.t))
                 elif isinstance(export, exports.VTXTemperatureExport):
                     export.writer.write(float(self.t))
+                elif isinstance(export, exports.VTXInterfaceResidualExport):
+                    export.write(float(self.t))
                 else:
                     raise NotImplementedError(
                         f"Export type {type(export)} not implemented"

From c79f043a649147c017c3887a4710342ae54f8626 Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Tue, 2 Jun 2026 09:54:45 -0600
Subject: [PATCH 04/14] evaluate residual

---
 src/festim/exports/vtx.py | 113 ++++++++++++++++++++++++++++++++++++--
 1 file changed, 107 insertions(+), 6 deletions(-)

diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py
index 26a68f386..8561fc31b 100644
--- a/src/festim/exports/vtx.py
+++ b/src/festim/exports/vtx.py
@@ -4,6 +4,8 @@
 from pathlib import Path
 from typing import Union
 
+import dolfinx
+import numpy as np
 import ufl
 from dolfinx import fem, io
 
@@ -346,11 +348,11 @@ def check_valid_inputs(self, kwargs: dict):
 
 
 class VTXInterfaceResidualExport(ExportBaseClass):
-    """Export the concentration jump (residual) at an interface to a VTX file.
+    """Export the interface condition residual ``c_0/K_S_0 - c_1/K_S_1`` to a VTX file.
 
-    The residual is defined as ``c_0 - c_1``, where ``c_0`` and ``c_1`` are the
-    concentrations of the species on each side of the interface, interpolated
-    onto a submesh built from the interface facets.
+    This quantity measures how well the penalty/Nitsche interface condition is
+    satisfied. It is zero when the condition holds exactly; a non-zero value
+    indicates the penalty term is too small.
 
     Args:
         field: The species whose interface residual is exported.
@@ -363,9 +365,9 @@ class VTXInterfaceResidualExport(ExportBaseClass):
         field: The species to export.
         interface: The interface between the two subdomains.
         function: The residual function on the interface submesh (set during
-            initialisation).
+            ``initialise``).
         writer: The VTXWriter object used to write the file (set during
-            initialisation).
+            ``initialise``).
     """
 
     def __init__(self, field, filename, interface, times=None):
@@ -373,6 +375,105 @@ def __init__(self, field, filename, interface, times=None):
         self.field = field
         self.interface = interface
 
+    def initialise(self, temperature_fenics):
+        """Create the interface submesh, interpolation data, and VTX writer.
+
+        Args:
+            temperature_fenics: The temperature field (``fem.Constant`` or
+                ``fem.Function`` on the parent mesh).
+        """
+        parent_mesh = self.interface.parent_mesh
+        fdim = parent_mesh.topology.dim - 1
+        interface_facets = self.interface.mt.find(self.interface.id)
+
+        interface_submesh, _, _, _ = dolfinx.mesh.create_submesh(
+            parent_mesh, fdim, interface_facets
+        )
+        V_interface = fem.functionspace(interface_submesh, ("CG", 1))
+
+        self._u_0 = fem.Function(V_interface)
+        self._u_1 = fem.Function(V_interface)
+        self.function = fem.Function(V_interface)
+        self.function.name = f"{self.field.name}_interface_residual"
+
+        imap = interface_submesh.topology.index_map(interface_submesh.topology.dim)
+        self._interface_cells = np.arange(
+            imap.size_local + imap.num_ghosts, dtype=np.int32
+        )
+
+        subdomain_0, subdomain_1 = self.interface.subdomains
+        V_0 = self.field.subdomain_to_post_processing_solution[
+            subdomain_0
+        ].function_space
+        V_1 = self.field.subdomain_to_post_processing_solution[
+            subdomain_1
+        ].function_space
+
+        self._interp_data_0 = fem.create_interpolation_data(
+            V_interface, V_0, self._interface_cells, padding=1e-11
+        )
+        self._interp_data_1 = fem.create_interpolation_data(
+            V_interface, V_1, self._interface_cells, padding=1e-11
+        )
+
+        self._K_S_0 = subdomain_0.material.get_K_S_0(self.field)
+        self._E_K_S_0 = subdomain_0.material.get_E_K_S(self.field)
+        self._K_S_1 = subdomain_1.material.get_K_S_0(self.field)
+        self._E_K_S_1 = subdomain_1.material.get_E_K_S(self.field)
+
+        self._temperature_fenics = temperature_fenics
+        if isinstance(temperature_fenics, fem.Constant):
+            self._T_func = None
+        else:
+            self._T_func = fem.Function(V_interface)
+            self._T_interp_data = fem.create_interpolation_data(
+                V_interface,
+                temperature_fenics.function_space,
+                self._interface_cells,
+                padding=1e-11,
+            )
+
+        self.writer = io.VTXWriter(
+            comm=interface_submesh.comm,
+            filename=self.filename,
+            output=self.function,
+            engine="BP5",
+        )
+
+    def write(self, t: float):
+        """Compute ``c_0/K_S_0 - c_1/K_S_1`` on the interface and write to file.
+
+        Args:
+            t: Current simulation time.
+        """
+        subdomain_0, subdomain_1 = self.interface.subdomains
+        u_0 = self.field.subdomain_to_post_processing_solution[subdomain_0]
+        u_1 = self.field.subdomain_to_post_processing_solution[subdomain_1]
+
+        self._u_0.interpolate_nonmatching(
+            u_0, self._interface_cells, interpolation_data=self._interp_data_0
+        )
+        self._u_1.interpolate_nonmatching(
+            u_1, self._interface_cells, interpolation_data=self._interp_data_1
+        )
+
+        if self._T_func is not None:
+            self._T_func.interpolate_nonmatching(
+                self._temperature_fenics,
+                self._interface_cells,
+                interpolation_data=self._T_interp_data,
+            )
+            T = self._T_func.x.array
+        else:
+            T = float(self._temperature_fenics)
+
+        K_S_0 = self._K_S_0 * np.exp(-self._E_K_S_0 / (_k_B * T))
+        K_S_1 = self._K_S_1 * np.exp(-self._E_K_S_1 / (_k_B * T))
+        self.function.x.array[:] = self._u_0.x.array / K_S_0 - self._u_1.x.array / K_S_1
+
+        print(np.max(np.absolute(self.function.x.array[:])))
+        self.writer.write(t)
+
 
 class ReactionRateExport(CustomFieldExport):
     """Export a reaction rate to a VTX file

From 1ec30863e9bc7357ece849e65732ce642f99a48d Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Tue, 2 Jun 2026 10:11:44 -0600
Subject: [PATCH 05/14] additional details in the doc strings

---
 src/festim/exports/vtx.py | 39 ++++++++++++++++++++++++++++++---------
 1 file changed, 30 insertions(+), 9 deletions(-)

diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py
index 8561fc31b..6a8bf589c 100644
--- a/src/festim/exports/vtx.py
+++ b/src/festim/exports/vtx.py
@@ -13,6 +13,7 @@
 from festim.helpers import get_interpolation_points
 from festim.reaction import Reaction
 from festim.species import ImplicitSpecies, Species
+from festim.subdomain.interface import Interface
 from festim.subdomain.volume_subdomain import VolumeSubdomain
 
 
@@ -364,23 +365,39 @@ class VTXInterfaceResidualExport(ExportBaseClass):
     Attributes:
         field: The species to export.
         interface: The interface between the two subdomains.
-        function: The residual function on the interface submesh (set during
-            ``initialise``).
-        writer: The VTXWriter object used to write the file (set during
-            ``initialise``).
+        function: Residual function on the interface submesh. Set by
+            ``initialise``.
+        writer: VTXWriter used to write the output file. Set by ``initialise``.
     """
 
-    def __init__(self, field, filename, interface, times=None):
+    field: Species
+    interface: Interface
+    function: fem.Function
+    writer: io.VTXWriter
+
+    def __init__(
+        self,
+        field: Species,
+        filename: str | Path,
+        interface: Interface,
+        times: list[float | int] | None = None,
+    ):
         super().__init__(filename, ".bp", times)
         self.field = field
         self.interface = interface
 
-    def initialise(self, temperature_fenics):
+    def initialise(self, temperature_fenics: fem.Constant | fem.Function) -> None:
         """Create the interface submesh, interpolation data, and VTX writer.
 
+        Called by the problem during ``initialise_exports``. Builds a CG1
+        function space on the interface submesh, pre-computes
+        ``create_interpolation_data`` for both subdomain concentrations and
+        (if space-dependent) temperature, and opens the VTXWriter.
+
         Args:
-            temperature_fenics: The temperature field (``fem.Constant`` or
-                ``fem.Function`` on the parent mesh).
+            temperature_fenics: Temperature field on the parent mesh. Either a
+                ``fem.Constant`` (uniform) or a ``fem.Function`` (spatially
+                varying).
         """
         parent_mesh = self.interface.parent_mesh
         fdim = parent_mesh.topology.dim - 1
@@ -440,9 +457,13 @@ def initialise(self, temperature_fenics):
             engine="BP5",
         )
 
-    def write(self, t: float):
+    def write(self, t: float) -> None:
         """Compute ``c_0/K_S_0 - c_1/K_S_1`` on the interface and write to file.
 
+        Interpolates the concentration from each subdomain onto the interface
+        submesh, evaluates K_S at the current temperature, computes the
+        residual, and writes the result via the VTXWriter.
+
         Args:
             t: Current simulation time.
         """

From c5e5f79f5d8d389f9486d0032ed7e55e3bc15df3 Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Tue, 2 Jun 2026 10:12:25 -0600
Subject: [PATCH 06/14] remove print statement

---
 src/festim/exports/vtx.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py
index 6a8bf589c..ba3688e16 100644
--- a/src/festim/exports/vtx.py
+++ b/src/festim/exports/vtx.py
@@ -492,7 +492,6 @@ def write(self, t: float) -> None:
         K_S_1 = self._K_S_1 * np.exp(-self._E_K_S_1 / (_k_B * T))
         self.function.x.array[:] = self._u_0.x.array / K_S_0 - self._u_1.x.array / K_S_1
 
-        print(np.max(np.absolute(self.function.x.array[:])))
         self.writer.write(t)
 
 

From eb9e2122e0dcc78262896d93e450d061f227a447 Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Tue, 2 Jun 2026 10:39:41 -0600
Subject: [PATCH 07/14] make function for getting interface term

---
 src/festim/subdomain/interface.py | 67 ++++++++++++++++++-------------
 1 file changed, 40 insertions(+), 27 deletions(-)

diff --git a/src/festim/subdomain/interface.py b/src/festim/subdomain/interface.py
index 8b568d40e..088e56a40 100644
--- a/src/festim/subdomain/interface.py
+++ b/src/festim/subdomain/interface.py
@@ -1,8 +1,10 @@
+from abc import ABC, abstractmethod
 from enum import Enum
 from typing import TYPE_CHECKING
 
 import dolfinx
 import numpy as np
+import numpy.typing as npt
 import ufl
 from dolfinx.cpp.fem import compute_integration_domains
 from packaging.version import Version
@@ -13,7 +15,38 @@
 if TYPE_CHECKING:
     from festim.species import Species
 
-from abc import ABC, abstractmethod
+
+def interface_condition_term(
+    c: ufl.core.expr.Expr | npt.NDArray[np.floating],
+    K_S: ufl.core.expr.Expr | npt.NDArray[np.floating],
+    law: SolubilityLaw,
+    other_law: SolubilityLaw,
+) -> ufl.core.expr.Expr | npt.NDArray[np.floating]:
+    """Return the interface condition expression for one side.
+
+    When both sides share the same solubility law the condition simplifies to
+    ``c/K_S`` on each side.  When the laws differ Henry is expressed as
+    ``c/K_S`` and Sievert as ``(c/K_S)^2``.
+
+    Args:
+        c: Concentration on this side.
+        K_S: Solubility coefficient on this side.
+        law: Solubility law for this side.
+        other_law: Solubility law for the other side.
+
+    Returns:
+        The interface condition term for this side.
+
+    Raises:
+        ValueError: If ``law`` is not ``HENRY`` or ``SIEVERT`` in the mixed-law
+            case.
+    """
+    if law == other_law or law == SolubilityLaw.HENRY:
+        return c / K_S
+    elif law == SolubilityLaw.SIEVERT:
+        return (c / K_S) ** 2
+    else:
+        raise ValueError(f"Unsupported solubility law: {law}")
 
 
 class InterfaceMethod(Enum):
@@ -393,32 +426,12 @@ def penalty_method(self, dS, species, temperature):
         u_0, u_1 = self.us(species)
         v_0, v_1 = self.vs(species)
         K_0, K_1 = self.Ks(species, temperature)
-        if subdomain_0.material.solubility_law == subdomain_1.material.solubility_law:
-            left = u_0 / K_0
-            right = u_1 / K_1
-        else:
-            match subdomain_0.material.solubility_law:
-                case SolubilityLaw.HENRY:
-                    left = u_0 / K_0
-                case SolubilityLaw.SIEVERT:
-                    left = (u_0 / K_0) ** 2
-                case _:
-                    raise ValueError(
-                        "Unsupported material law "
-                        + f"{subdomain_0.material.solubility_law}"
-                    )
-
-            match subdomain_1.material.solubility_law:
-                case SolubilityLaw.HENRY:
-                    right = u_1 / K_1
-                case SolubilityLaw.SIEVERT:
-                    right = (u_1 / K_1) ** 2
-                case _:
-                    raise ValueError(
-                        f"Unsupported material law "
-                        f"{subdomain_1.material.solubility_law}"
-                    )
-
+        law_0, law_1 = (
+            subdomain_0.material.solubility_law,
+            subdomain_1.material.solubility_law,
+        )
+        left = interface_condition_term(u_0, K_0, law_0, law_1)
+        right = interface_condition_term(u_1, K_1, law_1, law_0)
         equality = right - left
 
         F_0 = self.penalty_term * ufl.inner(equality, v_0) * dS(self.id)

From a00a9358f39aeaf72262ba7dd7821a994ed54735 Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Tue, 2 Jun 2026 10:39:51 -0600
Subject: [PATCH 08/14] use new interface function

---
 src/festim/exports/vtx.py | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py
index ba3688e16..faff90e7f 100644
--- a/src/festim/exports/vtx.py
+++ b/src/festim/exports/vtx.py
@@ -13,7 +13,7 @@
 from festim.helpers import get_interpolation_points
 from festim.reaction import Reaction
 from festim.species import ImplicitSpecies, Species
-from festim.subdomain.interface import Interface
+from festim.subdomain.interface import Interface, interface_condition_term
 from festim.subdomain.volume_subdomain import VolumeSubdomain
 
 
@@ -437,6 +437,8 @@ def initialise(self, temperature_fenics: fem.Constant | fem.Function) -> None:
         self._E_K_S_0 = subdomain_0.material.get_E_K_S(self.field)
         self._K_S_1 = subdomain_1.material.get_K_S_0(self.field)
         self._E_K_S_1 = subdomain_1.material.get_E_K_S(self.field)
+        self._law_0 = subdomain_0.material.solubility_law
+        self._law_1 = subdomain_1.material.solubility_law
 
         self._temperature_fenics = temperature_fenics
         if isinstance(temperature_fenics, fem.Constant):
@@ -490,7 +492,13 @@ def write(self, t: float) -> None:
 
         K_S_0 = self._K_S_0 * np.exp(-self._E_K_S_0 / (_k_B * T))
         K_S_1 = self._K_S_1 * np.exp(-self._E_K_S_1 / (_k_B * T))
-        self.function.x.array[:] = self._u_0.x.array / K_S_0 - self._u_1.x.array / K_S_1
+        left = interface_condition_term(
+            self._u_0.x.array, K_S_0, self._law_0, self._law_1
+        )
+        right = interface_condition_term(
+            self._u_1.x.array, K_S_1, self._law_1, self._law_0
+        )
+        self.function.x.array[:] = right - left
 
         self.writer.write(t)
 

From 359b16a49330cec9316233da93667fcdc9d81c6f Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Tue, 2 Jun 2026 11:13:40 -0600
Subject: [PATCH 09/14] update doc strings, showing residual depending on
 interface condition

---
 src/festim/exports/vtx.py | 31 +++++++++++++++++++++++++------
 1 file changed, 25 insertions(+), 6 deletions(-)

diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py
index faff90e7f..7d7207065 100644
--- a/src/festim/exports/vtx.py
+++ b/src/festim/exports/vtx.py
@@ -349,11 +349,21 @@ def check_valid_inputs(self, kwargs: dict):
 
 
 class VTXInterfaceResidualExport(ExportBaseClass):
-    """Export the interface condition residual ``c_0/K_S_0 - c_1/K_S_1`` to a VTX file.
+    """Export the interface condition residual to a VTX file.
 
     This quantity measures how well the penalty/Nitsche interface condition is
-    satisfied. It is zero when the condition holds exactly; a non-zero value
-    indicates the penalty term is too small.
+    satisfied at the interface. It is zero when the condition holds exactly; the lower
+    the value, the better.
+
+    The residual is ``right - left`` where each side's term depends on the
+    solubility laws of the two subdomains:
+
+    - **Same law on both sides** (Henry-Henry or Sievert-Sievert):
+      ``residual = c_1/K_S_1 - c_0/K_S_0``
+    - **Henry (side 0) - Sievert (side 1)**:
+      ``residual = (c_1/K_S_1)^2 - c_0/K_S_0``
+    - **Sievert (side 0) - Henry (side 1)**:
+      ``residual = c_1/K_S_1 - (c_0/K_S_0)^2``
 
     Args:
         field: The species whose interface residual is exported.
@@ -364,6 +374,7 @@ class VTXInterfaceResidualExport(ExportBaseClass):
 
     Attributes:
         field: The species to export.
+        filename: The name of the output file.
         interface: The interface between the two subdomains.
         function: Residual function on the interface submesh. Set by
             ``initialise``.
@@ -372,6 +383,8 @@ class VTXInterfaceResidualExport(ExportBaseClass):
 
     field: Species
     interface: Interface
+    times: list[float] | list[int] | None
+
     function: fem.Function
     writer: io.VTXWriter
 
@@ -460,11 +473,17 @@ def initialise(self, temperature_fenics: fem.Constant | fem.Function) -> None:
         )
 
     def write(self, t: float) -> None:
-        """Compute ``c_0/K_S_0 - c_1/K_S_1`` on the interface and write to file.
+        """Compute the interface condition residual and write to file.
 
         Interpolates the concentration from each subdomain onto the interface
-        submesh, evaluates K_S at the current temperature, computes the
-        residual, and writes the result via the VTXWriter.
+        submesh, evaluates ``K_S = K_S_0 * exp(-E_K_S / (k_B * T))`` at the
+        current temperature, then computes ``right - left`` via
+        :func:`interface_condition_term`:
+
+        - **Same law or Henry (side 0)**: ``left = c_0 / K_S_0``
+        - **Sievert (side 0, mixed only)**: ``left = (c_0 / K_S_0)^2``
+
+        and symmetrically for ``right``.
 
         Args:
             t: Current simulation time.

From fc52f6757951446565b81324b163797fcda9cc2b Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Mon, 22 Jun 2026 11:07:07 +0200
Subject: [PATCH 10/14] use dolin expression instead

---
 src/festim/exports/vtx.py                | 18 ++++++++----------
 src/festim/hydrogen_transport_problem.py |  4 +++-
 2 files changed, 11 insertions(+), 11 deletions(-)

diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py
index 7d7207065..5fac4c7a5 100644
--- a/src/festim/exports/vtx.py
+++ b/src/festim/exports/vtx.py
@@ -472,8 +472,8 @@ def initialise(self, temperature_fenics: fem.Constant | fem.Function) -> None:
             engine="BP5",
         )
 
-    def write(self, t: float) -> None:
-        """Compute the interface condition residual and write to file.
+    def set_dolfinx_expression(self) -> None:
+        """Compute the interface condition residual.
 
         Interpolates the concentration from each subdomain onto the interface
         submesh, evaluates ``K_S = K_S_0 * exp(-E_K_S / (k_B * T))`` at the
@@ -511,15 +511,13 @@ def write(self, t: float) -> None:
 
         K_S_0 = self._K_S_0 * np.exp(-self._E_K_S_0 / (_k_B * T))
         K_S_1 = self._K_S_1 * np.exp(-self._E_K_S_1 / (_k_B * T))
-        left = interface_condition_term(
-            self._u_0.x.array, K_S_0, self._law_0, self._law_1
-        )
-        right = interface_condition_term(
-            self._u_1.x.array, K_S_1, self._law_1, self._law_0
-        )
-        self.function.x.array[:] = right - left
+        left = interface_condition_term(self._u_0, K_S_0, self._law_0, self._law_1)
+        right = interface_condition_term(self._u_1, K_S_1, self._law_1, self._law_0)
 
-        self.writer.write(t)
+        self.dolfinx_expression = fem.Expression(
+            right - left,
+            get_interpolation_points(self.function.function_space.element),
+        )
 
 
 class ReactionRateExport(CustomFieldExport):
diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index 886a82b47..23e55b3bc 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -1857,7 +1857,9 @@ def post_processing(self):
                 elif isinstance(export, exports.VTXTemperatureExport):
                     export.writer.write(float(self.t))
                 elif isinstance(export, exports.VTXInterfaceResidualExport):
-                    export.write(float(self.t))
+                    export.set_dolfinx_expression()
+                    export.function.interpolate(export.dolfinx_expression)
+                    export.writer.write(float(self.t))
                 else:
                     raise NotImplementedError(
                         f"Export type {type(export)} not implemented"

From ac1186f06f42837a9a7ae7b3a985c1d66c27fd9c Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Mon, 22 Jun 2026 16:33:23 +0200
Subject: [PATCH 11/14] refactoring and export consituents of the residual

---
 src/festim/exports/vtx.py | 67 +++++++++++++++++++++++++++++----------
 1 file changed, 50 insertions(+), 17 deletions(-)

diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py
index 5fac4c7a5..f267f1c4f 100644
--- a/src/festim/exports/vtx.py
+++ b/src/festim/exports/vtx.py
@@ -371,6 +371,8 @@ class VTXInterfaceResidualExport(ExportBaseClass):
         interface: The interface between the two subdomains.
         times: if provided, the field will be exported at these timesteps.
             Otherwise exports at all timesteps. Defaults to None.
+        export_constituents: If True, also export the left and right terms of the
+            residual as separate fields. Defaults to False.
 
     Attributes:
         field: The species to export.
@@ -384,6 +386,7 @@ class VTXInterfaceResidualExport(ExportBaseClass):
     field: Species
     interface: Interface
     times: list[float] | list[int] | None
+    export_constituents: bool
 
     function: fem.Function
     writer: io.VTXWriter
@@ -394,10 +397,12 @@ def __init__(
         filename: str | Path,
         interface: Interface,
         times: list[float | int] | None = None,
+        export_constituents: bool = False,
     ):
         super().__init__(filename, ".bp", times)
         self.field = field
         self.interface = interface
+        self.export_constituents = export_constituents
 
     def initialise(self, temperature_fenics: fem.Constant | fem.Function) -> None:
         """Create the interface submesh, interpolation data, and VTX writer.
@@ -421,8 +426,8 @@ def initialise(self, temperature_fenics: fem.Constant | fem.Function) -> None:
         )
         V_interface = fem.functionspace(interface_submesh, ("CG", 1))
 
-        self._u_0 = fem.Function(V_interface)
-        self._u_1 = fem.Function(V_interface)
+        self._u_0_interface = fem.Function(V_interface)
+        self._u_1_interface = fem.Function(V_interface)
         self.function = fem.Function(V_interface)
         self.function.name = f"{self.field.name}_interface_residual"
 
@@ -446,13 +451,6 @@ def initialise(self, temperature_fenics: fem.Constant | fem.Function) -> None:
             V_interface, V_1, self._interface_cells, padding=1e-11
         )
 
-        self._K_S_0 = subdomain_0.material.get_K_S_0(self.field)
-        self._E_K_S_0 = subdomain_0.material.get_E_K_S(self.field)
-        self._K_S_1 = subdomain_1.material.get_K_S_0(self.field)
-        self._E_K_S_1 = subdomain_1.material.get_E_K_S(self.field)
-        self._law_0 = subdomain_0.material.solubility_law
-        self._law_1 = subdomain_1.material.solubility_law
-
         self._temperature_fenics = temperature_fenics
         if isinstance(temperature_fenics, fem.Constant):
             self._T_func = None
@@ -472,6 +470,26 @@ def initialise(self, temperature_fenics: fem.Constant | fem.Function) -> None:
             engine="BP5",
         )
 
+        if self.export_constituents:
+            u0_filename = self.filename.with_name(
+                self.filename.stem + "_u0" + self.filename.suffix
+            )
+            u1_filename = self.filename.with_name(
+                self.filename.stem + "_u1" + self.filename.suffix
+            )
+            self.writer_u0 = io.VTXWriter(
+                comm=interface_submesh.comm,
+                filename=u0_filename,
+                output=self._u_0_interface,
+                engine="BP5",
+            )
+            self.writer_u1 = io.VTXWriter(
+                comm=interface_submesh.comm,
+                filename=u1_filename,
+                output=self._u_1_interface,
+                engine="BP5",
+            )
+
     def set_dolfinx_expression(self) -> None:
         """Compute the interface condition residual.
 
@@ -492,10 +510,10 @@ def set_dolfinx_expression(self) -> None:
         u_0 = self.field.subdomain_to_post_processing_solution[subdomain_0]
         u_1 = self.field.subdomain_to_post_processing_solution[subdomain_1]
 
-        self._u_0.interpolate_nonmatching(
+        self._u_0_interface.interpolate_nonmatching(
             u_0, self._interface_cells, interpolation_data=self._interp_data_0
         )
-        self._u_1.interpolate_nonmatching(
+        self._u_1_interface.interpolate_nonmatching(
             u_1, self._interface_cells, interpolation_data=self._interp_data_1
         )
 
@@ -509,13 +527,28 @@ def set_dolfinx_expression(self) -> None:
         else:
             T = float(self._temperature_fenics)
 
-        K_S_0 = self._K_S_0 * np.exp(-self._E_K_S_0 / (_k_B * T))
-        K_S_1 = self._K_S_1 * np.exp(-self._E_K_S_1 / (_k_B * T))
-        left = interface_condition_term(self._u_0, K_S_0, self._law_0, self._law_1)
-        right = interface_condition_term(self._u_1, K_S_1, self._law_1, self._law_0)
+        K_S_0 = subdomain_0.material.get_K_S_0(self.field) * np.exp(
+            -subdomain_0.material.get_E_K_S(self.field) / (_k_B * T)
+        )
+        u0 = interface_condition_term(
+            self._u_0_interface,
+            K_S_0,
+            subdomain_0.material.solubility_law,
+            subdomain_1.material.solubility_law,
+        )
 
-        self.dolfinx_expression = fem.Expression(
-            right - left,
+        K_S_1 = subdomain_1.material.get_K_S_0(self.field) * np.exp(
+            -subdomain_1.material.get_E_K_S(self.field) / (_k_B * T)
+        )
+        u1 = interface_condition_term(
+            self._u_1_interface,
+            K_S_1,
+            subdomain_1.material.solubility_law,
+            subdomain_0.material.solubility_law,
+        )
+
+        self.residual_expr = fem.Expression(
+            u0 - u1,
             get_interpolation_points(self.function.function_space.element),
         )
 

From 757151518c89f7dbbd583a917fe424de9f4ad256 Mon Sep 17 00:00:00 2001
From: jhdark <darkj385@mit.edu>
Date: Mon, 22 Jun 2026 16:33:47 +0200
Subject: [PATCH 12/14] refactoring

---
 src/festim/hydrogen_transport_problem.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index 23e55b3bc..3fc3ffc50 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -1858,8 +1858,11 @@ def post_processing(self):
                     export.writer.write(float(self.t))
                 elif isinstance(export, exports.VTXInterfaceResidualExport):
                     export.set_dolfinx_expression()
-                    export.function.interpolate(export.dolfinx_expression)
+                    export.function.interpolate(export.residual_expr)
                     export.writer.write(float(self.t))
+                    if export.export_constituents:
+                        export.writer_u0.write(float(self.t))
+                        export.writer_u1.write(float(self.t))
                 else:
                     raise NotImplementedError(
                         f"Export type {type(export)} not implemented"

From 006e4c3531d7580ae5bf615395d380ac1c9bae0b Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Tue, 23 Jun 2026 09:45:00 -0400
Subject: [PATCH 13/14] all residuals + ufl instead of np

---
 src/festim/exports/vtx.py                | 47 ++++++++----------------
 src/festim/hydrogen_transport_problem.py |  3 --
 2 files changed, 15 insertions(+), 35 deletions(-)

diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py
index f267f1c4f..febe86cc6 100644
--- a/src/festim/exports/vtx.py
+++ b/src/festim/exports/vtx.py
@@ -371,8 +371,6 @@ class VTXInterfaceResidualExport(ExportBaseClass):
         interface: The interface between the two subdomains.
         times: if provided, the field will be exported at these timesteps.
             Otherwise exports at all timesteps. Defaults to None.
-        export_constituents: If True, also export the left and right terms of the
-            residual as separate fields. Defaults to False.
 
     Attributes:
         field: The species to export.
@@ -386,7 +384,6 @@ class VTXInterfaceResidualExport(ExportBaseClass):
     field: Species
     interface: Interface
     times: list[float] | list[int] | None
-    export_constituents: bool
 
     function: fem.Function
     writer: io.VTXWriter
@@ -397,12 +394,10 @@ def __init__(
         filename: str | Path,
         interface: Interface,
         times: list[float | int] | None = None,
-        export_constituents: bool = False,
     ):
         super().__init__(filename, ".bp", times)
         self.field = field
         self.interface = interface
-        self.export_constituents = export_constituents
 
     def initialise(self, temperature_fenics: fem.Constant | fem.Function) -> None:
         """Create the interface submesh, interpolation data, and VTX writer.
@@ -426,8 +421,12 @@ def initialise(self, temperature_fenics: fem.Constant | fem.Function) -> None:
         )
         V_interface = fem.functionspace(interface_submesh, ("CG", 1))
 
-        self._u_0_interface = fem.Function(V_interface)
-        self._u_1_interface = fem.Function(V_interface)
+        self._u_0_interface = fem.Function(
+            V_interface, name=f"{self.field.name}_interface_u0"
+        )
+        self._u_1_interface = fem.Function(
+            V_interface, name=f"{self.field.name}_interface_u1"
+        )
         self.function = fem.Function(V_interface)
         self.function.name = f"{self.field.name}_interface_residual"
 
@@ -466,30 +465,10 @@ def initialise(self, temperature_fenics: fem.Constant | fem.Function) -> None:
         self.writer = io.VTXWriter(
             comm=interface_submesh.comm,
             filename=self.filename,
-            output=self.function,
+            output=[self.function, self._u_0_interface, self._u_1_interface],
             engine="BP5",
         )
 
-        if self.export_constituents:
-            u0_filename = self.filename.with_name(
-                self.filename.stem + "_u0" + self.filename.suffix
-            )
-            u1_filename = self.filename.with_name(
-                self.filename.stem + "_u1" + self.filename.suffix
-            )
-            self.writer_u0 = io.VTXWriter(
-                comm=interface_submesh.comm,
-                filename=u0_filename,
-                output=self._u_0_interface,
-                engine="BP5",
-            )
-            self.writer_u1 = io.VTXWriter(
-                comm=interface_submesh.comm,
-                filename=u1_filename,
-                output=self._u_1_interface,
-                engine="BP5",
-            )
-
     def set_dolfinx_expression(self) -> None:
         """Compute the interface condition residual.
 
@@ -510,6 +489,10 @@ def set_dolfinx_expression(self) -> None:
         u_0 = self.field.subdomain_to_post_processing_solution[subdomain_0]
         u_1 = self.field.subdomain_to_post_processing_solution[subdomain_1]
 
+        # FIXME: once we support dolfinx 0.11 we should be able to mix parent and
+        # sub-meshes in the same expression
+
+        # NOTE: we need to do this multiple times to update the residual
         self._u_0_interface.interpolate_nonmatching(
             u_0, self._interface_cells, interpolation_data=self._interp_data_0
         )
@@ -523,11 +506,11 @@ def set_dolfinx_expression(self) -> None:
                 self._interface_cells,
                 interpolation_data=self._T_interp_data,
             )
-            T = self._T_func.x.array
+            T = self._T_func
         else:
-            T = float(self._temperature_fenics)
+            T = self._temperature_fenics
 
-        K_S_0 = subdomain_0.material.get_K_S_0(self.field) * np.exp(
+        K_S_0 = subdomain_0.material.get_K_S_0(self.field) * ufl.exp(
             -subdomain_0.material.get_E_K_S(self.field) / (_k_B * T)
         )
         u0 = interface_condition_term(
@@ -537,7 +520,7 @@ def set_dolfinx_expression(self) -> None:
             subdomain_1.material.solubility_law,
         )
 
-        K_S_1 = subdomain_1.material.get_K_S_0(self.field) * np.exp(
+        K_S_1 = subdomain_1.material.get_K_S_0(self.field) * ufl.exp(
             -subdomain_1.material.get_E_K_S(self.field) / (_k_B * T)
         )
         u1 = interface_condition_term(
diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index 3fc3ffc50..26bc365ed 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -1860,9 +1860,6 @@ def post_processing(self):
                     export.set_dolfinx_expression()
                     export.function.interpolate(export.residual_expr)
                     export.writer.write(float(self.t))
-                    if export.export_constituents:
-                        export.writer_u0.write(float(self.t))
-                        export.writer_u1.write(float(self.t))
                 else:
                     raise NotImplementedError(
                         f"Export type {type(export)} not implemented"

From 523372e50043bc00943d49c4a1db44a7db475ae8 Mon Sep 17 00:00:00 2001
From: RemDelaporteMathurin <rdelaportemathurin@gmail.com>
Date: Tue, 23 Jun 2026 10:08:56 -0400
Subject: [PATCH 14/14] export f0 and f1 instead of c0 and c1

---
 src/festim/exports/vtx.py                | 50 ++++++++++++++++--------
 src/festim/hydrogen_transport_problem.py |  4 ++
 2 files changed, 38 insertions(+), 16 deletions(-)

diff --git a/src/festim/exports/vtx.py b/src/festim/exports/vtx.py
index febe86cc6..66aa741a3 100644
--- a/src/festim/exports/vtx.py
+++ b/src/festim/exports/vtx.py
@@ -421,11 +421,18 @@ def initialise(self, temperature_fenics: fem.Constant | fem.Function) -> None:
         )
         V_interface = fem.functionspace(interface_submesh, ("CG", 1))
 
-        self._u_0_interface = fem.Function(
-            V_interface, name=f"{self.field.name}_interface_u0"
+        self._c_0_interface = fem.Function(
+            V_interface, name=f"{self.field.name}_interface_c0"
         )
-        self._u_1_interface = fem.Function(
-            V_interface, name=f"{self.field.name}_interface_u1"
+        self._c_1_interface = fem.Function(
+            V_interface, name=f"{self.field.name}_interface_c1"
+        )
+
+        self._f_0_interface = fem.Function(
+            V_interface, name=f"{self.field.name}_interface_f0"
+        )
+        self._f_1_interface = fem.Function(
+            V_interface, name=f"{self.field.name}_interface_f1"
         )
         self.function = fem.Function(V_interface)
         self.function.name = f"{self.field.name}_interface_residual"
@@ -465,7 +472,7 @@ def initialise(self, temperature_fenics: fem.Constant | fem.Function) -> None:
         self.writer = io.VTXWriter(
             comm=interface_submesh.comm,
             filename=self.filename,
-            output=[self.function, self._u_0_interface, self._u_1_interface],
+            output=[self.function, self._f_0_interface, self._f_1_interface],
             engine="BP5",
         )
 
@@ -486,20 +493,23 @@ def set_dolfinx_expression(self) -> None:
             t: Current simulation time.
         """
         subdomain_0, subdomain_1 = self.interface.subdomains
-        u_0 = self.field.subdomain_to_post_processing_solution[subdomain_0]
-        u_1 = self.field.subdomain_to_post_processing_solution[subdomain_1]
+
+        # Get the concentration from each subdomain on the interface submesh
+        c_0 = self.field.subdomain_to_post_processing_solution[subdomain_0]
+        c_1 = self.field.subdomain_to_post_processing_solution[subdomain_1]
 
         # FIXME: once we support dolfinx 0.11 we should be able to mix parent and
         # sub-meshes in the same expression
 
         # NOTE: we need to do this multiple times to update the residual
-        self._u_0_interface.interpolate_nonmatching(
-            u_0, self._interface_cells, interpolation_data=self._interp_data_0
+        self._c_0_interface.interpolate_nonmatching(
+            c_0, self._interface_cells, interpolation_data=self._interp_data_0
         )
-        self._u_1_interface.interpolate_nonmatching(
-            u_1, self._interface_cells, interpolation_data=self._interp_data_1
+        self._c_1_interface.interpolate_nonmatching(
+            c_1, self._interface_cells, interpolation_data=self._interp_data_1
         )
 
+        # Get the temperature on the interface submesh
         if self._T_func is not None:
             self._T_func.interpolate_nonmatching(
                 self._temperature_fenics,
@@ -510,11 +520,12 @@ def set_dolfinx_expression(self) -> None:
         else:
             T = self._temperature_fenics
 
+        # Compute f_i = (c_i / K_S_i) ^ {1, 2} on the interface submesh
         K_S_0 = subdomain_0.material.get_K_S_0(self.field) * ufl.exp(
             -subdomain_0.material.get_E_K_S(self.field) / (_k_B * T)
         )
-        u0 = interface_condition_term(
-            self._u_0_interface,
+        f0 = interface_condition_term(
+            self._c_0_interface,
             K_S_0,
             subdomain_0.material.solubility_law,
             subdomain_1.material.solubility_law,
@@ -523,15 +534,22 @@ def set_dolfinx_expression(self) -> None:
         K_S_1 = subdomain_1.material.get_K_S_0(self.field) * ufl.exp(
             -subdomain_1.material.get_E_K_S(self.field) / (_k_B * T)
         )
-        u1 = interface_condition_term(
-            self._u_1_interface,
+        f1 = interface_condition_term(
+            self._c_1_interface,
             K_S_1,
             subdomain_1.material.solubility_law,
             subdomain_0.material.solubility_law,
         )
 
+        self.f_0_expr = fem.Expression(
+            f0, get_interpolation_points(self.function.function_space.element)
+        )
+        self.f_1_expr = fem.Expression(
+            f1, get_interpolation_points(self.function.function_space.element)
+        )
+
         self.residual_expr = fem.Expression(
-            u0 - u1,
+            f0 - f1,
             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 26bc365ed..68705043e 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -1857,8 +1857,12 @@ def post_processing(self):
                 elif isinstance(export, exports.VTXTemperatureExport):
                     export.writer.write(float(self.t))
                 elif isinstance(export, exports.VTXInterfaceResidualExport):
+                    # FIXME: don't need to remake the whole expression everytime just
+                    # need to update "sub" functions
                     export.set_dolfinx_expression()
                     export.function.interpolate(export.residual_expr)
+                    export._f_0_interface.interpolate(export.f_0_expr)
+                    export._f_1_interface.interpolate(export.f_1_expr)
                     export.writer.write(float(self.t))
                 else:
                     raise NotImplementedError(