refactoring + documentation

Author: RemDelaporteMathurin

src/festim/boundary_conditions/dirichlet_bc.py

@@ -3,8 +3,10 @@
 import numpy as np
 import numpy.typing as npt
 import ufl
+import ufl.argument
 import ufl.core
 import ufl.core.expr
+import ufl.indexed
 from dolfinx import fem
 from dolfinx import mesh as _mesh
 
@@ -142,6 +144,47 @@ def update(self, t: float):
         elif self.bc_expr is not None:
             self.value_fenics.interpolate(self.bc_expr)
 
+    def weak_formulation(
+        self,
+        u: fem.Function | ufl.indexed.Indexed,
+        v: ufl.argument.Argument | ufl.indexed.Indexed,
+        ds: ufl.Measure,
+    ) -> ufl.core.expr.Expr:
+        """
+        Returns the Nitsche weak formulation for the BC
+        This follows the dolfinx tutorial
+        https://jsdokken.com/dolfinx-tutorial/chapter1/nitsche.html
+
+        Args:
+            u: the solution function associated to the species for which the BC
+                is applied
+            v: the test function
+            ds: the surface measure
+
+        Returns:
+            the weak formulation
+        """
+        mesh = ds.ufl_domain()
+        n = ufl.FacetNormal(mesh)
+        h = ufl.Circumradius(mesh)  # FIXME this doesn't work for rectangles
+        alpha = self.penalty
+        assert alpha is not None, (
+            "Penalty parameter must be given for weakly enforced Dirichlet BCs"
+        )
+        assert self.value_fenics is not None, (
+            "value_fenics must be defined for weakly enforced Dirichlet BCs"
+        )
+
+        form = -ufl.inner(n, ufl.grad(u)) * v * ds(self.subdomain.id)
+
+        form += (
+            +ufl.inner(n, ufl.grad(v)) * (u - self.value_fenics) * ds(self.subdomain.id)
+        )
+        form += (
+            -alpha / h * ufl.inner((u - self.value_fenics), v) * ds(self.subdomain.id)
+        )
+        return form
+
 
 class FixedConcentrationBC(DirichletBCBase):
     """

src/festim/hydrogen_transport_problem.py

@@ -898,34 +898,7 @@ def create_formulation(self):
                 if bc.enforce_weakly:
                     u = bc.species.solution
                     v = bc.species.test_function
-                    n = ufl.FacetNormal(self.mesh.mesh)
-                    h = ufl.Circumradius(
-                        self.mesh.mesh
-                    )  # FIXME this doesn't work for rectangles
-                    alpha = bc.penalty
-                    assert alpha is not None, (
-                        "Penalty parameter must be given for weakly enforced Dirichlet BCs"
-                    )
-                    assert bc.value_fenics is not None, (
-                        "value_fenics must be defined for weakly enforced Dirichlet BCs"
-                    )
-
-                    self.formulation += (
-                        -ufl.inner(n, ufl.grad(u)) * v * self.ds(bc.subdomain.id)
-                    )
-
-                    self.formulation += (
-                        +ufl.inner(n, ufl.grad(v))
-                        * (u - bc.value_fenics)
-                        * self.ds(bc.subdomain.id)
-                    )
-                    self.formulation += (
-                        -alpha
-                        / h
-                        * ufl.inner((u - bc.value_fenics), v)
-                        * self.ds(bc.subdomain.id)
-                    )
-
+                    self.formulation += bc.weak_formulation(u, v, self.ds)
         for adv_term in self.advection_terms:
             # create vector functionspace based on the elements in the mesh
 
@@ -1585,31 +1558,8 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain):
                 if bc.enforce_weakly:
                     u = bc.species.subdomain_to_solution[subdomain]
                     v = bc.species.subdomain_to_test_function[subdomain]
-                    n = ufl.FacetNormal(self.mesh.mesh)
-                    h = ufl.Circumradius(
-                        self.mesh.mesh
-                    )  # FIXME this doesn't work for rectangles
-                    alpha = bc.penalty
-                    assert alpha is not None, (
-                        "Penalty parameter must be given for weakly enforced Dirichlet BCs"
-                    )
-                    assert bc.value_fenics is not None, (
-                        "value_fenics must be defined for weakly enforced Dirichlet BCs"
-                    )
-
-                    form += -ufl.inner(n, ufl.grad(u)) * v * self.ds(bc.subdomain.id)
+                    form += bc.weak_formulation(u, v, self.ds)
 
-                    form += (
-                        +ufl.inner(n, ufl.grad(v))
-                        * (u - bc.value_fenics)
-                        * self.ds(bc.subdomain.id)
-                    )
-                    form += (
-                        -alpha
-                        / h
-                        * ufl.inner((u - bc.value_fenics), v)
-                        * self.ds(bc.subdomain.id)
-                    )
         # add volumetric sources
         for source in self.sources:
             v = source.species.subdomain_to_test_function[subdomain]