Zero internal Neumann boundary

1

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.
SubDomain2
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 :slight_smile:

Best Regards

Emanuel

2

Based on @MiroK’s answer here, the "+" cell for a facet is the one for which the subdomain MeshFunction has a higher value, so you want to use "+" instead of "-" to integrate flux using the function restricted to subdomain 1. Trying this, int_boundary_val is a small nonzero value that converges to zero as the mesh is refined. (You can only expect exact conservation over \Gamma\subset\partial\Omega if there exists a test function whose trace on \partial\Omega is the characteristic function on \Gamma. With continuous basis functions, this is only possible if \Gamma = \partial\Omega.)

2 Likes
3

Dear kamensky,

thank you for the reply. I actually tried both variants with ‘+’ and ‘-’. I’m always a bit unsure what to use. However good to know, that the general approach is fitting, I’ll read more into using a different function space. Thank you!

Best Regards

Emanuel

4

Hello all, I have a similar problem with dolfinx 0.5.0 and I can’t seem to apply the method above without throwing an AssertionError. I can’t get my form to compile.

The idea is to have 2D Navier-Stockes with a 1D wall inside the domain. I’d like the pressure to verify \partial_yp^+=\partial_yp^-=0 at the wall but allowing for discontinuities through, as in the picture below.

p_jump
p_jump761×834 9.37 KB

The relevant part of my current code is largely taken from dokken’s tutorial and looks like this :

boundaries = [(1, lambda x: np.logical_or(top(x),outlet(x))), (2, symmetry), (3, nozzle)]

facet_indices, facet_markers = [], []
fdim = spy.mesh.topology.dim - 1
for (marker, locator) in boundaries:
	facets = dfx.mesh.locate_entities(mesh, fdim, locator)
	facet_indices.append(facets)
	facet_markers.append(np.full_like(facets, marker))
	
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)

face_tag = dfx.mesh.meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure("ds", domain=mesh, subdomain_data=face_tag)
n = ufl.FacetNormal(mesh)
weak=ufl.inner(p.dx(1)('+')*spy.r,s('+')*n[1])*ds(3)

p is a DG type element, which seemed the simplest way to allow for such sharp jumps.

I have looked around this forum but not yet managed to find anything related to this specific problem - i.e. a Neumann or no-flux condition inside the domain.

5

Ok it seems that I made it work :

fdim = spy.mesh.topology.dim - 1
facet_indices = dfx.mesh.locate_entities(spy.mesh, fdim, nozzle).astype(np.int32)
facet_markers = np.full_like(facet_indices, 1)	.astype(np.int32)
sorted_facets = np.argsort(facet_indices)
	
face_tag = dfx.mesh.meshtags(spy.mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
dS =  ufl.Measure("dS", domain=spy.mesh, subdomain_data=face_tag) # dS for internal surfaces

return ufl.inner((spy.r*p.dx(1))('+'),s('+'))*dS(1)+ufl.inner((spy.r*p.dx(1))('-'),s('-'))*dS(1)

Key things are to use dS, DG elements and to be consistent with the ('+') everywhere.

6

Did you work out what was going on here, I am stuck with a similar problem to you. The neuman condtions seems to work when I use a non-zero flux, but a setting the flux to zero via a boundary term seems to do nothing.

Edit: To add to this, I suspect the answer is that DG elements are needed as suggested by the above reply. But I’m not sure so confirmation would be appreciated, as well as a link to a tutorial of some kind that discusses how to implement this.