Boundary condition handling

src/festim/boundary_conditions/flux_bc.py

@@ -187,6 +187,7 @@ def __init__(self, subdomain, value, species, species_dependent_value={}):
         super().__init__(subdomain=subdomain, value=value)
         self.species = species
         self.species_dependent_value = species_dependent_value
+        self._volume_subdomain = None
 
     def create_value_fenics(self, mesh, temperature, t: fem.Constant):
         """Creates the value of the boundary condition as a fenics object and sets it to
@@ -227,7 +228,12 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant):
                     kwargs["T"] = temperature
 
                 for name, species in self.species_dependent_value.items():
-                    kwargs[name] = species.concentration
+                    if species.concentration:
+                        kwargs[name] = species.concentration
+                    else:  # probably in discontinuous case
+                        kwargs[name] = species.subdomain_to_solution[
+                            self._volume_subdomain
+                        ]
 
                 self.value_fenics = self.value(**kwargs)
 

src/festim/boundary_conditions/surface_reaction.py

@@ -58,9 +58,22 @@ def create_value_fenics(self, mesh, temperature, t: fem.Constant):
         else:
             gas_pressure = self.gas_pressure
 
-        product_of_reactants = self.reactant[0].concentration
+        # initialise with the concentration of the first reactant
+        if self.reactant[0].concentration:
+            product_of_reactants = self.reactant[0].concentration
+        else:  # probably used in HydrogenTransportProblemDiscontinuous
+            product_of_reactants = self.reactant[0].subdomain_to_solution[
+                self._volume_subdomain
+            ]
+
+        # then compute the product of the concentrations of the other reactants
         for reactant in self.reactant[1:]:
-            product_of_reactants *= reactant.concentration
+            if reactant.concentration:
+                product_of_reactants *= reactant.concentration
+            else:  # probably used in HydrogenTransportProblemDiscontinuous
+                product_of_reactants *= reactant.subdomain_to_solution[
+                    self._volume_subdomain
+                ]
 
         self.value_fenics = kd * gas_pressure - kr * product_of_reactants
 

src/festim/hydrogen_transport_problem.py

@@ -296,6 +296,17 @@ def volume_meshtags(self, value):
         else:
             raise TypeError("value must be of type dolfinx.mesh.MeshTags")
 
+    @property
+    def _unpacked_bcs(self):
+        """Returns all boundary conditions, including fluxes from surface reactions"""
+        all_boundary_conditions = []
+        for bc in self.boundary_conditions:
+            if isinstance(bc, boundary_conditions.SurfaceReactionBC):
+                all_boundary_conditions.extend(bc.flux_bcs)
+            else:
+                all_boundary_conditions.append(bc)
+        return all_boundary_conditions
+
     def initialise(self):
         self.create_species_from_traps()
         self.define_function_spaces(element_degree=self.settings.element_degree)
@@ -641,17 +652,9 @@ def assign_functions_to_species(self):
             spe.test_function = sub_test_functions[idx]
 
     def define_boundary_conditions(self):
-        # @jhdark this all_bcs could be a property
-        # I just don't want to modify self.boundary_conditions
+        """Defines the boundary conditions of the model"""
 
-        # create all_bcs which includes all flux bcs from SurfaceReactionBC
-        all_bcs = self.boundary_conditions.copy()
-        for bc in self.boundary_conditions:
-            if isinstance(bc, boundary_conditions.SurfaceReactionBC):
-                all_bcs += bc.flux_bcs
-                all_bcs.remove(bc)
-
-        for bc in all_bcs:
+        for bc in self._unpacked_bcs:
             if isinstance(bc.species, str):
                 # if name of species is given then replace with species object
                 bc.species = _species.find_species_from_name(bc.species, self.species)
@@ -1322,6 +1325,13 @@ def convert_advection_term_to_fenics_objects(self):
                             function_space=V, t=self.t
                         )
 
+    def define_boundary_conditions(self):
+        for bc in self._unpacked_bcs:
+            if isinstance(bc, boundary_conditions.ParticleFluxBC):
+                bc._volume_subdomain = self.surface_to_volume[bc.subdomain]
+
+        super().define_boundary_conditions()
+
     def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain):
         """
         Creates the variational formulation for each subdomain and stores it in

test/test_complex_reaction.py

@@ -174,3 +174,43 @@ def test_reaction(k, p, c_A_0):
     model = model_test_reaction(stepsize=0.5, k=k, p=p, c_A_0=c_A_0)
     error = compute_error(model)
     assert error < 1e-2
+
+
+def test_surface_reaction_in_discontinuous_case():
+    """Test the SurfaceReactionBC in a discontinuous problem setup"""
+    my_model = F.HydrogenTransportProblemDiscontinuous()
+
+    my_model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=1000))
+
+    material = F.Material(D_0=1e-2, E_D=0, K_S_0=2.0, E_K_S=0)
+    subdomain = F.VolumeSubdomain1D(1, borders=[0, 1], material=material)
+    left_surface = F.SurfaceSubdomain1D(id=1, x=0)
+    right_surface = F.SurfaceSubdomain1D(id=2, x=1)
+
+    my_model.subdomains = [subdomain, left_surface, right_surface]
+    my_model.surface_to_volume = {
+        right_surface: subdomain,
+        left_surface: subdomain,
+    }
+
+    H = F.Species("H", mobile=True, subdomains=[subdomain])
+    my_model.species = [H]
+
+    my_model.boundary_conditions = [
+        F.FixedConcentrationBC(left_surface, value=100, species=H),
+        F.SurfaceReactionBC(
+            reactant=[H, H],
+            gas_pressure=1,
+            k_r0=1,
+            E_kr=0.2,
+            k_d0=0,
+            E_kd=0,
+            subdomain=right_surface,
+        ),
+    ]
+
+    my_model.temperature = 500
+
+    my_model.settings = F.Settings(atol=1e-5, rtol=1e-5, transient=False)
+
+    my_model.initialise()