Add clear static variables near the end of interpolate_any1

fenicstools/fem/interpolation.cpp

@@ -541,7 +541,7 @@ void interpolate_any1(const Function& u0, Function& u)
 
   // Create vector to hold all local values of u
   std::vector<double> local_u_vector(u.vector()->local_size());
-
+  
   // Get dofmap of u
   dolfin_assert(V1.dofmap());
   const GenericDofMap& dofmap = *V1.dofmap();
@@ -693,7 +693,7 @@ void interpolate_any1(const Function& u0, Function& u)
       // Store received result in map
       for (uint i = 0; i < u0.value_size(); i++)
           values[i] = vals[j*u0.value_size() + i];
-
+      
       coords_to_values.insert(std::make_pair(x, values));
       }
   }
@@ -751,7 +751,9 @@ void interpolate_any1(const Function& u0, Function& u)
               local_u_vector[d] = cell_coefficients[i];
       }
   }
-
+
+  coords_to_values.clear();
+  coords.clear();
   // Set and finalize vector
   u.vector()->set_local(local_u_vector);
   u.vector()->apply("insert");

tests/test_Interpolation.py

@@ -44,6 +44,13 @@ def test_functional2D():
     u1 = interpolate_nonmatching_mesh_any(u0, V1)
     assert round(assemble(dot(u0, u0)*dx) - assemble(dot(u1, u1)*dx), 4) == 0
 
+    # Test with another expression
+    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)
+    u1 = interpolate_nonmatching_mesh_any(u0, V1)
+    assert round(assemble(dot(u0, u0)*dx) - assemble(dot(u1, u1)*dx), 4) == 0
+
 
 def test_functional3D():
     """Test integration of function interpolated in non-matching meshes"""
@@ -76,6 +83,15 @@ def test_functional3D():
     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), 2) == 0
+
+    # Test with another expression
+    f = Expression(("2*(x[0]*x[0] + x[1]*x[1])",
+                    "2*(x[0]*x[0] + x[1]*x[1] + 1)",
+                    "2*(x[0]*x[0] + x[1]*x[1] + 2)"), degree=2)
+    u0 = interpolate(f, V0)
+    u1 = interpolate_nonmatching_mesh_any(u0, V1)
+    assert round(assemble(dot(u0, u0)*dx) - assemble(dot(u1, u1)*dx), 1) == 0
 
 if __name__ == '__main__':
     test_functional2D()
+    #test_functional3D()