diff --git a/CUDA_STRATEGY.md b/CUDA_STRATEGY.md
index fbfe99639..ba8fb9d3d 100644
--- a/CUDA_STRATEGY.md
+++ b/CUDA_STRATEGY.md
@@ -12,7 +12,7 @@ The work is split into small PRs that can be reviewed and merged one at a time.
 - [x] **PR 4: `Pusher` accepts `Kernel`** — the kernel for the active backend is chosen once, when the pusher is created; a plain `PyccelKernel` is wrapped, so the propagators do not change (no behaviour change on CPU).
 - [x] **PR 5: `Domain` on the GPU** — domain arguments are selected and stored at domain construction; CUDA arguments reference device arrays, and deepcopy/unpickling rebuilds the arguments from the copied or restored arrays.
 - [x] **PR 6: `Particles` on the GPU** — `Particles` can be created on the CuPy backend, and `Particles.args_markers` is selected as the CUDA or Pyccel argument bundle at construction.
-- [ ] **PR 7: `Derham` on the GPU** — `Derham` can be created on the CuPy backend, plus `Derham.cuda_args_derham`.
+- [x] **PR 7: `Derham` on the GPU** — `Derham` can be created on the CuPy backend, plus `Derham.cuda_args_derham`.
 - [ ] **PR 8: Shared CUDA headers for the argument classes** — one `.cuh` per argument class instead of long flat kernel signatures.
 - [ ] **PR 9: One folder per kernel, starting with `pic/pushing`** — pure refactor, no behaviour change.
 - [ ] **PR 10: Device versions of helper kernels** — B-spline evaluation, mapping evaluation (per domain), small linear algebra, as `__device__` functions in `.cuh` headers.
@@ -41,12 +41,12 @@ CUDA kernels can be added one by one. If the code runs on the GPU and needs a ke
 - **No silent CPU fallback on the GPU.** A kernel without a CUDA version raises an error on the GPU backend. Falling back would mean copying data to the host and back at every call.
 - **Small steps.** Every PR keeps the CPU code path working and tested.
 
-## Current state (PR 6)
+## Current state (PR 7)
 
 | File | Content |
 |---|---|
 | `src/struphy/utils/kernel_backends.py` | `is_cuda_backend()`, `CudaKernel` (wraps a `cupy.RawKernel`, compiled lazily; expands `Argument.get_cuda_args()` and takes `n_threads`), `Kernel` and `KernelCatalog` for backend selection and discovery |
-| `src/struphy/utils/cuda_arguments.py` | `Argument` contract plus `CudaMarkerArguments` and `CudaDomainArguments`; CUDA arrays are stored individually and returned in signature order by `get_cuda_args()` |
+| `src/struphy/utils/cuda_arguments.py` | `Argument` contract plus `CudaMarkerArguments`, `CudaDerhamArguments` and `CudaDomainArguments`; CUDA arrays are stored individually and returned in signature order by `get_cuda_args()` |
 | `src/struphy/geometry/base.py` | `Domain.args_domain` is selected once at construction; CUDA domains use device arrays, while direct Pyccel geometry calls retain a host argument bundle |
 | `src/struphy/pic/base.py` | `Particles` arrays and `args_markers` use the backend selected at construction; direct Pyccel methods retain a private host bundle |
 | `src/struphy/pic/tests/test_kernel_backends.py` | the demo kernel pair `push_eta_linear` (pyccel function compiled with `epyccel` at test time, CUDA source string) and tests on both backends |
@@ -126,7 +126,7 @@ kernel = catalog["push_eta_stage"]  # Kernel: pyccel or CUDA depending on the ba
 - The CUDA argument objects hold references. If an owner reallocates an array (today the markers are allocated once), it must rebuild its CUDA arguments at the same place, exactly like for the pyccel arguments.
 - First these classes must be creatable on the CuPy backend at all:
   - `Particles`: wrap Python lists in `xp.array` before reductions, and use host buffers for scalar MPI gathers (`pic/base.py`).
-  - `Derham`: NumPy arrays from feectools reach `cupy.ascontiguousarray`.
+  - `Derham`: feectools and struphy moved host data (knots, grids, collocation matrices) to CuPy before calling pyccel kernels (see PR 7 below).
   - `Domain`: deepcopy and unpickling on CuPy failed (see PR 5 below).
 
 ### PR 5: `Domain` on the GPU (complete)
@@ -142,7 +142,19 @@ kernel = catalog["push_eta_stage"]  # Kernel: pyccel or CUDA depending on the ba
 - Particle arrays, validity masks, and boundary-condition codes are allocated through `cunumpy`, so they live on CuPy when the CuPy backend is active.
 - `args_markers` is built as `CudaMarkerArguments` from device arrays on CuPy, or as `MarkerArguments` on NumPy. A private host bundle remains for direct Pyccel calls.
 - Domain decomposition now wraps the Python `nprocs` list with `xp.array` before calling `xp.prod`. Scalar MPI gathers use small NumPy buffers and copy the results back to the active array backend, avoiding unsupported CuPy buffers in MPI calls.
-- Full GPU particle pushing still depends on CUDA versions of the required kernels and on PR 7's `Derham` support.
+- Full GPU particle pushing still depends on CUDA versions of the required kernels.
+
+### PR 7: `Derham` on the GPU (complete)
+
+- feectools and the struphy code that builds `Derham` followed `xp` everywhere. So on CuPy, the knots, quadrature grids and decomposition metadata became device arrays and then reached pyccel kernels, SciPy or MPI, which only take host arrays.
+  - feectools: struphy-hub/feectools#85, the first of the feectools CUDA PRs, makes feectools run on the CuPy backend. It must be merged, and the submodule or the feectools version bumped, before `Derham` can be created on CuPy.
+  - struphy: data that describes the spline spaces is host data on every backend. Only the coefficients (`StencilVector` data) and stencil matrices live on the device. The projection and quadrature grids of `Derham` (`get_pts_and_wts`, ...) and `spline_types_pyccel` are NumPy. `domain_array`, `index_array(_N/_D)` and `neighbours` are gathered with NumPy MPI buffers and then converted with `xp.asarray`, so they are device arrays on CuPy like `Particles.domain_array`.
+- `Derham.args_derham` is built from the host knots, degrees and starts on both backends. `Derham.cuda_args_derham` lazily builds `CudaDerhamArguments` with one device copy of these small arrays. On NumPy it raises (host arrays are never copied to the device). The pyccel scratch arrays (`bn1`, ..., `bd3`) are not part of it; they become per-thread local arrays in CUDA (PR 10).
+- Not supported on CuPy yet:
+  - Local projectors (`DerhamOptions.local_projectors=True`) raise `NotImplementedError` when the `Derham` is created. `CommutingProjectorLocal` builds its data with `xp` and calls pyccel kernels on it, like `Derham` did.
+  - Polar splines need a spline mapping, which cannot be created on CuPy yet (see PR 5).
+  - Field evaluation (`SplineFunction.__call__`, ...) still calls pyccel kernels with the coefficients, which are device arrays on CuPy. It needs CUDA evaluation kernels (PR 10+).
+- Tests: `feec/tests/test_derham_gpu.py`. Without a GPU, a strict host stand-in for CuPy (rejects host/device mixing, cannot run kernels; not part of the repository) was used, on top of struphy-hub/feectools#85. With it, a `Derham` created on the "CuPy" backend matches the NumPy one on 1, 2 and 4 MPI processes.
 
 ### PR 8: Argument structs in shared headers
 
diff --git a/src/struphy/feec/psydac_derham.py b/src/struphy/feec/psydac_derham.py
index e73f4cdf2..e16af24c5 100644
--- a/src/struphy/feec/psydac_derham.py
+++ b/src/struphy/feec/psydac_derham.py
@@ -44,6 +44,8 @@
 from struphy.polar.extraction_operators import PolarExtractionBlocksC1
 from struphy.polar.linear_operators import PolarExtractionOperator, PolarLinearOperator
 from struphy.topology.grids import TensorProductGrid
+from struphy.utils.cuda_arguments import CudaDerhamArguments
+from struphy.utils.kernel_backends import is_cuda_backend
 
 NonTrivialBC = LiteralOptions.OptsNonTrivialBoundaryCondition
 space_to_form = {
@@ -57,14 +59,6 @@
 logger = logging.getLogger("struphy")
 
 
-def _to_numpy_for_kernel(value):
-    """Convert CuPy arrays to NumPy for compiled kernel calls."""
-    if hasattr(value, "get"):
-        # This is a CuPy array
-        return value.get()
-    return value
-
-
 class DiscreteDerham:
     """Discrete 3D de Rham sequence built from four FE spaces.
 
@@ -423,7 +417,8 @@ def __init__(
             self._quad_grid_spans[-1] = tuple(self.quad_grid_spans[-1])
             self._quad_grid_bases[-1] = tuple(self.quad_grid_bases[-1])
 
-            self._spline_types_pyccel[-1] = xp.array(
+            # inputs of pyccel kernels, kept on the host
+            self._spline_types_pyccel[-1] = np.array(
                 self._spline_types_pyccel[-1],
             )
 
@@ -485,7 +480,7 @@ def spline_types(self) -> tuple[tuple[str]]:
         return self._spline_types
 
     @property
-    def spline_types_pyccel(self) -> tuple[xp.ndarray]:
+    def spline_types_pyccel(self) -> tuple[np.ndarray]:
         """Tuple of spline types in each direction as integers (0 for 'B', 1 for 'M') for each component of the vector space."""
         return self._spline_types_pyccel
 
@@ -608,6 +603,10 @@ def __init__(
         polar_splines = options.polar_splines
         # local commuting projectors
         local_projectors = options.local_projectors
+        if local_projectors and is_cuda_backend():
+            raise NotImplementedError(
+                "Local projectors (DerhamOptions.local_projectors=True) are not supported on the CuPy backend yet."
+            )
 
         # number of elements and spline degrees in each direction
         assert len(num_elements) == 3
@@ -899,14 +898,20 @@ def __init__(
 
         self._neighbours = self._get_neighbours()
 
-        # collect arguments for kernels
-        self._args_derham = DerhamArguments(
-            _to_numpy_for_kernel(xp.array(self.degree)),
-            _to_numpy_for_kernel(self.V0fem.knots[0]),
-            _to_numpy_for_kernel(self.V0fem.knots[1]),
-            _to_numpy_for_kernel(self.V0fem.knots[2]),
-            _to_numpy_for_kernel(xp.array(self.V0.starts)),
+        # collect arguments for kernels (the knots of feectools are host arrays on every array backend)
+        self._pyccel_args_derham = DerhamArguments(
+            np.array(self.degree),
+            *self.V0fem.knots,
+            np.array(self.V0.starts),
         )
+        if is_cuda_backend():
+            self._args_derham = CudaDerhamArguments(
+                xp.asarray(self._pyccel_args_derham.pn),
+                *(xp.asarray(t) for t in self.V0fem.knots),
+                xp.asarray(self._pyccel_args_derham.starts),
+            )
+        else:
+            self._args_derham = self._pyccel_args_derham
 
         logger.debug("\nDERHAM:")
         logger.debug(f"{'number of elements:'.ljust(25)} {num_elements}")
@@ -1477,8 +1482,8 @@ def div_bcfree(self):
         return self._div_bcfree
 
     @property
-    def args_derham(self):
-        """Collection of mandatory arguments for pusher kernels."""
+    def args_derham(self) -> DerhamArguments | CudaDerhamArguments:
+        """Mandatory pusher kernel arguments for the backend used at initialization."""
         return self._args_derham
 
     # --------------------------
@@ -1778,7 +1783,7 @@ def _discretize_space(
         )
 
         # Create uniform grid
-        grids = [xp.linspace(xmin, xmax, num=ne + 1) for xmin, xmax, ne in zip(min_coords, max_coords, ncells)]
+        grids = [np.linspace(xmin, xmax, num=ne + 1) for xmin, xmax, ne in zip(min_coords, max_coords, ncells)]
 
         # Create 1D finite element spaces and precompute quadrature data
         spaces_1d = [
@@ -1925,11 +1930,11 @@ def _get_domain_array(self):
         else:
             nproc = 1
 
-        # send buffer
-        dom_arr_loc = xp.zeros(9, dtype=float)
+        # send buffer (MPI buffers are host arrays on every array backend)
+        dom_arr_loc = np.zeros(9, dtype=float)
 
         # main array (receive buffers)
-        dom_arr = xp.zeros(nproc * 9, dtype=float)
+        dom_arr = np.zeros(nproc * 9, dtype=float)
 
         # Get global starts and ends of domain decomposition
         gl_s = self.domain_decomposition.starts
@@ -1947,7 +1952,7 @@ def _get_domain_array(self):
         else:
             dom_arr[:] = dom_arr_loc
 
-        return dom_arr.reshape(nproc, 9)
+        return xp.asarray(dom_arr.reshape(nproc, 9))
 
     def _get_index_array(self, decomposition):
         """
@@ -1972,11 +1977,11 @@ def _get_index_array(self, decomposition):
         else:
             nproc = 1
 
-        # send buffer
-        ind_arr_loc = xp.zeros(6, dtype=int)
+        # send buffer (MPI buffers are host arrays on every array backend)
+        ind_arr_loc = np.zeros(6, dtype=int)
 
         # main array (receive buffers)
-        ind_arr = xp.zeros(nproc * 6, dtype=int)
+        ind_arr = np.zeros(nproc * 6, dtype=int)
 
         # Get global starts and ends of cart OR domain decomposition
         gl_s = decomposition.starts
@@ -1993,7 +1998,7 @@ def _get_index_array(self, decomposition):
         else:
             ind_arr[:] = ind_arr_loc
 
-        return ind_arr.reshape(nproc, 6)
+        return xp.asarray(ind_arr.reshape(nproc, 6))
 
     def _get_neighbours(self):
         """
@@ -2025,7 +2030,7 @@ def _get_neighbours(self):
             neighbours along the edges only have one 1, neighbours along the edges have no 1 in the index.
         """
 
-        neighs = xp.empty((3, 3, 3), dtype=int)
+        neighs = np.empty((3, 3, 3), dtype=int)
 
         for i in range(3):
             for j in range(3):
@@ -2034,7 +2039,7 @@ def _get_neighbours(self):
                     ind = tuple(comp)
                     neighs[ind] = self._get_neighbour_one_component(comp)
 
-        return neighs
+        return xp.asarray(neighs)
 
     def _get_neighbour_one_component(self, comp):
         """
@@ -2070,8 +2075,9 @@ def _get_neighbour_one_component(self, comp):
         if comp == [1, 1, 1]:
             return neigh_id
 
-        comp = xp.array(comp)
-        kinds = xp.array(kinds)
+        # computed on the host: the start/end indices below are an object array (with None entries)
+        comp = np.array(comp)
+        kinds = np.array(kinds)
 
         # if only one process: check if comp is neighbour in non-peridic directions, if this is not the case then return the rank as neighbour id
         if size == 1:
@@ -2084,12 +2090,13 @@ def _get_neighbour_one_component(self, comp):
             # elements with index 2n are the starts and 2n + 1 are the ends.
 
             neigh_inds = [None] * 6
+            index_array = xp.to_numpy(self.index_array)
 
             # in each direction find start/end index for neighbour
             for k, co in enumerate(comp):
                 if co == 1:
-                    neigh_inds[2 * k + 0] = self.index_array[rank, 2 * k + 0]
-                    neigh_inds[2 * k + 1] = self.index_array[rank, 2 * k + 1]
+                    neigh_inds[2 * k + 0] = index_array[rank, 2 * k + 0]
+                    neigh_inds[2 * k + 1] = index_array[rank, 2 * k + 1]
 
                 elif co == 0:
                     neigh_inds[2 * k + 1] = gl_s[k] - 1
@@ -2106,15 +2113,15 @@ def _get_neighbour_one_component(self, comp):
                         "Wrong value for component; must be 0 or 1 or 2 !",
                     )
 
-            neigh_inds = xp.array(neigh_inds)
+            neigh_inds = np.array(neigh_inds)
 
             # only use indices where information is present to find the neighbours rank
-            inds = xp.where(xp.not_equal(neigh_inds, None))
+            inds = np.where(np.not_equal(neigh_inds, None))
 
             # find ranks (row index of domain_array) which agree in start/end indices
-            index_temp = xp.squeeze(self.index_array[:, inds])
-            unique_ranks = xp.where(
-                xp.equal(index_temp, neigh_inds[inds]).all(1),
+            index_temp = np.squeeze(index_array[:, inds])
+            unique_ranks = np.where(
+                np.equal(index_temp, neigh_inds[inds]).all(1),
             )[0]
 
             # if any row satisfies condition, return its index (=rank of neighbour)
@@ -3493,18 +3500,16 @@ def get_pts_and_wts(space_1d, start, end, n_quad=None, polar_shift=False):
     histopol_loc = space_1d.histopolation_grid[start : end + 2].copy()
 
     # make sure that greville points used for interpolation are in [0, 1]
-    # Use numpy for comparison since greville points are NumPy arrays
-    greville_loc_np = greville_loc.get() if hasattr(greville_loc, "get") else greville_loc
-    assert np.all(np.logical_and(greville_loc_np >= 0.0, greville_loc_np <= 1.0))
+    assert np.all(np.logical_and(greville_loc >= 0.0, greville_loc <= 1.0))
 
     # interpolation
     if space_1d.basis == "B":
         x_grid = greville_loc
         pts = greville_loc[:, None]
-        wts = xp.ones(pts.shape, dtype=float)
+        wts = np.ones(pts.shape, dtype=float)
 
         # sub-interval index is always 0 for interpolation.
-        subs = xp.zeros(pts.shape[0], dtype=int)
+        subs = np.zeros(pts.shape[0], dtype=int)
 
         # !! shift away first interpolation point in eta_1 direction for polar domains !!
         if pts[0] == 0.0 and polar_shift:
@@ -3518,32 +3523,32 @@ def get_pts_and_wts(space_1d, start, end, n_quad=None, polar_shift=False):
             union_breaks = space_1d.breaks[:-1]
 
         # Make union of Greville and break points
-        # tmp = set(xp.round(space_1d.histopolation_grid, decimals=14)).union(
-        #     xp.round(union_breaks, decimals=14),
+        # tmp = set(np.round(space_1d.histopolation_grid, decimals=14)).union(
+        #     np.round(union_breaks, decimals=14),
         # )
         # tmp = list(tmp)
         # tmp.sort()
-        # tmp_a = xp.array(tmp)
+        # tmp_a = np.array(tmp)
 
-        tmp = set(xp.round(space_1d.histopolation_grid, decimals=14).tolist()).union(
-            xp.round(union_breaks, decimals=14).tolist()
+        tmp = set(np.round(space_1d.histopolation_grid, decimals=14).tolist()).union(
+            np.round(union_breaks, decimals=14).tolist()
         )
         tmp = sorted(tmp)
-        tmp_a = xp.array(tmp)
+        tmp_a = np.array(tmp)
 
         x_grid = tmp_a[
-            xp.logical_and(
+            np.logical_and(
                 tmp_a
-                >= xp.min(
+                >= np.min(
                     histopol_loc,
                 )
                 - 1e-14,
-                tmp_a <= xp.max(histopol_loc) + 1e-14,
+                tmp_a <= np.max(histopol_loc) + 1e-14,
             )
         ]
 
         # determine subinterval index (= 0 or 1):
-        subs = xp.zeros(x_grid[:-1].size, dtype=int)
+        subs = np.zeros(x_grid[:-1].size, dtype=int)
         for n, x_h in enumerate(x_grid[:-1]):
             add = 1
             for x_g in histopol_loc:
@@ -3558,12 +3563,6 @@ def get_pts_and_wts(space_1d, start, end, n_quad=None, polar_shift=False):
 
         pts_loc, wts_loc = np.polynomial.legendre.leggauss(n_quad)
 
-        if "cupy" in xp.__name__:
-            import cupy as cp
-
-            pts_loc = cp.array(pts_loc)
-            wts_loc = cp.array(wts_loc)
-
         x, wts = bsp.quadrature_grid(x_grid, pts_loc, wts_loc)
 
         pts = x % 1.0
@@ -3575,7 +3574,7 @@ def get_pts_and_wts_quasi(
     space_1d: SplineSpace,
     *,
     polar_shift: bool = False,
-) -> tuple[xp.ndarray, xp.ndarray]:
+) -> tuple[np.ndarray, np.ndarray]:
     r"""Obtain local projection point sets and weights in one grid direction for the quasi-interpolation method.
     The quasi-interpolation points are :math:`\nu - \mu +p` equidistant points :math:`\{ x^i_j \}_{0 \leq j < \nu - \mu +p}` in the sub-interval :math:`Q = [\eta_\mu , \eta_\nu]` given by:
 
@@ -3623,12 +3622,12 @@ def get_pts_and_wts_quasi(
         # interpolation
         if space_1d.basis == "B":
             if degree == 1 and h != 1.0:
-                x_grid = xp.linspace(-(degree - 1) * h, 1.0 - h + (h / 2.0), (N + degree - 1) * 2)
+                x_grid = np.linspace(-(degree - 1) * h, 1.0 - h + (h / 2.0), (N + degree - 1) * 2)
             else:
-                x_grid = xp.linspace(-(degree - 1) * h, 1.0 - h, (N + degree - 1) * 2 - 1)
+                x_grid = np.linspace(-(degree - 1) * h, 1.0 - h, (N + degree - 1) * 2 - 1)
 
             pts = x_grid[:, None] % 1.0
-            wts = xp.ones(pts.shape, dtype=float)
+            wts = np.ones(pts.shape, dtype=float)
 
             # !! shift away first interpolation point in eta_1 direction for polar domains !!
             if pts[0] == 0.0 and polar_shift:
@@ -3639,16 +3638,16 @@ def get_pts_and_wts_quasi(
             # The computation of histopolation points breaks in case we have num_elements=1 and periodic boundary conditions since we end up with only one x_grid point.
             # We need to build the histopolation points by hand in this scenario.
             if degree == 0 and h == 1.0:
-                x_grid = xp.array([0.0, 0.5, 1.0])
+                x_grid = np.array([0.0, 0.5, 1.0])
             elif degree == 0 and h != 1.0:
-                x_grid = xp.linspace(-degree * h, 1.0 - h + (h / 2.0), (N + degree) * 2)
+                x_grid = np.linspace(-degree * h, 1.0 - h + (h / 2.0), (N + degree) * 2)
             else:
-                x_grid = xp.linspace(-degree * h, 1.0 - h, (N + degree) * 2 - 1)
+                x_grid = np.linspace(-degree * h, 1.0 - h, (N + degree) * 2 - 1)
 
             n_quad = degree + 1
             # Gauss - Legendre quadrature points and weights
             # products of basis functions are integrated exactly
-            pts_loc, wts_loc = xp.polynomial.legendre.leggauss(n_quad)
+            pts_loc, wts_loc = np.polynomial.legendre.leggauss(n_quad)
 
             x, wts = bsp.quadrature_grid(x_grid, pts_loc, wts_loc)
             pts = x % 1.0
@@ -3662,26 +3661,26 @@ def get_pts_and_wts_quasi(
             N_b = N + degree
 
             # Filling the quasi-interpolation points for i=0 and i=1 (since they are equal)
-            x_grid = xp.linspace(0.0, knots[degree + 1], degree + 1)
-            x_aux = xp.linspace(0.0, knots[degree + 1], degree + 1)
-            x_grid = xp.append(x_grid, x_aux)
+            x_grid = np.linspace(0.0, knots[degree + 1], degree + 1)
+            x_aux = np.linspace(0.0, knots[degree + 1], degree + 1)
+            x_grid = np.append(x_grid, x_aux)
             # Now we append those for 1<i<degree-1
             for i in range(2, degree - 1):
-                x_aux = xp.linspace(knots[degree], knots[degree + i], degree + i)
-                x_grid = xp.append(x_grid, x_aux)
+                x_aux = np.linspace(knots[degree], knots[degree + i], degree + i)
+                x_grid = np.append(x_grid, x_aux)
 
             # Now we append the points for degree-1<= i <= N_b-degree
-            x_aux = xp.linspace(0.0, 1.0, 2 * N + 1)
-            x_grid = xp.append(x_grid, x_aux)
+            x_aux = np.linspace(0.0, 1.0, 2 * N + 1)
+            x_grid = np.append(x_grid, x_aux)
 
             # Now the points for N_b-degree < i < N_b-1
             for i in range(N_b - degree + 1, N_b - 1):
-                x_aux = xp.linspace(knots[i + 1], knots[N_b], N_b + degree - i - 1)
-                x_grid = xp.append(x_grid, x_aux)
+                x_aux = np.linspace(knots[i + 1], knots[N_b], N_b + degree - i - 1)
+                x_grid = np.append(x_grid, x_aux)
             # Finally we add the pointset for i = N_b-1, which is the same as the one for i = N_b-2
             i = N_b - 2
-            x_aux = xp.linspace(knots[i + 1], knots[N_b], N_b + degree - i - 1)
-            x_grid = xp.append(x_grid, x_aux)
+            x_aux = np.linspace(knots[i + 1], knots[N_b], N_b + degree - i - 1)
+            x_grid = np.append(x_grid, x_aux)
 
             if polar_shift:
                 for i in range(len(x_grid)):
@@ -3689,7 +3688,7 @@ def get_pts_and_wts_quasi(
                         x_grid[i] += 0.00001
 
             pts = x_grid[:, None]
-            wts = xp.ones(pts.shape, dtype=float)
+            wts = np.ones(pts.shape, dtype=float)
 
         # histopolation
         elif space_1d.basis == "M":
@@ -3706,31 +3705,31 @@ def get_pts_and_wts_quasi(
             # Thus, we must substract 1 to all the indices of the knots here to refere to the same point.
 
             # Filling the quasi-interpolation points for i=0 and i=1 (since they are equal)
-            x_grid = xp.linspace(0.0, knots[degree], degree + 1)
-            x_aux = xp.linspace(0.0, knots[degree], degree + 1)
-            x_grid = xp.append(x_grid, x_aux)
+            x_grid = np.linspace(0.0, knots[degree], degree + 1)
+            x_aux = np.linspace(0.0, knots[degree], degree + 1)
+            x_grid = np.append(x_grid, x_aux)
             # Now we append those for 1<i<degree-1
             for i in range(2, degree - 1):
-                x_aux = xp.linspace(knots[degree - 1], knots[degree + i - 1], degree + i)
-                x_grid = xp.append(x_grid, x_aux)
+                x_aux = np.linspace(knots[degree - 1], knots[degree + i - 1], degree + i)
+                x_grid = np.append(x_grid, x_aux)
 
             # Now we append the points for degree-1<= i <= N_b-degree
-            x_aux = xp.linspace(0.0, 1.0, 2 * N + 1)
-            x_grid = xp.append(x_grid, x_aux)
+            x_aux = np.linspace(0.0, 1.0, 2 * N + 1)
+            x_grid = np.append(x_grid, x_aux)
 
             # Now the points for N_b-degree < i < N_b-1
             for i in range(N_b - degree + 1, N_b - 1):
-                x_aux = xp.linspace(knots[i], knots[N_b - 1], N_b + degree - i - 1)
-                x_grid = xp.append(x_grid, x_aux)
+                x_aux = np.linspace(knots[i], knots[N_b - 1], N_b + degree - i - 1)
+                x_grid = np.append(x_grid, x_aux)
             # Finally we add the pointset for i = N_b-1, which is the same as the one for i = N_b-2
             i = N_b - 2
-            x_aux = xp.linspace(knots[i], knots[N_b - 1], N_b + degree - i - 1)
-            x_grid = xp.append(x_grid, x_aux)
+            x_aux = np.linspace(knots[i], knots[N_b - 1], N_b + degree - i - 1)
+            x_grid = np.append(x_grid, x_aux)
 
             # Gauss - Legendre quadrature points and weights
             # products of basis functions are integrated exactly
             n_quad = degree
-            pts_loc, wts_loc = xp.polynomial.legendre.leggauss(n_quad)
+            pts_loc, wts_loc = np.polynomial.legendre.leggauss(n_quad)
 
             x, wts = bsp.quadrature_grid(x_grid, pts_loc, wts_loc)
             pts = x
@@ -3743,7 +3742,7 @@ def get_span_and_basis(pts, space):
 
     Parameters
     ----------
-    pts : xp.array
+    pts : np.array
         2d array of points (ii, iq) = (interval, quadrature point).
 
     space : SplineSpace
@@ -3751,10 +3750,10 @@ def get_span_and_basis(pts, space):
 
     Returns
     -------
-    span : xp.array
+    span : np.array
         2d array indexed by (n, nq), where n is the interval and nq is the quadrature point in the interval.
 
-    basis : xp.array
+    basis : np.array
         3d array of values of basis functions indexed by (n, nq, basis function).
     """
 
@@ -3762,8 +3761,8 @@ def get_span_and_basis(pts, space):
     T = space.knots
     degree = space.degree
 
-    span = xp.zeros(pts.shape, dtype=int)
-    basis = xp.zeros((*pts.shape, degree + 1), dtype=float)
+    span = np.zeros(pts.shape, dtype=int)
+    basis = np.zeros((*pts.shape, degree + 1), dtype=float)
 
     for n in range(pts.shape[0]):
         for nq in range(pts.shape[1]):
@@ -3789,7 +3788,7 @@ def get_weights_local_projector(pts, fem_space):
 
     Parameters
     ----------
-    pts : xp.array
+    pts : np.array
         3d array of points. Contains the quasi-interpolation points in each direction.
 
     fem_space : SplineSpace
@@ -3797,10 +3796,10 @@ def get_weights_local_projector(pts, fem_space):
 
     Returns
     -------
-    wij : List of xp.array
+    wij : List of np.array
         List of 2d array indexed by (space_direction, i, j), where i determines for which FEEC coefficient this weights are needed. Used for interpolation.
 
-    whij : List of xp.array
+    whij : List of np.array
         List of 2d array indexed by (space_direction, i, j), where i determines for which FEEC coefficient this weights are needed. Used for histopolation.
     """
     wij = []
@@ -3816,7 +3815,7 @@ def get_weights_local_projector(pts, fem_space):
     #######
 
     # List with the degree of the B-splines in each spatial direction
-    plist = xp.zeros(3, dtype=int)
+    plist = np.zeros(3, dtype=int)
     # List with a bool that tell us if the B-splines in each spatial direction are periodic
     periodiclist = []
     # We iterate over each one of the spatial dimension of the 0 fem_space
@@ -3824,9 +3823,9 @@ def get_weights_local_projector(pts, fem_space):
         plist[d] = space.degree
         periodiclist.append(space.periodic)
 
-    periodiclist = xp.array(periodiclist)
+    periodiclist = np.array(periodiclist)
     # We get the maximum number of j entries for wij
-    lenj1, lenj2, lenj3 = get_local_problem_size(periodiclist, plist, xp.array([False, False, False], dtype=bool))
+    lenj1, lenj2, lenj3 = get_local_problem_size(periodiclist, plist, np.array([False, False, False], dtype=bool))
 
     maxjwij = [lenj1, lenj2, lenj3]
 
@@ -3837,7 +3836,7 @@ def get_weights_local_projector(pts, fem_space):
     #######
 
     # We get the maximum number of j entries for whij
-    lenj1, lenj2, lenj3 = get_local_problem_size(periodiclist, plist, xp.array([True, True, True], dtype=bool))
+    lenj1, lenj2, lenj3 = get_local_problem_size(periodiclist, plist, np.array([True, True, True], dtype=bool))
 
     maxjwhij = [lenj1, lenj2, lenj3]
 
@@ -3870,21 +3869,21 @@ def get_weights_local_projector(pts, fem_space):
             counter = 1
             minicol = colmatrix[xstart:xend, bstart]
             while counter < 2 * degree - 1:
-                minicol = xp.column_stack(
+                minicol = np.column_stack(
                     (minicol, colmatrix[xstart:xend, (bstart + counter) % Nbasis]),
                 )
                 counter += 1
 
             # We need to consider the case in which our minicollocation matrix ends up being just one number
-            if xp.shape(minicol)[0] == 1:
+            if np.shape(minicol)[0] == 1:
                 # There seems to be a bug with the bsp.collocation_matrix function for the case num_elements = 1, degree = 1 and periodic, when evaluating the only B-spline at 0 the answer should be 1 not 0.
                 if degree == 1 and Nbasis == 1:
                     minicol[0] = 1.0
                 invmini = 1.0 / minicol[0]
                 for i in range(Nbasis):
-                    wijaux.append(xp.array([invmini]))
+                    wijaux.append(np.array([invmini]))
             else:
-                invmini = xp.linalg.inv(minicol)
+                invmini = np.linalg.inv(minicol)
                 for i in range(Nbasis):
                     wijaux.append(invmini[degree - 1, :])
         else:
@@ -3896,7 +3895,7 @@ def get_weights_local_projector(pts, fem_space):
                 # We can finally build the minicollocation matrix necessary to obtain the weights wij
                 minicol = colmatrix[xstart:xend, bstart:bend]
                 # Now we get its inverse
-                invmini = xp.linalg.inv(minicol)
+                invmini = np.linalg.inv(minicol)
 
                 # Now we need to extract the row of invmini that corresponds to the ith histopolation coefficient.
                 if i == 0:
@@ -3916,9 +3915,9 @@ def get_weights_local_projector(pts, fem_space):
                 for j in range(len(auxiliar), maxjwij[d]):
                     auxiliar.append(0.0)
 
-                wijaux.append(xp.array(auxiliar))
+                wijaux.append(np.array(auxiliar))
 
-        wij.append(xp.array(wijaux))
+        wij.append(np.array(wijaux))
 
         # Now that we know the wij we must use them to compute the whij
         # We begin by adressing the special case degree=1 and periodic
@@ -3928,14 +3927,14 @@ def get_weights_local_projector(pts, fem_space):
             nD = Nbasis
             if degree == 1:
                 for i in range(nD):
-                    whijaux.append(xp.array([wijaux[i][0], wijaux[i][0]]))
+                    whijaux.append(np.array([wijaux[i][0], wijaux[i][0]]))
             else:
                 whats = [wijaux[0][0], wijaux[0][0] + wijaux[0][1]]
                 for j in range(2, 2 * degree - 1):
                     whats.append(wijaux[0][j - 1] + wijaux[0][j])
                 whats.append(wijaux[0][2 * degree - 2])
                 for i in range(nD):
-                    whijaux.append(xp.array(whats))
+                    whijaux.append(np.array(whats))
 
         else:
             # Number of D-splines
@@ -4013,9 +4012,9 @@ def get_weights_local_projector(pts, fem_space):
                         else:
                             whats.append(0.0)
 
-                whijaux.append(xp.array(whats))
+                whijaux.append(np.array(whats))
 
-        whij.append(xp.array(whijaux))
+        whij.append(np.array(whijaux))
 
     return wij, whij
 
diff --git a/src/struphy/feec/tests/test_derham_gpu.py b/src/struphy/feec/tests/test_derham_gpu.py
new file mode 100644
index 000000000..93a165b04
--- /dev/null
+++ b/src/struphy/feec/tests/test_derham_gpu.py
@@ -0,0 +1,75 @@
+"""Derham on the CuPy backend: creation and the CUDA Derham arguments (CUDA_STRATEGY.md, PR 7).
+
+The tests marked ``requires_cupy`` need a GPU; they compare the CuPy backend with the NumPy backend.
+"""
+
+import cunumpy
+import numpy as np
+import pytest
+
+requires_cupy = pytest.mark.skipif(not cunumpy.cupy_available(), reason="CuPy/GPU not available")
+
+DECOMPOSITION = ("domain_array", "index_array", "index_array_N", "index_array_D", "neighbours")
+
+
+def make_derham(bcs=(None, None, None), local_projectors=False):
+    """Derham on 8 x 6 x 4 elements with degrees (2, 3, 1), for the active cunumpy backend."""
+    from feectools.ddm.mpi import mpi as MPI
+
+    from struphy.feec.psydac_derham import Derham
+    from struphy.io.options import DerhamOptions
+    from struphy.topology.grids import TensorProductGrid
+
+    options = DerhamOptions(degree=(2, 3, 1), bcs=bcs, local_projectors=local_projectors)
+    return Derham(TensorProductGrid(num_elements=(8, 6, 4)), options, comm=MPI.COMM_WORLD)
+
+
+def test_args_derham_on_numpy():
+    """On the NumPy backend the general arguments are the Pyccel host arguments."""
+    from struphy.kernel_arguments.pusher_args_kernels import DerhamArguments
+
+    with cunumpy.use_backend("numpy"):
+        derham = make_derham()
+        assert isinstance(derham.args_derham, DerhamArguments)
+        assert derham.args_derham is derham._pyccel_args_derham
+        for name in ("pn", "tn1", "tn2", "tn3", "starts"):
+            assert isinstance(getattr(derham.args_derham, name), np.ndarray), name
+
+
+@requires_cupy
+@pytest.mark.parametrize("bcs", [(None, None, None), (("dirichlet", "free"), None, ("free", "dirichlet"))])
+def test_derham_on_cupy(bcs):
+    """Same decomposition and kernel arguments on both backends, and CUDA arguments holding device copies."""
+    from struphy.utils.cuda_arguments import CudaDerhamArguments
+
+    derhams = {}
+    for backend in ("numpy", "cupy"):
+        with cunumpy.use_backend(backend):
+            derhams[backend] = make_derham(bcs)
+    host, device = derhams["numpy"], derhams["cupy"]
+
+    for name in DECOMPOSITION:
+        assert cunumpy.is_gpu(getattr(device, name)), name
+        assert np.array_equal(cunumpy.to_numpy(getattr(device, name)), getattr(host, name)), name
+
+    # the pyccel arguments are host arrays on both backends (feectools knots are host arrays)
+    for name in ("pn", "tn1", "tn2", "tn3", "starts"):
+        assert isinstance(getattr(device._pyccel_args_derham, name), np.ndarray), name
+        assert np.array_equal(getattr(device._pyccel_args_derham, name), getattr(host._pyccel_args_derham, name)), name
+
+    # Arguments keep the construction backend even when accessed from the NumPy backend.
+    with cunumpy.use_backend("numpy"):
+        args = device.args_derham
+        assert isinstance(args, CudaDerhamArguments)
+        assert device.args_derham is args
+        expected = (host.args_derham.pn, *host.V0fem.knots, host.args_derham.starts)
+        for name, value in zip(("pn", "tn1", "tn2", "tn3", "starts"), expected):
+            assert cunumpy.is_gpu(getattr(args, name)), name
+            assert np.array_equal(cunumpy.to_numpy(getattr(args, name)), value), name
+
+
+@requires_cupy
+def test_local_projectors_not_supported_on_cupy():
+    """Local projectors have no device implementation yet; the Derham fails when it is created."""
+    with cunumpy.use_backend("cupy"), pytest.raises(NotImplementedError, match="Local projectors"):
+        make_derham(local_projectors=True)
diff --git a/src/struphy/pic/accumulation/particles_to_grid.py b/src/struphy/pic/accumulation/particles_to_grid.py
index 358fe9dea..571d9be9a 100644
--- a/src/struphy/pic/accumulation/particles_to_grid.py
+++ b/src/struphy/pic/accumulation/particles_to_grid.py
@@ -228,7 +228,7 @@ def _accumulate(self, *optional_args):
         with ProfileManager.profile_region("kernel: " + self.kernel.name):
             self.kernel(
                 self.particles._pyccel_args_markers,
-                self.derham.args_derham,
+                self.derham._pyccel_args_derham,
                 self.args_domain,
                 *self._args_data,
                 *optional_args,
@@ -556,7 +556,7 @@ def _accumulate(self, *optional_args):
         with ProfileManager.profile_region("kernel: " + self.kernel.name):
             self.kernel(
                 self.particles._pyccel_args_markers,
-                self.derham.args_derham,
+                self.derham._pyccel_args_derham,
                 self.args_domain,
                 *self._args_data,
                 *optional_args,
diff --git a/src/struphy/utils/cuda_arguments.py b/src/struphy/utils/cuda_arguments.py
index 104934f2b..70c54edaa 100644
--- a/src/struphy/utils/cuda_arguments.py
+++ b/src/struphy/utils/cuda_arguments.py
@@ -155,6 +155,35 @@ def get_cuda_args(self) -> tuple:
         )
 
 
+class CudaDerhamArguments(Argument):
+    """CUDA version of :class:`~struphy.kernel_arguments.pusher_args_kernels.DerhamArguments`.
+
+    CUDA signature of :meth:`get_cuda_args`: ``long long* pn, double* tn1, double* tn2, double* tn3, long long* starts``
+
+    The scratch arrays of the pyccel class (``bn1``, ..., ``bd3``) are not part of it; CUDA kernels use
+    per-thread local arrays instead.
+
+    Parameters
+    ----------
+    pn : cupy.ndarray[int]
+        Spline degrees of :class:`~struphy.feec.psydac_derham.Derham` (int64).
+
+    tn1, tn2, tn3 : cupy.ndarray[float]
+        Knot sequences of :class:`~struphy.feec.psydac_derham.Derham`.
+
+    starts : cupy.ndarray[int]
+        Start indices (current MPI process) of :class:`~struphy.feec.psydac_derham.Derham` (int64).
+    """
+
+    def __init__(self, pn, tn1, tn2, tn3, starts):
+        self.pn = _cupy_array("pn", pn, np.int64)
+        self.tn1, self.tn2, self.tn3 = (_cupy_array("tn", t, np.float64) for t in (tn1, tn2, tn3))
+        self.starts = _cupy_array("starts", starts, np.int64)
+
+    def get_cuda_args(self) -> tuple:
+        return (self.pn, self.tn1, self.tn2, self.tn3, self.starts)
+
+
 class CudaDomainArguments(Argument):
     """CUDA version of :class:`~struphy.kernel_arguments.pusher_args_kernels.DomainArguments`.