Calculating the E-Field resulting from the poisson equation

dolfinx
1

Hi everyone,

after i solved the poisson equation with a specific charge density, i wanted to calculate and simulate the E-field. i looked over other topic to maybe find a close problem to mine but somehow i’ve found very little about it.
i used the Gmsh Python API to define my mesh, then i followed the Tutorial for DOLFINX to solve the poisson equation and everything worked. I then defined the E function and the E_norm function as followed:

s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, s_cg1)
phi = Function(V)
...
...#Poisson Equation solution
...
Element = VectorElement("CG", mesh.ufl_cell(), 1)
Q = VectorFunctionSpace(mesh, Element)
E = Function(Q)
E_norm= Function(V)

E = - ufl.grad(phi)

my Kernel crashes everytime i run the code with the definition of VectorElement and VectorFunctionSpace. i would be very thankful if you could help me.

another this is that i know using DOLFIN the mathematical equation for E and E_norm are defined as followed:

E.assign(project(-grad(phi), Q))
E_norm.assign(project(sqrt(inner(E,E)), V))

I’m pretty sure that what i define in the DOLFINX code is not correct and will lead to an Error. could you please give me some hints on how to formulate my code right.

I look forward to your reply.

2

As phi is a scalar, grad(phi) is a vector, and you would have:

Element =FiniteElement("CG", mesh.ufl_cell(), 1)
Q = VectorFunctionSpace(mesh, Element)

See:
https://jorgensd.github.io/dolfinx-tutorial/chapter2/linearelasticity_code.html#stress-computation

You should also consider the choice of spaces. If phi is a first order Lagrange space, grad(phi) is a piecewise constant (DG 0), and the norm will also be discontinuous across element boundaries

1 Like
3

Dokken Thank you so much for your fast answer!
It did work as you told me to.
I have another question if you won’t mind, it’s about plotting this Function. I used the plotting method in the tutorial you showed me earlier and even when i install warped i keep getting this message(even after restarting the Kernel):

NameError: name ‘warped’ is not defined

I used the following code to save the solution to file but each time i try to run it in ParaView it crashes:

from dolfinx import io
with io.VTKFile(mesh.comm, "output2.pvd", "w") as vtk:
    vtk.write([E._cpp_object])
with io.XDMFFile(mesh.comm, "output2.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(E)

do you got any idea that could help me maybe?
I truly appreciate your help

4

Without any minimal code example reproducing the behavior, it is rather tricky to help you.

Seems like you have not defined the variable warped.

What errors does the following commands give?

5

as for the visualisation using Pyvista i used the following code:

import pyvista
from dolfinx import mesh, fem, plot, io
pyvista.set_jupyter_backend("pythreejs")

# Create plotter and pyvista grid
p = pyvista.Plotter()
topology, cell_types, geometry = plot.create_vtk_mesh(V)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

# Attach vector values to grid and warp grid by vector
grid["u"] = uh.x.array.real.reshape((geometry.shape[0], 1))
actor_0 = p.add_mesh(grid, style="wireframe", color="k")
warped = grid.warp_by_vector("u", factor=1.5)
actor_1 = p.add_mesh(warped, show_edges=True)
p.show_axes()
if not pyvista.OFF_SCREEN:
   p.show()
else:
   pyvista.start_xvfb()
   figure_as_array = p.screenshot("phi.png")

warped.cell_data["E"] = EE.vector.array.real
warped.set_active_scalars("E")
p = pyvista.Plotter()
p.add_mesh(warped)
p.show_axes()
if not pyvista.OFF_SCREEN:
   p.show()
else:
   pyvista.start_xvfb()
   Efield_figure = p.screenshot(f"EE.png")

i wrote this in 2 different cells and didn’t notice that i got an Error before as followed:

> ---------------------------------------------------------------------------
> ValueError                                Traceback (most recent call last)
> Input In [16], in <cell line: 13>()
>      11 grid["u"] = uh.x.array.real.reshape((geometry.shape[0], 1))
>      12 actor_0 = p.add_mesh(grid, style="wireframe", color="k")
> ---> 13 warped = grid.warp_by_vector("u", factor=1.5)
>      14 actor_1 = p.add_mesh(warped, show_edges=True)
>      15 p.show_axes()
> 
> File /usr/local/lib/python3.9/dist-packages/pyvista/core/filters/data_set.py:2299, in DataSetFilters.warp_by_vector(self, vectors, factor, inplace, progress_bar)
>    2297 # check that this is indeed a vector field
>    2298 if arr.ndim != 2 or arr.shape[1] != 3:
> -> 2299     raise ValueError(
>    2300         'Dataset can only by warped by a 3D vector point data array. '
>    2301         'The values you provided do not satisfy this requirement'
>    2302     )
>    2303 alg = _vtk.vtkWarpVector()
>    2304 alg.SetInputDataObject(self)
> 
> ValueError: Dataset can only by warped by a 3D vector point data array. The values you provided do not satisfy this requirement

i’m more interested in using Paraview for visualizing, therefor i used the following code:

from dolfinx import io
with io.VTKFile(mesh.comm, "output2.pvd", "w") as vtk:
    vtk.write([E._cpp_object])
with io.XDMFFile(mesh.comm, "output2.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(E)

these commands didn’t return any errors but as i put my pvd-file in Paraview for visualisation the ParaView software crashes.

6

For using the visualization on 2D vector fields, see:
https://jorgensd.github.io/dolfinx-tutorial/chapter3/component_bc.html#visualization
Without supplying information regarding what version of paraview and what error message you get, i cannot really help you any further.

As an additional note, the code you supplied in the previous post is not a minimal example, as you cannot simply copy paste it and reproduce the error message.