Integral of functions on non-matching meshes

General
1

It is about the implementation of fictitious domain method (FDM), which requires to assemble the following linear form

\int_{\Omega_s} {\lambda} v~d{x}

where {\lambda} \in V_h(\Omega_s) = H^1({\Omega}_s)\cap \mathcal{P}_2 and v\in V_h=H^1(\Omega)\cap \mathcal{P}_2 being finite element functions defined on different meshes.

Here \Omega is a background domain, and {\Omega}_s\subset \Omega is a subdomain suspending in \Omega.

Routine 1 in the following code is a way to assemble the form \int_{\Omega_s} {\lambda}{\color{red} (I_{h}^sv)}~d{x} by using an interpolation I_{h}^s : V_h(\Omega) \to V_h(\Omega_s). It is slow since it interpolates all of the basis functions.

How can I speed up the Routine 1 by using the local supported property of the basis function? Note that in my case, \Omega_s only overlap small portion of elements in the mesh of \Omega.

from dolfin import *
import time
import math

mesh = RectangleMesh(Point(0, 0), Point(1, 1), 4, 4)
meshs = RectangleMesh(Point(.25, .25), Point(.75, .75), 2, 2)

V = VectorFunctionSpace(mesh, "CG", 1)
Lambda = VectorFunctionSpace(meshs, "CG", 1)

# lamda is a given FE function defined on meshs
lmbda = Function(Lambda)
lmbda_ini = Constant((2.0,1.0)) # just for test, lmbda is a finite element function
LIP = LagrangeInterpolator
LIP.interpolate(lmbda, lmbda_ini)


v = TestFunction(V)
#we want to assemble $\int_{meshs} lmbda \cdot v dx$
##  use chi to eliminate the integral outside of meshs 
chi_ =  Expression("x[0]<.75 && x[0]>.25 && x[1]<.75 && x[1]>.25",degree=1)
DG = FunctionSpace(mesh, "DG", 0)
chi = interpolate(chi_, DG)
b_analytic = assemble( inner(lmbda_ini,v)*chi*dx )


# Routine 1: input lmbda and functin space V, output a vector
# slow if mesh becomes bigger
b_computed = Vector(b_analytic)
b_computed[:] = 0.0
for i in range(1, b_computed.size()):
    vv = Function(V)
    vv.vector()[i] = 1.0
    vv_s = Function(Lambda)
    LIP.interpolate(vv_s, vv)
    value = assemble(inner(lmbda, vv_s)*dx)
    b_computed[i] = value
# end of the routine

e = Vector(b_analytic)
e -= b_computed
info("Routine 1: norm of error {} \nnorm of b_analytic {} ".format(norm(e), norm(b_analytic)))