# Neumann BC for 2 conductors

**URL:** <https://fenicsproject.discourse.group/t/neumann-bc-for-2-conductors/9379>\
**Category:** dolfinx\
**Created:** [September 28, 2022, 8:55am UTC](https://fenicsproject.discourse.group/t/neumann-bc-for-2-conductors/9379 "2022-09-28T08:55:14Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Marryweather](https://avatars.discourse-cdn.com/v4/letter/m/b782af/32.png) [@Marryweather](https://fenicsproject.discourse.group/u/Marryweather)\
**Post date:** [September 28, 2022, 8:55am UTC](https://fenicsproject.discourse.group/t/neumann-bc-for-2-conductors/9379/1 "2022-09-28T08:55:14Z")

</div>

Hello,

I have a question about Neumann boundary conditions.  
I have a 2D plot in a room where the bottom is the ground with Dirichlet BC = 0 and 2 electrical conductors with a positive voltage and a negative voltage. The surface of the conductors should contain the Neumann boundary condition. I can’t get both conditions to be executed.

Here is the code. Unfortunately it is a bit big.  
But only the bottom part is important.

```auto
#Mesh Erzeugung:(mit GMSH)

from mpi4py import MPI
from dolfinx import mesh
import warnings

warnings.filterwarnings("ignore")

import gmsh

gmsh.initialize()

from dolfinx import fem
from dolfinx import nls
from dolfinx import log
from dolfinx import io
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx.fem import Constant, Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc
from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector, create_vector, set_bc
from dolfinx.graph import create_adjacencylist
from dolfinx.geometry import BoundingBoxTree, compute_collisions, compute_colliding_cells
from dolfinx.io import (XDMFFile, cell_perm_gmsh, distribute_entity_data, extract_gmsh_geometry,
                        extract_gmsh_topology_and_markers, ufl_mesh_from_gmsh)
from dolfinx import geometry
from dolfinx.mesh import create_mesh, meshtags_from_entities
import ufl
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
                 as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)
import numpy as np
import matplotlib.pyplot as plt
import tqdm.notebook
from mpi4py import MPI
from petsc4py import PETSc
import scipy.constants as c

#Domain-Definieren:

L =24.
H =18.
c_x1 = 8. 
c_y1 = 13.
r1 = 0.01
c_x2 = 16. 
c_y2 = 13.
r2 = 0.001
gdim = 2
rank = MPI.COMM_WORLD.rank

if rank == 0:
    rectangle = gmsh.model.occ.addRectangle(0,0,0, L, H, tag=1)
    obstacle1 = gmsh.model.occ.addDisk(c_x1, c_y1, 0, r1, r1)
    obstacle2 = gmsh.model.occ.addDisk(c_x2, c_y2, 0, r2, r2)

if rank == 0:
    Mesh_domain = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle1), (gdim, obstacle2)])
    gmsh.model.occ.synchronize()

#Verbinden aller Mesh-Punkte
Mesh_domain_marker = 1

if rank == 0:
    volumes = gmsh.model.getEntities(dim=gdim)
    gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], Mesh_domain_marker)
    gmsh.model.setPhysicalName(volumes[0][0], Mesh_domain_marker, "Mesh_domain")

inlet_marker, outlet_marker, wall_marker, obstacle1_marker,obstacle2_marker, wall_top = 2, 3, 4, 5, 6, 7
inflow, outflow, walls, obstacle1, obstacle2 , top = [], [], [], [], [], []

if rank == 0:
    boundaries = gmsh.model.getBoundary(volumes, oriented=False)
    for boundary in boundaries:
        center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
        if np.allclose(center_of_mass, [0, H/2, 0]):
            inflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L, H/2, 0]):
            outflow.append(boundary[1])
        elif np.allclose(center_of_mass, [L/2, 0, 0]):
            walls.append(boundary[1])
        elif np.allclose(center_of_mass, [L/2, H, 0]):
            top.append(boundary[1])
        elif np.allclose(center_of_mass, [c_x1, c_y1, 0]):
            obstacle1.append(boundary[1])
        elif np.allclose(center_of_mass, [c_x2, c_y2, 0]):
            obstacle2.append(boundary[1])
            
    gmsh.model.addPhysicalGroup(1, walls, wall_marker)
    gmsh.model.setPhysicalName(1, wall_marker, "Walls")
    gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
    gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
    gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
    gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
    gmsh.model.addPhysicalGroup(1, obstacle1, obstacle1_marker)
    gmsh.model.setPhysicalName(1, obstacle1_marker, "Obstacle1")
    gmsh.model.addPhysicalGroup(1, obstacle2, obstacle2_marker)
    gmsh.model.setPhysicalName(1, obstacle2_marker, "Obstacle2")
    gmsh.model.addPhysicalGroup(1, top, wall_top)
    gmsh.model.setPhysicalName(1, wall_top, "top")

if rank == 0:
    gmsh.option.setNumber('Mesh.MeshSizeMin', 2e-1)
    gmsh.option.setNumber('Mesh.MeshSizeMax', 2e-1)
    gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 8)
    gmsh.option.setNumber("Mesh.RecombineAll", 2)
    gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.optimize("Netgen") 

if MPI.COMM_WORLD.rank == 0:

    # Mesh Geometrie   
    x = extract_gmsh_geometry(gmsh.model)

    # Topologie für jedes Element vom Mesh
    topologies = extract_gmsh_topology_and_markers(gmsh.model)
    num_cell_types = len(topologies.keys())
    cell_information = {}
    cell_dimensions = np.zeros(num_cell_types, dtype=np.int32)
    
    for i, element in enumerate(topologies.keys()):
        properties = gmsh.model.mesh.getElementProperties(element)
        name, dim, order, num_nodes, local_coords, _ = properties
        cell_information[i] = {"id": element, "dim": dim, "num_nodes": num_nodes}
        cell_dimensions[i] = dim

    # Sort elements by ascending dimension

    perm_sort = np.argsort(cell_dimensions)

    # Broadcast cell type data and geometric dimension
    
    cell_id = cell_information[perm_sort[-1]]["id"]
    tdim = cell_information[perm_sort[-1]]["dim"]
    num_nodes = cell_information[perm_sort[-1]]["num_nodes"]
    cell_id, num_nodes = MPI.COMM_WORLD.bcast([cell_id, num_nodes], root=0)

    if tdim - 1 in cell_dimensions:
        num_facet_nodes = MPI.COMM_WORLD.bcast( cell_information[perm_sort[-2]]["num_nodes"], root=0)
        gmsh_facet_id = cell_information[perm_sort[-2]]["id"]
        marked_facets = np.asarray(topologies[gmsh_facet_id]["topology"], dtype=np.int64)
        facet_values = np.asarray(topologies[gmsh_facet_id]["cell_data"], dtype=np.int32)

    cells = np.asarray(topologies[cell_id]["topology"], dtype=np.int64)
    cell_values = np.asarray(topologies[cell_id]["cell_data"], dtype=np.int32)

else:
    cell_id, num_nodes = MPI.COMM_WORLD.bcast([None, None], root=0)
    cells, x = np.empty([0, num_nodes], np.int64), np.empty([0, gdim])
    cell_values = np.empty((0,), dtype=np.int32)
    num_facet_nodes = MPI.COMM_WORLD.bcast(None, root=0)
    marked_facets = np.empty((0, num_facet_nodes), dtype=np.int64)
    facet_values = np.empty((0,), dtype=np.int32)

    
# Create distributed mesh

ufl_domain = ufl_mesh_from_gmsh(cell_id, gdim)
gmsh_cell_perm = cell_perm_gmsh(to_type(str(ufl_domain.ufl_cell())), num_nodes)
cells = cells[:, gmsh_cell_perm]
mesh = create_mesh(MPI.COMM_WORLD, cells, x[:, :gdim], ufl_domain)
tdim = mesh.topology.dim
fdim = tdim - 1

# Permute facets from MSH to DOLFINx ordering

facet_type = cell_entity_type(to_type(str(ufl_domain.ufl_cell())), fdim, 0)
gmsh_facet_perm = cell_perm_gmsh(facet_type, num_facet_nodes)
marked_facets = np.asarray(marked_facets[:, gmsh_facet_perm], dtype=np.int64)

local_entities, local_values = distribute_entity_data(mesh, fdim, marked_facets, facet_values)
mesh.topology.create_connectivity(fdim, tdim)
adj = create_adjacencylist(local_entities)
bb_tree = geometry.BoundingBoxTree(mesh, mesh.topology.dim)

# Create DOLFINx MeshTags

ft = meshtags_from_entities(mesh, fdim, adj, np.int32(local_values))
ft.name = "Facet tags"

#Funktionspace und Vektorfunktionspace defenieren
s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, s_cg1)
V_E= fem.VectorFunctionSpace(mesh, ("DG", 0))
fdim = mesh.topology.dim - 1

E_cp = -16.5e5
E_cn = 16.5e5

#Boundary conditions definieren
# Walls
u1_nonslip = np.array((0,) *mesh.geometry.dim, dtype=PETSc.ScalarType)
u11_nonslip = np.array((0,) *mesh.geometry.dim, dtype=PETSc.ScalarType)
wall_facets = ft.indices[ft.values == wall_marker]
bc_walls = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(V, fdim, wall_facets), V)

#Boundary conditions for equations: 
bcwall = [bc_walls] 

#Laplace
u1 = ufl.TrialFunction(V) 
v1 = ufl.TestFunction(V)
dObt = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=obstacle1_marker)
dObs = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=obstacle2_marker)
a = ufl.inner(ufl.grad(u1), ufl.grad(v1)) * ufl.dx 
x = SpatialCoordinate(mesh)
L = ufl.inner(E_cp, v1) * dObs + ufl.inner(E_cn,v1) * dObt 
problem1 = fem.petsc.LinearProblem(a, L, bcwall, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
phi = problem1.solve()

print("phi_max: ", np.max(phi.x.array), " / phi_min: ", np.min(phi.x.array))
print(phi.x.array)

with io.VTKFile(mesh.comm, "Laplace.pvd", "w") as vtk:
    vtk.write([phi._cpp_object])

with io.XDMFFile(mesh.comm, "Laplace.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(phi)
```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 28, 2022, 12:04pm UTC](https://fenicsproject.discourse.group/t/neumann-bc-for-2-conductors/9379/2 "2022-09-28T12:04:54Z")

</div>

> [@Marryweather](#):
>
> ```auto
> dObt = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=obstacle1_marker)
> dObs = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=obstacle2_marker)
> 
> ```

Start by making sure that the facet markers are correct by visualizing `ft` in Paraview.  
then assemble simple forms such as `dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*dObt))` and `dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*dObs))` and check that they give reasonable results

---

<div class="post-metadata">

**Author:** ![Marryweather](https://avatars.discourse-cdn.com/v4/letter/m/b782af/32.png) [@Marryweather](https://fenicsproject.discourse.group/u/Marryweather)\
**Post date:** [September 29, 2022, 9:15am UTC](https://fenicsproject.discourse.group/t/neumann-bc-for-2-conductors/9379/3 "2022-09-29T09:15:59Z")

</div>

How do I visualize ft in ParaView?

Is that what you mean?

```auto
u1 = ufl.TrialFunction(V) 
v1 = ufl.TestFunction(V)
k = Constant(mesh, PETSc.ScalarType(0))
f = k/c.epsilon_0
dObt = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=obstacle1_marker)
dObs = Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=obstacle2_marker)
fem.assemble_scalar(fem.form(1*dObt))
fem.assemble_scalar(fem.form(1*dObs))
a = ufl.inner(ufl.grad(u1), ufl.grad(v1)) * ufl.dx 
x = SpatialCoordinate(mesh)
L = ufl.inner(E_cp, v1) * dObs + ufl.inner(E_cn,v1) * dObt 
u1 = fem.Function(V)
problem1 = fem.petsc.LinearProblem(a, L, bcwall, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
phi = problem1.solve()
```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 29, 2022, 10:25am UTC](https://fenicsproject.discourse.group/t/neumann-bc-for-2-conductors/9379/4 "2022-09-29T10:25:28Z")

</div>

`ft` can be written to an XDMFFile as shown in [dolfinx/demo\_gmsh.py at 2be366d20632377e0688df585fcd7e201d5dc86c · FEniCS/dolfinx · GitHub](https://github.com/FEniCS/dolfinx/blob/2be366d20632377e0688df585fcd7e201d5dc86c/python/demo/demo_gmsh.py#L67)

---

<div class="post-metadata">

**Author:** ![Marryweather](https://avatars.discourse-cdn.com/v4/letter/m/b782af/32.png) [@Marryweather](https://fenicsproject.discourse.group/u/Marryweather)\
**Post date:** [October 1, 2022, 10:59am UTC](https://fenicsproject.discourse.group/t/neumann-bc-for-2-conductors/9379/5 "2022-10-01T10:59:45Z")

</div>

I have checked the ft and it seems to be ok.

Do you have any idea why it still doesn’t work?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 3, 2022, 7:31am UTC](https://fenicsproject.discourse.group/t/neumann-bc-for-2-conductors/9379/6 "2022-10-03T07:31:00Z")

</div>

I do not have time to parse though a code of over 200 lines to find a potential typo. I would suggest starting with one of the tutorials, which as verification checks:  
[https://jorgensd.github.io/dolfinx-tutorial/chapter3/robin\_neumann\_dirichlet.html](https://jorgensd.github.io/dolfinx-tutorial/chapter3/robin_neumann_dirichlet.html)  
and do incremental adjustments to the code to get to your problem.

---

<div class="post-metadata">

**Author:** ![Marryweather](https://avatars.discourse-cdn.com/v4/letter/m/b782af/32.png) [@Marryweather](https://fenicsproject.discourse.group/u/Marryweather)\
**Post date:** [October 4, 2022, 10:01am UTC](https://fenicsproject.discourse.group/t/neumann-bc-for-2-conductors/9379/7 "2022-10-04T10:01:15Z")

</div>

I have changed the geometry of the 2d space and now it works with the 2 Neumann boundary conditions.

I just do not know how.

Thank you for the good support!