How to extract the DOFs of a subdomain?

Question

As a follow up of my recent question, I want to remesh a subdomain and to compute integrals on this new mesh on a subdomain.

If I have a function v on the actual domain, then interpolate(v, V_sub) gives the interpolation of v in the function space V_sub that relates to the mesh of the subdomain.

Since I want to assemble forms, in particular, I need to interpolate some basis functions of V in V_sub.

So the question is, how can I extract the DOFs of the basis functions, that are located within the subdomain?

I have a running example for this based on this answer.

However, this does not work for higher order elements, where the DOFs are also located at the facets. (because of vertex_to_dof_map does not work for these elements)

from dolfin import *
import numpy as np

def extract_dofs_subdomain(V, subd):
    """ extract dofs of subdomain"""

    dofs_of_V = V.dofmap().vertex_to_dof_map(V.mesh())
    coord_of_dofs = V.mesh().coordinates()

    ncords = coord_of_dofs.shape[0]

    subd_bools = np.zeros(V.dim())
    for inds in dofs_of_V:
        subd_bools[inds] = subd.inside(coord_of_dofs[np.mod(inds, ncords)],
                                       False)
    return np.arange(V.dim())[subd_bools.astype(bool)]

class ContDomain(dolfin.SubDomain):
    """define a subdomain"""
    def __init__(self, ddict):
        dolfin.SubDomain.__init__(self)
        self.minxy = [ddict['xmin'], ddict['ymin']]
        self.maxxy = [ddict['xmax'], ddict['ymax']]
    def inside(self, x, on_boundary):
        return (dolfin.between(x[0], (self.minxy[0], self.maxxy[0]))
                and
                dolfin.between(x[1], (self.minxy[1], self.maxxy[1])))

odcoo = dict(xmin=0.45,
             xmax=0.55,
             ymin=0.6,
             ymax=0.8)

mesh = UnitSquareMesh(2, 5)
V = VectorFunctionSpace(mesh, "CG", 1)

odom = ContDomain(odcoo)

ininds = extract_dofs_subdomain(V, odom)
print ininds

How to proceed then? Or is there a better approach to this problem?

Comments

You could loop over all basis functions $\phi_k$ and compute the integral of $\phi_k \chi_{\text{subdomain}}$ to see if $\phi_k$ vanishes on the subdomain or not.

Answer 1

Thanks to Nico, I have found the following solution to identify the DOFs that share their support with the subdomain.

• Test all basis functions against the characteristic function of the subdomain
• Then a nonzero value of the integral identifies basis functions that are nonzero on the subdomain

Nevertheless, there is probably a cleaner solution to this task.

from dolfin import *
import numpy as np

class ContDomain(dolfin.SubDomain):
    """define a subdomain"""
    def __init__(self, ddict):
        dolfin.SubDomain.__init__(self)
        self.minxy = [ddict['xmin'], ddict['ymin']]
        self.maxxy = [ddict['xmax'], ddict['ymax']]
    def inside(self, x, on_boundary):
        return (dolfin.between(x[0], (self.minxy[0], self.maxxy[0]))
                and
                dolfin.between(x[1], (self.minxy[1], self.maxxy[1])))

class CharactFun(dolfin.Expression):
    """ characteristic function of subdomain """
    def __init__(self, subdom):
        self.subdom = subdom
    def eval(self, value, x):
        if self.subdom.inside(x, False):
            value[:] = 1
        else:
            value[:] = 0
    def value_shape(self):
        return (2,)


odcoo = dict(xmin=0.45,
             xmax=0.55,
             ymin=0.6,
             ymax=0.8)

mesh = UnitSquareMesh(32, 32)
V = VectorFunctionSpace(mesh, "CG", 1)

odom = ContDomain(odcoo)

charfun = CharactFun(odom)
v = TestFunction(V)

checkf = assemble(inner(v, charfun) * dx)

dofs_on_subd = np.where(checkf.array() > 0)[0]

basfun = Function(V)
basfun.vector()[dofs_on_subd] = 0.2
basfun.vector()[0] = 1  # for scaling the others only 
plot(basfun)

interactive(True)

Comments

Alternatively you maybe could iterate all cells of the mesh and check which ones interesect with your subdomain. The dofs are then determined by DofMap.cell_dofs(cell). In case that the subdomain is aligned with the mesh, things become pretty simple.
Besides, you could restrict the measure dx to the subdomain. This would make the approach of above at least somewhat more efficient.

---

I agree with your first point!

I refrained from restricting the measure to the subdomain, as it will then be restricted to the cells that are completely inside the subdomain. And in the case that the cells are not aligned this will not capture particular dofs belonging to cells of partial intersection.