I am interested in solving the following equations:
-\nabla \times \boldsymbol{\eta} = \boldsymbol{\gamma} \quad \text{in} \; \Omega\\
\nabla \cdot \boldsymbol{\eta} = \mathbf{0} \quad \text{in} \;\Omega \\
\boldsymbol{\eta} \cdot \hat{\boldsymbol{\nu}} = \mathbf{0} \quad \text{on} \; \partial\Omega
where \boldsymbol{\eta} and \boldsymbol{\gamma} are second-order tensors (\boldsymbol{\gamma} is discontinuous) and \hat{\boldsymbol{\nu}} is the outward unit surface normal to \partial\Omega. In this case, \boldsymbol{\gamma} is given data, and I’m solving for \boldsymbol{\eta}.
The approach I’ve been using so far is the Least-Squares FEM, where I enforce the boundary condition weakly using a penalty term, like so:
\int_{\Omega} (\boldsymbol{\nabla} \times \boldsymbol{\eta}) : (\boldsymbol{\nabla} \times \boldsymbol{\delta\eta}) \, dV
+ \int_{\Omega} (\boldsymbol{\nabla} \cdot \boldsymbol{\eta}) \cdot (\boldsymbol{\nabla} \cdot \boldsymbol{\delta\eta}) \, dV \\
+ \frac{C_\eta}{h}\int_{\partial\Omega} (\boldsymbol{\eta} \cdot \hat{\boldsymbol{\nu}}) \cdot (\boldsymbol{\delta\eta} \cdot \hat{\boldsymbol{\nu}}) \, dS = -\int_{\Omega} \boldsymbol{\gamma} : (\boldsymbol{\nabla} \times \boldsymbol{\delta\eta}) \, dV
where C_\eta > 0 is the penalty factor (I set it to 100), and h = 2 * mesh.rmin()is related to the size of the smallest element.
Here is a MWE of how I have implemented it on FEniCS:
import dolfin as dlf
import ufl
# Compute curl of a tensor
def tcurl(A):
perm = ufl.PermutationSymbol(3)
i, j, k, l = ufl.indices(4)
return ufl.as_tensor(A[i, l].dx(k) * perm[j, k, l], (i,j,))
mesh = dlf.UnitCubeMesh(20, 20, 1)
h = 2 * mesh.rmin()
V_DG = dlf.TensorFunctionSpace(mesh, "DG", 1)
V_CG = dlf.TensorFunctionSpace(mesh, "CG", 1)
gamma_exp = dlf.Expression(
"fabs(x[0] - 0.5) <= 0.05 & fabs(x[1] - 0.5) <= 0.05 ? 1 : 0",
element=V_DG.ufl_element()
)
gamma_init = dlf.Expression(
[
['0', '0', 'val'],
['0', '0', '0'],
['0', '0', '0']
],
val=gamma_exp, element=V_DG.ufl_element()
)
gamma = dlf.project(gamma_init, V_DG)
eta_h = dlf.TrialFunction(V_CG)
deta = dlf.TestFunction(V_CG)
eta = dlf.Function(V_CG, name="eta")
dx_ = dlf.Measure("dx", domain=mesh)
ds_ = dlf.Measure("ds", domain=mesh)
nF = dlf.FacetNormal(mesh)
a = (
ufl.inner(tcurl(eta_h), tcurl(deta)) * dx_
+ ufl.dot(ufl.div(eta_h), ufl.div(deta)) * dx_
+ 100./h * ufl.dot(ufl.dot(eta_h, nF), ufl.dot(deta, nF)) * ds_ # The boundary condition eta \cdot n = 0 is enforced through the surface integral
)
L = ufl.inner(-gamma, tcurl(deta)) * dx_
solver = dlf.PETScKrylovSolver("bicgstab", "sor")
solver.parameters["monitor_convergence"] = False
solver.parameters["report"] = False
solver.parameters["absolute_tolerance"] = 1e-10
solver.parameters["relative_tolerance"] = 1e-6
# solver.ksp().setConvergenceHistory()
A = dlf.PETScMatrix()
b = dlf.PETScVector()
# dlf.assemble(_a, tensor=_A)
dlf.assemble_system(
a, L, bcs=[], A_tensor=A, b_tensor=b
)
solver.set_operator(A)
solver.solve(eta.vector(), b)
with dlf.XDMFFile("testing.xdmf") as outfile:
outfile.parameters["flush_output"] = True
outfile.parameters["functions_share_mesh"] = True
outfile.write(eta, 0.)
outfile.write(gamma, 0.)
I have some questions concerning this approach:
- Is this a consistent way of enforcing the boundary condition \boldsymbol{\eta} \cdot \hat{\boldsymbol{\nu}} = \mathbf{0}?
- Would it be better to use a different function space for \boldsymbol{\eta}? Which one and why?
