Neumann BCs when using mesh files

1

Hi,
I am having issues solving a Poisson equation with Neumann BCs using external mesh files.
The code doesn’t show any errors but Neumann BCs seem not to be applied. When solving Poisson equations using a FEnics mesh function (e.g. UnitSquareMesh) the simulation works fine by marking boundaries, instantiating DirichletBC() objects as well as generating expressions for NeumannBCs as described in tutorial 14.
(https://fenicsproject.org/pub/tutorial/html/._ftut1014.html)

When using an external mesh files, Neumann BCs don’t seem to work. As an example, I am considering a geometry (generated in gmsh and converted using dolfin-convert tool) that consists of a rectangle with a hole and 5 boundaries defined as follows; Upper part noted up, down part noted dp, left part noted lp, right part noted rp and the hole noted hp. Constant values are defined for each boundary parts and Neumann BC for the hole.

from dolfin import*

mesh = Mesh("mesh.xml")
subdomains = MeshFunction("size_t", mesh, "mesh_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "mesh_facet_region.xml")
V = FunctionSpace(mesh, 'P', 1)

C_up = Constant("20")
C_dp = Constant("200")
C_lp = Constant("60")
C_rp = Constant("120")
C_hp = Constant("150")

kappa = Constant("2")
f = Constant("0")

boundary_conditions = {0: {'Dirichlet': C_up},   
                       1: {'Dirichlet': C_dp},   
                       2: {'Dirichlet': C_lp}, 
                       3: {'Dirichlet': C_rp},
                       4: {'Neumann': C_hp}} 
                       
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

bcs = []
for i in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[i]:
        bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'],
                            boundaries, i)
        bcs.append(bc)

u = TrialFunction(V)
v = TestFunction(V)

integrals_N = []
for i in boundary_conditions:
    if 'Neumann' in boundary_conditions[i]:
        if boundary_conditions[i]['Neumann'] != 0:
            g = boundary_conditions[i]['Neumann']
            integrals_N.append(g*v*ds(subdomain_id=i, domain=mesh, subdomain_data=boundaries))


a = kappa*inner(grad(u), grad(v))*dx
L = f*v*dx - sum(integrals_N)

u = Function(V)
solve(a == L, u, bcs)

file = File("poisson.pvd")
file << u

Solving the problem with only Dirichlet BCs is working.

boundary_conditions = {0: {'Dirichlet': C_up},   
                       1: {'Dirichlet': C_dp},   
                       2: {'Dirichlet': C_lp}, 
                       3: {'Dirichlet': C_rp},
                       4: {'Dirichlet': C_hp}}

Is there anything missing when applying Neumann BCs?
Many thanks in advance.

2

I cant see anything wrong with your code. However, as i do not have the meshes, i cant really test it.
First suggestion is to not use dolfin-convert, as it is no longer maintained, see: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio
for the appropriate syntax using the latest version of meshio. How to read it into dolfin is found in earlier posts of that topic.
Secondly, check that you are integrering over the correct boundaries using a mesh file.
One way to do that is to visualize the mesh function in paraview. The second is to Calculate the area of each neumann segment, printing assemble(1*ds(subdomain_id=i, domain=mesh, subdomain_data=boundaries))

3

I’m experiencing a similar unclear issue. If I follow these steps my integrals with the measures ds over the “Neumann” boundaries are all just empty. My guess is that they are not initialized properly (as it probably happens in the DirichletBC maybe…)

4

Make sure that the mesh function used to mark boundary facets are correct.

Without a minimal reproducible example I cant Give any further guidance.

5

Thanks for having a look at my issue.

I have made this a separate topic – please check here: