Simple 3d electrostatic problem

1

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:

image
image766×380 10.6 KB

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:

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

You have not solved the linear variational problem.
I.e. you need to change the last 4 lines to:

uh = Function(V)

solve(a==L, uh, bcs=[bc])
vtkfile = File('solutionNew.pvd')
vtkfile << uh
3

Ah yeah, that was it. I must have lost the line during debugging.Thank you very much!
How about my second question? How do I properly convert to a numpy array that I can easily analyze with matplotlib?

4

You can simply call plot(uh) and plt.show() to visualize your solution with matplotlib (prior to writing the solution to pvd.
Iv you want to read the data in as a vtu file, I suggest using external libraries such as GitHub - marcomusy/vedo: A python module for scientific analysis of 3D data or GitHub - pyvista/pyvista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK) to visualize your data.

Note that matplotlib is not really suitable for 3D analysis.

5

I tried these solutions as well, but ran into some environment-related problems.
I solved my issue with the (admittedly expensive) solution of generating equidistant arrays of the coordinates I am interested in and evaluating my solution at these points, which is easy to do by just adding the following:

N = 10
x = np.linspace(0,1,N)
y = np.linspace(0,1,N)
z = np.linspace(0,1,N)
xg, yg ,zg = np.meshgrid(x, y, z, indexing='ij', sparse=False)
Ex_array = np.zeros((N,N,N))
Ey_array = np.zeros((N,N,N))
Ez_array = np.zeros((N,N,N))

for i in range(N):
    for j in range(N):
        for k in range(N):
            #print((xg[i,j,k], yg[i,j,k], zg[i,j,k]))
            Ex_array[i,j,k] = uh((xg[i,j,k], yg[i,j,k], zg[i,j,k]))[0]
            Ey_array[i,j,k] = uh((xg[i,j,k], yg[i,j,k], zg[i,j,k]))[1]
            Ez_array[i,j,k] = uh((xg[i,j,k], yg[i,j,k], zg[i,j,k]))[2]

plt.imshow(Ez_array[:,:,5])
plt.savefig("Ex.png")