Altering vectorfunction values

Question

Hello FEniCS experts!
I am very new to using fenics and actually also python (started using both for my master thesis because it was recommended to me)

I am in a situation where I need to solve some magnetostatic differential equations in 2D due to some domain assumptions I have a differential equation as following:

B = magnetic density field, J = current density field.
[dB_z/dx,dB_z/dy] = mu_0*[-Jy,Jx]

solving with my boundary conditions I can attain the potential and the current density J which is a function in the vector function space V_g = VectorFunctionSpace(mesh,'cg',1).

Now, is there an easy way to create a new vector function in V_g given by [-Jy,Jx].

I have tried something like: Jhat = Function(V_g)
Jhat.vector()[0->len/2] = -J_y.vector()[:] (from J_x,J_y = J.split(deepcopy=True))
Jhat.vector()[len/2->end] = J_x.vector()[:]

However, I am not sure if this is correct. I think I remember reading somewhere that x,y,z values of a vectorfunction are given after each other which then prompted my above solution...

So the question is two-part: 1. is my attempt correct, 2. if not - how can I do this? :)
edit:
I also tried
class Jnew(Expression):
def init(self,Jx,Jy):
self.Jx = Jx
self.Jy = Jy
def eval(self,value,x):
value[0] = (-self.Jy(x[0],x[1]),self.Jx(x[0],x[1]))
def value_shape(self):
return (2,1)

But this creates an rank-error when I try Jhat = interpolate(Jnew(Jx=J_x,Jy=J_y),V_g)
saying that my function is rank 2 and my function space is rank 1... which I somewhat understand...
best regards,
Rho

Answer 1

Hi, the first solution will not work because for function from VectorFunctionSpace entries in
its vector are aligned something like [0_x, 0_y, 1_x, 1_y, ...] and not [all_x_comps, all_y_comps] as you assume. In the second solution you need to return (2, ).

from dolfin import *
from numpy import zeros

mesh = UnitSquareMesh(10, 10)
V = VectorFunctionSpace(mesh, "CG", 1)

J = interpolate(Expression(('x[0]', 'x[1]')), V)
Jx, Jy = J.split(deepcopy=True)

class MyJ(Expression):
  def eval(self, values, x):
    values[0] = -Jy(x)
    values[1] = Jx(x)

  def value_shape(self):
    return (2, )

my_j = interpolate(MyJ(), V)
plot(my_j, interactive=True)

Jx_v = Jx.vector().array()
Jy_v = Jy.vector().array()
dim = V.dim()
my_j_v = zeros(dim)

for i in range(dim/2):
  my_j_v[2*i] = -Jy_v[i]
  my_j_v[2*i+1] = Jx_v[i]

my_j.vector()[:] = my_j_v
plot(my_j, interactive=True)

Comments

Thank you very much!