Restriction of DG space and measures to nodal patch

Question

I am having some issues with restricting a DG space and the assembly to a nodal patch. In the following example based on the dolfin DG Poisson demo, the patch of some arbitrary node (here with index 12) is determined. Subsequently, the inner and outer (patch boundary) facets are determined. With this information, the required measures are restricted and the DG Poisson problem is assembled on the restricted DG space. However, the solution does neither seem to be correct nor to be restricted to the selected patch.

Cheers, Martin

from dolfin import *
import itertools as iter

# Create mesh and define function space
mesh = UnitSquareMesh(4, 4)
mesh.init()

# determine vertex patch data
selected_vertex = Vertex(mesh,12)
cf = CellFunction('size_t', mesh)
patch_cells = set(selected_vertex.entities(2).astype(int).tolist())
for cid in patch_cells:
    cf[cid] = 1
patch_facets = [Cell(mesh,cid).entities(1).astype(int).tolist() for cid in patch_cells]
patch_facets = set(iter.chain(*patch_facets))
ff_inner = FacetFunction('size_t', mesh)
ff_outer = FacetFunction('size_t', mesh)
for fid in patch_facets:
    facet_cells = set(Facet(mesh,fid).entities(2).astype(int).tolist())
    if len(facet_cells)==1 or len(facet_cells-patch_cells)==1:
        ff_outer[fid] = 1
    else:
        ff_inner[fid] = 1
plot(cf)
plot(ff_inner, title='inner')
plot(ff_outer, title='outer', interactive=True)

# construct restricted function space
restriction = Restriction(cf,1)
V = FunctionSpace(restriction, "DG", 1)

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

# Define normal component, mesh size and right-hand side
n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2
f = Expression("500.0*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")

# Define parameters
alpha = 4.0
gamma = 8.0

# Define restricted measures
dx = Measure('dx')[cf]
ds = Measure('ds')[ff_outer]
dS = Measure('dS')[ff_inner]

# Define variational problem
a = dot(grad(v), grad(u))*dx(1) \
   - dot(avg(grad(v)), jump(u, n))*dS(1) \
   - dot(jump(v, n), avg(grad(u)))*dS(1) \
   + alpha/h_avg*dot(jump(v, n), jump(u, n))*dS(1) \
   - dot(grad(v), u*n)*ds(1) \
   - dot(v*n, grad(u))*ds(1) \
   + gamma/h*v*u*ds(1)
L = v*f*dx(1)

# Compute solution
u = Function(V)
solve(a == L, u)

# Plot solution
plot(u, interactive=True)

Answer 1

ds = Measure('ds')[ff_outer]

this obviously creates a measure of exterior (w. r. t. mesh) facets taking values ff_outter==1 which is an empty set. If you switch to

ds = Measure('dS')[ff_outer]

then it is not obvious how to restrict integrand. Could you file an issue so that it can be thought over and designed better?

Comments

So what you are saying is that

ds = Measure('ds')[ff_outer]

requires the marked facets in ff_outer to coincide with some exterior facets of mesh, otherwise ds(1) integrals never get evaluated? In that case, it should be possible to just use dS and define the exterior facet integrals of the patch as dS(2) (replacing ds(1)), right?
If not, do you have any suggestion how to discretise the described patch problem?

---

In that case, it should be possible to just use dS and define the exterior facet integrals of the patch as dS(2) (replacing ds(1)), right?

Yeah, I just said it. The problem is that integrand must be restricted and it is not obvious how.

If not, do you have any suggestion how to discretise the described patch problem?

I suggested that you file an issue.

---

The problem is that integrand must be restricted and it is not obvious how.

Could you elaborate? I don't get why the integrand is not restricted although it should only be defined on marked facets. In fact, I get the error message

ufl.log.UFLException: Form argument must be restricted.

which seems to support your point.

I suggested that you file an issue.

Sure, can do. My question was: Is there a way around the problem which is readily available?

---

Could you elaborate?

While integrating say foo over some interior facet, you need to say on which cell adjacent to the facet you take the integrand foo. You can do:

dS = Measure("dS")[facet_function]
foo('+') * dS(1)
foo('-') * dS(1)
jump(foo) * dS(1)
avg(foo) * dS(1)
...

but not

foo * dS(1)

Is there a way around the problem which is readily available?

avg may be a way around provided that integrand is defined on cells outside of restriction which probably does not hold. The issue itself is then characterized as following: exterior facet integral (as currently working) is not capable to integrate on actual exterior facets of restricted domain. Your example code is very suitable to demonstrate the issue.