Set DG function to zero on boundary elements

Question

Currently I use the following code which does what I need:

from dolfin import *

NUM_CELL = 33
mesh = UnitSquareMesh(NUM_CELL,NUM_CELL)
DG0 = FunctionSpace(mesh, "DG", 0)

class ZeroOnBoundary(Expression):
    def eval(self, value, x):
        if x[0] < 0.03:
            value[0] = 0.0
        elif x[0] > 0.97:
            value[0] = 0.0
        elif x[1] < 0.03:
            value[0] = 0.0
        elif x[1] > 0.97:
            value[0] = 0.0
        else:
            value[0] = 1.0
outcome = Function(DG0)
outcome.interpolate(ZeroOnBoundary())

plot(outcome)
interactive()

But as you see, it is not good, because it should be done from the topology of the mesh (without any interpolate()). I tried to use measures (I know how to mark boundary elements and define the measure -> Marking elements), but I was not succesful. I also tried bc.apply() but naturally, it does not help in case of DG spaces.

Comments

Hi, so you want to set all dofs that belong to a cell that has at least one facet on the boundary to zero? Also, do you plan to consider function spaces other than DG0?

---

Hello Miro! Thank you for your reply!
I want to set all dofs of boundary elements to zero. By "boundary elements" I mean elements which have either vertex or facet on the boundary (I will use both possibilities). I will consider only DG spaces of low dimension (DG0 or DG1).
I know already, how to mark elements, but I do not know how to produce what I need from this knowledge.

Answer 1

Hi, here's a possible solution. It should work in 2d and 3d as well in serial and parallel.

from dolfin import *
import numpy as np

domain = Rectangle(0, 0, 1, 1) - Circle(0.5, 0.5, 0.25)
mesh = Mesh(domain, 20)

# Get cell facet connectivity
tdim = mesh.topology().dim()
mesh.init(tdim, tdim-1)

# Mark facets on the boundary
bdry_facets = FacetFunction('bool', mesh, False)
BDRY_FACETS = bdry_facets.array()
DomainBoundary().mark(bdry_facets, True)

# Get all dofs which belong to cells with some facet on the boundary
order = 3
V = FunctionSpace(mesh, 'DG', order)
dofmap = V.dofmap()
first_dof, last_dof = dofmap.ownership_range()

bdry_dofs = np.concatenate([dofmap.cell_dofs(cell.index())
                            for cell in cells(mesh)
                            if any(BDRY_FACETS[cell.entities(tdim - 1)])])

bdry_dofs.sort()

# Consider parallel
bdry_dofs -= first_dof

# Some function to work with
u = Expression('2 + fabs(x[0]-0.5) + fabs(x[1]-0.5)')
u = interpolate(u, V)
U = u.vector()

# Set function's boundary dofs
values = U.get_local()
values[bdry_dofs] = 0
U.set_local(values)

# Check. Project higher order to 0 to make plots same
if order:
    u = project(u, FunctionSpace(mesh, 'DG', 0))
File('u.pvd') << u
plot(u)
interactive() 

If your version of DOLFIN is 1.3.0 then the plot in the parallel run will be different then the
Paraview plot. From the latter, I believe the solution is correct.

Comments

Great, thank you! I tried it and Iooked at fenicsproject pages for some, for me unknown, functions you used, e.g. dofmap.ownership_range() or cell.entities(). It seems really good. Maybe I will ask something here (as a comment) again later.