Integrating over an interior surface

variational formulation
4

Hi, consider using restriction using + and -. Here foo('+') restricts foo to the positive cell and FacetNormal('+') is the outer normal of a positive cell. Following the discussion here the +, - side can be changed by providing a form which has a dx measure term which has specified subdomain_data by some cell function g. Then a positive cell for a facet is the one for which value of g is higher. This allows for the following

from dolfin import *
from mshr import *

omega1 = Rectangle(Point(0, 0), Point(0.5, 2))
omega2 = Rectangle(Point(0.5, 0), Point(1, 2))

domain = omega1 + omega2
domain.set_subdomain(1, omega2)

mesh = generate_mesh(domain, 32)

interior_surface = 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_surface.mark(interior_surface_marker,1)

V = FunctionSpace(mesh, 'CG', 1)
# Make up something with discontinuous gradient
f = interpolate(Expression('std::max(0., x[0]-0.5)', degree=1), V)

dS = Measure('dS', domain=mesh, subdomain_data=interior_surface_marker, subdomain_id=1)

n = FacetNormal(mesh)

# +, - sides are by default assigned arbitrarily so assembling below
# will not yield meaningful result
L0 = dot(grad(f)('+'), n('+'))*dS
print(assemble(L0))

L1 = dot(grad(f)('-'), n('-'))*dS
print(assemble(L1))

mf = MeshFunction('size_t', mesh, 2, mesh.domains())
# Add dummy terms that will allow for controlling restrictions
# This is (0, 0) * (1, 0) * 2
L2 = dot(grad(f)('-'), n('-'))*dS + Constant(0)*dx(domain=mesh, subdomain_data=mf)
print(assemble(L2))

# This is (1, 0) . (-1, 0) * 2
L3 = dot(grad(f)('+'), n('+'))*dS + Constant(0)*dx(domain=mesh, subdomain_data=mf)
print(assemble(L3))
11 Likes