# Neumann Condition does not work after cutting the domain

**URL:** <https://fenicsproject.discourse.group/t/neumann-condition-does-not-work-after-cutting-the-domain/3651>\
**Category:** mesh\
**Created:** [June 26, 2020, 1:52pm UTC](https://fenicsproject.discourse.group/t/neumann-condition-does-not-work-after-cutting-the-domain/3651 "2020-06-26T13:52:14Z")\
**Posts on this page:** 1\
**Showing post:** 4

<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:** [June 30, 2020, 2:17pm UTC](https://fenicsproject.discourse.group/t/neumann-condition-does-not-work-after-cutting-the-domain/3651/4 "2020-06-30T14:17:08Z")

</div>

Dear all.

Yep, and here is the solution 😂, the issue was a mesh thing: if one has two or more ‘connected’ domains, which shall be created by ‘cutting’, then one has to look forward that the interface has only **one** well-defined facet (before I always had 2 !) To do so, one needs to use `geom.boolean_fragments`, see the code at the end.

Cheers all.

Here is an image …

 ![Image_8](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/4/4b1914c3da77a078ee802eb4d85ddfeda23c434e.png)  
… and here is the code. Note that one can now cut under a specific angle … 0.0° \< angle \< 90.0°. If the angle is larger, the tags change and need to be adapted.

```
import pygmsh
import meshio
from fenics import *

import matplotlib.pyplot as plt
import numpy as np
import scipy.constants as const
import os

# __________________________________________________________________ Parameters

z = 5.0
V_sphere = 5.0
V_plane = -2.0
eps_domain = 1.0
eps_dielec = 3.0
eps_domain = 1.0

FLAG_cut_vertic_4 = True
# The cutting angle, valid for: 0.0° < cut_angle < 90.0°
cut_angle = 15.0

# Domain
domain_size = 500.0

# Sphere
sphere_r = 15.0
sphere_x = 0.0
sphere_y = 0.0
sphere_z = 0.0

# Dielectric film
dielectric_size = 80.0
dielectric_height = 40.0
dielectric_x = 0.0
dielectric_y = 0.0
dielectric_z = 0.0
dielectric_vx = 0.0
dielectric_vy = 0.0
dielectric_vz = 0.0

# Conducting plane
plane_size = 100.0
plane_x = 0.0
plane_y = 0.0
plane_z = 0.0
plane_vx = 0.0
plane_vy = 0.0
plane_vz = 10.0

# ___________________________________________________________________ Directory

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

# ____________________________________________________ Mesh creation via pygmsh

sphere_z = sphere_r + z
plane_z = -(plane_vz + dielectric_height + 8.0)

# The height of the dielectric film
dielectric_vz = -dielectric_height

# We use opencascade for the 3D objects.
geom = pygmsh.opencascade.Geometry()

domain = geom.add_ball([0.0, 0.0, 0.0], domain_size)

sphere = geom.add_ball([sphere_x, sphere_y, sphere_z], sphere_r)

plane = geom.add_cylinder([plane_x, plane_y, plane_z],
                            [plane_vx, plane_vy, plane_vz],
                            plane_size, 2.0*const.pi)

dielect = geom.add_cylinder([dielectric_x, dielectric_y, dielectric_z],
                             [dielectric_vx, dielectric_vy, dielectric_vz],
                              dielectric_size, 2.0*const.pi)

# ... take the difference (boolean).
space_1 = geom.boolean_difference([domain],[sphere],
                                delete_first=True,
                                delete_other=True)

space_2 = geom.boolean_difference([space_1],[plane],
                                delete_first=True,
                                delete_other=True)

# We vertically cut out a 4th of the scene.
if FLAG_cut_vertic_4:
    d = domain_size*1.5
    box_cut_1 = geom.add_box([0.0 , -d, -d],
                             [2.0*d, 2.0*d, 2.0*d], char_length=1.0)
    box_cut_2 = geom.add_box([0.0 , -d, -d],
                             [2.0*d, 2.0*d, 2.0*d], char_length=1.0)
    # Rotate box No. 2
    geom.rotate(box_cut_2,
                point=(0.0, 0.0, 0.0),
                axis =(0.0, 0.0, 1.0),
                angle=(180.0-cut_angle)/360.0 * (2.0*const.pi))

    # First, we cut the dielectric film
    dielect_final = geom.boolean_difference([dielect],[box_cut_1, box_cut_2],
                                            delete_first=True,
                                            delete_other=False)

    # Now we cut the domain.
    space_3 = geom.boolean_difference([space_2],[box_cut_1, box_cut_2],
                                      delete_first=True,
                                      delete_other=True)

    # And here it comes: we have to use fractions such that we don't have
    # multiple interfaces. This was the problem before!
    geom.boolean_fragments([space_3],[dielect_final],
                           delete_first=True,
                           delete_other=True)
else:
    # See just above.
    geom.boolean_fragments([space_2],[dielect],
                           delete_first=True,
                           delete_other=True)

## This is for finding the identifiers ...
#mesh = pygmsh.generate_mesh(geom, geo_filename=os.path.join(dir_base, "base.geo"))
#exit(-1)

if FLAG_cut_vertic_4:
    # Sphere
    geom.add_raw_code("Physical Surface(\"1\") = {2};")
    # Plane
    geom.add_raw_code("Physical Surface(\"2\") = {7,8,9};")
    # Domain
    geom.add_raw_code("Physical Surface(\"3\") = {10};") # Sides: 1,3
    # Dielectricum
    geom.add_raw_code("Physical Surface(\"4\") = {4,5,6};") # Sides: 11,12
    # Dielectricum interface
    geom.add_raw_code("Physical Surface(\"5\") = {4,5,6};")
    # Volume space
    geom.add_raw_code("Physical Volume(\"6\") = {33};") # This needs to be determined manually!
    # Volume dielectricum
    geom.add_raw_code("Physical Volume(\"7\") = {32};") # This needs to be determined manually!

else:
    # Sphere
    geom.add_raw_code("Physical Surface(\"1\") = {10};")
    # Plane
    geom.add_raw_code("Physical Surface(\"2\") = {11,12,13};")
    # Domain
    geom.add_raw_code("Physical Surface(\"3\") = {9};")
    # Dielectricum
    geom.add_raw_code("Physical Surface(\"4\") = {6,7,8};")
    # Dielectricum interface
    geom.add_raw_code("Physical Surface(\"5\") = {6,7,8};")
    # Volume space
    geom.add_raw_code("Physical Volume(\"6\") = {11};")
    # Volume dielectricum
    geom.add_raw_code("Physical Volume(\"7\") = {10};")

# Some quick mesh refinement:
geom.add_raw_code("Field[1] = Box;")
geom.add_raw_code("Field[1].VIn = " + str(5.0) + ";")
geom.add_raw_code("Field[1].VOut = " + str(50.0) + ";")
geom.add_raw_code("Field[1].XMax = " + str(dielectric_size) + ";")
geom.add_raw_code("Field[1].XMin = " + str(-dielectric_size) + ";")
geom.add_raw_code("Field[1].YMax = " + str(dielectric_size) + ";")
geom.add_raw_code("Field[1].YMin = " + str(-dielectric_size) + ";")
geom.add_raw_code("Field[1].ZMax = " + str(sphere_r*3.0) + ";")
geom.add_raw_code("Field[1].ZMin = " + str(-dielectric_height) + ";")
geom.add_raw_code("Background Field = 1;")

# Build the .geo file.
file_geo = os.path.join(dir_base, "Geometry.geo")
mesh = pygmsh.generate_mesh(geom, geo_filename=file_geo)

# ______________________________________________ Mesh conversion and re-reading

tetra_cells = None
tetra_data = None
triangle_cells = None
triangle_data = None

for key in mesh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = mesh.cell_data_dict["gmsh:physical"][key]
    elif key == "tetra":
        tetra_data = mesh.cell_data_dict["gmsh:physical"][key]

for cell in mesh.cells:
    if cell.type == "tetra":
        if tetra_cells is None:
            tetra_cells = cell.data
        else:
            tetra_cells = np.vstack([tetra_cells, cell.data])
            print("aha")
    elif cell.type == "triangle":
        if triangle_cells is None:
            triangle_cells = cell.data
        else:
            triangle_cells = np.vstack([triangle_cells, cell.data])

tetra_mesh = meshio.Mesh(points=mesh.points,
                         cells={"tetra": tetra_cells},
                         cell_data={"name_to_read":[tetra_data]})

triangle_mesh =meshio.Mesh(points=mesh.points,
                           cells=[("triangle", triangle_cells)],
                           cell_data={"name_to_read":[triangle_data]})
meshio.write(os.path.join(dir_base,"tetra_mesh.xdmf"), tetra_mesh)
meshio.write(os.path.join(dir_base,"mf.xdmf") , triangle_mesh)

# All the latter .xdmf files are read such that they can be used in FEniCS.
mesh = Mesh()
with XDMFFile(os.path.join(dir_base,"tetra_mesh.xdmf")) as infile:
    infile.read(mesh)

mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile(os.path.join(dir_base,"mf.xdmf")) as infile:
    infile.read(mvc, "name_to_read")
boundaries = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc2 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile(os.path.join(dir_base,"tetra_mesh.xdmf")) as infile:
    infile.read(mvc2, "name_to_read")
subdomains = cpp.mesh.MeshFunctionSizet(mesh, mvc2)

dx = Measure("dx", domain=mesh, subdomain_data=subdomains)
ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
dS = Measure("dS", domain=mesh, subdomain_data=boundaries)

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

# The epsilon in divers considered media (vacuum, dielectric, etc.)
def get_epsilon(subdomains):

    # Here, we define the eps for the thin film (dielectric_eps) and the vacuum.
    class permittivity(UserExpression):
        def __init__ (self, markers, **kwargs):
            super(permittivity, self). __init__ (**kwargs)
            self.markers = markers

        def eval_cell(self, values, x, cell):
            if self.markers[cell.index] == 7:
                values[0] = eps_dielec
            else:
                values[0] = eps_domain

    return permittivity(subdomains, degree=1)

# Permittivity
eps = get_epsilon(subdomains)

# Mesh function
V = FunctionSpace(mesh, 'P', 2)

# The Dirichlet boundary conditions for the sphere and plane:
b_cond = []
b_cond.append(DirichletBC(V, Constant(V_sphere), boundaries, 1))
b_cond.append(DirichletBC(V, Constant(V_plane), boundaries, 2))
# We put the potential at the domain boundary onto 0V:
b_cond.append(DirichletBC(V, Constant(0.0), boundaries, 3))

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

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

u = Function(V)
solve(a == L, u, b_cond,
    solver_parameters = {"linear_solver" : "gmres",
                        "preconditioner" : "amg"})
if FLAG_cut_vertic_4:
    file_pot = os.path.join(dir_base, "Potential_cut.pvd")
else:
    file_pot = os.path.join(dir_base, "Potential_uncut.pvd")

potentialFile = File(file_pot)
potentialFile << u
```

---

_[View the full topic](https://fenicsproject.discourse.group/t/neumann-condition-does-not-work-after-cutting-the-domain/3651)._