# Solve a tensor PDE using the Least-Squares FEM with BC applied weakly

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:

1. Is this a consistent way of enforcing the boundary condition \boldsymbol{\eta} \cdot \hat{\boldsymbol{\nu}} = \mathbf{0}?
2. Would it be better to use a different function space for \boldsymbol{\eta}? Which one and why?

Note that most of the mathematics in this post is not rendering

That’s weird, it seems fine to me:

Hm, Seems to be working now:)