Hello everyone, I was able to solve the problem by myself. The code is not the cleanest and has a lot of things that could be polished but I’ll post it in case anyone encounters a similar problem:

```
import dolfin as df
import numpy as np
import matplotlib.pyplot as plt
#Mesh
mesh = df.RectangleMesh(df.Point(-10,-10), df.Point(10,10), 50, 50, 'left/right')
#Function Spaces for vector potential and B-Field
V = df.FunctionSpace(mesh, 'CG', 2)
Vdim = df.VectorFunctionSpace(mesh, 'CG', 2)
#Dirichlet BC in all boundaries.
def boundary(x, on_boundary):
return on_boundary
bc = df.DirichletBC(V, df.Constant(0), boundary)
#Get node coordinates
bmesh = df.BoundaryMesh(mesh, "exterior", True)
b_coord = bmesh.coordinates()
b_coord_x = np.sort(np.unique(b_coord[:,0]))
b_coord_z = np.sort(np.unique(b_coord[:,1]))
N_layers = 3
#Define meshfunction and set all markers to 0
mshfunc = df.MeshFunction("size_t", mesh, mesh.topology().dim())
mshfunc.set_all(0)
#Overwrite markers of last N_layers cells
for i in reversed(range(N_layers+1)):
class Layer(df.SubDomain):
def inside(self, x, on_boundary):
return (df.between(x[0], (b_coord_z[-i-1], b_coord_z[-i])) or df.between(x[1], (b_coord_x[-i-1], b_coord_x[-i]))) or (df.between(x[0], (b_coord_z[i-1], b_coord_z[i])) or df.between(x[1], (b_coord_x[i-1], b_coord_x[i])))
Layer().mark(mshfunc, N_layers+1-i)
# Define integration measure
dx = df.Measure('dx', domain=mesh, subdomain_data=mshfunc)
#Assign permeability to each domain
class Permeability(df.UserExpression):
def __init__(self, markers, **kwargs):
super(Permeability, self).__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
for i in range(N_layers+1):
if self.markers[cell.index] == i:
values[0] = (i+1)**2
def value_shape(self):
return ()
mu = df.interpolate(Permeability(mshfunc), V)
#Currents
class Currents(df.UserExpression):
def eval(self, value, x):
if np.sqrt(x[0]**2 + (x[1]-1)**2)< 1:
value[0] = 1
elif np.sqrt(x[0]**2 + (x[1]+1)**2)< 1:
value[0] = -1
else:
value[0] = 0
def value_shape(self):
return ()
J = df.interpolate(Currents(), V)
#Weak form
A_y = df.TrialFunction(V)
v = df.TestFunction(V)
a = (1/mu)*df.dot(df.grad(A_y), df.grad(v))*dx
L = J*v*dx
A_y = df.Function(V)
df.solve(a == L, A_y, bc)
#############################
#Postprocessing
x = np.linspace(-10,10, 100)
z = np.linspace(-10,10,100)
X,Z = np.meshgrid(z,x)
psi_grid = [A_y(j,i) for i in x for j in z]
psi_grid = np.reshape(psi_grid, np.shape(X))
mu_grid = [mu(j,i) for i in x for j in z]
mu_grid = np.reshape(mu_grid, np.shape(X))
fig = plt.figure()
ax = fig.add_axes([0,0,1,1])
pcol = ax.pcolor(X, Z, mu_grid)
plt.colorbar(pcol)
cont = ax.contour(X, Z, psi_grid, np.linspace(-0.3, 0.3, 20), colors='cyan', linestyles='solid', linewidths=0.8)
```

And the outcome is: