Computing surface integral over interior boundary

Question

Hi all,

In the fenics tutorial robin boundary condition are desribed as:

− ∂ u / ∂ n = p ( u − q ) on Γ R

Applying these on an outer boundary seems to be alright (I checked by computing the surface integral of both the left and right hand side of the equation).

However, as can be seen in the simple example below, computing the interior surface integrals does not seem to be correct here. How should I correctly compute these?

from dolfin import *

# Create mesh and define function space

mesh = UnitCubeMesh(25, 25, 25)

V = FunctionSpace(mesh, 'Lagrange', 1)

Define boundary conditions

U0_BOUNDARY_MARKER = 1

U1_BOUNDARY_MARKER = 2

boundary_func = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

boundary_func.set_all(0)

# Dirichlet boundary

class u0_boundary_class(SubDomain):

    def inside(self, x, on_boundary):

        return on_boundary and x[0]     <= DOLFIN_EPS

u0_boundary = u0_boundary_class()

u0_boundary.mark(boundary_func, U0_BOUNDARY_MARKER)

u0 = Expression('1')

bc = DirichletBC(V, u0, boundary_func, U0_BOUNDARY_MARKER)

class u1_boundary_class(SubDomain):

        def inside(self, x, on_boundary):

                return x[0] == 0.52

u1_boundary = u1_boundary_class()

u1_boundary.mark(boundary_func, U1_BOUNDARY_MARKER)

dSs = dS(subdomain_data = boundary_func)

dx = dx(subdomain_data=boundary_func)

Define variational problem

u = TrialFunction(V)

v = TestFunction(V)

f = Constant(0.)        # charge density

p = Constant(1.)        # robin condition

q = Constant(0.75)      # robin condition

a = inner(nabla_grad(u), nabla_grad(v)) dx + pavg(u)avg(v)dSs(U1_BOUNDARY_MARKER)

L = fvdx + pqavg(v)*dSs(U1_BOUNDARY_MARKER)

Compute solution

A = assemble(a)

b = assemble(L)

bc.apply(A, b)

u = Function(V)

solve(A, u.vector(), b, 'lu')

n = FacetNormal(mesh)

uproj = project(grad(u), VectorFunctionSpace(mesh, 'Lagrange', 1))

# check − ∂ u ∂ n = p ( u − q )  on  Γ R

f = project(Constant(1.), V)

surface = assemble(f*dSs(U1_BOUNDARY_MARKER))

# compute p ( u − q )  on  Γ R

print p(0)(assemble(udSs(U1_BOUNDARY_MARKER))/surface-q(0))

 # compute − ∂ u ∂ n on  Γ R

print assemble(dot(uproj,n)('+')*dSs(U1_BOUNDARY_MARKER))

Answer 1

The problem is in your projection step

uproj = project(grad(u), VectorFunctionSpace(mesh, 'Lagrange', 1))

where you erroneously project the discontinuous gradient to a function having a continuous (piecewise linear) basis.
In order to solve this you can either project to a DG space with order 0:

uproj = project(grad(u), VectorFunctionSpace(mesh, 'DG', 0))
print assemble(p*(u-q)*dSs(U1_BOUNDARY_MARKER))
print -assemble(dot(uproj,n)('+')*dSs(U1_BOUNDARY_MARKER)) 

Or simply:

print assemble(p*(u-q)*dSs(U1_BOUNDARY_MARKER))
print -assemble(dot(grad(u),n)('+')*dSs(U1_BOUNDARY_MARKER)) 

Both approaches give the same (and correct!) result.

Comments

I'm curious--with this approach also work with discontinuities across boundaries? I'm wondering if I can project the solution from one subdomain onto a shared interface with another subdomain and use that as a boundary condition. Would that be possible in the FEniCS framework?