diff --git a/Operators/Mass/MassCUDA.hpp b/Operators/Mass/MassCUDA.hpp
index d48b50c5a5c36c27065911e8a26f65d4a6adec1d..8575a961dbf55afd673776eaf05c2eda31f8191f 100644
--- a/Operators/Mass/MassCUDA.hpp
+++ b/Operators/Mass/MassCUDA.hpp
@@ -18,8 +18,6 @@ public:
               Field<TData, FieldState::Phys>::template create<MemoryRegionCUDA>(
                   GetBlockAttributes(FieldState::Phys, expansionList)))
     {
-        auto nCoord = this->m_expansionList->GetCoordim(0);
-
         m_BwdTransOp = BwdTrans<>::create(this->m_expansionList, "CUDA");
         m_IProductWRTBaseOp =
             IProductWRTBase<>::create(this->m_expansionList, "CUDA");
@@ -44,12 +42,10 @@ public:
 
     static std::string className;
 
-private:
+protected:
     std::shared_ptr<OperatorBwdTrans<TData>> m_BwdTransOp;
     std::shared_ptr<OperatorIProductWRTBase<TData>> m_IProductWRTBaseOp;
     Field<TData, FieldState::Phys> m_field;
-    size_t m_blockSize = 32;
-    size_t m_gridSize;
 };
 
 } // namespace Nektar::Operators::detail
diff --git a/Operators/Matrix/MatrixCUDA.hpp b/Operators/Matrix/MatrixCUDA.hpp
index 40dfd5f8b30283b11efeabf631b1b57602e67b03..40910be81397e9940b1648fc4f59ffbda1eadeb0 100644
--- a/Operators/Matrix/MatrixCUDA.hpp
+++ b/Operators/Matrix/MatrixCUDA.hpp
@@ -1,3 +1,5 @@
+#pragma once
+
 #include "MemoryRegionCUDA.hpp"
 #include "Operators/Matrix/MatrixCUDAKernels.cuh"
 #include "Operators/OperatorHelper.cuh"
@@ -6,11 +8,6 @@
 namespace Nektar::Operators::detail
 {
 
-template <typename TData>
-void MatrixKernel(const size_t gridSize, const size_t blockSize,
-                  const size_t numPts, const size_t nElmts, const size_t size,
-                  const TData *mat, const TData *in, TData *out);
-
 // Matrix implementation
 template <typename TData, FieldState TFieldState>
 class OperatorMatrixImpl<TData, TFieldState, ImplCUDA>
@@ -34,7 +31,6 @@ public:
     ~OperatorMatrixImpl()
     {
         cudaFree(m_matrix);
-        m_matrix = nullptr;
     }
 
     size_t size()
@@ -99,8 +95,8 @@ public:
             m_gridSize = nElmts / m_blockSize;
             m_gridSize += (nElmts % m_blockSize == 0) ? 0 : 1;
 
-            MatrixKernel(m_gridSize, m_blockSize, numPts, nElmts, m_size,
-                         m_matrix, inptr, outptr);
+            MatrixKernel<<<m_gridSize, m_blockSize>>>(numPts, nElmts, m_size,
+                                                      m_matrix, inptr, outptr);
 
             // Increment pointer and index for next element type.
             inptr += numPts * nElmts;
@@ -126,13 +122,4 @@ private:
     size_t m_blockSize = 32;
 };
 
-template <typename TData>
-void MatrixKernel(const size_t gridSize, const size_t blockSize,
-                  const size_t numPts, const size_t nElmts, const size_t size,
-                  const TData *mat, const TData *in, TData *out)
-{
-    MatrixKernel<TData>
-        <<<gridSize, blockSize>>>(numPts, nElmts, size, mat, in, out);
-}
-
 } // namespace Nektar::Operators::detail
diff --git a/Operators/Matrix/MatrixStdMat.hpp b/Operators/Matrix/MatrixStdMat.hpp
index 10c10f776ad07b4b6a9a2978532dc268c32af7d2..fd045a139de3342725294bd4e18bcc86e6f44407 100644
--- a/Operators/Matrix/MatrixStdMat.hpp
+++ b/Operators/Matrix/MatrixStdMat.hpp
@@ -1,8 +1,6 @@
 #pragma once
 
 #include <numeric>
-
-#include "Field.hpp"
 #include "Operators/OperatorMatrix.hpp"
 
 namespace Nektar::Operators::detail