diff --git a/feectools/core/bsplines.py b/feectools/core/bsplines.py
index 25e8c3ca8..189d95f0e 100644
--- a/feectools/core/bsplines.py
+++ b/feectools/core/bsplines.py
@@ -17,6 +17,7 @@
 """
 import cunumpy as xp
 from cunumpy.xp import array_backend
+import numpy as np
 
 from feectools.core.bsplines_kernels import (find_span_p,
                                           find_spans_p,
@@ -352,9 +353,9 @@ def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None,
         if periodic:
             nb -= degree + 1 - multiplicity
 
-        out = xp.zeros((int(xgrid.shape[0]), int(nb)), dtype=float)
+        out = np.zeros((int(xgrid.shape[0]), int(nb)), dtype=float)
     else:
-        assert out.shape == ((int(xgrid.shape[0]), int(nb))) and out.dtype == xp.dtype('float')
+        assert out.shape == ((int(xgrid.shape[0]), int(nb))) and out.dtype == np.dtype('float')
 
     bool_normalization = normalization == "M"
     multiplicity = int(multiplicity)
@@ -437,9 +438,9 @@ def histopolation_matrix(knots, degree, periodic, normalization, xgrid, multipli
 
     if out is None:
         if periodic:
-            out = xp.zeros((len(xgrid), len(knots) - 2 * int(degree) - 2 + int(multiplicity)), dtype=float)
+            out = np.zeros((len(xgrid), len(knots) - 2 * int(degree) - 2 + int(multiplicity)), dtype=float)
         else:
-            out = xp.zeros((len(xgrid) - 1, len(elevated_knots) - (int(degree) + 1) - 1 - 1), dtype=float)
+            out = np.zeros((len(xgrid) - 1, len(elevated_knots) - (int(degree) + 1) - 1 - 1), dtype=float)
     else:
         if periodic:
             assert out.shape == (len(xgrid), len(knots) - 2 * degree - 2 + multiplicity)
@@ -514,10 +515,15 @@ def greville(knots, degree, periodic, out=None, multiplicity=1):
         Abscissas of all Greville points.
 
     """
-    knots = xp.ascontiguousarray(knots, dtype=float)
+    # Greville points are index arrays, keep on NumPy
+    if isinstance(knots, (list, tuple)):
+        knots = np.asarray(knots, dtype=float)
+    if hasattr(knots, 'get'):
+        knots = knots.get()  # Convert CuPy to NumPy
+    knots = np.ascontiguousarray(knots, dtype=float)
     if out is None:
         n = len(knots) - 2 * degree - 2 + multiplicity if periodic else len(knots) - degree - 1
-        out = xp.zeros(int(n))
+        out = np.zeros(int(n))
     multiplicity = int(multiplicity)
     greville_p(knots, degree, periodic, out, multiplicity)
     return out
@@ -568,7 +574,7 @@ def elements_spans(knots, degree, out=None):
     """
     knots = xp.ascontiguousarray(knots, dtype=float)
     if out is None:
-        out = xp.zeros(len(knots), dtype=xp.int64)
+        out = np.zeros(len(knots), dtype=xp.int64)
     else:
         assert out.shape == knots.shape and out.dtype == xp.dtype('int64')
     i_final = elements_spans_p(knots, degree, out)
@@ -617,7 +623,11 @@ def make_knots(breaks, degree, periodic, multiplicity=1, out=None):
 
     # Consistency checks
     assert len(breaks) > 1
-    assert all( xp.diff(breaks) > 0 )
+    # Convert to numpy for comparison since assertion needs Python bool
+    breaks_np = breaks.get() if hasattr(breaks, 'get') else breaks
+    if isinstance(breaks_np, (list, tuple)):
+        breaks_np = np.asarray(breaks_np)
+    assert all( np.diff(breaks_np) > 0 )
     assert degree >= 0
     assert 1 <= multiplicity and multiplicity <= degree + 1
     # Cast potential numpy.int64 into python native int
@@ -626,9 +636,14 @@ def make_knots(breaks, degree, periodic, multiplicity=1, out=None):
     if periodic:
         assert len(breaks) > degree
 
-    breaks = xp.ascontiguousarray(breaks, dtype=float)
+    # Keep breaks on NumPy for initialization - knots are index arrays needed for CPU operations
+    breaks = np.asarray(breaks, dtype=float) if isinstance(breaks, (list, tuple)) else breaks
+    if hasattr(breaks, 'get'):
+        breaks = breaks.get()  # Convert CuPy to NumPy
+    breaks = np.ascontiguousarray(breaks, dtype=float)
     if out is None:
-        out = xp.zeros(multiplicity * len(breaks[1:-1]) + 2 + 2 * degree)
+        # Knots are index arrays, keep them on NumPy
+        out = np.zeros(multiplicity * len(breaks[1:-1]) + 2 + 2 * degree)
     else:
         assert out.shape == (multiplicity * len(breaks[1:-1]) + 2 + 2 * degree,) \
             and out.dtype == xp.dtype('float')
@@ -676,25 +691,29 @@ def elevate_knots(knots, degree, periodic, multiplicity=1, tol=1e-15, out=None):
         Knots sequence of spline space of degree p+1.
     """
     multiplicity = int(multiplicity)
-    knots = xp.ascontiguousarray(knots, dtype=float)
+    if isinstance(knots, (list, tuple)):
+        knots = np.asarray(knots, dtype=float)
+    if hasattr(knots, 'get'):
+        knots = knots.get()  # Convert CuPy to NumPy
+    knots = np.ascontiguousarray(knots, dtype=float)
     if out is None:
         if periodic:
-            out = xp.zeros(knots.shape[0] + 2, dtype=float)
+            out = np.zeros(knots.shape[0] + 2, dtype=float)
         else:
             shape = 2*(degree + 2)
             if len(knots) - 2 * (degree + 1) > 0:
-                uniques = (xp.diff(knots[degree + 1:-degree - 1]) > tol).nonzero()
+                uniques = (np.diff(knots[degree + 1:-degree - 1]) > tol).nonzero()
                 shape += multiplicity * (1 + uniques[0].shape[0])
-            out = xp.zeros(shape, dtype=float)
+            out = np.zeros(shape, dtype=float)
     else:
         if periodic:
-            assert out.shape == (knots.shape[0] + 2,) and out.dtype == xp.dtype('float')
+            assert out.shape == (knots.shape[0] + 2,) and out.dtype == np.dtype('float')
         else:
             shape = 2*(degree + 2)
             if len(knots) - 2 * (degree + 1) > 0:
-                uniques = (xp.diff(knots[degree + 1:-degree - 1]) > tol).nonzero()
+                uniques = (np.diff(knots[degree + 1:-degree - 1]) > tol).nonzero()
                 shape += multiplicity * (1 + uniques[0].shape[0])
-            assert out.shape == shape and out.dtype == xp.dtype('float')
+            assert out.shape == shape and out.dtype == np.dtype('float')
 
     elevate_knots_p(knots, degree, periodic, out, multiplicity, tol)
     return out
@@ -751,14 +770,18 @@ def quadrature_grid(breaks, quad_rule_x, quad_rule_w):
     assert min(quad_rule_x) >= -1
     assert max(quad_rule_x) <= +1
 
-    breaks = xp.ascontiguousarray(breaks, dtype=float)
+    # Convert breaks to numpy if CuPy (breaks/grids should stay on CPU)
+    if hasattr(breaks, 'get'):
+        breaks = breaks.get()
+    breaks = np.ascontiguousarray(breaks, dtype=float)
 
     if array_backend.backend == "cupy":
-        quad_rule_x = xp.ascontiguousarray(xp.array(quad_rule_x), dtype=float)
-        quad_rule_w = xp.ascontiguousarray( xp.array(quad_rule_w), dtype=float )
-    else:
-        quad_rule_x = xp.ascontiguousarray(quad_rule_x, dtype=float)
-        quad_rule_w = xp.ascontiguousarray( quad_rule_w, dtype=float )
+        # Convert CuPy arrays to NumPy
+        quad_rule_x = quad_rule_x.get() if hasattr(quad_rule_x, 'get') else quad_rule_x
+        quad_rule_w = quad_rule_w.get() if hasattr(quad_rule_w, 'get') else quad_rule_w
+    
+    quad_rule_x = np.ascontiguousarray(quad_rule_x, dtype=float)
+    quad_rule_w = np.ascontiguousarray(quad_rule_w, dtype=float)
     
 
     out1 = xp.zeros((len(breaks) - 1, len(quad_rule_x)))
@@ -914,7 +937,7 @@ def cell_index(breaks, i_grid, tol=1e-15, out=None):
     breaks = xp.ascontiguousarray(breaks, dtype=float)
     i_grid = xp.ascontiguousarray(i_grid, dtype=float)
     if out is None:
-        out = xp.zeros_like(i_grid, dtype=xp.int64)
+        out = np.zeros_like(i_grid, dtype=xp.int64)
     else:
         assert out.shape == i_grid.shape and out.dtype == xp.dtype('int64')
     status = cell_index_p(breaks, i_grid, tol, out)
diff --git a/feectools/core/bsplines_kernels.py b/feectools/core/bsplines_kernels.py
index e019e59f8..f59f0673c 100644
--- a/feectools/core/bsplines_kernels.py
+++ b/feectools/core/bsplines_kernels.py
@@ -7,6 +7,8 @@
 # This file holds the pyccelisable versions of the functions in bsplines.py
 # This will be changed once pyccel can return arrays and can get out=None arguments
 # like Numpy functions.
+# NOTE: This file must use ONLY numpy for pyccel compilation compatibility.
+# Backend conversion (NumPy/CuPy) happens at the Python wrapper level.
 
 from pyccel.decorators import pure
 from numpy import shape, abs
@@ -401,6 +403,9 @@ def basis_funs_all_ders_p(knots: 'float[:]', degree: int, x: float, span: int, n
     .. [1] L. Piegl and W. Tiller. The NURBS Book, 2nd ed.,
         Springer-Verlag Berlin Heidelberg GmbH, 1997.
     """
+    # Detect backend from output array
+    # Backend array operations removed - always use numpy
+    
     sh_a  = np.empty(2)
     sh_b  = np.empty(2)
     left  = np.empty(degree)
@@ -556,6 +561,9 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali
         multiplicity applies to each interior knot.
         
     """
+    # Detect backend from output array
+    # Backend array operations removed - always use numpy
+    
     # Number of basis functions (in periodic case remove degree repeated elements)
     nb = len(knots)-degree-1
     if periodic:
@@ -565,7 +573,7 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali
     nx = len(xgrid)
 
     basis = np.zeros((nx, degree + 1))
-    spans = np.zeros(nx, dtype=int)
+    spans = np.zeros(nx, dtype=int)  # Keep indices on CPU
     find_spans_p(knots, degree, xgrid, spans)
     basis_funs_array_p(knots, degree, xgrid, spans, basis)
 
@@ -644,6 +652,9 @@ def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, norma
     contains the integrals of each B-spline basis function :math:`B_j` between
     two successive grid points.
     """
+    # Detect backend from output array
+    # Backend array operations removed - always use numpy
+    
     nb = len(knots) - degree - 1
     if periodic:
         nb -= degree + 1 - multiplicity
@@ -700,7 +711,7 @@ def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, norma
     m = colloc.shape[0] - 1
     n = colloc.shape[1] - 1
 
-    spans = np.zeros(colloc.shape[0], dtype=int)
+    spans = np.zeros(colloc.shape[0], dtype=int)  # Keep indices on CPU
     for i in range(colloc.shape[0]):
         local_span = 0
         for j in range(colloc.shape[1]):
@@ -767,7 +778,10 @@ def merge_sort(a: 'float[:]'):
     """
     if len(a) != 1 and len(a) != 0:
         n = len(a)
-
+        
+        # Detect backend and use the appropriate array module
+        # Backend array operations removed - always use numpy
+        
         a1 = np.zeros(n // 2)
         a1[:] = a[:n // 2]
         a2 = np.zeros(n - n // 2)
@@ -1167,12 +1181,16 @@ def basis_ders_on_quad_grid_p(knots: 'float[:]', degree: int, quad_grid: 'float[
     """
     ne = quad_grid.shape[0]
     nq = quad_grid.shape[1]
+    
+    # Detect backend from output array
+    # Backend array operations removed - always use numpy
+    
     if normalization:
         integrals = np.zeros(knots.shape[0] - degree - 1)
         basis_integrals_p(knots, degree, integrals)
         scaling = 1.0 /integrals
 
-    temp_spans = np.zeros(len(knots), dtype=int)
+    temp_spans = np.zeros(len(knots), dtype=int)  # Keep indices on CPU
     actual_index = elements_spans_p(knots, degree, temp_spans)
     spans = temp_spans[:actual_index]
 
diff --git a/feectools/ddm/blocking_data_exchanger.py b/feectools/ddm/blocking_data_exchanger.py
index 7c506fa32..b8be40bf3 100644
--- a/feectools/ddm/blocking_data_exchanger.py
+++ b/feectools/ddm/blocking_data_exchanger.py
@@ -1,8 +1,5 @@
-#---------------------------------------------------------------------------#
-# This file is part of PSYDAC which is released under MIT License. See the  #
-# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
-# for full license details.                                                 #
-#---------------------------------------------------------------------------#
+# coding: utf-8
+
 import cunumpy as xp
 import numpy as np
 from feectools.ddm.mpi import mpi as MPI
diff --git a/feectools/ddm/cart.py b/feectools/ddm/cart.py
index f00efff90..2b2b58b41 100644
--- a/feectools/ddm/cart.py
+++ b/feectools/ddm/cart.py
@@ -1,9 +1,20 @@
 # coding: utf-8
 
 import os
+import numpy as np
 import cunumpy as xp
+from cunumpy.xp import array_backend
 from itertools import product
 
+# Initialize CUDA context before MPI if using CuPy backend
+if array_backend.backend == "cupy":
+    try:
+        import cupy as cp
+        cp.cuda.Device(0).use()
+        cp.cuda.Stream.null.synchronize()
+    except Exception:
+        pass
+
 from feectools.ddm.mpi import mpi as MPI
 from feectools.ddm.mpi import MockMPI
 from feectools.ddm.partition import compute_dims, partition_procs_per_patch
@@ -43,6 +54,17 @@ def find_mpi_type( dtype ):
 
     return mpi_type
 
+#===============================================================================
+def _cuda_sync_before_mpi():
+    """Synchronize CUDA before MPI operations to avoid conflicts."""
+    if array_backend.backend == "cupy":
+        try:
+            import cupy as cp
+            cp.cuda.Stream.null.synchronize()
+        except Exception:
+            pass
+
+#===============================================================================
 class MultiPatchDomainDecomposition:
     """
     Cartesian decomposition of multiple N-Cube grids.
@@ -260,11 +282,14 @@ def __init__(self, ncells, periods, comm=None, global_comm=None, num_threads=Non
 
         if comm is None:
             # compute the coords for all processes
-            self._global_coords = xp.array([xp.unravel_index(xp.int64(rank), nprocs) for rank in range(self._size)])
+            self._global_coords = np.array([np.unravel_index(int(rank), nprocs) for rank in range(self._size)])
             self._coords        = self._global_coords[self._rank]
             self._rank_in_topo  = 0
             self._ranks_in_topo = xp.array([0])
         else:
+            # Synchronize CUDA before MPI operations
+            _cuda_sync_before_mpi()
+            
             # Create a MPI cart
             self._comm_cart = comm.Create_cart(
                 dims    = self._nprocs,
@@ -288,6 +313,7 @@ def __init__(self, ncells, periods, comm=None, global_comm=None, num_threads=Non
         # Create (N-1)-dimensional communicators within the Cartesian topology
         self._subcomm = [None]*self._ndims
         for i in range(self._ndims):
+            _cuda_sync_before_mpi()  # Synchronize before each Sub() call
             remain_dims     = [i==j for j in range( self._ndims )]
             self._subcomm[i] = self._comm_cart.Sub( remain_dims )
 
@@ -467,8 +493,9 @@ def __init__( self, domain_decomposition, npts, global_starts, global_ends, pads
         # Store input arguments
         self._domain_decomposition = domain_decomposition
         self._npts          = tuple( npts    )
-        self._global_starts = tuple( [ xp.asarray(gs) for gs in global_starts]  )
-        self._global_ends   = tuple( [ xp.asarray(ge) for ge in global_ends]    )
+        # Convert to NumPy arrays for MPI compatibility (MPI can't handle CuPy arrays)
+        self._global_starts = tuple( [ np.asarray(gs.get() if hasattr(gs, 'get') else gs) for gs in global_starts]  )
+        self._global_ends   = tuple( [ np.asarray(ge.get() if hasattr(ge, 'get') else ge) for ge in global_ends]    )
         self._pads          = tuple( pads    )
         self._shifts        = tuple( shifts  )
         self._periods       = domain_decomposition.periods
@@ -494,10 +521,12 @@ def __init__( self, domain_decomposition, npts, global_starts, global_ends, pads
 
         # Know my coordinates in the topology
         self._coords = domain_decomposition.coords
+        # Convert coords to NumPy for indexing (MPI coords should be on CPU)
+        coords_np = [c.get() if hasattr(c, 'get') else c for c in self._coords]
 
         # Start/end values of global indices (without ghost regions)
-        self._starts = tuple( self._global_starts[axis][c] for axis,c in zip(range(self._ndims), self._coords) )
-        self._ends   = tuple( self._global_ends  [axis][c] for axis,c in zip(range(self._ndims), self._coords) )
+        self._starts = tuple( self._global_starts[axis][c] for axis,c in zip(range(self._ndims), coords_np) )
+        self._ends   = tuple( self._global_ends  [axis][c] for axis,c in zip(range(self._ndims), coords_np) )
 
         # List of 1D global indices (without ghost regions)
         # self._grids = tuple( range(s,e+1) for s,e in zip( self._starts, self._ends ) )
@@ -906,7 +935,7 @@ def _compute_shift_info_non_blocking( self, shift ):
         if len([i for i in shift if i==0]) == 2 and rank_dest != MPI.PROC_NULL:
             direction = [i for i,s in enumerate(shift) if s != 0][0]
             comm = self._subcomm[direction]
-            # local_dest_rank = self._comm_cart.group.Translate_ranks(xp.array([rank_dest]), comm.group)[0]
+            # local_dest_rank = self._comm_cart.group.Translate_ranks(np.array([rank_dest]), comm.group)[0]
             local_dest_rank = self._comm_cart.group.Translate_ranks([int(rank_dest)], comm.group)[0]
 
         else:
@@ -921,7 +950,7 @@ def _compute_shift_info_non_blocking( self, shift ):
         if len([i for i in shift if i==0]) == 2 and rank_source != MPI.PROC_NULL:
             direction = [i for i,s in enumerate(shift) if s != 0][0]
             comm = self._subcomm[direction]
-            # local_source_rank = self._comm_cart.group.Translate_ranks(xp.array([rank_source]), comm.group)[0]
+            # local_source_rank = self._comm_cart.group.Translate_ranks(np.array([rank_source]), comm.group)[0]
             local_source_rank = self._comm_cart.group.Translate_ranks([int(rank_source)], comm.group)[0]
         else:
             local_source_rank = rank_source
diff --git a/feectools/ddm/mpi.py b/feectools/ddm/mpi.py
index d79d37a1c..9b6caf23d 100644
--- a/feectools/ddm/mpi.py
+++ b/feectools/ddm/mpi.py
@@ -81,6 +81,11 @@ def COMM_WORLD(self):
 
 
 try:
+    # Disable MPI when using CuPy due to known segfault issues with OpenMPI + CUDA
+    import os
+    if os.environ.get('ARRAY_BACKEND') == 'cupy':
+        raise ImportError("MPI disabled when using CuPy backend")
+    
     from mpi4py import MPI
 
     _comm = MPI.COMM_WORLD
diff --git a/feectools/ddm/partition.py b/feectools/ddm/partition.py
index ef0e821c8..8b2b0d3b7 100644
--- a/feectools/ddm/partition.py
+++ b/feectools/ddm/partition.py
@@ -1,9 +1,5 @@
-#---------------------------------------------------------------------------#
-# This file is part of PSYDAC which is released under MIT License. See the  #
-# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
-# for full license details.                                                 #
-#---------------------------------------------------------------------------#
 import cunumpy as xp
+import numpy as np
 import numpy.ma as ma
 
 from sympy.ntheory import factorint
diff --git a/feectools/ddm/petsc.py b/feectools/ddm/petsc.py
index 54719062b..4f3a8b9c3 100644
--- a/feectools/ddm/petsc.py
+++ b/feectools/ddm/petsc.py
@@ -1,8 +1,6 @@
-#---------------------------------------------------------------------------#
-# This file is part of PSYDAC which is released under MIT License. See the  #
-# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
-# for full license details.                                                 #
-#---------------------------------------------------------------------------#
+# coding: utf-8
+
+import cunumpy as xp
 from itertools import product
 
 import cunumpy as xp
@@ -113,4 +111,3 @@ def create_g2n(self, gvec, natural):
         to_is   = self.petsc.IS().createGeneral(indices, comm=cart.comm)
 
         return self.petsc.Scatter().create(gvec, from_is, natural, to_is)
-
diff --git a/feectools/feec/global_geometric_projectors.py b/feectools/feec/global_geometric_projectors.py
index 1a44ac0fb..3dc7752be 100644
--- a/feectools/feec/global_geometric_projectors.py
+++ b/feectools/feec/global_geometric_projectors.py
@@ -2,6 +2,7 @@
 
 import cunumpy as xp
 from cunumpy.xp import array_backend
+import numpy as np
 
 from feectools.linalg.kron           import KroneckerLinearSolver, KroneckerStencilMatrix
 from feectools.linalg.stencil        import StencilMatrix, StencilVectorSpace
@@ -25,6 +26,18 @@
            'evaluate_dofs_2d_0form', 'evaluate_dofs_2d_1form_hcurl', 'evaluate_dofs_2d_1form_hdiv', 'evaluate_dofs_2d_2form',
            'evaluate_dofs_3d_0form', 'evaluate_dofs_3d_1form', 'evaluate_dofs_3d_2form', 'evaluate_dofs_3d_3form')
 
+
+def _to_numpy_for_kernel(*args):
+    """Convert CuPy arrays to NumPy for compiled kernel calls."""
+    result = []
+    for arg in args:
+        if hasattr(arg, 'get'):  # CuPy array
+            result.append(arg.get())
+        else:
+            result.append(arg)
+    return result if len(result) > 1 else result[0]
+
+
 #==============================================================================
 class GlobalGeometricProjector(metaclass=ABCMeta):
     """
@@ -207,8 +220,9 @@ def __init__(self, space, nquads = None):
                     solvercells += [V._histopolator]
                     
                     # make 1D collocation matrix in stencil format
+                    # Always use NumPy for indices since they're used for indexing/comparison
                     if array_backend.backend == "cupy":
-                        row_indices, col_indices = xp.nonzero(xp.array(V.hmat))
+                        row_indices, col_indices = np.nonzero(np.asarray(V.hmat))
                     else:
                         row_indices, col_indices = xp.nonzero(V.hmat)
 
@@ -800,7 +814,10 @@ def evaluate_dofs_1d_0form(
     
     F_temp = xp.zeros_like(F, order='C')
     
-    dof_kernels.evaluate_dofs_1d_0form(F_temp, f_pts)
+    F_temp_np, f_pts_np = _to_numpy_for_kernel(F_temp, f_pts)
+    dof_kernels.evaluate_dofs_1d_0form(F_temp_np, f_pts_np)
+    if hasattr(F_temp, 'get'):
+        F_temp[:] = xp.asarray(F_temp_np)
     
     F[:] = F_temp
         
@@ -819,7 +836,10 @@ def evaluate_dofs_1d_1form(
     # call kernel
     F_temp = xp.zeros_like(F, order='C')
     
-    dof_kernels.evaluate_dofs_1d_1form(quad_w1, F_temp, f_pts)
+    quad_w1_np, F_temp_np, f_pts_np = _to_numpy_for_kernel(quad_w1, F_temp, f_pts)
+    dof_kernels.evaluate_dofs_1d_1form(quad_w1_np, F_temp_np, f_pts_np)
+    if hasattr(F_temp, 'get'):
+        F_temp[:] = xp.asarray(F_temp_np)
     
     F[:] = F_temp
 
@@ -842,7 +862,10 @@ def evaluate_dofs_2d_0form(
     
     F_temp = xp.zeros_like(F, order='C')
     
-    dof_kernels.evaluate_dofs_2d_0form(F_temp, f_pts)
+    F_temp_np, f_pts_np = _to_numpy_for_kernel(F_temp, f_pts)
+    dof_kernels.evaluate_dofs_2d_0form(F_temp_np, f_pts_np)
+    if hasattr(F_temp, 'get'):
+        F_temp[:] = xp.asarray(F_temp_np)
     
     F[:, :] = F_temp
 
@@ -869,7 +892,12 @@ def evaluate_dofs_2d_1form_hcurl(
     F1_temp = xp.zeros_like(F1, order='C')
     F2_temp = xp.zeros_like(F2, order='C')
     
-    dof_kernels.evaluate_dofs_2d_1form_hcurl(quad_w1, quad_w2, F1_temp, F2_temp, f1_pts, f2_pts)
+    quad_w1_np, quad_w2_np, F1_temp_np, F2_temp_np, f1_pts_np, f2_pts_np = _to_numpy_for_kernel(quad_w1, quad_w2, F1_temp, F2_temp, f1_pts, f2_pts)
+    dof_kernels.evaluate_dofs_2d_1form_hcurl(quad_w1_np, quad_w2_np, F1_temp_np, F2_temp_np, f1_pts_np, f2_pts_np)
+    if hasattr(F1_temp, 'get'):
+        F1_temp[:] = xp.asarray(F1_temp_np)
+    if hasattr(F2_temp, 'get'):
+        F2_temp[:] = xp.asarray(F2_temp_np)
     
     F1[:, :] = F1_temp
     F2[:, :] = F2_temp
@@ -897,7 +925,12 @@ def evaluate_dofs_2d_1form_hdiv(
     F1_temp = xp.zeros_like(F1, order='C')
     F2_temp = xp.zeros_like(F2, order='C')
     
-    dof_kernels.evaluate_dofs_2d_1form_hdiv(quad_w1, quad_w2, F1_temp, F2_temp, f1_pts, f2_pts)
+    quad_w1_np, quad_w2_np, F1_temp_np, F2_temp_np, f1_pts_np, f2_pts_np = _to_numpy_for_kernel(quad_w1, quad_w2, F1_temp, F2_temp, f1_pts, f2_pts)
+    dof_kernels.evaluate_dofs_2d_1form_hdiv(quad_w1_np, quad_w2_np, F1_temp_np, F2_temp_np, f1_pts_np, f2_pts_np)
+    if hasattr(F1_temp, 'get'):
+        F1_temp[:] = xp.asarray(F1_temp_np)
+    if hasattr(F2_temp, 'get'):
+        F2_temp[:] = xp.asarray(F2_temp_np)
     
     F1[:, :] = F1_temp
     F2[:, :] = F2_temp
@@ -917,7 +950,10 @@ def evaluate_dofs_2d_2form(
     # call kernel
     F_temp = xp.zeros_like(F, order='C')
     
-    dof_kernels.evaluate_dofs_2d_2form(quad_w1, quad_w2, F_temp, f_pts)
+    quad_w1_np, quad_w2_np, F_temp_np, f_pts_np = _to_numpy_for_kernel(quad_w1, quad_w2, F_temp, f_pts)
+    dof_kernels.evaluate_dofs_2d_2form(quad_w1_np, quad_w2_np, F_temp_np, f_pts_np)
+    if hasattr(F_temp, 'get'):
+        F_temp[:] = xp.asarray(F_temp_np)
     
     F[:, :] = F_temp
 
@@ -940,7 +976,12 @@ def evaluate_dofs_2d_vec(
     F1_temp = xp.zeros_like(F1, order='C')
     F2_temp = xp.zeros_like(F2, order='C')
     
-    dof_kernels.evaluate_dofs_2d_vec(F1_temp, F2_temp, f1_pts, f2_pts)
+    F1_temp_np, F2_temp_np, f1_pts_np, f2_pts_np = _to_numpy_for_kernel(F1_temp, F2_temp, f1_pts, f2_pts)
+    dof_kernels.evaluate_dofs_2d_vec(F1_temp_np, F2_temp_np, f1_pts_np, f2_pts_np)
+    if hasattr(F1_temp, 'get'):
+        F1_temp[:] = xp.asarray(F1_temp_np)
+    if hasattr(F2_temp, 'get'):
+        F2_temp[:] = xp.asarray(F2_temp_np)
     
     F1[:, :] = F1_temp
     F2[:, :] = F2_temp
@@ -965,7 +1006,10 @@ def evaluate_dofs_3d_0form(
     
     F_temp = xp.zeros_like(F, order='C')
     
-    dof_kernels.evaluate_dofs_3d_0form(F_temp, f_pts)
+    F_temp_np, f_pts_np = _to_numpy_for_kernel(F_temp, f_pts)
+    dof_kernels.evaluate_dofs_3d_0form(F_temp_np, f_pts_np)
+    if hasattr(F_temp, 'get'):
+        F_temp[:] = xp.asarray(F_temp_np)
     
     F[:, :, :] = F_temp
 
@@ -997,7 +1041,14 @@ def evaluate_dofs_3d_1form(
     F2_temp = xp.zeros_like(F2, order='C')
     F3_temp = xp.zeros_like(F3, order='C')
     
-    dof_kernels.evaluate_dofs_3d_1form(quad_w1, quad_w2, quad_w3, F1_temp, F2_temp, F3_temp, f1_pts, f2_pts, f3_pts)
+    quad_w1_np, quad_w2_np, quad_w3_np, F1_temp_np, F2_temp_np, F3_temp_np, f1_pts_np, f2_pts_np, f3_pts_np = _to_numpy_for_kernel(quad_w1, quad_w2, quad_w3, F1_temp, F2_temp, F3_temp, f1_pts, f2_pts, f3_pts)
+    dof_kernels.evaluate_dofs_3d_1form(quad_w1_np, quad_w2_np, quad_w3_np, F1_temp_np, F2_temp_np, F3_temp_np, f1_pts_np, f2_pts_np, f3_pts_np)
+    if hasattr(F1_temp, 'get'):
+        F1_temp[:] = xp.asarray(F1_temp_np)
+    if hasattr(F2_temp, 'get'):
+        F2_temp[:] = xp.asarray(F2_temp_np)
+    if hasattr(F3_temp, 'get'):
+        F3_temp[:] = xp.asarray(F3_temp_np)
     
     F1[:, :, :] = F1_temp
     F2[:, :, :] = F2_temp
@@ -1031,7 +1082,14 @@ def evaluate_dofs_3d_2form(
     F2_temp = xp.zeros_like(F2, order='C')
     F3_temp = xp.zeros_like(F3, order='C')
     
-    dof_kernels.evaluate_dofs_3d_2form(quad_w1, quad_w2, quad_w3, F1_temp, F2_temp, F3_temp, f1_pts, f2_pts, f3_pts)
+    quad_w1_np, quad_w2_np, quad_w3_np, F1_temp_np, F2_temp_np, F3_temp_np, f1_pts_np, f2_pts_np, f3_pts_np = _to_numpy_for_kernel(quad_w1, quad_w2, quad_w3, F1_temp, F2_temp, F3_temp, f1_pts, f2_pts, f3_pts)
+    dof_kernels.evaluate_dofs_3d_2form(quad_w1_np, quad_w2_np, quad_w3_np, F1_temp_np, F2_temp_np, F3_temp_np, f1_pts_np, f2_pts_np, f3_pts_np)
+    if hasattr(F1_temp, 'get'):
+        F1_temp[:] = xp.asarray(F1_temp_np)
+    if hasattr(F2_temp, 'get'):
+        F2_temp[:] = xp.asarray(F2_temp_np)
+    if hasattr(F3_temp, 'get'):
+        F3_temp[:] = xp.asarray(F3_temp_np)
     
     F1[:, :, :] = F1_temp
     F2[:, :, :] = F2_temp
@@ -1052,7 +1110,10 @@ def evaluate_dofs_3d_3form(
     # call kernel
     F_temp = xp.zeros_like(F, order='C')
     
-    dof_kernels.evaluate_dofs_3d_3form(quad_w1, quad_w2, quad_w3, F_temp, f_pts)
+    quad_w1_np, quad_w2_np, quad_w3_np, F_temp_np, f_pts_np = _to_numpy_for_kernel(quad_w1, quad_w2, quad_w3, F_temp, f_pts)
+    dof_kernels.evaluate_dofs_3d_3form(quad_w1_np, quad_w2_np, quad_w3_np, F_temp_np, f_pts_np)
+    if hasattr(F_temp, 'get'):
+        F_temp[:] = xp.asarray(F_temp_np)
     
     F[:, :, :] = F_temp
 
@@ -1078,7 +1139,14 @@ def evaluate_dofs_3d_vec(
     F2_temp = xp.zeros_like(F2, order='C')
     F3_temp = xp.zeros_like(F3, order='C')
     
-    dof_kernels.evaluate_dofs_3d_vec(F1_temp, F2_temp, F3_temp, f1_pts, f2_pts, f3_pts)
+    F1_temp_np, F2_temp_np, F3_temp_np, f1_pts_np, f2_pts_np, f3_pts_np = _to_numpy_for_kernel(F1_temp, F2_temp, F3_temp, f1_pts, f2_pts, f3_pts)
+    dof_kernels.evaluate_dofs_3d_vec(F1_temp_np, F2_temp_np, F3_temp_np, f1_pts_np, f2_pts_np, f3_pts_np)
+    if hasattr(F1_temp, 'get'):
+        F1_temp[:] = xp.asarray(F1_temp_np)
+    if hasattr(F2_temp, 'get'):
+        F2_temp[:] = xp.asarray(F2_temp_np)
+    if hasattr(F3_temp, 'get'):
+        F3_temp[:] = xp.asarray(F3_temp_np)
     
     F1[:, :, :] = F1_temp
     F2[:, :, :] = F2_temp
diff --git a/feectools/fem/grid.py b/feectools/fem/grid.py
index 7dbf8f5a6..eea96c146 100644
--- a/feectools/fem/grid.py
+++ b/feectools/fem/grid.py
@@ -1,8 +1,7 @@
-#---------------------------------------------------------------------------#
-# This file is part of PSYDAC which is released under MIT License. See the  #
-# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
-# for full license details.                                                 #
-#---------------------------------------------------------------------------#
+# coding: utf-8
+#
+# Copyright 2018 Yaman Güçlü
+
 import cunumpy as xp
 
 from feectools.core.bsplines         import elements_spans
diff --git a/feectools/fem/partitioning.py b/feectools/fem/partitioning.py
index 104b506ea..4df8bae0b 100644
--- a/feectools/fem/partitioning.py
+++ b/feectools/fem/partitioning.py
@@ -1,6 +1,7 @@
 # -*- coding: UTF-8 -*-
 import os
 
+import numpy as np
 import cunumpy as xp
 
 from feectools.ddm.cart       import CartDecomposition, InterfaceCartDecomposition, create_interfaces_cart
@@ -56,6 +57,10 @@ def partition_coefficients(domain_decomposition, spaces, min_blocks=None):
         ee = domain_decomposition.global_element_ends  [axis]
         m  = multiplicity[axis]
 
+        # Convert to numpy if CuPy (needed for MPI operations later)
+        if hasattr(ee, 'get'):
+            ee = ee.get()
+        
         global_ends  [axis]     = m*(ee+1)-1
         global_ends  [axis][-1] = npts[axis]-1
         global_starts[axis]     = xp.array([0] + (global_ends[axis][:-1]+1).tolist())
@@ -64,14 +69,18 @@ def partition_coefficients(domain_decomposition, spaces, min_blocks=None):
         min_blocks = [None] * ndims
 
     for s, e, V, mb in zip(global_starts, global_ends, spaces, min_blocks):
+        s_host = s.get() if hasattr(s, 'get') else np.asarray(s)
+        e_host = e.get() if hasattr(e, 'get') else np.asarray(e)
+        local_sizes = e_host - s_host + 1
+
         if V.periodic or mb is None:
-            assert all(e-s+1 >= V.degree), f"Local number of elements (after domain decomposition) is to small for spline degree p={V.degree}: {e-s+1} is not >= {V.degree} everywhere.\n \
+            assert all(local_sizes >= V.degree), f"Local number of elements (after domain decomposition) is to small for spline degree p={V.degree}: {local_sizes} is not >= {V.degree} everywhere.\n \
             You can:\n \
                 1. increase Nel\n \
                 2. lower p\n \
                 3. decrease the MPI size."
         else:
-            assert all(e-s+1 >= mb)
+            assert all(local_sizes >= mb)
 
     return global_starts, global_ends
 
diff --git a/feectools/fem/splines.py b/feectools/fem/splines.py
index 41c41c4a8..ac62f9fa7 100644
--- a/feectools/fem/splines.py
+++ b/feectools/fem/splines.py
@@ -140,7 +140,7 @@ def __init__(self, degree, knots=None, grid=None, multiplicity=None, parent_mult
 
         # Create space of spline coefficients
         domain_decomposition = DomainDecomposition([self._ncells], [periodic])
-        cart     = CartDecomposition(domain_decomposition, [nbasis], [xp.array([0])],[xp.array([nbasis-1])], [self._pads], [multiplicity])
+        cart     = CartDecomposition(domain_decomposition, [nbasis], [_np.array([0])],[_np.array([nbasis-1])], [self._pads], [multiplicity])
         self._coeff_space = StencilVectorSpace(cart)
 
         # Store flag: object NOT YET prepared for interpolation / histopolation
@@ -189,7 +189,7 @@ def init_interpolation( self, dtype=float ):
             # Convert to CSC format and compute sparse LU decomposition
             
             # Convert to LAPACK banded format (see DGBTRF function)
-            if array_backend.backend == "cupy":
+            if hasattr(imat, 'get'):
                 imat = imat.get()
             else:
                 imat = _np.asanyarray(imat)
@@ -231,7 +231,7 @@ def init_histopolation( self, dtype=float):
             xgrid    = self.ext_greville,
             multiplicity = self._multiplicity
         )
-        if array_backend.backend == "cupy":
+        if hasattr(imat, 'get'):
             imat = imat.get()
         else:
             imat = _np.asanyarray(imat)
diff --git a/feectools/fem/tensor.py b/feectools/fem/tensor.py
index 61d774fe2..513635c52 100644
--- a/feectools/fem/tensor.py
+++ b/feectools/fem/tensor.py
@@ -8,6 +8,7 @@
 from feectools.ddm.mpi import mpi as MPI
     
 import cunumpy as xp
+import numpy as np
 import itertools
 import h5py
 import os
diff --git a/feectools/fem/tests/utilities.py b/feectools/fem/tests/utilities.py
index 1685a585f..fa1cd0d26 100644
--- a/feectools/fem/tests/utilities.py
+++ b/feectools/fem/tests/utilities.py
@@ -1,8 +1,6 @@
-#---------------------------------------------------------------------------#
-# This file is part of PSYDAC which is released under MIT License. See the  #
-# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
-# for full license details.                                                 #
-#---------------------------------------------------------------------------#
+# coding: utf-8
+# Copyright 2018 Yaman Güçlü
+
 import cunumpy as xp
 
 #===============================================================================
diff --git a/feectools/linalg/stencil.py b/feectools/linalg/stencil.py
index 58df349a3..4848595c3 100644
--- a/feectools/linalg/stencil.py
+++ b/feectools/linalg/stencil.py
@@ -9,7 +9,6 @@
 
 import cunumpy as xp
 from cunumpy.xp import array_backend
-from types        import MappingProxyType
 from scipy.sparse import coo_matrix, diags as sp_diags
 
 from feectools.ddm.mpi import mpi as MPI
@@ -35,8 +34,23 @@
     'StencilInterfaceMatrix'
 )
 
-#===============================================================================
-# Dictionary used to select correct kernel functions based on dimensionality
+#========================================================================
+def _to_numpy_int64(val):
+    """Convert CuPy or NumPy scalar/array to numpy int64."""
+    import numpy as _np
+    if hasattr(val, 'get'):
+        # CuPy array - convert to NumPy first
+        val = val.get()
+    return _np.int64(val)
+
+def _to_numpy_array(val):
+    """Convert CuPy array to NumPy array, preserving dtype. Return as-is if already NumPy."""
+    if hasattr(val, 'get'):
+        # CuPy array - convert to NumPy
+        return val.get()
+    return val
+
+#========================================================================# Dictionary used to select correct kernel functions based on dimensionality
 kernels = {
     'axpy'  : (None,   axpy_1d,   axpy_2d,   axpy_3d),
     'inner' : (None,  inner_1d,  inner_2d,  inner_3d),
@@ -47,7 +61,7 @@
                     'C': (None, stencil2coo_1d_C, stencil2coo_2d_C, stencil2coo_3d_C)}
 }
 
-#===============================================================================
+#========================================================================
 def compute_diag_len(pads, shifts_domain, shifts_codomain, return_padding=False):
     """
     Compute the diagonal length and the padding of the stencil matrix for each direction,
@@ -83,7 +97,7 @@ def compute_diag_len(pads, shifts_domain, shifts_codomain, return_padding=False)
     else:
         return n.astype('int')
 
-#===============================================================================
+#========================================================================
 class StencilVectorSpace(VectorSpace):
     """
     Vector space for n-dimensional stencil format. Two different initializations
@@ -301,10 +315,22 @@ def axpy(self, a, x, y):
             else:
                 a = float(a)
 
-        self._axpy_func(a, x._data, y._data)
+        x_data_np = _to_numpy_array(x._data)
+        y_data_np = _to_numpy_array(y._data)
+        self._axpy_func(a, x_data_np, y_data_np)
+        # Copy result back if CuPy
+        if hasattr(y._data, 'get'):
+            import cupy as cp
+            y._data[:] = cp.asarray(y_data_np)
 
         for axis, ext in self.interfaces:
-            self._axpy_func(a, x._interface_data[axis, ext], y._interface_data[axis, ext])
+            x_int_np = _to_numpy_array(x._interface_data[axis, ext])
+            y_int_np = _to_numpy_array(y._interface_data[axis, ext])
+            self._axpy_func(a, x_int_np, y_int_np)
+            # Copy result back if CuPy
+            if hasattr(y._interface_data[axis, ext], 'get'):
+                import cupy as cp
+                y._interface_data[axis, ext][:] = cp.asarray(y_int_np)
 
         x._sync = x._sync and y._sync
 
@@ -431,7 +457,7 @@ def set_interface(self, axis, ext, cart):
 
             self._interfaces[axis, ext] = space
 
-#===============================================================================
+#========================================================================
 class StencilVector(Vector):
     """
     Vector in n-dimensional stencil format.
@@ -638,7 +664,7 @@ def toarray_local(self , *, order='C'):
 
     # ...
     def _toarray_parallel_no_pads(self, order='C'):
-        a         = xp.zeros( self.space.npts, dtype=self.dtype )
+        a         = xp.zeros(self.space.npts, dtype=self.dtype)
         idx_from  = tuple( slice(m*p,-m*p) if p != 0 else slice(0, None) for p,m in zip(self.pads, self.space.shifts) )
         idx_to    = tuple( slice(s,e+1) for s,e in zip(self.starts,self.ends) )
         a[idx_to] = self._data[idx_from]
@@ -650,7 +676,7 @@ def _toarray_parallel_with_pads(self, order='C'):
         pads = [m*p for m,p in zip(self.space.shifts, self.pads)]
         # Step 0: create extended n-dimensional array with zero values
         shape = tuple( n+2*p for n,p in zip( self.space.npts, pads ) )
-        a = xp.zeros( shape, dtype=self.dtype )
+        a = xp.zeros(shape, dtype=self.dtype)
 
         # Step 1: write extended data chunk (local to process) onto array
         idx = tuple( slice(s,e+2*p+1) for s,e,p in
@@ -882,7 +908,7 @@ def _getindex(self, key):
             index.append(l)
         return tuple(index)
 
-#===============================================================================
+#========================================================================
 class StencilMatrix(LinearOperator):
     """
     Matrix in n-dimensional stencil format.
@@ -1073,7 +1099,28 @@ def dot(self, v, out=None):
         if not v.ghost_regions_in_sync:
             v.update_ghost_regions()
 
-        self._func(self._data, v._data, out._data, **self._args)
+        # Convert arrays for compiled kernel - create NumPy output
+        import numpy as _np
+        self_data_np = _to_numpy_array(self._data)
+        v_data_np = _to_numpy_array(v._data)
+        # zeros, not empty: the compiled kernel only writes the interior
+        # (non-padding) region, so padding must be pre-initialized to avoid
+        # leaking uninitialized memory into ghost regions of the output.
+        out_data_np = _np.zeros(out._data.shape, dtype=out._data.dtype)
+
+        # Convert args that might be CuPy arrays
+        args_np = {}
+        for key, val in self._args.items():
+            args_np[key] = _to_numpy_array(val)
+
+        self._func(self_data_np, v_data_np, out_data_np, **args_np)
+        
+        # Copy result back to CuPy array if needed
+        if hasattr(out._data, 'get'):
+            import cupy as cp
+            out._data[:] = cp.asarray(out_data_np)
+        else:
+            out._data[:] = out_data_np
 
         # IMPORTANT: flag that ghost regions are not up-to-date
         out.ghost_regions_in_sync = False
@@ -1112,7 +1159,28 @@ def vdot( self, v, out=None):
         if not v.ghost_regions_in_sync:
             v.update_ghost_regions()
 
+        # Convert arrays for compiled kernel - create NumPy output
+        import numpy as _np
+        self_data_np = _to_numpy_array(self._data)
+        v_data_conj_np = _to_numpy_array(xp.conjugate(v._data))
+        # zeros, not empty: see comment in dot() above.
+        out_data_np = _np.zeros(out._data.shape, dtype=out._data.dtype)
+        
+        # Convert args that might be CuPy arrays
+        args_np = {}
+        for key, val in self._args.items():
+            args_np[key] = _to_numpy_array(val)
+        
         # Instead of computing A_*x, this function computes (A*x_)_
+        self._func(self_data_np, v_data_conj_np, out_data_np, **args_np)
+        
+        # Copy result back to CuPy array if needed
+        if hasattr(out._data, 'get'):
+            import cupy as cp
+            out_data_conj = cp.conjugate(cp.asarray(out_data_np))
+            out._data[:] = out_data_conj
+        else:
+            out._data[:] = _np.conjugate(out_data_np)
         self._func(self._data, xp.conjugate(v._data), out._data, **self._args)
         xp.conjugate(out._data, out=out._data)
 
@@ -1150,10 +1218,21 @@ def transpose(self, conjugate=False, out=None):
             out = StencilMatrix(M.codomain, M.domain, pads=self._pads, backend=self._backend, precompiled=self._precompiled)
 
         # Call low-level '_transpose' function (works on Numpy arrays directly)
+        # Convert CuPy arrays to NumPy for compiled kernels
+        M_data_np = _to_numpy_array(M._data)
+        out_data_np = _to_numpy_array(out._data)
+        
         if conjugate:
+            self._transpose_func(_to_numpy_array(xp.conjugate(M_data_np)), out_data_np, **self._transpose_args)
             self._transpose_func(xp.conjugate(M._data), out._data, **self._transpose_args)
         else:
-            self._transpose_func(M._data, out._data, **self._transpose_args)
+            self._transpose_func(M_data_np, out_data_np, **self._transpose_args)
+        
+        # Copy results back to CuPy if needed
+        if array_backend.backend == "cupy":
+            import cupy as cp
+            out._data[:] = cp.asarray(out_data_np)
+        
         return out
 
     # ...
@@ -1630,6 +1709,11 @@ def _tocoo_no_pads(self , order='C'):
 
         import numpy as _np
         # Pyccel kernels require explicit numpy.int64 type arguments
+        # Handle CuPy arrays by explicitly converting to NumPy
+        pp = []
+        for p, mi, mj in zip(self._pads, cm, dm):
+            diag_len = compute_diag_len(p, mj, mi) - (p + 1)
+            pp.append(_to_numpy_int64(diag_len))
         pp = [_np.int64(compute_diag_len(p,mj,mi)-(p+1)) for p,mi,mj in zip(self._pads, cm, dm)]
 
         # Range of data owned by local process (no ghost regions)
@@ -1640,6 +1724,31 @@ def _tocoo_no_pads(self , order='C'):
         rows = xp.zeros(size, dtype='int64')
         cols = xp.zeros(size, dtype='int64')
         data = xp.zeros(size, dtype=self.dtype)
+        nrl = [_to_numpy_int64(e-s+1) for s,e in zip(self.codomain.starts, self.codomain.ends)]
+        ncl = [_to_numpy_int64(i) for i in self._data.shape[nd:]]
+        ss = [_to_numpy_int64(i) for i in ss]
+        nr = [_to_numpy_int64(i) for i in nr]
+        nc = [_to_numpy_int64(i) for i in nc]
+        dm = [_to_numpy_int64(i) for i in dm]
+        cm = [_to_numpy_int64(i) for i in cm]
+        cpads = [_to_numpy_int64(i) for i in cpads]
+        pp = [_to_numpy_int64(i) for i in pp]
+
+        stencil2coo = kernels['stencil2coo'][order][nd]
+        # Convert CuPy arrays to NumPy for the compiled kernel
+        self_data_np = _to_numpy_array(self._data)
+        data_np = _to_numpy_array(data)
+        rows_np = _to_numpy_array(rows)
+        cols_np = _to_numpy_array(cols)
+        
+        ind = stencil2coo(self_data_np, data_np, rows_np, cols_np, *nrl, *ncl, *ss, *nr, *nc, *dm, *cm, *cpads, *pp)
+        
+        # Copy results back to CuPy arrays if needed
+        if array_backend.backend == "cupy":
+            import cupy as cp
+            data[:ind] = cp.asarray(data_np[:ind])
+            rows[:ind] = cp.asarray(rows_np[:ind])
+            cols[:ind] = cp.asarray(cols_np[:ind])
         nrl = [_np.int64(e-s+1) for s,e in zip(self.codomain.starts, self.codomain.ends)]
         ncl = [_np.int64(i) for i in self._data.shape[nd:]]
         ss = [_np.int64(i) for i in ss]
@@ -2071,7 +2180,7 @@ def _get_diagonal_indices(self):
 
         return self._diag_indices
 
-#===============================================================================
+#========================================================================
 class StencilDiagonalMatrix(LinearOperator):
     """
     Linear operator which operates between stencil vector spaces, and which can
@@ -2246,6 +2355,7 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None):
         # Calculate entries, or set `out=self` in default case
         if inverse:
             data = xp.divide(1, diag, out=data)
+        elif out:
             if sqrt:
                 data = xp.sqrt(data, out=data)
         elif sqrt:
@@ -2261,8 +2371,7 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None):
 
         return out
 
-#===============================================================================
-# TODO [YG, 28.01.2021]:
+#========================================================================# TODO [YG, 28.01.2021]:
 # - Check if StencilMatrix should be subclassed
 # - Reimplement magic methods (some are simply copied from StencilMatrix)
 def flip_axis(index, n):
@@ -3011,5 +3120,6 @@ def set_backend(self, backend, precompiled=False):
 
             self._func = dot.func
 
-#===============================================================================
+#========================================================================
+
 del VectorSpace, Vector
diff --git a/feectools/utilities/quadratures.py b/feectools/utilities/quadratures.py
index 241f8f61c..3c9708589 100644
--- a/feectools/utilities/quadratures.py
+++ b/feectools/utilities/quadratures.py
@@ -9,9 +9,8 @@
 with weights equal to 1
 """
 
-from math import cos, pi
-
 import cunumpy as xp
+from math import cos, pi
 
 
 __all__ = ('gauss_legendre', 'gauss_lobatto', 'quadrature')
diff --git a/feectools/utilities/utils.py b/feectools/utilities/utils.py
index ed8e16c1e..b0043b171 100644
--- a/feectools/utilities/utils.py
+++ b/feectools/utilities/utils.py
@@ -1,8 +1,9 @@
-#---------------------------------------------------------------------------#
-# This file is part of PSYDAC which is released under MIT License. See the  #
-# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
-# for full license details.                                                 #
-#---------------------------------------------------------------------------#
+# coding: utf-8
+#
+# Copyright 2018 Yaman Güçlü
+
+import cunumpy as xp
+import numpy as np
 from numbers import Number
 
 import cunumpy as xp
@@ -71,7 +72,13 @@ def unroll_edges(domain, xgrid):
 
     xA, xB = domain
 
-    assert all(xp.diff(xgrid) >= 0)
+    # Convert to numpy if needed (grid arrays should be on CPU)
+    if hasattr(xgrid, 'get'):
+        xgrid = xgrid.get()
+    xgrid = np.asarray(xgrid)
+
+    # Convert to numpy for comparison
+    assert all(np.diff(xgrid) >= 0)
     assert xA < xB
     assert xA <= xgrid[0]
     assert xgrid[-1] <= xB
@@ -80,10 +87,14 @@ def unroll_edges(domain, xgrid):
         return xgrid
 
     elif xgrid[0] != xA:
-        return xp.array([xgrid[-1] - (xB-xA), *xgrid])
+        # Make sure scalars are converted to Python float
+        new_point = float(xgrid[-1]) - float(xB-xA)
+        return np.concatenate([[new_point], xgrid])
 
     elif xgrid[-1] != xB:
-        return xp.array([*xgrid, xgrid[0] + (xB-xA)])
+        # Make sure scalars are converted to Python float
+        new_point = float(xgrid[0]) + float(xB-xA)
+        return np.concatenate([xgrid, [new_point]])
 
 #===============================================================================
 def roll_edges(domain, points):
@@ -92,9 +103,27 @@ def roll_edges(domain, points):
     """
     xA, xB = domain
     assert xA < xB
-    points -=xA
-    points %=(xB-xA)
-    points +=xA
+    
+    # Convert domain bounds to same backend as points to ensure compatibility
+    # First, normalize xA and xB to Python float or correct backend
+    if hasattr(xA, 'get'):
+        xA = float(xA.get())
+    elif hasattr(xA, '__array__'):
+        xA = float(xA)
+    
+    if hasattr(xB, 'get'):
+        xB = float(xB.get())
+    elif hasattr(xB, '__array__'):
+        xB = float(xB)
+    
+    # Now convert to backend of points if needed
+    if hasattr(points, 'get'):  # CuPy array
+        xA = xp.asarray(xA)
+        xB = xp.asarray(xB)
+    
+    points -= xA
+    points %= (xB - xA)
+    points += xA
 
 #===============================================================================
 def split_field(uh, spaces, out=None):
diff --git a/pyproject.toml b/pyproject.toml
index 71cfd2299..ae8f33f18 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name            = "feectools"
-version         = "0.1.7"
+version         = "0.1.8"
 description     = "Slimmed-down fork of Psydac (https://github.com/pyccel/psydac) with less functionality and fewer dependencies."
 readme          = "README.md"
 requires-python = ">= 3.10"