identify and assemble on boundary elements

Question

1) What would be a good way to identify all cells which have at least one vertex on the (exterior) boundary of a domain?
2) Is it possible to solve a problem on this restricted set of cells? I know that assembly works for marked cells by employing the respective measures. Is it also possible to solve on a restricted set of dofs or would I have to extract the sub-matrix (or set the remaining dofs to identity)?

Thanks, Martin

Answer 1

1)

#mesh.init(j, k) # this may be needed for some j, k
boundary_adjacent_cells = [c for c in cells(mesh) 
                             if any(f.exterior() for f in facets(v)
                                                 for v in vertices(c))]

2)

cell_domains = CellFunction('size_t', mesh)
#cell_domains.set_all(0) # needed in older versions of DOLFIN
for c in boundary_adjacent_cells:
    cell_domains[c] = 42

# for example
foo = assemble(qux*dx(42), cell_domains=cell_domains)

# or any other DOLFIN function doing assembly,
# taking cell_domains argument
# e.g. solve(), FooSolver.solve(), ...

# or we can define measure
dxb = Measure('dx')[cell_domains]

If you're asking how to avoid restricted problem being singular, one possibility is to eliminate null space using DirichletBC which is little tricky. There is also probably more efficient and cleaner Restriction functionality but is not polished yet - see restriction demo.

Comments

Thanks, this seems promising. Will also have a look at the restriction demo. Apart from that, how would I determine the dofs belonging to the restricted problem?

---

Currently DOLFIN user is being encouraged to avoid direct indexing of system tensors. But if you want to determine a correspondence of DOF indices with mesh entities you can use helper functions in FunctionSpace.dofmap(). It depends on a particular space used.

Also DirichletBC approach for regularizing restricted problem depends on a particular set-up and you would better ask in a separate, specific question.

On the other hand with Restriction approach you don't need list of restricted DOFs. Restricted problem is defined on a an already restricted function space lacking redundant DOFs. This is why this should be a preferred approach (at least in a future).

---

A follow-up re #1: How can one identify all cells that are adjacent to an interior boundary defined by some SubDomain? Facet.exterior() unfortunately won't work in that case.

---

Just play with FacetFunctions, CellFunctions and entity iterators (cells(mesh), cells(cell), facets(cell), ...). For example

class MySubdomain (SubDomain):
    ....
subdomain = MySubdomain()

ff = FacetFunction('size_t', mesh)
# ff.set_all(0) # in older DOLFIN
subdomain.mark(ff, 42)

cf = CellFunction('size_t', mesh)
# cf.set_all(0) # in older DOLFIN
for f in facets(mesh):
    if ff[f]==42:
        for c in cells(f):
            cf[c] = 999

Just note that SubDomain.mark() method marks entity iff SubDomain.inside() returns true all its vertices plus midpoint. It is source of many confusions...

---

ok, will do, thanks.

Answer 2

for the first one, I would extract the facet and then point on the exterior boundary using the BoundaryDomain. Then with the connectivity, find out the cells you want.
For the second one, if you have the boundary condition for the restricted domains, then you only need to mark the vertices and do not need to extract the sub-matrix to solve the restrict set of dofs, which will have exact solution as you solve the problem after extracting and solving the sub-problem.