Hello all,

after my initial installation hurdles I managed to use the dockerfile to install fenics on the cluster.

Now, as a first test, I want to calculate a simple electrostatic problem with a cylindrical current source. My code is below:

```
from fenics import *
from ufl import nabla_div
import numpy as np
mesh = UnitCubeMesh(10, 10, 10)
V = VectorFunctionSpace(mesh, family='CG', degree=1, dim=3)
def boundary(x):
tol = 1E-14
return near(x[0], 0, tol) or near(x[1], 0, tol) or near(x[2], 0, tol) \
or near(x[0], 1, tol) or near(x[1], 1, tol) or near(x[2], 1, tol) \
materials = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
r = np.sqrt(x[0]**2+x[1]**2)
return r <= 0.2 + tol
class Omega_1(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
tol = 1E-14
r = np.sqrt(x[0]**2+x[1]**2)
return r >= 0.2 + tol
class F(UserExpression):
def __init__(self, materials, f_0, f_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.f_0 = f_0
self.f_1 = f_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = self.f_0
else:
values[0] = self.f_1
f_value = F(materials, 0, 1, degree=0)
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(materials, 0)
subdomain_1.mark(materials, 1)
omega = 1
bc = DirichletBC(V,[0, 0, 0],boundary)
u = TrialFunction(V)
v = TestFunction(V)
f = as_vector((0,0,omega*f_value))
a = (div(u)*div(v))*dx
L = dot(f,v)*dx
u = Function(V)
vtkfile = File('solutionNew.pvd')
vtkfile << u
```

Aside from one warning (WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.), it works.

The only problem is: When I open the generated file with pyvista via jupyter-notebook, the output looks like this:

Which does not seem right. Sadly, so far, I did not figure out how to plot the solution; But I suspect there is an error before that.

So can you help me answer either (or preferably both) of these questions:

- What am I doing wrong in setting up the problem?
- How can I convert the solution to a numpy array which facillitates correct ploting?