Hey,

I have some issues with a zero Neumann boundary inside my domain. Basically I have two variables u_0 and u_1. u_0 is defined on the whole domain, u_1 is only defined on a subdomain (the yellow one). There should be no flow of u_1 in between the subdomains.

I use the following code:

```
from __future__ import print_function
from mshr import *
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
tol = 1E-14
mesh = UnitSquareMesh(64, 64)
omega1 = Rectangle(Point(0, 0), Point(0.5, 1))
omega2 = Rectangle(Point(0.5, 0), Point(1, 1))
domain = omega1 + omega2
domain.set_subdomain(1, omega2)
mesh = generate_mesh(domain, 32)
# Create Subdomain Omega_0
subdomain1 = CompiledSubDomain('x[0] <= 0.5 + DOLFIN_EPS')
interior_boundary = CompiledSubDomain('near(x[0], 0.5, DOLFIN_EPS)')
interior_surface_marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
interior_surface_marker.set_all(0)
interior_boundary.mark(interior_surface_marker,1)
materials = MeshFunction("size_t", mesh, 2)
materials.set_all(0)
subdomain1.mark(materials, 1)
# Introduce FunctionSpaces for different variables
V0 = FiniteElement('CG', mesh.ufl_cell(), 1)
V1 = FiniteElement('CG', mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, MixedElement([V0, V1]))
## Init Trial and Test Function
u0, u1 = TrialFunctions(V)
v0, v1 = TestFunctions(V)
# dx and ds for the integrals
dx = Measure("dx", domain=mesh, subdomain_data=materials)
dS_1 = Measure('dS', domain=mesh, subdomain_data=interior_surface_marker, subdomain_id=1)
# Boundary conditions
def boundary_Left(x, on_boundary):
return on_boundary and near(x[0], 0., tol)
bc_0 = DirichletBC(V.sub(0), Constant(1.), boundary_Left)
bc_1 = DirichletBC(V.sub(1), Constant(1.), boundary_Left)
bcs = [bc_0, bc_1]
### Set up the variational problem
n1 = FacetNormal(V.sub(1).mesh())
a_0 = dot(grad(u0), grad(v0)) * dx(tuple([0,1]))
a_1 = dot(grad(u1), grad(v1)) * dx(1) - Constant(0.)('-')*dot(grad(u1)('-'), n1('-'))*v1('-')*dS_1
a = a_0 + a_1
L_0 = Constant(1.) * v0 * dx(tuple([0,1]))
L_1 = Constant(1.) * v1 * dx(1)
L = L_0 + L_1
u = Function(V)
A = assemble(a, keep_diagonal=True)
b = assemble(L)
for bc in bcs:
bc.apply(A,b)
A.ident_zeros()
solver = KrylovSolver('gmres', 'ilu')
solver.solve(A, u.vector(), b)
u0_sol, u1_sol = split(u)
# value of the flux on the internal boundary
int_boundary_val = assemble( dot(grad(u1_sol)('-'), n1('-'))*dS_1 + Constant(0.)*dx )
print('Internal Neumann Boundary: ' + str(int_boundary_val))
```

When I evaluate the flux at the bottom, it is non-zero.

I suspect it has something to do with my FunctionSpace (as I use Lagrange for u1), however I don’t get the example working when using ‘DG’.

Any help is appreciated

Best Regards

Emanuel