removed the need for override_solution

Author: RemDelaporteMathurin

src/festim/hydrogen_transport_problem.py

@@ -1467,9 +1467,6 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain):
             if reaction.volume != subdomain:
                 continue
 
-            # TODO remove
-            # temporarily overide the solution to the one of the subdomain
-            self.override_solution_attributes(reaction)
             # reactant
             for reactant in reaction.reactant:
                 if isinstance(reactant, _species.Species):
@@ -1592,27 +1589,6 @@ def create_formulation(self):
             },
         )
 
-    def override_solution_attributes(self, reaction: _reaction.Reaction):
-        """
-        Reaction.reaction_term() relies on the .solution attribute of the species
-        however, in the discontinuous class, this attribute doesn't really make sense
-        since there is one solution per subdomain.
-        Therefore we temporarily override the .solution attribute based on the reactants,
-        products, and `others` if there are implicit species
-        """
-        list_of_species_to_override = reaction.reactant + reaction.product
-
-        # check if we have implicit species:
-        for reactant in reaction.reactant:
-            if isinstance(reactant, _species.ImplicitSpecies):
-                for other_spe in reactant.others:
-                    if other_spe not in list_of_species_to_override:
-                        list_of_species_to_override.append(other_spe)
-
-        for species in list_of_species_to_override:
-            if isinstance(species, _species.Species):
-                species.solution = species.subdomain_to_solution[reaction.volume]
-
     def create_solver(self):
         if Version(dolfinx.__version__) == Version("0.9.0"):
             self.solver = BlockedNewtonSolver(

src/festim/reaction.py

@@ -134,6 +134,24 @@ def reaction_term(
             The reaction term to be used in a formulation.
         """
 
+        # make sure products is a list
+        products = self.product if isinstance(self.product, list) else [self.product]
+
+        # detect if mixed_domain
+        mixed_domain = any(
+            isinstance(reactant, _Species) and reactant.subdomain_to_solution != {}
+            for reactant in self.reactant
+        ) or any(
+            isinstance(product, _Species) and product.subdomain_to_solution != {}
+            for product in products
+        )
+
+        def get_concentration(species):
+            if mixed_domain:
+                return species.concentration_submesh(self.volume)
+            else:
+                return species.concentration
+
         if self.product == []:
             if self.p_0 is not None:
                 raise ValueError(
@@ -155,8 +173,6 @@ def reaction_term(
                     "E_p cannot be None when reaction products are present."
                 )
 
-        products = self.product if isinstance(self.product, list) else [self.product]
-
         # reaction rates
         k = self.k_0 * exp(-self.E_k / (_k_B * temperature))
 
@@ -173,18 +189,22 @@ def reaction_term(
             assert len(reactant_concentrations) == len(reactants)
             for i, reactant in enumerate(reactants):
                 if reactant_concentrations[i] is None:
-                    reactant_concentrations[i] = reactant.concentration
+                    reactant_concentrations[i] = get_concentration(reactant)
         else:
-            reactant_concentrations = [reactant.concentration for reactant in reactants]
+            reactant_concentrations = [
+                get_concentration(reactant) for reactant in reactants
+            ]
 
         # if product_concentrations is provided, use these concentrations
         if product_concentrations is not None:
             assert len(product_concentrations) == len(products)
             for i, product in enumerate(products):
                 if product_concentrations[i] is None:
-                    product_concentrations[i] = product.concentration
+                    product_concentrations[i] = get_concentration(product)
         else:
-            product_concentrations = [product.concentration for product in products]
+            product_concentrations = [
+                get_concentration(product) for product in products
+            ]
 
         # multiply all concentrations to be used in the term
         product_of_reactants = reactant_concentrations[0]

src/festim/species.py

@@ -111,6 +111,12 @@ def __str__(self) -> str:
     def concentration(self):
         return self.solution
 
+    def concentration_submesh(self, subdomain: _VolumeSubdomain):
+        assert subdomain in self.subdomains, (
+            f"Species {self.name} has no solution on subdomain {subdomain}."
+        )
+        return self.subdomain_to_solution[subdomain]
+
     @property
     def legacy(self) -> bool:
         """
@@ -174,6 +180,17 @@ def concentration(self):
                     )
         return self.value_fenics - sum([other.solution for other in self.others])
 
+    def concentration_submesh(self, subdomain: _VolumeSubdomain):
+        if len(self.others) > 0:
+            for other in self.others:
+                assert other.subdomain_to_solution[subdomain], (
+                    f"Cannot compute concentration of {self.name} because {other.name}"
+                    + f" has no solution on subdomain {subdomain}."
+                )
+        return self.value_fenics - sum(
+            [other.subdomain_to_solution[subdomain] for other in self.others]
+        )
+
     def create_value_fenics(self, mesh, t: fem.Constant):
         """Creates the value of the density as a fenics object and sets it to
         self.value_fenics.