Integration on an internal boundary dS

Question

Hello,

I have a 2D domain divided into 2 subdomains (a rectangle with a smaller circle in the middle).
I want to compute an interior facet integral over the interface between the two domains, i.e. the circle.
Having marked the edges of that circle with label 4, what I want to do is:

int_{interface} v*dS(4)

where v is the TestFunction of the FunctionSpace defined on the mesh.

I get an error doing the assemble of (v*dS(4)) :

Calling FFC just-in-time (JIT) compiler, this may take some time.
Discontinuous type Argument must be restricted.
Traceback (most recent call last):
File "test_forum.py", line 86, in
s0 = assemble(s0)\
File "/usr/local/lib/python2.7/dist-packages/dolfin/fem/assembling.py", line 192, in assemble
dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
File "/usr/local/lib/python2.7/dist-packages/dolfin/fem/assembling.py", line 66, in _create_dolfin_form
function_spaces=function_spaces)
File "/usr/local/lib/python2.7/dist-packages/dolfin/fem/form.py", line 94, in init
mpi_comm=mesh.mpi_comm())
File "/usr/local/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 65, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/local/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 124, in jit
result = ffc.jit(ufl_object, parameters=p)
File "/usr/local/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 201, in jit
module = jit_build_with_instant(ufl_object, module_name, parameters)
File "/usr/local/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 98, in jit_build_with_instant
code_h, code_c = jit_generate(ufl_object, module_name, parameters)
File "/usr/local/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 62, in jit_generate
parameters=parameters, jit=True)
File "/usr/local/lib/python2.7/dist-packages/ffc/compiler.py", line 154, in compile_form
analysis = analyze_forms(forms, parameters)
File "/usr/local/lib/python2.7/dist-packages/ffc/analysis.py", line 60, in analyze_forms
parameters) for form in forms)
File "/usr/local/lib/python2.7/dist-packages/ffc/analysis.py", line 60, in
parameters) for form in forms)
File "/usr/local/lib/python2.7/dist-packages/ffc/analysis.py", line 142, in _analyze_form
form_data = compute_form_data(form)
File "/usr/local/lib/python2.7/dist-packages/ufl/algorithms/compute_form_data.py", line 285, in compute_form_data
form = apply_restrictions(form)
File "/usr/local/lib/python2.7/dist-packages/ufl/algorithms/apply_restrictions.py", line 201, in apply_restrictions
return map_integrand_dags(rules, expression, only_integral_type=integral_types)
File "/usr/local/lib/python2.7/dist-packages/ufl/algorithms/map_integrands.py", line 60, in map_integrand_dags
return map_integrands(lambda expr: map_expr_dag(function, expr, compress), form, only_integral_type)
File "/usr/local/lib/python2.7/dist-packages/ufl/algorithms/map_integrands.py", line 39, in map_integrands
mapped_integrals = [map_integrands(function, itg, only_integral_type) for itg in form.integrals()]
File "/usr/local/lib/python2.7/dist-packages/ufl/algorithms/map_integrands.py", line 46, in map_integrands
return itg.reconstruct(function(itg.integrand()))
File "/usr/local/lib/python2.7/dist-packages/ufl/algorithms/map_integrands.py", line 60, in
return map_integrands(lambda expr: map_expr_dag(function, expr, compress), form, only_integral_type)
File "/usr/local/lib/python2.7/dist-packages/ufl/corealg/map_dag.py", line 35, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress)
File "/usr/local/lib/python2.7/dist-packages/ufl/corealg/map_dag.py", line 82, in map_expr_dags
r = handlersv._ufl_typecode_
File "/usr/local/lib/python2.7/dist-packages/ufl/algorithms/apply_restrictions.py", line 58, in _require_restriction
ufl_assert(self.current_restriction is not None, "Discontinuous type %s must be restricted." % o._ufl_class_.__name__)
File "/usr/local/lib/python2.7/dist-packages/ufl/assertions.py", line 45, in ufl_assert
error(*message)File "/usr/local/lib/python2.7/dist-packages/ufl/log.py", line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Discontinuous type Argument must be restricted.

This doesn't happen for integrals over the external boundary, i.e. int_bottom v*ds(1).

Has anyone come across this problem before?

---

Here is the relevant code:

from dolfin import *
from mshr import *
from math import pi

# -------------------- #
#     Create  mesh     #
# -------------------- #
rectangle = Rectangle(Point(-1., -1.), Point(1., 1.))
R = 0.5
origin = Point(0.,0.)
circle = Circle(origin, R, segments=32)
domain = rectangle
domain.set_subdomain(1, circle)
domain.set_subdomain(2, rectangle - circle)
mesh = generate_mesh(domain, 15) 

# Create subdomains
subdomains = MeshFunction("size_t", mesh, 2, mesh.domains())


# -------------------- #
#   Create boundaries  #
# -------------------- #
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[1] + 1.) < DOLFIN_EPS

boundaries = FacetFunction("size_t", mesh, 0 ) 
Bottom().mark(boundaries, 1)


# -------------------- #
# Add interface marker #
# -------------------- #
for f in facets(mesh): # For all edges in the mesh
    p0 = Vertex(mesh, f.entities(0)[0]) # save the two ending points p0 and p1
    p1 = Vertex(mesh, f.entities(0)[1])
    if near(p0.x(0)*p0.x(0) + p0.x(1)*p0.x(1), R*R, 0.05) and near(p1.x(0)*p1.x(0) + p1.x(1)*p1.x(1), R*R, 0.05): # check if the points lie on the circle - if yes, put a label on this     edge
        boundaries[f] = 4


# --------------------------------- #
# Check that everything is correct  #
# --------------------------------- #
# Define a function of ones on the mesh:
V = FunctionSpace(mesh, "CG", 1)
f = Function(V)
f.vector()[:] = 1.

dS = dS(subdomain_data=boundaries)
form4 = f*dS(4) # integration on the whole facet marked with number 4
ans = assemble(form4)   # compute the length of the circle
print ans, "vs.", 2*pi*R # dS on the inteface = 2*pi*R


# -------------------------------------------- #
# Mesh / Function Space / Measure / Integrals  #
# -------------------------------------------- #
# 1. Save and load mesh/subdomains/boundaries 
File("physical_region.xml") << subdomains
File("facet_region.xml") << boundaries

mesh = Mesh(mesh)
#subd = MeshFunction("size_t", mesh, 2)
#bound = MeshFunction("size_t", mesh, boundaries) # why is this not working?
subd = MeshFunction("size_t", mesh, "physical_region.xml")
bound = MeshFunction("size_t", mesh, "facet_region.xml")


# 2. Create Finite Element space (Lagrange P1) + Dirichlet BC 
V = FunctionSpace(mesh, "Lagrange", 1)

# Define measures:
dx = Measure("dx")(subdomain_data=subd)
ds = Measure("ds")(subdomain_data=bound)
dS = Measure("dS")(subdomain_data=bound)

# Define trial/test functions
u = TrialFunction(V)
v = TestFunction(V)

# Integral on the external boundary
f0 = v*ds(1) + 1e-15*v*dx
F0 = assemble(f0)\

# Integral on the internal boundary (interface)
s0 = v*dS(4) + 1e-15*v*dx
s0 = assemble(s0)\

Answer 1

When integrating over facets, you have to specify on which side of the facet you want to integrate ('restricted integration') as follows (respectively for the '+' side, the '-' side and both sides):

s0 = v('+')*dS(4) + 1e-15*v*dx
s0 = v('-')*dS(4) + 1e-15*v*dx
s0 = v('+')*dS(4) + v('-')*dS(4) + 1e-15*v*dx

Comments

OK Thank you!!

---

Hi, I am also doing some surface integral over internal boundary. I wonder is there any reason why you include 1e-15vdx at the end? Does it help on removing some of the instability along the boundary?