FEniCSx, Gmsh, electrostatics: how to model point charges?

variational formulation
1

Dear all.

In a capacitor, I would like to model point charges (some electron charges strong), which are on the surface of an dielectric insulating film (cylinder) that is sandwiched between a conducting support (cylinder) and an AFM tip (sphere). I use FEniCSx v 0.6.0 and Gmsh 4.11.0 for describing the 3D objects, see the MWE at the end of this message.

How can I do this?

I found a description of the PointLoad class but it seems that this is not yet implemented in FEniCSx (see here and here). And to be honest, I also don’t know how to deal with this (sorry). :wink:

(EDIT: I modified the text and removed the code, which was just below, because they both were quite misleading …) Because I work with object tags and Gmsh, I have to create a point in the Gmsh mesh and a corresponding tag such that I can somewhat assign a charge at this point, is this right?

Thanks for some comments in advance.

MWE

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

cut_angle          =   15.0

# Domain
domain_size        =  500.0   # in nm

# Sphere
sphere_size        =   15.0   # in nm

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

# Conducting plane
plane_size         =  400.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)


# 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_unrolled"))
gmsh.write(os.path.join(dir_base, "Mesh.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)
1 Like
2

Re all.

So, I have created a charge by just inserting a very small sphere into the Gmsh mesh. Then, I assign a large eps value for this small sphere (1e10) and a charge (e.g., 5 electrons). With this, I can put the small charged sphere into the roh data, see MWE.

If the sphere is very small, I get indeed Coulomb’s law (as can be seen in my much larger program, not shown here). However, I ask myself if one cannot build a ‘real point source’ into the FEM code, by using something similar like this … .

Anyway, so far, what I have is sufficient … . Cheers.

import os
import numpy           as np
import scipy.constants as const

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 locate_dofs_geometrical
from dolfinx.fem       import Expression
from dolfinx.fem       import VectorFunctionSpace
from dolfinx.fem       import assemble_scalar
from dolfinx.fem       import form
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io        import VTXWriter
from dolfinx.mesh      import locate_entities
from ufl               import Measure
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           =    1.0   # in V
V_plane            =    0.0   # in V
eps_domain         =    1.0
eps_dielec         =    7.0
eps_charge         =    1e10

cut_angle          =   15.0

# Domain
domain_size        =  500.0   # in nm

# Sphere
sphere_size        =   15.0   # in nm

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

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

# Charge
charge_size        =  1.0         # in nm
charge_x           =  0.0         # Position in x
charge_y           =  0.0         # Position in y
charge_z           =  charge_size # Position in z, charge is on the film surface
charge_q           = -5           # Nr. electrons

# ___________________________________________________________________ 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)]

# Create charge.
ID_charge = gmsh.model.occ.addSphere(charge_x,
                                     charge_y,
                                     charge_z,
                                     charge_size)
bb_charge = gmsh.model.occ.getBoundingBox(3, ID_charge)
bb_charge = [x+y for x, y in zip(bb_charge, 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_3 = gmsh.model.occ.fragment(cut_domain_2[0], [(3, ID_film)],
                                       removeObject=True, removeTool=True)

# Cut domain with charge.
cut_domain_4 = gmsh.model.occ.fragment(cut_domain_3[0], [(3, ID_charge)],
                                       removeObject=True, removeTool=True)

# 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)

# Cut the charge.
cut_charge_1 = gmsh.model.occ.cut([(3, ID_charge)],
                                  [(3, box_cut_1), (3, box_cut_2)],
                                  removeObject=True, removeTool=False)

# Cut the domain, delete the 2 box objects.
cut_domain_5 = gmsh.model.occ.cut(cut_domain_4[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]
entities  = gmsh.model.occ.getEntitiesInBoundingBox(*bb_charge, dim=2)
enties_surf_charge  = [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]
enties_surf_domain = [item for item in enties_surf_domain if item not in enties_surf_charge]
# 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]
entities   = gmsh.model.occ.getEntitiesInBoundingBox(*bb_charge, dim=3)
enties_vol_charge   = [x[1] for x in entities]
# Domain: remove the volume of the film and charge.
enties_vol_domain = [item for item in enties_vol_domain if item not in enties_vol_film]
enties_vol_domain = [item for item in enties_vol_domain if item not in enties_vol_charge]

# 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")

ID_group_charge = gmsh.model.addPhysicalGroup(2, enties_surf_charge,  50)
gmsh.model.setPhysicalName(2, ID_group_charge, "Surface charge")

# 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")

ID_group_Vol_charge = gmsh.model.addPhysicalGroup(3, enties_vol_charge, 300)
gmsh.model.setPhysicalName(3, ID_group_Vol_charge, "Volume charge")

# Some mesh refinement:
gmsh.model.mesh.field.add("Box", 1)
gmsh.model.mesh.field.setNumber(1, "VIn",    1.0)
gmsh.model.mesh.field.setNumber(1, "VOut",  80.0)
gmsh.model.mesh.field.setNumber(1, "XMax",  plane_size  * 0.07)
gmsh.model.mesh.field.setNumber(1, "XMin", -plane_size  * 0.07)
gmsh.model.mesh.field.setNumber(1, "YMax",  plane_size  * 0.07)
gmsh.model.mesh.field.setNumber(1, "YMin", -plane_size  * 0.07)
gmsh.model.mesh.field.setNumber(1, "ZMax",  sphere_size * 4.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_unrolled"))
gmsh.write(os.path.join(dir_base, "Mesh.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()

dx = Measure("dx", domain=mesh, subdomain_data=cell_markers)

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)

charge_V = (360/cut_angle) * assemble_scalar(form(ScalarType(1.0) * dx(300)))
print(charge_V, 4/3*np.pi*charge_size**3)

# The epsilon values
eps = Function(Q)
for tag in material_tags:
    cells = cell_markers.find(tag)
    # Set values for eps
    if   tag == 300:
        eps_ = eps_charge
    elif 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 == 300:
        roh_ = charge_q * const.e / (charge_V)
    elif tag == 200:
        roh_ = 0.0
    else:
        roh_ = 0.0

    roh_ = roh_ * 10.0**(9.0) / const.epsilon_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.

#aha = dot(PC,v)*dx

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)