What happens with MeshFunction after mesh.move ?

Question

Hi,

I am interested in accessing a displacement field on a subdomain of a mesh indexed by 1 in a MeshFunction "subdomains" :

mesh = UnitCubeMesh(10, 10, 10)

TOL = 0.00001  
class Crest(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0],0., TOL) and near(x[1], 0., TOL)

crest = Crest()
subdomains = MeshFunction("size_t", mesh, 1)
subdomains.set_all(0)
crest.mark(subdomains, 1)

I move the mesh before launching computations to introduce a small geometrical perturbation :

lp = 0.01
V = VectorFunctionSpace(mesh, "Lagrange", 1)

class Perturbation(Expression):
    def eval(self, value, x):
        value[0] = lp*exp(-(x[0]**2.+ x[1]**2.+(x[2]-0.5)**2.)/lp**2)
        value[1] = -lp*exp(-(x[0]**2.+ x[1]**2.+(x[2]-0.5)**2.)/lp**2)
        value[2] = 0
    def value_shape(self):
        return (3,)

u3D = interpolate(Perturbation(),V)
if dolfin_version()=='1.7.0dev':
    ALE.move(mesh,u3D)
else:
    mesh.move(u3D)

And then I try to access displacement field along the subdomain marked with 1 :

u  = Function(V)
marked_cells = SubsetIterator(subdomains, 1)

for cell in marked_cells:
    x = cell.midpoint().x()
    y = cell.midpoint().y()
    z = cell.midpoint().z()
    print(z)
    print(u(x,y,z))

But I have the impression that the cells in marked_cells are no more inside the domain, I get the following error message :

The point is not inside the domain. 
Consider calling "Function::set_allow_extrapolation(true)" on this Function to allow extrapolation.

So what happens with the MeshFunction "subdomains" after applying mesh.move?

Many thanks in advance,

Claire

Answer 1

Hi, MeshFunction does not hold the mesh entities themselves but only their identifiers/indices. The identifiers, unlike the entities, do not change during deformation. So what you see is due to mesh being modified. Consider

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(4, 4)
f = EdgeFunction('size_t', mesh, 0)
DomainBoundary().mark(f, 1)

# Get ids for 1 marked edges
before = [edge.index() for edge in SubsetIterator(f, 1)]
# and physical coordinates of their midpoints
before_x = [Edge(mesh, edge).midpoint() for edge in before]
# Deform
ALE.move(mesh, Constant((1, 1)))
# The mesh has moved but ids have not changed
after = set(np.where(f.array() == 1)[0])
assert len(before) == len(after)
assert all(edge in after for edge in before)
# Only the physical coordinates of the midpoints did
after_x = [Edge(mesh, edge).midpoint() for edge in before]
print min(a.distance(b) for a, b in zip(after_x, before_x))