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

**URL:** <https://fenicsproject.discourse.group/t/fenicsx-gmsh-electrostatics-how-to-model-point-charges/11166>\
**Category:** variational formulation\
**Created:** [May 12, 2023, 2:57pm UTC](https://fenicsproject.discourse.group/t/fenicsx-gmsh-electrostatics-how-to-model-point-charges/11166 "2023-05-12T14:57:14Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Nexius](https://avatars.discourse-cdn.com/v4/letter/n/50afbb/32.png) [@Nexius](https://fenicsproject.discourse.group/u/Nexius)\
**Post date:** [May 12, 2023, 2:57pm UTC](https://fenicsproject.discourse.group/t/fenicsx-gmsh-electrostatics-how-to-model-point-charges/11166/1 "2023-05-12T14:57:14Z")

</div>

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](https://fenicsproject.org/olddocs/dolfin/1.5.0/python/programmers-reference/cpp/fem/PointSource.html) but it seems that this is not yet implemented in FEniCSx (see [here](https://fenicsproject.discourse.group/t/pointsource-in-dolfinx/8337) and [here](https://github.com/FEniCS/dolfinx/issues/28)). And to be honest, I also don’t know how to deal with this (sorry). 😉

(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

```auto
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)

```

---

<div class="post-metadata">

**Author:** ![Nexius](https://avatars.discourse-cdn.com/v4/letter/n/50afbb/32.png) [@Nexius](https://fenicsproject.discourse.group/u/Nexius)\
**Post date:** [May 15, 2023, 9:48am UTC](https://fenicsproject.discourse.group/t/fenicsx-gmsh-electrostatics-how-to-model-point-charges/11166/2 "2023-05-15T09:48:48Z")

</div>

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](https://fenicsproject.org/olddocs/dolfin/1.4.0/python/programmers-reference/cpp/fem/PointSource.html) … .

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

```auto
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)

```