problem marking interior facets

[Chaffra Affouda]

The following example proposed by @blechta is not marking interior facets as expected. It's also marking exterior facets to zero when it should not. What  am I doing wrong here? Can anyone reproduce? Only 2 corners remain at 1 where the whole boundary should stay at 1. Any suggestion is appreciated.

from dolfin import *
import dolfin
import pylab as pb

mesh = UnitSquareMesh(10, 10)
plot(mesh)

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

contact_facets = MeshFunction('size_t', mesh, 1, 0)

for f in facets(mesh):
    if any(ff.exterior() for ff in facets(f)):
        contact_facets[f] = 1

anti_bc = DirichletBC(V, Constant(0.0), contact_facets, 0)

f_contact_facets = Function(V)
f_contact_facets.vector()[:] = 1.0

anti_bc.apply(f_contact_facets.vector());

pb.figure()
plt = plot(f_contact_facets, mode='warp')
pb.colorbar(plt)

[Jan Blechta]

There is a bug, should be

anti_bc = DirichletBC(V, Constant(0.0), contact_facets, 1)

[Jan Blechta]

Also note that behaviour of facets(f) might have changed. Any d-d connectivities are now just identity, i.e., facet 6 is connected just to facet 6. You might want d-0-d or d-tdim-d - choose one of those and implement that along the lines of

ff.exterior() for ff in facets(v) for v in vertices(f)
ff.exterior() for ff in facets(c) for v in cells(f)

[Chaffra Affouda]

I am a little confused. I am trying to mark facets with marker 0 (interior) so it should be anti_bc = DirichletBC(V, Constant(0.0), contact_facets, 0), correct? I will look into your other suggestion.

[Chaffra Affouda]

Thanks @blechta I had to write the marking as

for f in facets(mesh):
    for v in vertices(f):
        if any(ff.exterior() for ff in facets(v)):
            contact_facets[f] = 1