Repeated call to interpolate_nonmatching_mesh_any seems not to work

[sma-user]

Hi!

The code below should produce two arrays of all 0s and two arrays of all 1s, but I get
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]

The code is
from dolfin import *
from fenicstools import interpolate_nonmatching_mesh_any
import numpy

mesh = UnitSquareMesh(2, 2)

V = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 1)

u0=Function(V)
u0.vector().set_local( u0.vector().get_local()*0.0 )
u0.vector().apply('insert')
print(u0.vector().get_local()) # Should be all 0
i0 = interpolate_nonmatching_mesh_any(u0, V)
print(i0.vector().get_local()) # Should be all 0

u1=Function(V)
u1.vector().set_local( u1.vector().get_local()*0.0+1.0 )
u1.vector().apply('insert')
print(u1.vector().get_local()) # Should be all 1
i1 = interpolate_nonmatching_mesh_any(u1, V)
print(i1.vector().get_local()) # Should be all 1

Is this a bug or am I doing something wrong?

Best Regards,
Søren

PS This function is a really nice and useful addition to fenics. Thanks for sharing it!

[mikaem]

Hi Søren,

I haven't really touched this function in quite a few years so it's a bit difficult for me to say what's wrong. Are you saying it works if you apply interpolate_nonmatching_mesh_any once, but not twice in the same namespace? That could easily be a bug, but it does not seem to be captured by the tests (that are also calling the same function several times). You understand that the function is supposed to be used to interpolate from one space to another different space, right? So there's actually no test for interpolating on the same space.

[sma-user]

Yes, when I call it twice in a row on the same functions but with modified values, it simply spits out the values from the first time it was called. It does that also when using two different function spaces (u0 \in V1 and i0 \in V2), but I just kept the code above as simple as possible.

In test_Interpolation.py it seems also to be called only once pr. dimension pr. function-type (1 time for an expression and 1 time for a function), so maybe this is not really tested?

If I change 'test_functional2D' to use the function twice on two different expressions, i.e.
def test_functional2D():
"""Test integration of function interpolated in non-matching meshes"""

f = Quadratic2D(degree=2)

# Interpolate quadratic function on course mesh
mesh0 = UnitSquareMesh(8, 8)
V0 = FunctionSpace(mesh0, "Lagrange", 2)
u0 = interpolate_nonmatching_mesh(f, V0)

# Interpolate FE function on finer mesh
mesh1 = UnitSquareMesh(31, 31)
V1 = FunctionSpace(mesh1, "Lagrange", 2)
u1 = interpolate_nonmatching_mesh(u0, V1)
assert round(assemble(u0*dx) - assemble(u1*dx), 10) == 0

mesh1 = UnitSquareMesh(30, 30)
V1 = FunctionSpace(mesh1, "Lagrange", 2)
u1 = interpolate_nonmatching_mesh(u0, V1)
assert round(assemble(u0*dx) - assemble(u1*dx), 10) == 0

f = Expression(("x[0]*x[0] + x[1]*x[1]",
                "x[0]*x[0] + x[1]*x[1] + 1"), degree=2)
V0 = FunctionSpace(mesh0, "Nedelec 1st kind H(curl)", 2)
u0 = interpolate(f, V0)

# Interpolate FE function on finer mesh
V1 = FunctionSpace(mesh1, "Nedelec 1st kind H(curl)", 2)
u1 = interpolate_nonmatching_mesh_any(u0, V1)
assert round(assemble(dot(u0, u0)*dx) - assemble(dot(u1, u1)*dx), 4) == 0

f = Expression(("2*(x[0]*x[0] + x[1]*x[1])",
                "2*(x[0]*x[0] + x[1]*x[1] + 1)"), degree=2)
u0 = interpolate_nonmatching_mesh_any(f, V0)

# Interpolate FE function on finer mesh
u1 = interpolate_nonmatching_mesh_any(u0, V1)
assert round(assemble(dot(u0, u0)*dx) - assemble(dot(u1, u1)*dx), 4) == 0

then I get an assertion error when running 'test_Interpolation.py':
Traceback (most recent call last):
File "test_Interpolation.py", line 89, in
test_functional2D()
File "test_Interpolation.py", line 53, in test_functional2D
assert round(assemble(dot(u0, u0)*dx) - assemble(dot(u1, u1)*dx), 4) == 0
AssertionError

[mikaem]

Ok, so then it's highly likely a bug. Thanks for reporting and feel free to dig deeper:-) I'd be happy to accept a PR that fixes the issue:-)

[sma-user]

It seems to me that the fault is due to the map coords_to_values is being declared as static, and the next time the method is called it already contains all the old values from the first call. Adding

coords_to_values.clear();

near the end of interpolate_any1 (in interpolation.cpp) seems to fix the problem:
....
// Place result in local vector
for (uint i = 0; i < cell_dofs.size(); i++)
{
uint d = cell_dofs[i];
if (d < local_size)
local_u_vector[d] = cell_coefficients[i];
}
}

coords_to_values.clear();
// Set and finalize vector
u.vector()->set_local(local_u_vector);
u.vector()->apply("insert");
}

Thanks again for this great functionality. It should IMO be included directly in fenics :)

[sma-user]

Maybe there should also be a

coords.clear();

at the same place, since coords is also declared static.

[mikaem]

Great work, will you submit a PR?

[sma-user]

Yes, just did. Let me know if it works.