Interior Neumann B.C. in BDM/DG mixed Poisson

variational formulation
1

Hello,

I’m trying to solve the mixed Poisson in the subdomain. Firstly, I set up a reference case with submesh created with MeshView():

from dolfin import *
# Create mesh
mesh_all = UnitSquareMesh(MPI.comm_world, 20, 20, 'crossed')

class TopBottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1]*(1-x[1]), 0.0)
    
class Wall(SubDomain):
    def inside(self, x, on_boundary):
        return near((x[0]-0.2)*(0.8-x[0]), 0.0)

class Block(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0.2, 0.8))
                       
subdomains = MeshFunction("size_t",mesh_all, mesh_all.topology().dim())
subdomains.set_all(0)
Block().mark(subdomains, 1)

mesh = Mesh(MeshView.create(subdomains, 1))

boundaries = MeshFunction("size_t",mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
TopBottom().mark(boundaries, 1)
Wall().mark(boundaries, 2)
ds = Measure('ds', domain = mesh, subdomain_data = boundaries)
n = FacetNormal(mesh)               
# Define function spaces and mixed (product) space
BDM = FiniteElement("BDM", mesh.ufl_cell(),1)
DG = FiniteElement("DG", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, BDM * DG)
D = FunctionSpace(mesh, DG)

# Define trial and test functions
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)

# Define source function
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", element=D.ufl_element())

# Define variational form
a = (dot(sigma, tau)+ div(tau)*u + div(sigma)*v)*dx
L = - f*v*dx + Constant(1.0)*dot(tau, n)*ds(2)

# Define function G such that G \cdot n = g
class BoundarySource(UserExpression):
    def __init__(self, mesh, **kwargs):
        self.mesh = mesh
        super().__init__(**kwargs)
    def eval_cell(self, values, x, ufc_cell):
        cell = Cell(self.mesh, ufc_cell.index)
        n = cell.normal(ufc_cell.local_facet)
        g = sin(5*x[0])
        values[0] = g*n[0]
        values[1] = g*n[1]
    def value_shape(self):
        return (2,)

G = BoundarySource(mesh, degree=2)

# Define essential boundary
def boundary(x):
    return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

bc = DirichletBC(W.sub(0), G, boundaries, 1)

# Compute solution
w = Function(W)
solve(a == L, w, bc)
(sigma, u) = w.split(True)

print("u_max =", u.vector().max(), 'u_min =', u.vector().min())
plot(u)

Here the Dirichlet B.C. u=1.0 is impelemented weakly by adding the term Constant(1.0)*dot(tau, n)*ds(2) on the RHS.

The results give
u_max = 1.1637949250619102 u_min = 0.991970223504017.

Then I switch to solve the problem on the subdomain and appyling A.ident_zeros() to the domain out of interest, and my question is how to apply the Neumann B.C. above now. Since the functional u is defined on element-wise constent space DG0, the facet integral after the integrate by parts reads \sum _{e \in \mathcal{F}_i} \int_{e} [[u \tau]]ds, where \mathcal{F}_i is the set of interior facets. I need your help to convert the interior B.C. u=1.0 into the interior facet integral, below is the WME that i tried but failed to match the refernce case above with submesh

from dolfin import *

# Create mesh
mesh = UnitSquareMesh(MPI.comm_world, 20, 20, 'crossed')

class TopBottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1]*(1-x[1]), 0.0)
    
class Wall(SubDomain):
    def inside(self, x, on_boundary):
        return near((x[0]-0.2)*(0.8-x[0]), 0.0)

class Block(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0.2, 0.8))
                       
subdomains = MeshFunction("size_t",mesh, mesh.topology().dim())
subdomains.set_all(0)
Block().mark(subdomains, 1)
dx = Measure('dx', domain = mesh, subdomain_data = subdomains)

boundaries = MeshFunction("size_t",mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
TopBottom().mark(boundaries, 1)
Wall().mark(boundaries, 2)
dS = Measure('dS', domain = mesh, subdomain_data = boundaries)
n = FacetNormal(mesh)               
# Define function spaces and mixed (product) space
BDM = FiniteElement("BDM", mesh.ufl_cell(),1)
DG = FiniteElement("DG", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, BDM * DG)
D = FunctionSpace(mesh, DG)

# Define trial and test functions
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)

# Define source function
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", element=D.ufl_element())

# Define variational form
a = (dot(sigma, tau)+ div(tau)*u + div(sigma)*v)*dx(1)
L = - f*v*dx(1) + Constant(1.0)*jump(tau, n)*dS(2)

# Define function G such that G \cdot n = g
class BoundarySource(UserExpression):
    def __init__(self, mesh, **kwargs):
        self.mesh = mesh
        super().__init__(**kwargs)
    def eval_cell(self, values, x, ufc_cell):
        cell = Cell(self.mesh, ufc_cell.index)
        n = cell.normal(ufc_cell.local_facet)
        g = sin(5*x[0])
        values[0] = g*n[0]
        values[1] = g*n[1]
    def value_shape(self):
        return (2,)

G = BoundarySource(mesh, degree=2)

# Define essential boundary
def boundary(x):
    return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

bc = DirichletBC(W.sub(0), G, boundaries, 1)

A, b = assemble_system(a, L, keep_diagonal=True)
bc.apply(A, b)
A.ident_zeros()

# Compute solution
w = Function(W)
solve(A, w.vector(), b)
(sigma, u) = w.split(True)

print("u_max =", u.vector().max(), 'u_min =', u.vector().min())
plot(u)

with results
u_max = 0.15264490930389638 u_min = -0.005560735505760749.

Any help will be really appreciated!

2

Hello everyone,

unfortunately after a month, I`m still struggling with this problem, anyone has some idea? Many thanks!

3

Looks like you already know about the principle of Mixed function spaces. Make sure you’ve checked against examples at https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/
e.g. https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/mixed-poisson/demo_mixed-poisson.py.rst

I guess you might handle your inner boundary condition with a lagrange multiplier.
Dokken gives more references (published papers) to Cécile Daversin-Catty’s work on mixed spaces at Prescribe discontinuity on internal boundary and also mentions Nitsche’s method as another approach (for discontinuous systems).

1 Like