From 586e7553b0bec89299b3c60df2cd61bc2070056c Mon Sep 17 00:00:00 2001
From: Gabriele Bozzola <sbozzolator@gmail.com>
Date: Fri, 20 Oct 2023 16:46:52 -0700
Subject: [PATCH] Add support for GPUs in interpolate_array

---
 src/Remapping/interpolate_array.jl | 23 ++++++++++++-----------
 1 file changed, 12 insertions(+), 11 deletions(-)

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