diff --git a/src/struphy/feec/banded_to_stencil_kernels.py b/src/struphy/feec/banded_to_stencil_kernels.py
index 6fa244be1..118caecda 100644
--- a/src/struphy/feec/banded_to_stencil_kernels.py
+++ b/src/struphy/feec/banded_to_stencil_kernels.py
@@ -1,27 +1,29 @@
+import numpy as np
 from numpy import mod, shape
 from pyccel.decorators import pure
 
 
 @pure
-def band_to_stencil_1d(arr: "float[:, :]", out: "float[:, :]"):
+def band_to_stencil_1d(arr: "float[:, :]", p: "int") -> "float[:, :]":
     """Converts the 2d banded matrix arr, of shape (n, m) and with 2*p + 1 < m bands centered around the diagonal,
-    into the array out, of shape (n, 2*p + 1), which corresponds to the StencilMatrix format.
+    into an array of shape (n, 2*p + 1), which corresponds to the StencilMatrix format.
     """
     s = shape(arr)
-    p = shape(out)[1] // 2
+    out = np.zeros((s[0], 2 * p + 1), dtype=float)
 
     for i in range(s[0]):
         for j in range(2 * p + 1):
             jj = mod(i - p + j, s[1])
             out[i, j] = arr[i, jj]
 
+    return out
+
 
 @pure
-def band_to_stencil_2d(arr: "float[:, :, :, :]", out: "float[:, :, :, :]"):
+def band_to_stencil_2d(arr: "float[:, :, :, :]", p1: "int", p2: "int") -> "float[:, :, :, :]":
     """Converts a 4d banded matrix to StencilMatrix format (see band_to_stencil_1d)."""
     s = shape(arr)
-    p1 = shape(out)[2] // 2
-    p2 = shape(out)[3] // 2
+    out = np.zeros((s[0], s[1], 2 * p1 + 1, 2 * p2 + 1), dtype=float)
 
     for i1 in range(s[0]):
         for j1 in range(2 * p1 + 1):
@@ -32,14 +34,14 @@ def band_to_stencil_2d(arr: "float[:, :, :, :]", out: "float[:, :, :, :]"):
                     jj2 = mod(i2 - p2 + j2, s[3])
                     out[i1, i2, j1, j2] = arr[i1, i2, jj1, jj2]
 
+    return out
+
 
 @pure
-def band_to_stencil_3d(arr: "float[:, :, :, :, :, :]", out: "float[:, :, :, :, :, :]"):
+def band_to_stencil_3d(arr: "float[:, :, :, :, :, :]", p1: "int", p2: "int", p3: "int") -> "float[:, :, :, :, :, :]":
     """Converts a 6d banded matrix to StencilMatrix format (see band_to_stencil_1d)."""
     s = shape(arr)
-    p1 = shape(out)[3] // 2
-    p2 = shape(out)[4] // 2
-    p3 = shape(out)[5] // 2
+    out = np.zeros((s[0], s[1], s[2], 2 * p1 + 1, 2 * p2 + 1, 2 * p3 + 1), dtype=float)
 
     for i1 in range(s[0]):
         for j1 in range(2 * p1 + 1):
@@ -53,3 +55,5 @@ def band_to_stencil_3d(arr: "float[:, :, :, :, :, :]", out: "float[:, :, :, :, :
                         for j3 in range(2 * p3 + 1):
                             jj3 = mod(i3 - p3 + j3, s[5])
                             out[i1, i2, i3, j1, j2, j3] = arr[i1, i2, i3, jj1, jj2, jj3]
+
+    return out
diff --git a/src/struphy/feec/local_projectors_kernels.py b/src/struphy/feec/local_projectors_kernels.py
index 8d8dd5f34..75a45963c 100644
--- a/src/struphy/feec/local_projectors_kernels.py
+++ b/src/struphy/feec/local_projectors_kernels.py
@@ -5,7 +5,7 @@
 from struphy.kernel_arguments.local_projectors_args_kernels import LocalProjectorsArguments
 
 
-def compute_shifts(IoH: "bool[:]", p: "int[:]", B_nbasis: "int[:]", shift: "int[:]"):
+def compute_shifts(IoH: "bool[:]", p: "int[:]", B_nbasis: "int[:]") -> "int[:]":
     """This function computes by how much we must shift the indices in case we loop over the evaluation points.
 
     Parameters
@@ -19,9 +19,13 @@ def compute_shifts(IoH: "bool[:]", p: "int[:]", B_nbasis: "int[:]", shift: "int[
     B_nbasis: 1d int array
         Array with the number of B-splines in each direction.
 
-    shifts : 1d int array
+    Returns
+    -------
+    shift : 1d int array
         array of 3 ints, each one denotes the amout by which we must shift the indices to loop around the quasi-points for each spatial direction.
     """
+    shift = zeros(3, dtype=int)
+
     for i, ioh in enumerate(IoH):
         # Histopolation
         if ioh:
@@ -33,6 +37,8 @@ def compute_shifts(IoH: "bool[:]", p: "int[:]", B_nbasis: "int[:]", shift: "int[
         else:
             shift[i] = -2 * B_nbasis[i]
 
+    return shift
+
 
 def get_local_problem_size(periodic: "bool[:]", p: "int[:]", IoH: "bool[:]"):
     """Determines the number of interpolation or histopolation weights present for a fixed index i.
@@ -82,9 +88,11 @@ def get_dofs_local_1_form_ec_component_weighted(
     basis1: "float[:]",
     basis2: "float[:]",
     arezeroc: "int[:]",
-    f_eval_aux: "float[:,:,:]",
+    n1: int,
+    n2: int,
+    n3: int,
     c: int,
-):
+) -> "float[:,:,:]":
     """Kernel for evaluating the degrees of freedom for the c-th component of 1-forms. This function is for local commuting projetors.
 
     Parameters
@@ -107,14 +115,21 @@ def get_dofs_local_1_form_ec_component_weighted(
         arezeroc : 1d int array
             Array of zeros or ones. A one means that for this particular set of quadrature points, in the c-th direction, the basis function is not zero for at least one of them.
 
-        f_eval_aux : 3d float array
-            Output array where the evaluated degrees of freedom are stored. It is passed to this function with zeros in each entry.
+        n1, n2, n3 : int
+            Shape of the output array of evaluated degrees of freedom.
 
         c : int
             This int tell us whichone of the three components of the 1-form vector we are dealing with. It must be 0, 1 or 2.
+
+    Returns
+    -------
+    f_eval_aux : 3d float array
+        The evaluated degrees of freedom.
     """
     p = args_solve.degree[c]
 
+    f_eval_aux = zeros((n1, n2, n3), dtype=float)
+
     for i in range(shape(f_eval_aux)[0]):
         if c == 0:
             computei = arezeroc[i] != 0
@@ -165,14 +180,18 @@ def get_dofs_local_1_form_ec_component_weighted(
                                         * args_solve.wts2[args_solve.inv_index_translation2[k], ii]
                                     )
 
+    return f_eval_aux
+
 
 @stack_array("shp")
 def get_dofs_local_1_form_ec_component(
     args_solve: LocalProjectorsArguments,
     f3: "float[:,:,:]",
-    f_eval_aux: "float[:,:,:]",
+    n1: int,
+    n2: int,
+    n3: int,
     c: int,
-):
+) -> "float[:,:,:]":
     """Kernel for evaluating the degrees of freedom for the c-th component of 1-forms.  This function is for local commuting projetors.
 
     Parameters
@@ -183,13 +202,20 @@ def get_dofs_local_1_form_ec_component(
         f3 : 3d float array
             Evaluation for the c-th component of the 1-form function over all the interpolation points in e_a and e_b (a!=b!=c), as well as all the Gauss-Legendre quadrature point in e_c.
 
-        f_eval_aux : 3d float array
-            Output array where the evaluated degrees of freedom are stored. It is passed to this function with zeros in each entry.
+        n1, n2, n3 : int
+            Shape of the output array of evaluated degrees of freedom.
 
         c : int
             This integer determines which of the three components of the 1-form vector we are working on. Must be 0,1 or 2.
+
+    Returns
+    -------
+    f_eval_aux : 3d float array
+        The evaluated degrees of freedom.
     """
 
+    f_eval_aux = zeros((n1, n2, n3), dtype=float)
+
     shp = zeros(3, dtype=int)
     shp[:] = shape(f_eval_aux)
 
@@ -220,14 +246,18 @@ def get_dofs_local_1_form_ec_component(
                     for ii in range(p):
                         f_eval_aux[i, j, k] += f3[i, j, in_start + ii] * wts[inv_index_translation[k], ii]
 
+    return f_eval_aux
+
 
 @stack_array("shp")
 def get_dofs_local_2_form_ec_component(
     args_solve: LocalProjectorsArguments,
     fc: "float[:,:,:]",
-    f_eval_aux: "float[:,:,:]",
+    n1: int,
+    n2: int,
+    n3: int,
     c: int,
-):
+) -> "float[:,:,:]":
     """Kernel for evaluating the degrees of freedom for the c-th component of 2-forms.  This function is for local commuting projetors.
 
     Parameters
@@ -239,13 +269,20 @@ def get_dofs_local_2_form_ec_component(
             Evaluation for the c-th component of the 2-form function over all the interpolation points in e_c, as well as all the Gauss-Legendre quadrature point in e_a and e_b.
             Of the two spatial directions different from e_c, e_a is the one with the smaller index, and e_b is the one with the larger index.
 
-        f_eval_aux : 3d float array
-            Output array where the evaluated degrees of freedom are stored. It is passed to this function with zeros in each entry.
+        n1, n2, n3 : int
+            Shape of the output array of evaluated degrees of freedom.
 
         c : int
             This integer determines which of the three components of the 2-form vector we are working on. Must be 0, 1 or 2.
+
+    Returns
+    -------
+    f_eval_aux : 3d float array
+        The evaluated degrees of freedom.
     """
 
+    f_eval_aux = zeros((n1, n2, n3), dtype=float)
+
     shp = zeros(3, dtype=int)
     shp[:] = shape(f_eval_aux)
 
@@ -283,6 +320,8 @@ def get_dofs_local_2_form_ec_component(
                                 * args_solve.wts1[args_solve.inv_index_translation1[j], kk]
                             )
 
+    return f_eval_aux
+
 
 def get_dofs_local_2_form_ec_component_weighted(
     args_solve: LocalProjectorsArguments,
@@ -292,9 +331,11 @@ def get_dofs_local_2_form_ec_component_weighted(
     basis2: "float[:]",
     arezero_a: "int[:]",
     arezero_b: "int[:]",
-    f_eval_aux: "float[:,:,:]",
+    n1: int,
+    n2: int,
+    n3: int,
     c: int,
-):
+) -> "float[:,:,:]":
     """Kernel for evaluating the degrees of freedom for the c-th component of 2-forms.  This function is for local commuting projetors.
 
     Parameters
@@ -321,13 +362,20 @@ def get_dofs_local_2_form_ec_component_weighted(
         arezero_b : 1d int array
             Array zeros or ones. A one means that for this particular set of quadrature points, in the e_b direction, the basis function is not zero for at least one of them.
 
-        f_eval_aux : 3d float array
-            Output array where the evaluated degrees of freedom are stored. It is passed to this function with zeros in each entry.
+        n1, n2, n3 : int
+            Shape of the output array of evaluated degrees of freedom.
 
         c : int
             This int tell us whichone of the three components of the 2-form vector we are dealing with. It must be 0, 1 or 2.
+
+    Returns
+    -------
+    f_eval_aux : 3d float array
+        The evaluated degrees of freedom.
     """
 
+    f_eval_aux = zeros((n1, n2, n3), dtype=float)
+
     for i in range(shape(f_eval_aux)[0]):
         if c == 0:
             computei = abs(basis0[i]) >= 10.0 ** (-16)
@@ -390,9 +438,17 @@ def get_dofs_local_2_form_ec_component_weighted(
                                             * args_solve.wts1[args_solve.inv_index_translation1[j], kk]
                                         )
 
+    return f_eval_aux
+
 
 @stack_array("shp")
-def get_dofs_local_3_form(args_solve: LocalProjectorsArguments, faux: "float[:,:,:]", f_eval: "float[:,:,:]"):
+def get_dofs_local_3_form(
+    args_solve: LocalProjectorsArguments,
+    faux: "float[:,:,:]",
+    n1: int,
+    n2: int,
+    n3: int,
+) -> "float[:,:,:]":
     """Kernel for evaluating the degrees of freedom for 3-forms.  This function is for local commuting projetors.
 
     Parameters
@@ -403,9 +459,16 @@ def get_dofs_local_3_form(args_solve: LocalProjectorsArguments, faux: "float[:,:
         faux : 3d float array
             Evaluation for the 3-form function over all the Gauss-Legendre quadrature point in e1, e2 and e3.
 
-        f_eval : 3d float array
-            Output array where the evaluated degrees of freedom are stored. It is passed to this function with zeros in each entry.
+        n1, n2, n3 : int
+            Shape of the output array of evaluated degrees of freedom.
+
+    Returns
+    -------
+    f_eval : 3d float array
+        The evaluated degrees of freedom.
     """
+    f_eval = zeros((n1, n2, n3), dtype=float)
+
     shp = zeros(3, dtype=int)
     shp[:] = shape(f_eval)
 
@@ -429,6 +492,8 @@ def get_dofs_local_3_form(args_solve: LocalProjectorsArguments, faux: "float[:,:
                                 * args_solve.wts2[args_solve.inv_index_translation2[k], kk]
                             )
 
+    return f_eval
+
 
 def get_dofs_local_3_form_weighted(
     args_solve: LocalProjectorsArguments,
@@ -439,8 +504,10 @@ def get_dofs_local_3_form_weighted(
     arezero0: "int[:]",
     arezero1: "int[:]",
     arezero2: "int[:]",
-    f_eval: "float[:,:,:]",
-):
+    n1: int,
+    n2: int,
+    n3: int,
+) -> "float[:,:,:]":
     """Kernel for evaluating the degrees of freedom for 3-forms.  This function is for local commuting projetors.
 
     Parameters
@@ -469,10 +536,17 @@ def get_dofs_local_3_form_weighted(
         arezero2 : 1d int array
             Array of zeros or ones. A one means that for this particular set of quadrature points, in the third direction, the basis function is not zero for at least one of them.
 
-        f_eval : 3d float array
-            Output array where the evaluated degrees of freedom are stored. It is passed to this function with zeros in each entry.
+        n1, n2, n3 : int
+            Shape of the output array of evaluated degrees of freedom.
+
+    Returns
+    -------
+    f_eval : 3d float array
+        The evaluated degrees of freedom.
     """
 
+    f_eval = zeros((n1, n2, n3), dtype=float)
+
     for i in range(shape(f_eval)[0]):
         if arezero0[i] != 0:
             for j in range(shape(f_eval)[1]):
@@ -499,6 +573,8 @@ def get_dofs_local_3_form_weighted(
                                             * args_solve.wts2[args_solve.inv_index_translation2[k], kk]
                                         )
 
+    return f_eval
+
 
 # We need a functions that tell us which of the quasi-interpolation points to take for a any given i
 def select_quasi_points(i: int, p: int, Nbasis: int, periodic: bool):
diff --git a/src/struphy/feec/projectors.py b/src/struphy/feec/projectors.py
index 6b015b34f..e795697af 100644
--- a/src/struphy/feec/projectors.py
+++ b/src/struphy/feec/projectors.py
@@ -1030,8 +1030,7 @@ def __init__(
 
             lenj = [lenj1, lenj2, lenj3]
 
-            self._shift = xp.array([0, 0, 0], dtype=int)
-            compute_shifts(self._IoH, self._degree, self._B_nbasis, self._shift)
+            self._shift = compute_shifts(self._IoH, self._degree, self._B_nbasis)
 
             split_points(
                 IoH_for_indices,
@@ -1252,7 +1251,7 @@ def __init__(
 
                 lenj = [lenj1, lenj2, lenj3]
 
-                compute_shifts(self._IoH[h], self._degree, self._B_nbasis, self._shift[h])
+                self._shift[h] = compute_shifts(self._IoH[h], self._degree, self._B_nbasis)
 
                 split_points(
                     IoH_for_indices[h],
@@ -1693,23 +1692,23 @@ def get_dofs(self, fun, dofs=None):
                     # Evaluation of the function to compute the h component
                     fh = fun[h](*self._meshgrid[h])
 
-                # Array into which we will write the Dofs.
-                f_eval_aux = xp.zeros(tuple(xp.shape(dim)[0] for dim in self._localpts[h]))
+                # Shape of the array into which we will write the Dofs.
+                n1, n2, n3 = (xp.shape(dim)[0] for dim in self._localpts[h])
 
                 # For 1-forms
                 if self._space_key == "1":
-                    get_dofs_local_1_form_ec_component(self._solve_args[h], fh, f_eval_aux, h)
+                    f_eval_aux = get_dofs_local_1_form_ec_component(self._solve_args[h], fh, n1, n2, n3, h)
                 # For 2-forms
                 else:
-                    get_dofs_local_2_form_ec_component(self._solve_args[h], fh, f_eval_aux, h)
+                    f_eval_aux = get_dofs_local_2_form_ec_component(self._solve_args[h], fh, n1, n2, n3, h)
 
                 f_eval.append(f_eval_aux)
 
         elif self._space_key == "3":
-            f_eval = xp.zeros(tuple(xp.shape(dim)[0] for dim in self._localpts))
+            n1, n2, n3 = (xp.shape(dim)[0] for dim in self._localpts)
             # Evaluation of the function at all Gauss-Legendre quadrature points
             faux = fun(*self._meshgrid)
-            get_dofs_local_3_form(self._solve_args, faux, f_eval)
+            f_eval = get_dofs_local_3_form(self._solve_args, faux, n1, n2, n3)
 
         elif self._space_key == "v":
             f_eval = []
@@ -1759,8 +1758,8 @@ def get_dofs_weighted(self, fun, dofs=None, first_go=True, pre_computed_dofs=Non
                 if first_go:
                     pre_computed_dofs.append(fun[h](*self._meshgrid[h]))
 
-                # Array into which we will write the Dofs.
-                f_eval_aux = xp.zeros(tuple(xp.shape(dim)[0] for dim in self._localpts[h]))
+                # Shape of the array into which we will write the Dofs.
+                n1, n2, n3 = (xp.shape(dim)[0] for dim in self._localpts[h])
 
                 # We check if the current set of basis functions is not one of those we have to compute in the current MPI rank.
                 if (
@@ -1770,8 +1769,9 @@ def get_dofs_weighted(self, fun, dofs=None, first_go=True, pre_computed_dofs=Non
                 ):
                     # We should do nothing here
                     self._do_nothing[h] = 1
+                    f_eval_aux = xp.zeros((n1, n2, n3))
                 elif self._space_key == "1":
-                    get_dofs_local_1_form_ec_component_weighted(
+                    f_eval_aux = get_dofs_local_1_form_ec_component_weighted(
                         self._solve_args[h],
                         pre_computed_dofs[h],
                         self.get_values(
@@ -1781,13 +1781,15 @@ def get_dofs_weighted(self, fun, dofs=None, first_go=True, pre_computed_dofs=Non
                         self.get_values(1, h),
                         self.get_values(2, h),
                         self.get_are_zero(h, h),
-                        f_eval_aux,
+                        n1,
+                        n2,
+                        n3,
                         h,
                     )
                 else:
                     # ind1 and ind2 are the indices of the two directions with histopolation, ind1 must be smaller than ind2.
                     (ind1, ind2) = [(1, 2), (0, 2), (0, 1)][h]
-                    get_dofs_local_2_form_ec_component_weighted(
+                    f_eval_aux = get_dofs_local_2_form_ec_component_weighted(
                         self._solve_args[h],
                         pre_computed_dofs[h],
                         self.get_values(
@@ -1798,19 +1800,21 @@ def get_dofs_weighted(self, fun, dofs=None, first_go=True, pre_computed_dofs=Non
                         self.get_values(2, h),
                         self.get_are_zero(ind1, h),
                         self.get_are_zero(ind2, h),
-                        f_eval_aux,
+                        n1,
+                        n2,
+                        n3,
                         h,
                     )
 
                 f_eval.append(f_eval_aux)
 
         elif self._space_key == "3":
-            f_eval = xp.zeros(tuple(xp.shape(dim)[0] for dim in self._localpts))
+            n1, n2, n3 = (xp.shape(dim)[0] for dim in self._localpts)
             # Evaluation of the function at all Gauss-Legendre quadrature points
             if first_go:
                 pre_computed_dofs = [fun(*self._meshgrid)]
 
-            get_dofs_local_3_form_weighted(
+            f_eval = get_dofs_local_3_form_weighted(
                 self._solve_args,
                 pre_computed_dofs[0],
                 self.get_values(0),
@@ -1821,7 +1825,9 @@ def get_dofs_weighted(self, fun, dofs=None, first_go=True, pre_computed_dofs=Non
                 self.get_are_zero(0),
                 self.get_are_zero(1),
                 self.get_are_zero(2),
-                f_eval,
+                n1,
+                n2,
+                n3,
             )
 
         elif self._space_key == "v":
diff --git a/src/struphy/feec/utilities.py b/src/struphy/feec/utilities.py
index 88784f85a..b63915307 100644
--- a/src/struphy/feec/utilities.py
+++ b/src/struphy/feec/utilities.py
@@ -294,18 +294,7 @@ def compare_arrays(arr_psy, arr, rank, atol=1e-14):
         if tmp_arr.shape == tmp1.shape:
             tmp2 = tmp_arr
         else:
-            tmp2 = xp.zeros(
-                (
-                    e[0] + 1 - s[0],
-                    e[1] + 1 - s[1],
-                    e[2] + 1 - s[2],
-                    2 * degree[0] + 1,
-                    2 * degree[1] + 1,
-                    2 * degree[2] + 1,
-                ),
-                dtype=float,
-            )
-            bts.band_to_stencil_3d(tmp_arr, tmp2)
+            tmp2 = bts.band_to_stencil_3d(tmp_arr, degree[0], degree[1], degree[2])
 
         assert xp.allclose(tmp1, tmp2, atol=atol)
 
@@ -338,18 +327,7 @@ def compare_arrays(arr_psy, arr, rank, atol=1e-14):
                 if tmp_mat.shape == tmp1.shape:
                     tmp2 = tmp_mat
                 else:
-                    tmp2 = xp.zeros(
-                        (
-                            e[0] + 1 - s[0],
-                            e[1] + 1 - s[1],
-                            e[2] + 1 - s[2],
-                            2 * degree[0] + 1,
-                            2 * degree[1] + 1,
-                            2 * degree[2] + 1,
-                        ),
-                        dtype=float,
-                    )
-                    bts.band_to_stencil_3d(tmp_mat, tmp2)
+                    tmp2 = bts.band_to_stencil_3d(tmp_mat, degree[0], degree[1], degree[2])
 
                 assert xp.allclose(tmp1, tmp2, atol=atol)
 
diff --git a/src/struphy/geometry/base.py b/src/struphy/geometry/base.py
index 557bc54e3..3e5d0db12 100644
--- a/src/struphy/geometry/base.py
+++ b/src/struphy/geometry/base.py
@@ -977,18 +977,12 @@ def _evaluate_metric_coefficient(self, *etas, which=0, **kwargs):
                 flat_eval=False,
             )
 
-            # to keep C-ordering the (3, 3)-part is in the last indices
-            out = xp.empty(
-                (E1.shape[0], E2.shape[1], E3.shape[2], 3, 3),
-                dtype=float,
-            )
-            evaluation_kernels.kernel_evaluate(
+            out = evaluation_kernels.kernel_evaluate(
                 E1,
                 E2,
                 E3,
                 which,
                 self.args_domain,
-                out,
                 is_sparse_meshgrid,
                 avoid_round_off,
             )
@@ -1180,12 +1174,7 @@ def _pull_push_transform(self, which, a, kind_fun, *etas, flat_eval=False, **kwa
                 X = self(E1, E2, E3)
                 A = Domain.prepare_arg(a, X[0], X[1], X[2], a_kwargs=a_kwargs)
 
-            # call evaluation kernel
-            out = xp.empty(
-                (E1.shape[0], E2.shape[1], E3.shape[2], 3),
-                dtype=float,
-            )
-            transform_kernels.kernel_pullpush(
+            out = transform_kernels.kernel_pullpush(
                 A,
                 E1,
                 E2,
@@ -1194,7 +1183,6 @@ def _pull_push_transform(self, which, a, kind_fun, *etas, flat_eval=False, **kwa
                 kind_int,
                 self.args_domain,
                 is_sparse_meshgrid,
-                out,
             )
 
             # move the (3, 3)-part to front
diff --git a/src/struphy/geometry/evaluation_kernels.py b/src/struphy/geometry/evaluation_kernels.py
index 4a2ae9bbe..12db8b3e6 100644
--- a/src/struphy/geometry/evaluation_kernels.py
+++ b/src/struphy/geometry/evaluation_kernels.py
@@ -828,10 +828,9 @@ def kernel_evaluate(
     eta3: "float[:,:,:]",
     kind_coeff: int,
     args: "DomainArguments",
-    mat_f: "float[:,:,:,:,:]",
     is_sparse_meshgrid: bool,
     avoid_round_off: bool,
-):
+) -> "float[:,:,:,:,:]":
     """
     Evaluation of metric coefficients on a given 3d grid of evaluation points.
 
@@ -839,6 +838,11 @@ def kernel_evaluate(
     ----------
     is_sparse_meshgrid : bool
         Whether the 3d evaluation points were obtained from a sparse meshgrid.
+
+    Returns
+    -------
+    mat_f : xp.ndarray
+        Evaluated metric coefficients, of shape (n1, n2, n3, 3, 3).
     """
     tmp0 = zeros(3, dtype=float)
     tmp1 = zeros((3, 3), dtype=float)
@@ -850,6 +854,8 @@ def kernel_evaluate(
     n2 = shape(eta2)[1]
     n3 = shape(eta3)[2]
 
+    mat_f = zeros((n1, n2, n3, 3, 3), dtype=float)
+
     if is_sparse_meshgrid:
         sparse_factor = 0
     else:
@@ -862,8 +868,6 @@ def kernel_evaluate(
                 e2 = eta2[i1 * sparse_factor, i2, i3 * sparse_factor]
                 e3 = eta3[i1 * sparse_factor, i2 * sparse_factor, i3]
 
-                out[:] = mat_f[i1, i2, i3, :, :]
-
                 select_metric_coeff(
                     e1,
                     e2,
@@ -880,6 +884,8 @@ def kernel_evaluate(
 
                 mat_f[i1, i2, i3, :, :] = out
 
+    return mat_f
+
 
 @stack_array("tmp0", "tmp1", "tmp2", "tmp3", "out")
 def kernel_evaluate_pic(
diff --git a/src/struphy/geometry/transform_kernels.py b/src/struphy/geometry/transform_kernels.py
index f9e6d8077..d97ddf88f 100644
--- a/src/struphy/geometry/transform_kernels.py
+++ b/src/struphy/geometry/transform_kernels.py
@@ -301,8 +301,7 @@ def kernel_pullpush(
     kind_fun: int,
     args_domain: "DomainArguments",
     is_sparse_meshgrid: bool,
-    out: "float[:,:,:,:]",
-):
+) -> "float[:,:,:,:]":
     """
     Pull-backs, pushforwards and transformations on a given 3d grid of evaluation points.
 
@@ -326,19 +325,21 @@ def kernel_pullpush(
     is_sparse_meshgrid : bool
         Whether the evaluation points were obtained from a sparse meshgrid.
 
+    Returns
+    -------
     out : float[:,:,:,:]
         Output values.
     """
 
     tmp1 = zeros(shape(a)[-1], dtype=float)
-    tmp2 = zeros(shape(out)[-1], dtype=float)
-    # tmp1 = zeros(3, dtype=float)
-    # tmp2 = zeros(3, dtype=float)
+    tmp2 = zeros(3, dtype=float)
 
     n1 = shape(eta1)[0]
     n2 = shape(eta2)[1]
     n3 = shape(eta3)[2]
 
+    out = zeros((n1, n2, n3, 3), dtype=float)
+
     if is_sparse_meshgrid:
         sparse_factor = 0
     else:
@@ -352,7 +353,6 @@ def kernel_pullpush(
                 e3 = eta3[i1 * sparse_factor, i2 * sparse_factor, i3]
 
                 tmp1[:] = a[i1, i2, i3, :]
-                tmp2[:] = out[i1, i2, i3, :]
 
                 if kind_transform == 0:
                     pull(tmp1, e1, e2, e3, kind_fun, args_domain, tmp2)
@@ -363,6 +363,8 @@ def kernel_pullpush(
 
                 out[i1, i2, i3, :] = tmp2
 
+    return out
+
 
 @stack_array("tmp1", "tmp2")
 def kernel_pullpush_pic(
@@ -403,8 +405,6 @@ def kernel_pullpush_pic(
 
     tmp1 = zeros(shape(a)[1], dtype=float)
     tmp2 = zeros(shape(out)[1], dtype=float)
-    # tmp1 = zeros((3,), dtype=float)
-    # tmp2 = zeros((3,), dtype=float)
 
     np = shape(markers)[0]
 
diff --git a/src/struphy/geometry/utilities.py b/src/struphy/geometry/utilities.py
index 625b384d2..753375230 100644
--- a/src/struphy/geometry/utilities.py
+++ b/src/struphy/geometry/utilities.py
@@ -290,8 +290,7 @@ def f_angles(xis, s_val):
             gp = xp.sqrt(psi(_R, _Z, dR=1) ** 2 + psi(_R, _Z, dZ=1) ** 2)
 
             # compute weighted arc_lengths between two successive points in xis_extended array
-            dl = xp.zeros(xis_extended.size - 1, dtype=float)
-            weighted_arc_lengths_flux_surface(_R, _Z, gp, dl, xi_param_dict[xi_param])
+            dl = weighted_arc_lengths_flux_surface(_R, _Z, gp, xi_param_dict[xi_param])
 
             # total length of the flux surface
             l = xp.sum(dl)
diff --git a/src/struphy/geometry/utilities_kernels.py b/src/struphy/geometry/utilities_kernels.py
index 1c26b25c5..6214be74b 100644
--- a/src/struphy/geometry/utilities_kernels.py
+++ b/src/struphy/geometry/utilities_kernels.py
@@ -1,7 +1,8 @@
+import numpy as np
 from numpy import shape, sqrt
 
 
-def weighted_arc_lengths_flux_surface(r: "float[:]", z: "float[:]", grad_psi: "float[:]", dwls: "float[:]", kind: int):
+def weighted_arc_lengths_flux_surface(r: "float[:]", z: "float[:]", grad_psi: "float[:]", kind: int) -> "float[:]":
     """
     Computes the weighted arc lengths
 
@@ -27,12 +28,14 @@ def weighted_arc_lengths_flux_surface(r: "float[:]", z: "float[:]", grad_psi: "f
     grad_psi : xp.ndarray
         Absolute values of the flux function gradient on the flux surface: |grad(psi)| = sqrt[ (d_R psi)**2 + (d_Z psi)**2 ].
 
-    dwls : xp.ndarray
-        The weighted arc lengths will be written into this array. Length must be one smaller than lengths of r, z and grad_psi.
-
     kind : int
         Which weight to use (see above table: 1: equal arc length, 2: straight field line, etc.)
 
+    Returns
+    -------
+    dwls : xp.ndarray
+        The weighted arc lengths. Length is one smaller than lengths of r, z and grad_psi.
+
     References
     ----------
     [1] Jardin Stephen, Computational Methods in Plasma Physics, Taylor and Francis Group 2010.
@@ -41,6 +44,8 @@ def weighted_arc_lengths_flux_surface(r: "float[:]", z: "float[:]", grad_psi: "f
     # number of angle boundaries
     n_th = r.size
 
+    dwls = np.zeros(n_th - 1, dtype=float)
+
     for j in range(n_th - 1):
         # local orthonormal coordinate system at line segment (j --> j+1)
         er_1 = r[j + 1] - r[j]
@@ -96,3 +101,5 @@ def weighted_arc_lengths_flux_surface(r: "float[:]", z: "float[:]", grad_psi: "f
         # h = 1 (constant volume)
         elif kind == 4:
             dwls[j] = dls * 1 / 2 * (r[j] / grad_psi[j] + r[j + 1] / grad_psi[j + 1])
+
+    return dwls