diff --git a/src/Remapping/interpolate_array.jl b/src/Remapping/interpolate_array.jl
index 6a2467b30f..dfe1f924d6 100644
--- a/src/Remapping/interpolate_array.jl
+++ b/src/Remapping/interpolate_array.jl
@@ -4,15 +4,15 @@ function interpolate_slab(
     (I1,)::Tuple{<:AbstractArray},
 )
     space = axes(field)
-    x = zero(eltype(field))
     QS = Spaces.quadrature_style(space)
     Nq = Spaces.Quadratures.degrees_of_freedom(QS)
 
-    for i in 1:Nq
-        ij = CartesianIndex((i,))
-        x += I1[i] * Operators.get_node(space, field, ij, slabidx)
-    end
-    return x
+    nodes = [
+        Operators.get_node(space, field, CartesianIndex((i,)), slabidx) for
+        i in 1:Nq
+    ]
+
+    return I1 ⋅ nodes
 end
 
 function interpolate_slab(
@@ -21,14 +21,15 @@ function interpolate_slab(
     (I1, I2)::Tuple{<:AbstractArray, <:AbstractArray},
 )
     space = axes(field)
-    x = zero(eltype(field))
     QS = Spaces.quadrature_style(space)
     Nq = Spaces.Quadratures.degrees_of_freedom(QS)
 
-    for j in 1:Nq, i in 1:Nq
-        ij = CartesianIndex((i, j))
-        x += I1[i] * I2[j] * Operators.get_node(space, field, ij, slabidx)
-    end
+    nodes = [
+        Operators.get_node(space, field, CartesianIndex((i, j)), slabidx)
+        for i in 1:Nq, j in 1:Nq
+    ]
+
+    x = (I1 * nodes) ⋅ I2
     return x
 end