FEniCSx, electrostatic field: how to 'project' or how to show the field lines?

dolfinx
1

Dear all.

I’m currently porting my FEniCS project to FEniCSx. The project basically deals with the electrostatics of conductors and insulators in atomic force microscopy (tip, surface, insulating films, etc.) and solves the standard Poisson equation. Everything is fine: I get what I want, like the electrostatic potential u between all conductors and inside insulating media.

What I would like to do is to obtain the field lines and to make them visible in ParaView. In FEniCS, I used a code similar to:

# u is known ...

# E-field.
E = -grad(u)

# Fact normals.
n = FacetNormal(mesh)

E_field = project(E,
                  VectorFunctionSpace(mesh,'P',1),
                  solver_type="mumps",
                  form_compiler_parameters={"cpp_optimize": True,
                                            "representation": "uflacs",
                                            "quadrature_degree": 2} )

How can I do this in FEniCSx? Thanks for some ideas … .

2

I would strongly suggest interpolating

Into a P-th order vector DG space (assuming u is a Pth order scalar Lagrange function).

This operation is a lot faster than doing a projection.
See for instance:
Electromagnetics example — FEniCSx tutorial,
In particular

W = VectorFunctionSpace(mesh, ("DG", 0))
B = Function(W)
B_expr = Expression(as_vector((A_z.dx(1), -A_z.dx(0))), W.element.interpolation_points())
B.interpolate(B_expr)

and use VTXWriter to write a BP file you can open in paraview

1 Like
3

Thanks a lot dokken.

Although I think I do a similar thing like you (see WE), FEniCSx outputs an error at the end

with VTXWriter(mesh.comm, os.path.join(dir_base, "E_field.bp"), B) as vtx:
    vtx.write(0)

RuntimeError: VTK does not support cell-wise fields. See Adios VTX reader (VTU) doesn't support CellData (#18458) · Issues · VTK / VTK · GitLab.

Here is the WE (forget about the Gmsh stuff, is quite straight forward …):

import os
import numpy as np

import gmsh
from dolfinx.io        import gmshio
from dolfinx.io        import XDMFFile
from dolfinx.fem       import dirichletbc
from dolfinx.fem       import Function
from dolfinx.fem       import FunctionSpace
from dolfinx.fem       import locate_dofs_topological
from dolfinx.fem       import Expression
from dolfinx.fem       import VectorFunctionSpace
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io        import VTXWriter
from ufl               import TestFunction
from ufl               import TrialFunction
from ufl               import dot
from ufl               import dx
from ufl               import grad
from ufl               import as_vector
from petsc4py.PETSc    import ScalarType
from mpi4py            import MPI


# __________________________________________________________________ Parameters

z                  =    5.0   # in nm
V_sphere           =    5.0   # in V
V_plane            =   -2.0   # in V
eps_domain         =    1.0
eps_dielec         =    7.0
eps_domain         =    1.0

FLAG_cut_vertic    =    True
cut_angle          =   15.0

# Domain
domain_size        =  500.0   # in nm

# Sphere
sphere_size        =   15.0   # in nm

# Dielectric film
film_size          =   80.0   # in nm
film_height        =   40.0   # in nm

# Conducting plane
plane_size         =  100.0   # in nm
plane_height       =   10.0   # in nm
plane_vx           =    0.0
plane_vy           =    0.0
plane_vz           =   10.0   # in nm


# ___________________________________________________________________ Directory

text_spaces = "    "
dir_base = "./" + os.path.basename(__file__)[:-3]
if not os.path.isdir(dir_base):
    os.mkdir(dir_base)


# ________________________________________________________________________ Gmsh

# Adjust distances.
sphere_z      =  sphere_size + z
plane_z       = -plane_vz - film_height
film_z        = -film_height

# Increase a tiny bit the bounding boxes, see below.
add     = 0.01
bb_corr = [-add, -add, -add, add, add, add]

# Initialize gmsh.
gmsh.initialize()

# Create domain (sphere).
ID_domain = gmsh.model.occ.addSphere(0.0,
                                     0.0,
                                     0.0,
                                     domain_size)
bb_domain = gmsh.model.occ.getBoundingBox(3, ID_domain)
bb_domain = [x+y for x, y in zip(bb_domain, bb_corr)]


# Create sphere.
ID_sphere = gmsh.model.occ.addSphere(0.0,
                                     0.0,
                                     sphere_z,
                                     sphere_size)
bb_sphere = gmsh.model.occ.getBoundingBox(3, ID_sphere)
bb_sphere = [x+y for x, y in zip(bb_sphere, bb_corr)]

# Create dielectric film.
ID_film   = gmsh.model.occ.addCylinder(0.0,
                                       0.0,
                                       film_z,
                                       0.0,
                                       0.0,
                                       film_height,
                                       film_size,
                                       angle = 2.0 * np.pi)
bb_film = gmsh.model.occ.getBoundingBox(3, ID_film)
bb_film = [x+y for x, y in zip(bb_film, bb_corr)]

# Create plane.
ID_plane   = gmsh.model.occ.addCylinder(0.0,
                                        0.0,
                                        plane_z,
                                        0.0,
                                        0.0,
                                        plane_height,
                                        plane_size,
                                        angle = 2.0 * np.pi)
bb_plane = gmsh.model.occ.getBoundingBox(3, ID_plane)
bb_plane = [x+y for x, y in zip(bb_plane, bb_corr)]



# Cut domain with sphere.
cut_domain_1 = gmsh.model.occ.cut([(3, ID_domain)], [(3, ID_sphere)],
                                  removeObject=True, removeTool=True)

# Cut domain with plane.
cut_domain_2 = gmsh.model.occ.cut(cut_domain_1[0], [(3, ID_plane)],
                                  removeObject=True, removeTool=True)

# Cut domain with film.
cut_domain_2 = gmsh.model.occ.fragment(cut_domain_2[0], [(3, ID_film)],
                                       removeObject=True, removeTool=True)

if FLAG_cut_vertic:

    # Create two cutting boxes.
    d = domain_size * 1.5
    box_cut_1 = gmsh.model.occ.addBox(-0.0,-d,-d, 2.0*d,2.0*d,2.0*d)
    box_cut_2 = gmsh.model.occ.addBox(-0.0,-d,-d, 2.0*d,2.0*d,2.0*d)

    # Rotate box No. 2.
    gmsh.model.occ.rotate([(3, box_cut_2)], 0.0, 0.0, 0.0,  0.0, 0.0, 1.0,
                           (180.0-cut_angle)/360 * (2.0*np.pi))

    # Cut the film.
    cut_film_1 = gmsh.model.occ.cut([(3, ID_film)],
                                    [(3, box_cut_1), (3, box_cut_2)],
                                    removeObject=True, removeTool=False)

    # Domain, delete the 2 box objects.
    cut_domain_3 = gmsh.model.occ.cut(cut_domain_2[0],
                                      [(3, box_cut_1), (3, box_cut_2)],
                                      removeObject=True, removeTool=True)

# Synchronize.
gmsh.model.occ.synchronize()

# Entities for surfaces
entities = gmsh.model.occ.getEntitiesInBoundingBox(*bb_domain, dim=2)
enties_surf_domain = [x[1] for x in entities]
entities = gmsh.model.occ.getEntitiesInBoundingBox(*bb_sphere, dim=2)
enties_surf_sphere = [x[1] for x in entities]
entities   = gmsh.model.occ.getEntitiesInBoundingBox(*bb_film, dim=2)
enties_surf_film   = [x[1] for x in entities]
entities  = gmsh.model.occ.getEntitiesInBoundingBox(*bb_plane, dim=2)
enties_surf_plane  = [x[1] for x in entities]
# Domain: remove all surfaces from the other objects.
enties_surf_domain = [item for item in enties_surf_domain if item not in enties_surf_sphere]
enties_surf_domain = [item for item in enties_surf_domain if item not in enties_surf_film]
enties_surf_domain = [item for item in enties_surf_domain if item not in enties_surf_plane]
# Domain: only sphere-like objects
enties_surf_domain = [x for x in enties_surf_domain if gmsh.model.getType(2, x) == 'Sphere']
# Film: remove side planes
planes_remove = []
for plane in enties_surf_film:
    if gmsh.model.getType(2, plane) == 'Plane':
        normal = gmsh.model.getNormal(plane, [0.0,0.0])
        # If the z component is zero then we found the plane to be removed.
        if abs(normal[2]) < 0.1:
            planes_remove.append(plane)
enties_surf_film = [item for item in enties_surf_film if item not in planes_remove + [25]]
enties_surf_film = [item for item in enties_surf_film if item not in enties_surf_plane]

# Entities for volumes.
entities = gmsh.model.occ.getEntitiesInBoundingBox(*bb_domain, dim=3)
enties_vol_domain = [x[1] for x in entities]
entities   = gmsh.model.occ.getEntitiesInBoundingBox(*bb_film, dim=3)
enties_vol_film   = [x[1] for x in entities]
# Domain: remove the volume of the film.
enties_vol_domain = [item for item in enties_vol_domain if item not in enties_vol_film]

# Physical surfaces
ID_group_domain = gmsh.model.addPhysicalGroup(2, enties_surf_domain, 10)
gmsh.model.setPhysicalName(2, ID_group_domain, "Surface domain")

ID_group_sphere = gmsh.model.addPhysicalGroup(2, enties_surf_sphere, 20)
gmsh.model.setPhysicalName(2, ID_group_sphere, "Surface sphere")

ID_group_film = gmsh.model.addPhysicalGroup(2, enties_surf_film,     30)
gmsh.model.setPhysicalName(2, ID_group_film, "Surface film")

ID_group_plane = gmsh.model.addPhysicalGroup(2, enties_surf_plane,   40)
gmsh.model.setPhysicalName(2, ID_group_plane, "Surface plane")

# Physical volumes
ID_group_Vol_domain = gmsh.model.addPhysicalGroup(3, enties_vol_domain, 100)
gmsh.model.setPhysicalName(3, ID_group_Vol_domain, "Volume domain")

ID_group_Vol_film = gmsh.model.addPhysicalGroup(3, enties_vol_film,     200)
gmsh.model.setPhysicalName(3, ID_group_Vol_film, "Volume film")

# Some mesh refinement:
gmsh.model.mesh.field.add("Box", 1)
gmsh.model.mesh.field.setNumber(1, "VIn",    2.0)
gmsh.model.mesh.field.setNumber(1, "VOut",  80.0)
gmsh.model.mesh.field.setNumber(1, "XMax",  plane_size  * 0.8)
gmsh.model.mesh.field.setNumber(1, "XMin", -plane_size  * 0.8)
gmsh.model.mesh.field.setNumber(1, "YMax",  plane_size  * 0.8)
gmsh.model.mesh.field.setNumber(1, "YMin", -plane_size  * 0.8)
gmsh.model.mesh.field.setNumber(1, "ZMax",  sphere_size * 3.0)
gmsh.model.mesh.field.setNumber(1, "ZMin", -film_height)
gmsh.model.mesh.field.setAsBackgroundMesh(1)

# Generate the mesh.
gmsh.model.mesh.generate(3)

# Write the files.
gmsh.write(os.path.join(dir_base, "Mesh_geo.geo_unrolled"))
gmsh.write(os.path.join(dir_base, "Mesh_geo.msh"))


# ______________________________________________ Mesh conversion and re-reading

gmsh_model_rank = 0
mesh_comm       = MPI.COMM_WORLD

mesh, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model,
                                                         mesh_comm,
                                                         gmsh_model_rank,
                                                         gdim=3)
gmsh.finalize()


with XDMFFile(MPI.COMM_WORLD, os.path.join(dir_base, "mt.xdmf"), "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_markers)
    xdmf.write_meshtags(facet_markers)


# ______________________________________ (Sub)domains, boundary conditions, etc.

Q = FunctionSpace(mesh, ("DG", 0))
material_tags = np.unique(cell_markers.values)

# The epsilon values
eps = Function(Q)
for tag in material_tags:
    cells = cell_markers.find(tag)
    # Set values for eps
    if tag == 200:
        eps_ = eps_dielec
    else:
        eps_ = eps_domain

    eps.x.array[cells] = np.full_like(cells, eps_, dtype=ScalarType)

# The roh values
roh = Function(Q)
for tag in material_tags:
    cells = cell_markers.find(tag)
    # Set values for eps
    if tag == 200:
        roh_ = 0.0
    else:
        roh_ = 0.0

    roh.x.array[cells] = np.full_like(cells, roh_, dtype=ScalarType)


# Mesh function
V = FunctionSpace(mesh, ("CG", 1))

list_bcs = []
dofs   = locate_dofs_topological(V, 2, facet_markers.find(10))
list_bcs.append(dirichletbc(ScalarType(0.0), dofs, V))
dofs   = locate_dofs_topological(V, 2, facet_markers.find(20))
list_bcs.append(dirichletbc(ScalarType(V_sphere), dofs, V))
dofs   = locate_dofs_topological(V, 2, facet_markers.find(40))
list_bcs.append(dirichletbc(ScalarType(V_plane), dofs, V))

# Test and trial functions
u = TrialFunction(V)
v = TestFunction(V)

# Create a and L.
a = dot(grad(u), grad(v)) * eps * dx
L = roh * v * dx

A_z = Function(V)
problem = LinearProblem(a, L, u=A_z, bcs=list_bcs)
problem.solve()

xdmf = XDMFFile(mesh.comm, os.path.join(dir_base, "Potential.xdmf"), "w")
xdmf.write_mesh(mesh)
xdmf.write_function(A_z, 0)

E = -grad(A_z)

W = VectorFunctionSpace(mesh, ("DG", 0))
B = Function(W)
B_expr = Expression(E, W.element.interpolation_points())
B.interpolate(B_expr)

with VTXWriter(mesh.comm, os.path.join(dir_base, "E_field.bp"), B) as vtx:
    vtx.write(0)

EDIT:

  1. There must be E inside Expression, I guess.
  2. I now have a problem to write the bp file …
4

Just interpolate into DG-1 for now:

W = VectorFunctionSpace(mesh, ("DG", 1))
B = Function(W)
B_expr = Expression(E, W.element.interpolation_points())
B.interpolate(B_expr)

with VTXWriter(mesh.comm, os.path.join(dir_base, "E_field.bp"), B) as vtx:
    vtx.write(0)
1 Like
5

Yes, true, I found the same yesterday evening. Thanks a lot.

So here the full WE:

import os
import numpy as np

import gmsh
from dolfinx.io        import gmshio
from dolfinx.io        import XDMFFile
from dolfinx.fem       import dirichletbc
from dolfinx.fem       import Function
from dolfinx.fem       import FunctionSpace
from dolfinx.fem       import locate_dofs_topological
from dolfinx.fem       import Expression
from dolfinx.fem       import VectorFunctionSpace
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io        import VTXWriter
from ufl               import TestFunction
from ufl               import TrialFunction
from ufl               import dot
from ufl               import dx
from ufl               import grad
from ufl               import as_vector
from petsc4py.PETSc    import ScalarType
from mpi4py            import MPI


# __________________________________________________________________ Parameters

z                  =    5.0   # in nm
V_sphere           =    5.0   # in V
V_plane            =   -2.0   # in V
eps_domain         =    1.0
eps_dielec         =    7.0
eps_domain         =    1.0

FLAG_cut_vertic    =    True
cut_angle          =   15.0

# Domain
domain_size        =  500.0   # in nm

# Sphere
sphere_size        =   15.0   # in nm

# Dielectric film
film_size          =   80.0   # in nm
film_height        =   40.0   # in nm

# Conducting plane
plane_size         =  100.0   # in nm
plane_height       =   10.0   # in nm
plane_vx           =    0.0
plane_vy           =    0.0
plane_vz           =   10.0   # in nm


# ___________________________________________________________________ Directory

text_spaces = "    "
dir_base = "./" + os.path.basename(__file__)[:-3]
if not os.path.isdir(dir_base):
    os.mkdir(dir_base)


# ________________________________________________________________________ Gmsh

# Adjust distances.
sphere_z      =  sphere_size + z
plane_z       = -plane_vz - film_height
film_z        = -film_height

# Increase a tiny bit the bounding boxes, see below.
add     = 0.01
bb_corr = [-add, -add, -add, add, add, add]

# Initialize gmsh.
gmsh.initialize()

# Create domain (sphere).
ID_domain = gmsh.model.occ.addSphere(0.0,
                                     0.0,
                                     0.0,
                                     domain_size)
bb_domain = gmsh.model.occ.getBoundingBox(3, ID_domain)
bb_domain = [x+y for x, y in zip(bb_domain, bb_corr)]


# Create sphere.
ID_sphere = gmsh.model.occ.addSphere(0.0,
                                     0.0,
                                     sphere_z,
                                     sphere_size)
bb_sphere = gmsh.model.occ.getBoundingBox(3, ID_sphere)
bb_sphere = [x+y for x, y in zip(bb_sphere, bb_corr)]

# Create dielectric film.
ID_film   = gmsh.model.occ.addCylinder(0.0,
                                       0.0,
                                       film_z,
                                       0.0,
                                       0.0,
                                       film_height,
                                       film_size,
                                       angle = 2.0 * np.pi)
bb_film = gmsh.model.occ.getBoundingBox(3, ID_film)
bb_film = [x+y for x, y in zip(bb_film, bb_corr)]

# Create plane.
ID_plane   = gmsh.model.occ.addCylinder(0.0,
                                        0.0,
                                        plane_z,
                                        0.0,
                                        0.0,
                                        plane_height,
                                        plane_size,
                                        angle = 2.0 * np.pi)
bb_plane = gmsh.model.occ.getBoundingBox(3, ID_plane)
bb_plane = [x+y for x, y in zip(bb_plane, bb_corr)]



# Cut domain with sphere.
cut_domain_1 = gmsh.model.occ.cut([(3, ID_domain)], [(3, ID_sphere)],
                                  removeObject=True, removeTool=True)

# Cut domain with plane.
cut_domain_2 = gmsh.model.occ.cut(cut_domain_1[0], [(3, ID_plane)],
                                  removeObject=True, removeTool=True)

# Cut domain with film.
cut_domain_2 = gmsh.model.occ.fragment(cut_domain_2[0], [(3, ID_film)],
                                       removeObject=True, removeTool=True)

if FLAG_cut_vertic:

    # Create two cutting boxes.
    d = domain_size * 1.5
    box_cut_1 = gmsh.model.occ.addBox(-0.0,-d,-d, 2.0*d,2.0*d,2.0*d)
    box_cut_2 = gmsh.model.occ.addBox(-0.0,-d,-d, 2.0*d,2.0*d,2.0*d)

    # Rotate box No. 2.
    gmsh.model.occ.rotate([(3, box_cut_2)], 0.0, 0.0, 0.0,  0.0, 0.0, 1.0,
                           (180.0-cut_angle)/360 * (2.0*np.pi))

    # Cut the film.
    cut_film_1 = gmsh.model.occ.cut([(3, ID_film)],
                                    [(3, box_cut_1), (3, box_cut_2)],
                                    removeObject=True, removeTool=False)

    # Domain, delete the 2 box objects.
    cut_domain_3 = gmsh.model.occ.cut(cut_domain_2[0],
                                      [(3, box_cut_1), (3, box_cut_2)],
                                      removeObject=True, removeTool=True)

# Synchronize.
gmsh.model.occ.synchronize()

# Entities for surfaces
entities = gmsh.model.occ.getEntitiesInBoundingBox(*bb_domain, dim=2)
enties_surf_domain = [x[1] for x in entities]
entities = gmsh.model.occ.getEntitiesInBoundingBox(*bb_sphere, dim=2)
enties_surf_sphere = [x[1] for x in entities]
entities   = gmsh.model.occ.getEntitiesInBoundingBox(*bb_film, dim=2)
enties_surf_film   = [x[1] for x in entities]
entities  = gmsh.model.occ.getEntitiesInBoundingBox(*bb_plane, dim=2)
enties_surf_plane  = [x[1] for x in entities]
# Domain: remove all surfaces from the other objects.
enties_surf_domain = [item for item in enties_surf_domain if item not in enties_surf_sphere]
enties_surf_domain = [item for item in enties_surf_domain if item not in enties_surf_film]
enties_surf_domain = [item for item in enties_surf_domain if item not in enties_surf_plane]
# Domain: only sphere-like objects
enties_surf_domain = [x for x in enties_surf_domain if gmsh.model.getType(2, x) == 'Sphere']
# Film: remove side planes
planes_remove = []
for plane in enties_surf_film:
    if gmsh.model.getType(2, plane) == 'Plane':
        normal = gmsh.model.getNormal(plane, [0.0,0.0])
        # If the z component is zero then we found the plane to be removed.
        if abs(normal[2]) < 0.1:
            planes_remove.append(plane)
enties_surf_film = [item for item in enties_surf_film if item not in planes_remove + [25]]
enties_surf_film = [item for item in enties_surf_film if item not in enties_surf_plane]

# Entities for volumes.
entities = gmsh.model.occ.getEntitiesInBoundingBox(*bb_domain, dim=3)
enties_vol_domain = [x[1] for x in entities]
entities   = gmsh.model.occ.getEntitiesInBoundingBox(*bb_film, dim=3)
enties_vol_film   = [x[1] for x in entities]
# Domain: remove the volume of the film.
enties_vol_domain = [item for item in enties_vol_domain if item not in enties_vol_film]

# Physical surfaces
ID_group_domain = gmsh.model.addPhysicalGroup(2, enties_surf_domain, 10)
gmsh.model.setPhysicalName(2, ID_group_domain, "Surface domain")

ID_group_sphere = gmsh.model.addPhysicalGroup(2, enties_surf_sphere, 20)
gmsh.model.setPhysicalName(2, ID_group_sphere, "Surface sphere")

ID_group_film = gmsh.model.addPhysicalGroup(2, enties_surf_film,     30)
gmsh.model.setPhysicalName(2, ID_group_film, "Surface film")

ID_group_plane = gmsh.model.addPhysicalGroup(2, enties_surf_plane,   40)
gmsh.model.setPhysicalName(2, ID_group_plane, "Surface plane")

# Physical volumes
ID_group_Vol_domain = gmsh.model.addPhysicalGroup(3, enties_vol_domain, 100)
gmsh.model.setPhysicalName(3, ID_group_Vol_domain, "Volume domain")

ID_group_Vol_film = gmsh.model.addPhysicalGroup(3, enties_vol_film,     200)
gmsh.model.setPhysicalName(3, ID_group_Vol_film, "Volume film")

# Some mesh refinement:
gmsh.model.mesh.field.add("Box", 1)
gmsh.model.mesh.field.setNumber(1, "VIn",    2.0)
gmsh.model.mesh.field.setNumber(1, "VOut",  80.0)
gmsh.model.mesh.field.setNumber(1, "XMax",  plane_size  * 0.8)
gmsh.model.mesh.field.setNumber(1, "XMin", -plane_size  * 0.8)
gmsh.model.mesh.field.setNumber(1, "YMax",  plane_size  * 0.8)
gmsh.model.mesh.field.setNumber(1, "YMin", -plane_size  * 0.8)
gmsh.model.mesh.field.setNumber(1, "ZMax",  sphere_size * 3.0)
gmsh.model.mesh.field.setNumber(1, "ZMin", -film_height)
gmsh.model.mesh.field.setAsBackgroundMesh(1)

#gmsh.option.setNumber("Mesh.Algorithm", 7)
gmsh.model.mesh.optimize("Netgen")

# Generate the mesh.
gmsh.model.mesh.generate(3)


# Write the files.
gmsh.write(os.path.join(dir_base, "Mesh_geo.geo_unrolled"))
gmsh.write(os.path.join(dir_base, "Mesh_geo.msh"))


# ______________________________________________ Mesh conversion and re-reading

gmsh_model_rank = 0
mesh_comm       = MPI.COMM_WORLD

mesh, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model,
                                                         mesh_comm,
                                                         gmsh_model_rank,
                                                         gdim=3)
gmsh.finalize()


with XDMFFile(MPI.COMM_WORLD, os.path.join(dir_base, "mt.xdmf"), "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_markers)
    xdmf.write_meshtags(facet_markers)


# ______________________________________ (Sub)domains, boundary conditions, etc.

Q = FunctionSpace(mesh, ("DG", 0))
material_tags = np.unique(cell_markers.values)

# The epsilon values
eps = Function(Q)
for tag in material_tags:
    cells = cell_markers.find(tag)
    # Set values for eps
    if tag == 200:
        eps_ = eps_dielec
    else:
        eps_ = eps_domain

    eps.x.array[cells] = np.full_like(cells, eps_, dtype=ScalarType)

# The roh values
roh = Function(Q)
for tag in material_tags:
    cells = cell_markers.find(tag)
    # Set values for eps
    if tag == 200:
        roh_ = 0.0
    else:
        roh_ = 0.0

    roh.x.array[cells] = np.full_like(cells, roh_, dtype=ScalarType)


# Mesh function
V = FunctionSpace(mesh, ("CG", 1))

list_bcs = []
dofs   = locate_dofs_topological(V, 2, facet_markers.find(10))
list_bcs.append(dirichletbc(ScalarType(0.0), dofs, V))
dofs   = locate_dofs_topological(V, 2, facet_markers.find(20))
list_bcs.append(dirichletbc(ScalarType(V_sphere), dofs, V))
dofs   = locate_dofs_topological(V, 2, facet_markers.find(40))
list_bcs.append(dirichletbc(ScalarType(V_plane), dofs, V))

# Test and trial functions
u = TrialFunction(V)
v = TestFunction(V)

# Create a and L.
a = dot(grad(u), grad(v)) * eps * dx
L = roh * v * dx

A_z = Function(V)
problem = LinearProblem(a, L, u=A_z, bcs=list_bcs)
problem.solve()

xdmf = XDMFFile(mesh.comm, os.path.join(dir_base, "Potential.xdmf"), "w")
xdmf.write_mesh(mesh)
xdmf.write_function(A_z, 0)


W = VectorFunctionSpace(mesh, ("DG", 1))
E = Function(W)
E_expr = Expression(-grad(A_z), W.element.interpolation_points())
E.interpolate(E_expr)

with VTXWriter(mesh.comm, os.path.join(dir_base, "E_field.bp"), E) as vtx:
    vtx.write(0)