Fix interpolation into vector serendipity (#2316)

* Throw better runtime error

* block size

* block size

* Add test of interpolation in vector elements

* _

* Add test of identity-mapped vector-valued element that isn't a VectorElement

* Add test of vector-valued identity-mapped element that isn't a VectorElement

* Fix implementation into vector-valued identity-mapped element without a block size

* remove debugging

* runtime error

Co-authored-by: Garth N. Wells <gnw20@cam.ac.uk>

cpp/dolfinx/fem/interpolate.h

@@ -476,14 +476,11 @@ void interpolate(Function<T>& u, std::span<const T> f,
 
   // This assumes that any element with an identity interpolation matrix
   // is a point evaluation
-  if (element->interpolation_ident())
+  if (element->map_ident() && element->interpolation_ident())
   {
     // Point evaluation element *and* the geometric map is the identity,
     // e.g. not Piola mapped
 
-    if (!element->map_ident())
-      throw std::runtime_error("Element does not have identity map.");
-
     auto apply_inv_transpose_dof_transformation
         = element->get_dof_transformation_function<T>(true, true, true);
 
@@ -513,8 +510,13 @@ void interpolate(Function<T>& u, std::span<const T> f,
     // Not a point evaluation, but the geometric map is the identity,
     // e.g. not Piola mapped
 
-    if (_f.extent(0) != 1)
-      throw std::runtime_error("Interpolation data has the wrong shape.");
+    const int element_vs = element->value_size() / element_bs;
+
+    if (element_vs > 1 && element_bs > 1)
+    {
+      throw std::runtime_error(
+          "Interpolation into this element not supported.");
+    }
 
     // Get interpolation operator
     const auto [_Pi, pi_shape] = element->interpolation_operator();
@@ -535,9 +537,16 @@ void interpolate(Function<T>& u, std::span<const T> f,
       std::span<const std::int32_t> dofs = dofmap->cell_dofs(cell);
       for (int k = 0; k < element_bs; ++k)
       {
-        std::copy_n(std::next(f.begin(), k * f_shape1 + c * num_interp_points),
-                    num_interp_points, ref_data_b.begin());
-        impl::interpolation_apply(Pi, ref_data, std::span(_coeffs), element_bs);
+        for (int i = 0; i < element_vs; ++i)
+        {
+          std::copy_n(
+              std::next(f.begin(), (i + k) * f_shape1
+                                       + c * num_interp_points / element_vs),
+              num_interp_points / element_vs,
+              std::next(ref_data_b.begin(),
+                        i * num_interp_points / element_vs));
+        }
+        impl::interpolation_apply(Pi, ref_data, std::span(_coeffs), 1);
         apply_inv_transpose_dof_transformation(_coeffs, cell_info, cell, 1);
         for (int i = 0; i < num_scalar_dofs; ++i)
         {