extension of restricted function to full space

Question

Assume $u_r$ to be a function on some restricted space $V_r \subset V$. What would be a good way to extend $u_r$ by zero to $u \in V$? Interpolation does not work and I have not found a way to associate restricted dofs with global dofs (which would be a useful information to have in any case).
Here is an example to play with.

from dolfin import *

# create mesh and select patch of some cell -> C
mesh = UnitSquareMesh(4,4)
arbitrary_selected_index = 3
arbitrary_selected_cell = Cell(mesh, arbitrary_selected_index)
C = [c.index() for c in cells(arbitrary_selected_cell)]
C.append(arbitrary_selected_cell.index())

cf = CellFunction('size_t', mesh)
cf.set_all(1)
for i in C:
    cf[i] = 0

# create patch measure
dx = dx[cf]

# define restricted space
restriction = Restriction(cf, 0)
Vr = FunctionSpace(restriction, "CG", 1)

# define functions
u = TrialFunction(Vr)
v = TestFunction(Vr)
f = Constant(1.0)

# define forms
a = dot(grad(u),grad(v))*dx()
L = f*v*dx()

# define boundary condition
ff = FacetFunction("size_t", mesh, 0)
mesh.init()
# set value of ff to 1 on boundary of restriction
for cell in cells(mesh):
    if cf[cell.index()] == 0:
        for facet in facets(cell):
            if facet.num_entities(2) == 1: # exterior boundary
                ff[facet.index()] = 1
            else:
                c0, c1 = facet.entities(2)
                if cf[int(c0)] != cf[int(c1)]:   # interior boundary
                    ff[facet.index()] = 1

zero = Constant(0.0)
bc = DirichletBC(Vr, zero, ff, 1)

# solve system
ur = Function(Vr)
solve(a == L, ur, bcs=bc)

# transfer ur to V --- but HOW?
V = FunctionSpace(mesh, "CG", 1)
u = interpolate(ur, V)

# plot cell patch and solution
plot(cf, title="patch")
plot(ur, title="u restricted")
plot(u, title="u extended")
interactive()

Comments

This seems to me as a bug. Evaluation of ur seems to be correct - ur is automatically extended by zero outside of the restriction:

>>> for c in cells(mesh):
...   if cf[c] == 1:
...     for v in vertices(c):
...       print ur(v.point())
... 
(zeros or small numbers here)

Function::compute_vertex_values seems to work well also

>>> ur.compute_vertex_values(mesh)
array([ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.02083333,  0.02083333,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ])

You should file a bug probably. Alternatively we could ask developers for making FunctionAssigner work also for restrictions.

On the other hand I'm also going to solve for a quantity defined on a whole mesh given a plenty of local problems on patches. My idea is rather to assemble small local problems into one big and then take advantage of specialized routines in PETSc for solving block-diagonal problems. This should be scalable and also take advantage of a locality of a problem on each patch.

---

Ok, I see. I can file a bug report for this. However, I'd in particular be interested in a work-around. The local patch problems I am looking at are in some DG space. Thus, a simple evaluation at the vertices and mapping to global dofs won't do the trick. Are there any other suggestions or is it somehow possible to get the dof correspondence w.r.t. the full space?

The other approach you mention probably is more efficient. Could you elaborate on how you setup and assemble the global matrix?

---

Could you elaborate on how you setup and assemble the global matrix?

I was thinking that I could assemble a global problem by just summing up local ones. But now I see that in my case this does not make a resulting problem block-diagonal. So now I'm just thinking about patch-wise assembling another big matrix having diagonal blocks given by local restricted problems. Then there is some mapping between this auxiliary problem and a global problem which is probably represented just by (rectangular) matrix multiplication. But this will require some C++ coding and will be tricky in parallel.

---

In any case you would also need a mapping between restricted and global dofs. Is there currently no way to obtain this information?
I also have tried to reset the dofs that should be zero explicitly. Unfortunately, this does not work reliably since even on the patch the interpolation is incorrect.

# "postprocess" after creation of "interpolated" u
V_dofmap = V.dofmap()
V_dofs = dict([(cid, V_dofmap.cell_dofs(cid)) for cid in C])
V_dofs = [v.tolist() for v in V_dofs.values()]
import itertools as iter
V_dofs = set(iter.chain(*V_dofs))
zero_idx = set(range(V.dim())).difference(V_dofs)
u.vector()[list(zero_idx)] = 0

---

Please, file a bug report, so that relevant developer can take care of this issue.

---

ok, done now.

Answer 1

Since FEniCS has no problem interpolating from $V$ to $V_r$ (the other way around from what is desired), we can construct the inverse map as follows:

Start with the function in $V$ whose value for each degree of freedom is actually the index for that degree of freedom. That is, the function has the vector of dof values [0,1,2,3,..]. Then interpolate it into the smaller space $V_r$. After rounding to integers, the vector of values for the interpolated function in the smaller space will be the corresponding indices from the larger space.

Here is a minimal example:

import numpy as np
from dolfin import *

def make_submesh_local_to_global_map(V, Vr):
    # V is function space on whole mesh, 
    # Vr is restriction of V to submesh
    # u.vector()[local_to_global_map] = ur.vector()
    u = Function(V)
    u.vector()[:] = np.array(range(V.dim()))
    ur = interpolate(u, Vr)
    local_to_global_map = np.round(ur.vector().array()).astype(int)
    return local_to_global_map

# Test it out
mesh = UnitSquareMesh(30,30)

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

class My_domain(SubDomain):
    def inside(self,x,on_boundary):
        return x[0]<=0.5

subdomain = CellFunction('size_t',mesh)
subdomain.set_all(0)
my_domain = My_domain()
my_domain.mark(subdomain,1 )

submesh = SubMesh(mesh, subdomain, 1)
Vr = FunctionSpace(submesh,'Lagrange',1)

local_to_global_map = make_submesh_local_to_global_map(V, Vr)

wr = interpolate(Expression('sin(pi*x[0])*x[1]'), Vr)

w = Function(V)
w.vector()[local_to_global_map] = wr.vector().array()

plot(wr, title='restricted expression')
plot(w, title='prolongation')
interactive()`

Comments

Thanks for the suggestion of this work-around. I had the impression that the recent or planned developments in FEniCS would allow for a nicer way but haven't checked lately.

---

Maybe there is a more direct way, I don't know. If there is, google searches were not of much help finding it (that's how I ended up at this question).