Why curl-curl with imposed current not working (only small difference to original example)

Question

Hello Fenics-pros,

I am trying to run a simple example to simulate magnetic field around a current.
The current is supposedly running on a straight line (say, a piece of wire) in the middle of the domain. One would expect that the magnetic field would form around the wire solenoidally.
I used the curl-curl-example here

When I try and put my own current into the variable dbdt (actually i am imposing a current, but the equation to solve should stay the same) that is the current in the wire, I get total crap as solution. How could I fix that?

For comparison in the line dbdt delete "my_j" and uncomment what is there behind my_j.

from dolfin import *
from numpy import zeros

mesh=UnitCubeMesh(16,16,16)

class MyJ(Expression):
  def eval(self, values, x):
    if x[0]<= 0.6 and x[0]>=0.4 and x[1]>=0.1 and x[1]<=0.9 and x[2]<= 0.6 and x[2]>=0.4 :
           values[0] = 0.0
           values[1] = 1.0
           values[2] = 0.0
    else :
       values[0]=0.0
       values[1]=0
       values[2] = 0.0

  def value_shape(self):
    return (3, )


# Define function spaces
P1 = VectorFunctionSpace(mesh, "CG", 1)
PN = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 1)

my_j = interpolate(MyJ(), P1)
plot(my_j)
# Define test and trial functions
v0 = TestFunction(PN)
u0 = TrialFunction(PN)
v1 = TestFunction(P1)
u1 = TrialFunction(P1)

# Define functions
dbdt = my_j#Expression(("0.0", "0.0", "1.0"), degree=1)
zero = Expression(("0.0", "0.0", "0.0"), degree=1)
T = Function(PN)
J = Function(P1)

# Dirichlet boundary
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

# Boundary condition
bc = DirichletBC(PN, zero, DirichletBoundary())

# Solve currents equation (using potential T)
solve(inner(curl(u0), curl(v0))*dx == inner(dbdt, v0)*dx, T, bc)

# Solve density equation
solve(inner(v1, u1)*dx == dot(v1, curl(T))*dx, J)

# Plot solution
plot(curl(T))

file=File("current_density.pvd")
file << J

# Hold plot
interactive()

Answer 1

The current density you impose does not fulfill div(J)=0. It 'starts' at
x[1]=0.1 and 'ends' at x[1] = 0.9, which means, that you have a source and
a sink inside your domain.
Thus, neither curl(cur(T))) = J nor its discrete version can hold, because
the right hand side of the equation is not contained in the image of the operator
on the left.
For the linear system Ax=b you get after assembly, b is not contained
in the range of A (which is singular). Obviously, the linear solver is then
producing strange results.
If you change the line

if x[0]<= 0.6 and x[0]>=0.4 and x[1]>=0.1 and x[1]<=0.9 and x[2]<= 0.6 and x[2]>=0.4 :

to

if x[0]<= 0.6 and x[0]>=0.4 and x[2]<= 0.6 and x[2]>=0.4 :

the system will be consistent and the results are ok.
A second approach would be to project your current density into
the space curl(PN) to make it consistent. This is described in some
detail in

Z. Ren, Influence of the RHS on the convergence behaviour of the curl-curl equation,
IEEE Trans. Magn., vol. 32, no. 3, pp. 655–658, May 1996

where it also shown, that CG can be used to solve the consistent
singular linear equations.