# Solving Vectorial Poisson Equation in FENICS

I am trying to solve the following, "test problem" involving a vectorial Poisson equation:

$$-\nabla^2 \vec{A}=\vec{J} \quad \forall x\in\Omega=[-1,1]^3$$ $$\vec{A}=\vec{0} \quad \forall x\in\partial\Omega$$

The inhomogenity $$\vec{J}$$ is given as:

$$\vec{J}(x)=\vec{e_{x}} \quad \mathrm{if} \quad \mathrm{max}(\mathrm{abs}(x_1,x_2,x_3))<0.1$$

E.g. $$\vec{J}$$ is a vector pointing in the x-Direction in a cube of length 0.2 units centered around the origin.

I tried to solve this problem using the Python code listed below. Unfortunately, the code always produces an error message (also listed below), saying that $$a(\vec{A},\vec{v})$$, defined as a=inner(grad(A),grad(v))*dx, is not a bilinear form. Does anyone have an idea on how I could resolve this issue?

Thanks

Python-Code:

from fenics import *
from mshr import *

#define the domain
domain=Box(Point(-1,-1,-1),Point(1,1,1))
mesh = generate_mesh(domain, 32)
V=VectorFunctionSpace(mesh,"P",1)
#plot(mesh)

#define the boundary conditions
def boundary(x,on_boundary):
return on_boundary
A_D=Constant((0,0,0))
bc=DirichletBC(V,A_D,boundary)

#define the inhomogenity J
codex="x>-ext&&x<ext&&x>-ext&&x<ext&&x>-ext&&x<ext?1:0"
codex=codex.replace("ext","0.1")
codey=codez="0"

#define variational problem
A=TrialFunction(V)
v=TestFunction(V)
J=Expression((codex,codey,codez),degree=1)
L=dot(J,v)*dx

#solve the system
A=Function(V)
solve(L==a,A,bc)
plot(A)


Error Message:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to define linear variational problem a(u, v) == L(v) for all v.
*** Reason:  Expecting the left-hand side to be a bilinear form (not rank 1).
*** Where:   This error was encountered inside LinearVariationalProblem.cpp.
*** Process: 0
***
*** DOLFIN version: 2018.1.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------


The problem is in the line solve(L==a,A,bc) which needs to be replaced by solve(a==L,A,bc). The two versions seem to non interchangeable. Doing so will result in the magnetic field B=project(curl(A),V) which is shown in the picture below: 