# How to constraint the tangential gradient on a boundary in FEniCS?

The problem I'm considering is a 2D scalar PDE. The domain $$\Omega$$ is a disk with two holes $$\partial\Omega_1$$ and $$\partial\Omega_2$$ and an external boundary $$\partial\Omega_0$$. The PDE and boundary conditions are:

$$\left(\nabla^2-1\right)\eta=0$$ on $$\Omega$$

$$\nabla \eta\cdot\mathbf{n}=a_i$$ on $$\partial\Omega_i, {i\in\{1,2\}}$$

$$\nabla \eta\cdot\mathbf{n}=-\eta$$ on $$\partial\Omega_0$$

The weak form is

$$\int_\Omega \nabla \eta \cdot \nabla \varphi d\Omega + \int_\Omega \eta\varphi d\Omega + \int_{\partial\Omega_0}\eta\varphi d\Gamma = \sum_{i\in\{1,2\}} a_i \int_{\partial\Omega_i} \varphi d\Gamma$$

I need to add the following condition to the model:

$$\nabla \eta \cdot \mathbf{t} = 0$$ on $$\partial\Omega_i, {i\in\{1,2\}}$$

which is equivalent to

$$\exists b_i, \eta = b_i$$ on $$\partial\Omega_i, {i\in\{1,2\}}$$

How can I add this constraint in FEniCSx?

• I would ask this question on the FEniCS forum (fenicsproject.discourse.group) for FEniCSx implementational questions Commented May 24 at 18:13
• This is the electrode model in electrical impedance tomography. You can probably find more details about implementations (in FEniCS or other systems) if you search for that term. Commented May 25 at 2:20
• I managed to solve my problem, here is the FEniCS forum thread fenicsproject.discourse.group/t/… Commented Jun 3 at 13:17

Not really a solution for FeniCS, but this is how I would solve this problem in GetFEM

import getfem as gf
import numpy as np

# Parameters for the mesh
h = 0.05  # Mesh size

# Define the mesher objects for the circular domain and the holes
mo_outer = gf.MesherObject('ball', [0.0, 0.0], radius_outer)
mo_hole1 = gf.MesherObject('ball', [-0.5, 0.0], radius_inner)
mo_hole2 = gf.MesherObject('ball', [0.5, 0.0], radius_inner)
mo_holes = gf.MesherObject('union', mo_hole1, mo_hole2)
mo_domain = gf.MesherObject('set minus', mo_outer, mo_holes)

# Generate the mesh
mesh = gf.Mesh('generate', mo_domain, h, 2)
mesh.export_to_vtu("mesh.vtu")

# Identify the boundary regions
outer_boundary = mesh.outer_faces()
mesh.set_region(99, outer_boundary)

hole1_faces = mesh.outer_faces_in_ball([-0.5, 0.0], radius_inner + h)
hole2_faces = mesh.outer_faces_in_ball([0.5, 0.0], radius_inner + h)
mesh.set_region(101, hole1_faces)
mesh.set_region(102, hole2_faces)
mesh.region_subtract(100, 101)
mesh.region_subtract(100, 102)

# Define the finite element method for the solution (u) and test function (Test_u)
mf = gf.MeshFem(mesh, 1)
mf.set_fem(gf.Fem('FEM_PK(2,1)'))
# reduce the finite element space to constant values on boundaries 101 and 102
kept_dofs = list(set(range(mf.nb_basic_dof()))-set(mf.basic_dof_on_region(101)[1:])
-set(mf.basic_dof_on_region(102)[1:]))
fullsize = mf.nb_basic_dof()
newdofs = np.nan * np.ones(fullsize, np.int_)
newsize = 0
for dof in kept_dofs:
newdofs[dof] = newsize
newsize += 1
RR = gf.Spmat("empty", newsize, fullsize)
EE = gf.Spmat("empty", fullsize, newsize)
for dof in kept_dofs:
RR.assign(newdofs[dof], dof, 1.)
EE.assign(dof, newdofs[dof], 1.)
master_dof = newdofs[mf.basic_dof_on_region(101)[0]]
for dof in mf.basic_dof_on_region(101)[1:]:
EE.assign(dof, master_dof, 1.)
master_dof = newdofs[mf.basic_dof_on_region(102)[0]]
for dof in mf.basic_dof_on_region(102)[1:]:
EE.assign(dof, master_dof, 1.)
mf.reduction_matrices(RR, EE)

# Integration methods
mim = gf.MeshIm(mesh, gf.Integ('IM_TRIANGLE(6)'))

# Model
md = gf.Model('real')

# PDE terms

# Boundary conditions
a1 = 1.0  # Example value, set as needed
a2 = -1.0  # Example value, set as needed

# Normal derivative conditions on the holes (Neumann)

# Solve the model
md.solve('noisy')

# Extract the solution
U = md.variable('u')

# Visualization
mf.export_to_vtu('solution.vtu', U, 'Solution')

# Post-processing or further analysis as needed
print("Solution computed and exported to solution.vtu")


(code based on ChatGPT-4o suggestion, heavily adapted)