Weakly enforced Dirichlet BCs

.github/workflows/ci_docker.yml

@@ -20,7 +20,8 @@ jobs:
 
       - name: Install local package and dependencies
         run: |
-          pip install .[test]
+          python -m pip install scipy  # needed in scifem, can be removed when scifem adds it to their dependencies
+          python -m pip install .[test]
 
       - name: Overload adios4dolfinx
         if: ${{ matrix.container_version == 'nightly' }}

pyproject.toml

@@ -1,5 +1,9 @@
 [build-system]
-requires = ["setuptools >= 61", "wheel", "setuptools-scm[toml] >= 7.0.5"]
+requires = [
+    "setuptools>=77.0.3", # Ensuring license field is correctly parsed
+    "wheel",
+    "setuptools-scm[toml] >= 7.0.5",
+]
 build-backend = "setuptools.build_meta"
 
 [project]

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
 
@@ -19,13 +21,20 @@ class DirichletBCBase:
     Args:
         subdomain: The surface subdomain where the boundary condition is applied
         value: The value of the boundary condition
+        enforce_weakly: Whether to enforce the boundary condition weakly using Nitsche's
+            method. Defaults to False.
+        penalty: The penalty parameter to use if ``enforce_weakly`` is True.
+            Defaults to None
 
     Attributes:
         subdomain: The surface subdomain where the boundary condition is applied
         value: The value of the boundary condition
         value_fenics: The value of the boundary condition in fenics format
         bc_expr: The expression of the boundary condition that is used to
             update the `value_fenics`
+        enforce_weakly: Whether to enforce the boundary condition weakly using Nitsche's
+            method.
+        penalty: The penalty parameter to use if ``enforce_weakly`` is True.
     """
 
     subdomain: _subdomain.SurfaceSubdomain
@@ -45,9 +54,13 @@ def __init__(
         self,
         subdomain: _subdomain.SurfaceSubdomain,
         value: np.ndarray | fem.Constant | int | float | Callable,
+        enforce_weakly: bool = False,
+        penalty: float | None = None,
     ):
         self.subdomain = subdomain
         self.value = value
+        self.enforce_weakly = enforce_weakly
+        self.penalty = penalty
 
         self.value_fenics = None
         self.bc_expr = None
@@ -138,6 +151,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):
     """
@@ -147,6 +201,10 @@ class FixedConcentrationBC(DirichletBCBase):
         value: The value of the boundary condition. It can be a function of
         space and/or time
         species: The name of the species
+        enforce_weakly: Whether to enforce the boundary condition weakly using Nitsche's
+            method. Defaults to False.
+        penalty: The penalty parameter to use if ``enforce_weakly`` is True.
+            Defaults to None
 
     Attributes:
         temperature_dependent (bool): True if the value of the bc is
@@ -180,9 +238,13 @@ def __init__(
         subdomain: _subdomain.SurfaceSubdomain,
         value: np.ndarray | fem.Constant | int | float | Callable,
         species: Species,
+        enforce_weakly: bool = False,
+        penalty: float | None = None,
     ):
         self.species = species
-        super().__init__(subdomain, value)
+        super().__init__(
+            subdomain, value, enforce_weakly=enforce_weakly, penalty=penalty
+        )
 
     @property
     def temperature_dependent(self):

src/festim/hydrogen_transport_problem.py

@@ -878,7 +878,7 @@ def create_formulation(self):
                 * self.dx(source.volume.id)
             )
 
-        # add fluxes
+        # add boundary conditions (fluxes and weak dirichlet)
         for bc in self.boundary_conditions:
             if isinstance(bc, boundary_conditions.ParticleFluxBC):
                 self.formulation -= (
@@ -894,6 +894,11 @@ def create_formulation(self):
                         * self.ds(flux_bc.subdomain.id)
                     )
 
+            if isinstance(bc, boundary_conditions.FixedConcentrationBC):
+                if bc.enforce_weakly:
+                    u = bc.species.solution
+                    v = bc.species.test_function
+                    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
 
@@ -1549,6 +1554,11 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain):
                 if subdomain == self.surface_to_volume[bc.subdomain]:
                     v = bc.species.subdomain_to_test_function[subdomain]
                     form -= bc.value_fenics * v * self.ds(bc.subdomain.id)
+            if isinstance(bc, boundary_conditions.FixedConcentrationBC):
+                if bc.enforce_weakly:
+                    u = bc.species.subdomain_to_solution[subdomain]
+                    v = bc.species.subdomain_to_test_function[subdomain]
+                    form += bc.weak_formulation(u, v, self.ds)
 
         # add volumetric sources
         for source in self.sources:

src/festim/problem.py

@@ -124,7 +124,8 @@ def define_boundary_conditions(self):
         for bc in self.boundary_conditions:
             if isinstance(bc, F.DirichletBCBase):
                 form = self.create_dirichletbc_form(bc)
-                self.bc_forms.append(form)
+                if not bc.enforce_weakly:
+                    self.bc_forms.append(form)
 
     def get_petsc_options(self) -> dict[str, Any]:
         """Gets the PETSc options to pass to the NewtonProblem solver. Default options

test/system_tests/test_misc.py

@@ -3,10 +3,12 @@
 import dolfinx
 import numpy as np
 import pytest
+import ufl
 
 import festim as F
 
 from .test_multi_mat_penalty import generate_mesh
+from .tools import error_L2
 
 
 def test_petsc_options():
@@ -389,3 +391,90 @@ def test_sub_temperature_as_function_mixed_domain_not_updated():
     my_model.run()
 
     assert not np.allclose(avg_surf.data, 4.0)
+
+
+@pytest.mark.parametrize(
+    "model_class", [F.HydrogenTransportProblem, F.HydrogenTransportProblemDiscontinuous]
+)
+def test_MMS_weak_dirichlet(
+    model_class: type[F.HydrogenTransportProblem]
+    | type[F.HydrogenTransportProblemDiscontinuous],
+):
+    """
+    MMS test with one mobile species at steady state, with Dirichlet BCs enforced
+    weakly.
+    """
+
+    def u_exact(mod):
+        return lambda x: 1 + mod.sin(2 * mod.pi * x[0]) + mod.cos(2 * mod.pi * x[1])
+
+    H_analytical_ufl = u_exact(ufl)
+    H_analytical_np = u_exact(np)
+
+    mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 50, 50)
+
+    V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
+    T = dolfinx.fem.Function(V)
+
+    D_0 = 1
+    E_D = 0.1
+
+    def T_expr(x):
+        return 500 + 100 * x[0]
+
+    T.interpolate(T_expr)
+    D = D_0 * ufl.exp(-E_D / (F.k_B * T))
+
+    my_model = model_class()
+    my_model.mesh = F.Mesh(mesh=mesh)
+
+    my_mat = F.Material(name="mat", D_0=D_0, E_D=E_D)
+    vol = F.VolumeSubdomain(id=1, material=my_mat)
+    left = F.SurfaceSubdomain(id=1, locator=lambda x: np.isclose(x[0], 0))
+    right = F.SurfaceSubdomain(id=2, locator=lambda x: np.isclose(x[0], 1))
+    my_model.subdomains = [vol, left, right]
+
+    H = F.Species("H")
+    H2 = F.Species("D")
+    my_model.species = [H, H2]
+
+    if isinstance(my_model, F.HydrogenTransportProblemDiscontinuous):
+        H.subdomains = [vol]
+        H2.subdomains = [vol]
+
+    my_model.temperature = T_expr
+
+    my_model.boundary_conditions = [
+        F.FixedConcentrationBC(
+            subdomain=left,
+            value=H_analytical_ufl,
+            species=H,
+            enforce_weakly=True,
+            penalty=300,
+        ),
+        F.FixedConcentrationBC(
+            subdomain=right,
+            value=H_analytical_ufl,
+            species=H,
+            enforce_weakly=True,
+            penalty=300,
+        ),
+    ]
+
+    x = ufl.SpatialCoordinate(mesh)
+    f = -ufl.div(D * ufl.grad(H_analytical_ufl(x)))
+    my_model.sources = [F.ParticleSource(value=f, volume=vol, species=H)]
+
+    my_model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)
+
+    my_model.initialise()
+    my_model.run()
+
+    if isinstance(my_model, F.HydrogenTransportProblemDiscontinuous):
+        H_computed = H.subdomain_to_post_processing_solution[vol]
+    else:
+        H_computed = H.post_processing_solution
+
+    L2_error = error_L2(H_computed, H_analytical_np)
+
+    assert L2_error < 2e-3