Fix bug in interpolation of vector fields onto Cartesian grids (#617)

* Added a test case for this situation
* Fixed an unnecessarily strict assert in GridBase._vector_to_cartesian

pde/fields/vectorial.py

@@ -606,7 +606,19 @@ def interpolate_to_grid(
             )
 
         # determine the points at which data needs to be calculated
-        if isinstance(grid, CartesianGrid):
+        if (
+            self.grid.num_axes == grid.num_axes
+            and self.grid.__class__ is grid.__class__
+            or (  # UnitGrid and CartesianGrid should be treated the same here
+                isinstance(self.grid, CartesianGrid) and isinstance(grid, CartesianGrid)
+            )
+        ):
+            # convert within the same grid class
+            points = grid.cell_coords
+            # vectors are already given in the correct basis
+            data = self.interpolate(points, bc=bc, fill=fill)
+
+        elif isinstance(grid, CartesianGrid):
             # convert Cartesian coordinates to coordinates in current grid
             points = self.grid.c.pos_from_cart(grid.cell_coords)
             points_grid_sym = self.grid._coords_symmetric(points)
@@ -615,15 +627,6 @@ def interpolate_to_grid(
             # convert the vector to the cartesian basis
             data = self.grid._vector_to_cartesian(points, data_grid)
 
-        elif (
-            self.grid.__class__ is grid.__class__
-            and self.grid.num_axes == grid.num_axes
-        ):
-            # convert within the same grid class
-            points = grid.cell_coords
-            # vectors are already given in the correct basis
-            data = self.interpolate(points, bc=bc, fill=fill)
-
         else:
             # this type of interpolation is not supported
             grid_in = self.grid.__class__.__name__

pde/grids/base.py

@@ -761,7 +761,10 @@ def _vector_to_cartesian(
 
         # convert the basis of the vectors to Cartesian
         rot_mat = self.c.basis_rotation(points)
-        assert rot_mat.shape == (self.dim, self.dim) + shape
+        assert (
+            rot_mat.shape == (self.dim, self.dim)
+            or rot_mat.shape == (self.dim, self.dim) + shape
+        )
         return np.einsum("j...,ji...->i...", components, rot_mat)  # type: ignore
 
     def normalize_point(

tests/fields/test_vectorial_fields.py

@@ -300,6 +300,25 @@ def test_vector_bcs():
     np.testing.assert_allclose(s1, s2)
 
 
+def test_interpolation_vector_fields_cartesian():
+    """Test interpolation of a vector field on Cartesian coordinates."""
+    # upscale
+    grid = UnitGrid([10, 10])
+    vf = VectorField.from_expression(grid, ["x", "y**2"])
+    grid2 = CartesianGrid([[0, 10], [0, 10]], 20)
+    vf2 = vf.interpolate_to_grid(grid2)
+    vf_expect = VectorField.from_expression(grid2, ["x", "y**2"])
+    np.testing.assert_allclose(vf2.data[1:-1, 1:-1], vf_expect.data[1:-1, 1:-1])
+
+    # downscale
+    grid = CartesianGrid([[0, 10], [0, 10]], 20)
+    vf = VectorField.from_expression(grid, ["x", "y**2"])
+    grid2 = UnitGrid([10, 10])
+    vf2 = vf.interpolate_to_grid(grid2)
+    vf_expect = VectorField.from_expression(grid2, ["x", "y**2"])
+    np.testing.assert_allclose(vf2.data[1:-1, 1:-1], vf_expect.data[1:-1, 1:-1])
+
+
 def test_interpolation_vector_fields_cylindrical():
     """Test interpolation of a vector field on cylindrical coordinates."""
     grid = CylindricalSymGrid(5, [-2, 3], 10)