Nedelec Curl-Conforming Elements Do Not Converge for a Simple Test Problem

Question

The following minimum example tries to solve the problem

curl(curl(A)) = f
A subject to boundary conditions

using curl-conforming Nedelec elements.
For test purposes I set the vector field A to be constant, yielding f = 0.

from fenics import *


def boundary(_, on_boundary):
    return on_boundary

mesh = UnitSquareMesh(3, 3)
V = FunctionSpace(mesh, 'Nedelec 1st kind H(curl)', 1)

A = TrialFunction(V)
N = TestFunction(V)

A_D = Constant((0, 1))
bc = DirichletBC(V, A_D, boundary)
f = Constant((0, 0))

a = inner(curl(A), curl(N)) * dx
b = inner(f, N) * dx

A_ = Function(V)
solve(a == b, A_, bc)

print errornorm(A_D, A_, 'Hcurl')
plot(A_, interactive=True)

The result however is entirely wrong.
The norm of the error is between one and two, depending on whether I run FEniCS from a docker image or a self-compiled version, but it is nowhere near zero.
Plotting shows that the result is totally wrong.
The vectors point in all possible directions and their norm varies between zero and hundred on a mesh with 10x10 triangles (the error is then 24.3).

Note, that the test function I use for A is just a constant and thus the error should be zero with respect to machine precision, since it is contained in the space of basis functions.
Surprisingly, the error drops to magnitude 1e-16, when I use a mesh with only 1x1 or 2x2 triangles.
Increasing the number of triangles destroys everything.
Increasing the order of the function space increases the error even further.

What is wrong here?

Replacing the line

V = FunctionSpace(mesh, 'Nedelec 1st kind H(curl)', 1)

with

V = VectorFunctionSpace(mesh, 'Lagrange', 1)

solves the problem exactly, but is not what we want / need.

print dolfin.__version__
2016.1.0

Thanks in advance.

Answer 1

The code is doing the correct thing, but the mathematics is wrong :)

The problem

$$ curl(curl(A)) = f, $$

has a large nullspace. More specifically, if $A$ is a solution also $A + \nabla \phi$ is a solution for any $\phi \in H^1$.

Cosider solving the following saddle point system where you constain $A$ to be L^2 orthogonal to the space $\nabla H^1$, i.e.

Find $(A, phi) \in H(curl) \times H^1$
$(curl A, curl N) + (\nabla \phi, N) = (f, N) \forall N \in H(curl)$
$(A, \nabla \psi) = 0 \forall \psi \in H^1$

Best,

U.

Comments

If I'm not mistaken, this is equivalent to setting the gauge, right?
In this case, it's the Coulomb gauge, leading to $\nabla\cdot A=0$.

When testing with some more interesting functions, I can see, that the curl of the solution is actually correct.

I always thought, the Coulomb gauge condition is already intrinsically embedded in the Nedelec elements, where the ansatz functions have zero divergence per construction?

When I plot the divergence of the solution of the code posted above, I can see that it is nearly zero everywhere.
Correct me, if I'm wrong, but there should be only one unique vector field that fulfills
$ \nabla\times\nabla\times A = f $
$ \nabla\cdot A = 0 $
and the boundary conditions at the same time.
In my example above, this is the vector field $A^\ast = (0, 1)$, which is divergence free.
So, I don't understand, how the code above can produce a (nearly) divergence free solution, which differs from $A^\ast$ by several orders of magnitude.

When I incorporate your suggestions, the result gets even worse. The error is of order 1e15 then.
Here is my code:

from fenics import *


def boundary(_, on_boundary):
    return on_boundary

mesh = UnitSquareMesh(3, 3)
NE = FiniteElement('N1curl', triangle, 1)
LA = FiniteElement('Lagrange', triangle, 1)
V = FunctionSpace(mesh, NE * LA)

(A, phi) = TrialFunctions(V)
(N, psi) = TestFunctions(V)

A_D = Constant((0, 1))
bc = DirichletBC(V.sub(0), A_D, boundary)
f = Constant((0, 0))

a = inner(curl(A), curl(N)) * dx + inner(grad(phi), N) * dx \
  + inner(A, grad(psi)) * dx
b = inner(f, N) * dx

w = Function(V)
solve(a == b, w, bc)
(A_, _) = w.split()

print errornorm(A_D, A_, 'Hcurl')
plot(A_, title='A')
interactive()

Moreover, when incorporating your suggestions, the curl of the solution does not converge anymore, as it used to do before.
Here is a more complicated example:

from fenics import *


def boundary(_, on_boundary):
    return on_boundary

mesh = UnitSquareMesh(3, 3)
NE = FiniteElement('N1curl', triangle, 1)
LA = FiniteElement('Lagrange', triangle, 1)
V = FunctionSpace(mesh, NE * LA)

(A, phi) = TrialFunctions(V)
(N, psi) = TestFunctions(V)

A_D = Expression(('-x[1]*x[1]', 'x[0]*x[0]'))
B_D = Expression('2*x[0] + 2*x[1]')
bc = DirichletBC(V.sub(0), A_D, boundary)
f = Constant((2, -2))

a = inner(curl(A), curl(N)) * dx + inner(grad(phi), N) * dx \
  + inner(A, grad(psi)) * dx
b = inner(f, N) * dx

w = Function(V)
solve(a == b, w, bc)
(A_, _) = w.split()
B = project(curl(A_), FunctionSpace(mesh, 'Lagrange', 1))

print errornorm(A_D, A_, 'Hcurl')
print errornorm(B_D, B, 'L2')
plot(A_, title='A')
plot(B, title='B = curl(A)')

interactive()

It worked before for $B$, but not for $A$.
Now it works for none of the two.

Is there anywhere an example, that shows how to properly solve the curl-curl problem?
The examples I found so far immediately take the curl and don't care for setting the gauge.
However, I really need the vector potential.

Regards

---

I think I got it now.
When I add homogeneous boundary conditions for the second subspace, everything works as expected.

Is it correct, that I have to add boundary conditions for the second subspace?

All the literature I found characterizes the saddle point problem as: Find $(A, \phi)\in (H(\operatorname{curl})\times H^1 / \mathbb{R})$ such that
$\int (\nabla\times A)\cdot(\nabla\times v)\,\mathrm{d}\Omega + \int \nabla\phi\cdot v\,\mathrm{d}\Omega = \int f\cdot v\,\mathrm{d}\Omega\qquad \forall v\in H(\operatorname{curl})$

$\int A\cdot \nabla\psi\,\mathrm{d}\Omega = 0 \qquad \forall \psi\in H^1 / \mathbb{R}$

I have no idea, what the space $H^1 / \mathbb{R}$ actually is, but it seems to me, it is $H^1$ without constant valued functions, which can be modeled as $H^1_0$, is that correct?

Before switching to FEniCS, I had a Diffpack implementation that solved the same problem, but without any problems and without specifying $\nabla\cdot A = 0$ in any way.
I wonder why that worked...

---

Hi,

I am glad you figured out the issues of the boundary conditions. That was a very important point I overlooked in my previous answer.

Just to clarify the definition of the functional spaces:

$$ H^1/\mathbb{R} := { v \in H^1(\Omega) \, | \, \int_{\Omega} v d\Omega = 0 }$$
and
$$ H^1_0 := { v \in H^1(\Omega) \, | \, v = 0 \quad \text{on} \, \partial \Omega }$$

Regarding your Diffpack implementation, it is true that one can solve directly the curl-curl problem. However one needs to be extremely careful.

In fact, you can solve directly the problem
$$ curl \, curl A = f $$
using an iterative method (such as conjugate gradient) provided that
1) $f \perp \nabla \phi$ for all $\phi \in H^1$ (or $\phi \in H^1_0$ depending on your bc) .
2) You are using a preconditioner that does not see the singularity of the operator (i.e. Gauss-Seidel or Jacobi will work fine but slowly, incomplete factorization will most likely fail.)

Finally, in HYPRE there is a specific solver called HYPRE_AMS to solve the singular curl-curl problem.

see:
http://www.jstor.org/stable/43693530
and
http://epubs.siam.org/doi/ref/10.1137/110859361

Unfortunately, I don't think that FEniCS is able to support HYPRE_AMS yet.

MFEM (http://mfem.org/) is a finite element library that does support HYPRE_AMS, but requires you to get your hands dirty with c++ and the finite element assembly.

Best,

U.

---

Is it correct, that I have to add boundary conditions for the second
subspace?

Yes, the orthogonal decomposition of $H_0(curl)$ is given for example in

Toselli-Widlund, Domain Decomposition Methods – Algorithms and Theory
page 276, section 10.1.1 equation 10.7

And it uses the Grad H^1_0 (i.e. the space with dirichlet bc)

Best,

U.

---

Could you please tell me how to impliment the boundary condition, which if the problem's boundary condition is A in H0(curl)?

I can do Dirichlet boundary condition, Neumann boundary condition, Robin boundary condition, but not H0(curl).

I will appreciate your any help. Thanks!!!

---

Hi,

I simply used Dirichlet boundary conditions for the second subspace and thus used $H^1_0$ rather than $H^1 / \mathbb{R}$.

The details are given in "Toselli-Widlund, Domain Decomposition Methods – Algorithms and Theory", as said above by umberto.

Hope this helps you.