diff --git a/tests/init_assmbscatrfields.hpp b/tests/init_assmbscatrfields.hpp
index 9d38110a11230faaabcc340cea65df44aa0b2147..46aa1edd82c1f72ac5227a0c2096f3d113887d53 100644
--- a/tests/init_assmbscatrfields.hpp
+++ b/tests/init_assmbscatrfields.hpp
@@ -21,7 +21,7 @@ public:
     void SetTestCase(const std::vector<BlockAttributes> &blocks, double *inptr,
                      bool padding = true)
     {
-        for (auto const &block : fixt_in->GetBlocks())
+        for (auto const &block : blocks)
         {
             for (size_t el = 0; el < block.num_elements; ++el)
             {
@@ -50,7 +50,7 @@ public:
         Array<OneD, NekDouble> outcoeff(fixt_explist->GetNcoeffs());
 
         // Set test case
-        SetTestCase(blocks, incoeff.get(), false);
+        SetTestCase(fixt_in->GetBlocks(), incoeff.get(), false);
 
         // Calculate expected result from Nektar++
         auto map = fixt_explist->GetLocalToGlobalMap();
@@ -58,7 +58,7 @@ public:
         // Vmath::Zero(map->GetNumGlobalDirBndCoeffs(), outcoeff, 1);
         map->GlobalToLocal(outcoeff, outcoeff);
 
-        // Copy expected result from Array to fixt_expected
+        // Copy expected result from Array to pointer
         double *ptr = outcoeff.get();
         for (auto const &block : blocks)
         {
diff --git a/tests/init_bwdtransfields.hpp b/tests/init_bwdtransfields.hpp
index d18432858afe9f327c2b0c84adfae215c9f5ecde..b7c270409dbd7d1433f1b980c7ab123c12208268 100644
--- a/tests/init_bwdtransfields.hpp
+++ b/tests/init_bwdtransfields.hpp
@@ -16,7 +16,7 @@ public:
     void SetTestCase(const std::vector<BlockAttributes> &blocks, double *inptr,
                      bool padding = true)
     {
-        for (auto const &block : fixt_in->GetBlocks())
+        for (auto const &block : blocks)
         {
             for (size_t el = 0; el < block.num_elements; ++el)
             {
@@ -45,12 +45,12 @@ public:
         Array<OneD, NekDouble> outphys(fixt_explist->GetTotPoints());
 
         // Set test case
-        SetTestCase(blocks, incoeffs.get(), false);
+        SetTestCase(fixt_in->GetBlocks(), incoeffs.get(), false);
 
         // Calculate expected result from Nektar++
         fixt_explist->BwdTrans(incoeffs, outphys);
 
-        // Copy expected result from Array to fixt_expected
+        // Copy expected result from Array to pointer
         double *ptr = outphys.get();
         for (auto const &block : blocks)
         {
diff --git a/tests/init_diagpreconfields.hpp b/tests/init_diagpreconfields.hpp
index 5de28e9f60c3db355f34d4b5927f293fd73a8986..6a96e8d1bb9d60793ae2589cda8d643e3739f90b 100644
--- a/tests/init_diagpreconfields.hpp
+++ b/tests/init_diagpreconfields.hpp
@@ -22,7 +22,7 @@ public:
     void SetTestCase(const std::vector<BlockAttributes> &blocks, double *inptr,
                      bool padding = true)
     {
-        for (auto const &block : fixt_in->GetBlocks())
+        for (auto const &block : blocks)
         {
             for (size_t el = 0; el < block.num_elements; ++el)
             {
@@ -51,7 +51,7 @@ public:
         Array<OneD, NekDouble> outcoeffs(fixt_explist->GetNcoeffs(), 0.0);
 
         // Set test case
-        SetTestCase(blocks, incoeffs.get(), false);
+        SetTestCase(fixt_in->GetBlocks(), incoeffs.get(), false);
 
         // Calculate expected result from Nektar++
         StdRegions::ConstFactorMap factors;
@@ -68,7 +68,7 @@ public:
         precond->BuildPreconditioner();
         precond->DoPreconditioner(incoeffs, outcoeffs, true);
 
-        // Copy expected result from Array to fixt_expected
+        // Copy expected result from Array to pointer
         double *coeffptr = outcoeffs.get();
         for (auto const &block : blocks)
         {
diff --git a/tests/init_dirichletfields.hpp b/tests/init_dirichletfields.hpp
index c6fa17311442e8e0fe931808a5c8cfecc9a9e4fa..e1ea503d33536d55156d5ee883787be67db42b55 100644
--- a/tests/init_dirichletfields.hpp
+++ b/tests/init_dirichletfields.hpp
@@ -24,7 +24,7 @@ public:
         // Calculate expected result from Nektar++
         fixt_explist->ImposeDirichletConditions(outcoeffs);
 
-        // Copy expected result from Array to fixt_expected
+        // Copy expected result from Array to pointer
         double *coeffptr = outcoeffs.get();
         for (auto const &block : blocks)
         {
diff --git a/tests/init_fwdtransfields.hpp b/tests/init_fwdtransfields.hpp
index 7145a952c1fd2ddc71b33add50dc544f9618d7e7..ba151f2bf1c15556489b3338b0c5012a31f61cf8 100644
--- a/tests/init_fwdtransfields.hpp
+++ b/tests/init_fwdtransfields.hpp
@@ -34,7 +34,7 @@ public:
 
         double *xptr = x.get(), *yptr = y.get(), *zptr = z.get(),
                *fceptr = fce.get();
-        for (auto const &block : fixt_in->GetBlocks())
+        for (auto const &block : blocks)
         {
             for (size_t el = 0; el < block.num_elements; ++el)
             {
@@ -103,12 +103,12 @@ public:
         Array<OneD, NekDouble> outcoeffs(fixt_explist->GetNcoeffs(), 0.0);
 
         // Set test case
-        SetTestCase(blocks, inphys.get(), false);
+        SetTestCase(fixt_in->GetBlocks(), inphys.get(), false);
 
         // Calculate expected result from Nektar++
         fixt_explist->FwdTrans(inphys, outcoeffs);
 
-        // Copy expected result from Array to fixt_expected
+        // Copy expected result from Array to pointer
         double *coeffptr = outcoeffs.get();
         for (auto const &block : blocks)
         {
diff --git a/tests/init_helmholtzfields.hpp b/tests/init_helmholtzfields.hpp
index 8a5644c384ddda6be42247f9d06143ecd82c4480..ba78e21c30cfad3fff5e132366ce78e7c2bbe723 100644
--- a/tests/init_helmholtzfields.hpp
+++ b/tests/init_helmholtzfields.hpp
@@ -17,7 +17,7 @@ public:
     void SetTestCase(const std::vector<BlockAttributes> &blocks, double *inptr,
                      bool padding = true)
     {
-        for (auto const &block : fixt_in->GetBlocks())
+        for (auto const &block : blocks)
         {
             for (size_t el = 0; el < block.num_elements; ++el)
             {
@@ -47,7 +47,7 @@ public:
         Array<OneD, NekDouble> tmp;
 
         // Set test case
-        SetTestCase(blocks, incoeffs.get(), false);
+        SetTestCase(fixt_in->GetBlocks(), incoeffs.get(), false);
 
         // Calculate expected result from Nektar++
         auto e = 0, offset = 0;
@@ -69,7 +69,7 @@ public:
             }
         }
 
-        // Copy expected result from Array to fixt_expected
+        // Copy expected result from Array to pointer
         double *coeffptr = outcoeffs.get();
         for (auto const &block : blocks)
         {
diff --git a/tests/init_helmsolvefields.hpp b/tests/init_helmsolvefields.hpp
index ec37467344e835d00e1569aa4f8058ccadc625ee..33034087eddee5f0f7e3bf1d618fd213ec72d74f 100644
--- a/tests/init_helmsolvefields.hpp
+++ b/tests/init_helmsolvefields.hpp
@@ -33,7 +33,7 @@ public:
 
         double *xptr = x.get(), *yptr = y.get(), *zptr = z.get(),
                *fceptr = fce.get();
-        for (auto const &block : fixt_in->GetBlocks())
+        for (auto const &block : blocks)
         {
             for (size_t el = 0; el < block.num_elements; ++el)
             {
@@ -102,7 +102,7 @@ public:
         Array<OneD, NekDouble> outcoeffs(fixt_explist->GetNcoeffs(), 0.0);
 
         // Set test case
-        SetTestCase(blocks, inphys.get(), false);
+        SetTestCase(fixt_in->GetBlocks(), inphys.get(), false);
 
         // Calculate expected result from Nektar++
         StdRegions::ConstFactorMap factors;
@@ -112,7 +112,7 @@ public:
                 : 1.0;
         fixt_explist->HelmSolve(inphys, outcoeffs, factors);
 
-        // Copy expected result from Array to fixt_expected
+        // Copy expected result from Array to pointer
         double *coeffptr = outcoeffs.get();
         for (auto const &block : blocks)
         {
diff --git a/tests/init_identitycoefffields.hpp b/tests/init_identitycoefffields.hpp
index d61981f099bc6482dfac240ee79bf5983f17a6f5..704ffe7866a8d98be62b358a3a75bcba199e00c6 100644
--- a/tests/init_identitycoefffields.hpp
+++ b/tests/init_identitycoefffields.hpp
@@ -16,7 +16,7 @@ public:
     void SetTestCase(const std::vector<BlockAttributes> &blocks, double *inptr,
                      bool padding = true)
     {
-        for (auto const &block : fixt_in->GetBlocks())
+        for (auto const &block : blocks)
         {
             for (size_t el = 0; el < block.num_elements; ++el)
             {
@@ -41,8 +41,8 @@ public:
     void ExpectedSolution(const std::vector<BlockAttributes> &blocks,
                           double *inptr)
     {
-        // Copy expected result to fixt_expected
-        SetTestCase(blocks, fixt_expected->GetStorage().GetCPUPtr());
+        // Copy expected result to pointer
+        SetTestCase(blocks, inptr);
     }
 };
 
diff --git a/tests/init_identityphysfields.hpp b/tests/init_identityphysfields.hpp
index 59262fdc47b49770569178fba0edffea7fad5b70..03d9d786f34f4050f87764e46c1cc96c14e6e23a 100644
--- a/tests/init_identityphysfields.hpp
+++ b/tests/init_identityphysfields.hpp
@@ -16,7 +16,7 @@ public:
     void SetTestCase(const std::vector<BlockAttributes> &blocks, double *inptr,
                      bool padding = true)
     {
-        for (auto const &block : fixt_in->GetBlocks())
+        for (auto const &block : blocks)
         {
             for (size_t el = 0; el < block.num_elements; ++el)
             {
@@ -41,8 +41,8 @@ public:
     void ExpectedSolution(const std::vector<BlockAttributes> &blocks,
                           double *inptr)
     {
-        // Copy expected result to fixt_expected
-        SetTestCase(blocks, fixt_expected->GetStorage().GetCPUPtr());
+        // Copy expected result to pojnter
+        SetTestCase(blocks, inptr);
     }
 };
 
diff --git a/tests/init_ipwrtbasefields.hpp b/tests/init_ipwrtbasefields.hpp
index 71de7bc889a0bbc95305d7b814551e86cf5b151e..1524174eff2784cb96d3dbac1834298bc3e10649 100644
--- a/tests/init_ipwrtbasefields.hpp
+++ b/tests/init_ipwrtbasefields.hpp
@@ -17,7 +17,7 @@ public:
     void SetTestCase(const std::vector<BlockAttributes> &blocks, double *inptr,
                      bool padding = true)
     {
-        for (auto const &block : fixt_in->GetBlocks())
+        for (auto const &block : blocks)
         {
             for (size_t el = 0; el < block.num_elements; ++el)
             {
@@ -46,12 +46,12 @@ public:
         Array<OneD, NekDouble> outcoeffs(fixt_explist->GetNcoeffs());
 
         // Set test case
-        SetTestCase(blocks, inphys.get(), false);
+        SetTestCase(fixt_in->GetBlocks(), inphys.get(), false);
 
         // Calculate expected result from Nektar++
         fixt_explist->IProductWRTBase(inphys, outcoeffs);
 
-        // Copy expected result from Array to fixt_expected
+        // Copy expected result from Array to pointer
         double *ptr = outcoeffs.get();
         for (auto const &block : blocks)
         {
diff --git a/tests/init_ipwrtderivbasefields.hpp b/tests/init_ipwrtderivbasefields.hpp
index 4c2971766b134385943118c4d213445f4b265c67..3ee889ceda348e183daa6f3724037f13be82a42c 100644
--- a/tests/init_ipwrtderivbasefields.hpp
+++ b/tests/init_ipwrtderivbasefields.hpp
@@ -19,7 +19,7 @@ public:
     {
         for (size_t k = 0; k < fixt_explist->GetCoordim(0); k++)
         {
-            for (auto const &block : fixt_in->GetBlocks())
+            for (auto const &block : blocks)
             {
                 for (size_t el = 0; el < block.num_elements; ++el)
                 {
@@ -50,7 +50,7 @@ public:
         Array<OneD, NekDouble> outcoeffs(fixt_explist->GetNcoeffs(), 0.0);
 
         // Set test case
-        SetTestCase(blocks, inphys.get(), false);
+        SetTestCase(fixt_in->GetBlocks(), inphys.get(), false);
 
         // Calculate expected result from Nektar++
         Array<OneD, Array<OneD, NekDouble>> inphysarray(
@@ -69,7 +69,7 @@ public:
         }
         fixt_explist->IProductWRTDerivBase(inphysarray, outcoeffs);
 
-        // Copy expected result from Array to fixt_expected
+        // Copy expected result from Array to pointer
         double *ptr = outcoeffs.get();
         for (auto const &block : blocks)
         {
diff --git a/tests/init_massfields.hpp b/tests/init_massfields.hpp
index 0d8bb8e5a26e7d545cef5bbeb6a25c37d5058008..ae6cb7219b8f1abba7d3aeb3f414fc65d0e2861b 100644
--- a/tests/init_massfields.hpp
+++ b/tests/init_massfields.hpp
@@ -16,7 +16,7 @@ public:
     void SetTestCase(const std::vector<BlockAttributes> &blocks, double *inptr,
                      bool padding = true)
     {
-        for (auto const &block : fixt_in->GetBlocks())
+        for (auto const &block : blocks)
         {
             for (size_t el = 0; el < block.num_elements; ++el)
             {
@@ -46,13 +46,13 @@ public:
         Array<OneD, NekDouble> outcoeffs(fixt_explist->GetNcoeffs());
 
         // Set test case
-        SetTestCase(blocks, incoeffs.get(), false);
+        SetTestCase(fixt_in->GetBlocks(), incoeffs.get(), false);
 
         // Calculate expected result from Nektar++
         fixt_explist->BwdTrans(incoeffs, bwdtrans);
         fixt_explist->IProductWRTBase(bwdtrans, outcoeffs);
 
-        // Copy expected result from Array to fixt_expected
+        // Copy expected result from Array to pointer
         double *coeffptr = outcoeffs.get();
         for (auto const &block : blocks)
         {
diff --git a/tests/init_nullpreconfields.hpp b/tests/init_nullpreconfields.hpp
index 20ff840c06e82dca757a14ceb842dd082970aded..80a6633daf9e7d957b1252e9c83081cbdfb0dac8 100644
--- a/tests/init_nullpreconfields.hpp
+++ b/tests/init_nullpreconfields.hpp
@@ -22,7 +22,7 @@ public:
     void SetTestCase(const std::vector<BlockAttributes> &blocks, double *inptr,
                      bool padding = true)
     {
-        for (auto const &block : fixt_in->GetBlocks())
+        for (auto const &block : blocks)
         {
             for (size_t el = 0; el < block.num_elements; ++el)
             {
@@ -51,7 +51,7 @@ public:
         Array<OneD, NekDouble> outcoeffs(fixt_explist->GetNcoeffs(), 0.0);
 
         // Set test case
-        SetTestCase(blocks, incoeffs.get(), false);
+        SetTestCase(fixt_in->GetBlocks(), incoeffs.get(), false);
 
         // Calculate expected result from Nektar++
         auto map = fixt_explist->GetLocalToGlobalMap();
@@ -62,7 +62,7 @@ public:
             GetPreconFactory().CreateInstance("Null", globalSys, map);
         precond->DoPreconditioner(incoeffs, outcoeffs, true);
 
-        // Copy expected result from Array to fixt_expected
+        // Copy expected result from Array to pointer
         double *coeffptr = outcoeffs.get();
         for (auto const &block : blocks)
         {
diff --git a/tests/init_physderivfields.hpp b/tests/init_physderivfields.hpp
index f59c2e6dc717ea19f60d91cdfb4c4e0a2cd1d758..5adeb60a0135d1d0fa33c88f8372c5a49aaaf5a9 100644
--- a/tests/init_physderivfields.hpp
+++ b/tests/init_physderivfields.hpp
@@ -66,12 +66,12 @@ public:
         Array<OneD, NekDouble> outphys2 = NullNekDouble1DArray;
 
         // Set test case
-        SetTestCase(blocks, inphys.get(), false);
+        SetTestCase(fixt_in->GetBlocks(), inphys.get(), false);
 
         // Calculate expected result from Nektar++
         fixt_explist->PhysDeriv(inphys, outphys0, outphys1, outphys2);
 
-        // Copy expected result from Array to fixt_expected
+        // Copy expected result from Array to pointer
         double *ptr = outphys.get();
         for (auto const &block : blocks)
         {
@@ -182,7 +182,7 @@ public:
         Array<OneD, NekDouble> outphys2 = NullNekDouble1DArray;
 
         // Set test case
-        SetTestCase(blocks, inphys.get(), false);
+        SetTestCase(fixt_in->GetBlocks(), inphys.get(), false);
 
         // Calculate expected result from Nektar++
         fixt_explist->PhysDeriv(inphys, outphys0, outphys1, outphys2);
@@ -326,7 +326,7 @@ public:
             outphys1 + fixt_explist->GetTotPoints();
 
         // Set test case
-        SetTestCase(blocks, inphys.get(), false);
+        SetTestCase(fixt_in->GetBlocks(), inphys.get(), false);
 
         // Calculate expected result from Nektar++
         fixt_explist->PhysDeriv(inphys, outphys0, outphys1, outphys2);
diff --git a/tests/test_assmbscatrcuda.cpp b/tests/test_assmbscatrcuda.cpp
index 13768c41bda37318762d59bafd7093bb031c75d5..42f5f6c5000a8b2ac104c2977255d7b94963e8e6 100644
--- a/tests/test_assmbscatrcuda.cpp
+++ b/tests/test_assmbscatrcuda.cpp
@@ -24,8 +24,9 @@ BOOST_FIXTURE_TEST_CASE(assmbscatrcuda_seg, Seg)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -42,8 +43,9 @@ BOOST_FIXTURE_TEST_CASE(assmbscatrcuda_quad, Quad)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -60,8 +62,9 @@ BOOST_FIXTURE_TEST_CASE(assmbscatrcuda_tri, Tri)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -78,8 +81,9 @@ BOOST_FIXTURE_TEST_CASE(assmbscatrcuda_square_all_elements, SquareAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -96,8 +100,9 @@ BOOST_FIXTURE_TEST_CASE(assmbscatrcuda_hex, Hex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -114,8 +119,9 @@ BOOST_FIXTURE_TEST_CASE(assmbscatrcuda_prism, Prism)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -132,8 +138,9 @@ BOOST_FIXTURE_TEST_CASE(assmbscatrcuda_pyr, Pyr)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -150,8 +157,9 @@ BOOST_FIXTURE_TEST_CASE(assmbscatrcuda_tet, Tet)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -168,8 +176,9 @@ BOOST_FIXTURE_TEST_CASE(assmbscatrcuda_cube_prism_hex, CubePrismHex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -186,8 +195,9 @@ BOOST_FIXTURE_TEST_CASE(assmbscatrcuda_cube_all_elements, CubeAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
diff --git a/tests/test_bwdtranscuda.cpp b/tests/test_bwdtranscuda.cpp
index 60a10b91d7fa1b42e9cd9a83e64512b719aa1eb4..da83b92532e07b4223c9bca78de404a2e995d508 100644
--- a/tests/test_bwdtranscuda.cpp
+++ b/tests/test_bwdtranscuda.cpp
@@ -24,8 +24,9 @@ BOOST_FIXTURE_TEST_CASE(bwdtranscuda_seg, Seg)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -42,8 +43,9 @@ BOOST_FIXTURE_TEST_CASE(bwdtranscuda_quad, Quad)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -60,8 +62,9 @@ BOOST_FIXTURE_TEST_CASE(bwdtranscuda_tri, Tri)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -78,8 +81,9 @@ BOOST_FIXTURE_TEST_CASE(bwdtranscuda_square_all_elements, SquareAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -96,8 +100,9 @@ BOOST_FIXTURE_TEST_CASE(bwdtranscuda_hex, Hex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -114,8 +119,9 @@ BOOST_FIXTURE_TEST_CASE(bwdtranscuda_prism, Prism)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -132,8 +138,9 @@ BOOST_FIXTURE_TEST_CASE(bwdtranscuda_pyr, Pyr)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -150,8 +157,9 @@ BOOST_FIXTURE_TEST_CASE(bwdtranscuda_tet, Tet)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -168,8 +176,9 @@ BOOST_FIXTURE_TEST_CASE(bwdtranscuda_cube_prism_hex, CubePrismHex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -186,8 +195,9 @@ BOOST_FIXTURE_TEST_CASE(bwdtranscuda_cube_all_elements, CubeAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
diff --git a/tests/test_diagpreconcuda.cpp b/tests/test_diagpreconcuda.cpp
index df15cfa8dd7baa769854a6d23645411475ad3401..f48d74876a4ac5f623d505f7eaa6039cd2ff934a 100644
--- a/tests/test_diagpreconcuda.cpp
+++ b/tests/test_diagpreconcuda.cpp
@@ -15,8 +15,10 @@ BOOST_AUTO_TEST_SUITE(TestDiagPreconCUDA)
 BOOST_FIXTURE_TEST_CASE(diagpreconcuda_seg, Helmholtz1D_Seg)
 {
     Configure();
-    SetTestCase(fixtcuda_in->GetBlocks(), fixtcuda_in->GetStorage().GetCPUPtr());
-    auto HelmholtzOp  = Helmholtz<>::create(fixt_explist, "StdMat");
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    auto HelmholtzOp  = Helmholtz<>::create(fixt_explist, "CUDA");
     auto DiagPreconOp = DiagPrecon<>::create(fixt_explist, "CUDA");
     HelmholtzOp->setLambda(1.0);
     DiagPreconOp->configure(HelmholtzOp);
@@ -26,16 +28,19 @@ BOOST_FIXTURE_TEST_CASE(diagpreconcuda_seg, Helmholtz1D_Seg)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
 BOOST_FIXTURE_TEST_CASE(diagpreconcuda_tri_quad, Helmholtz2D_Tri_Quad)
 {
     Configure();
-    SetTestCase(fixtcuda_in->GetBlocks(), fixtcuda_in->GetStorage().GetCPUPtr());
-    auto HelmholtzOp  = Helmholtz<>::create(fixt_explist, "StdMat");
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    auto HelmholtzOp  = Helmholtz<>::create(fixt_explist, "CUDA");
     auto DiagPreconOp = DiagPrecon<>::create(fixt_explist, "CUDA");
     HelmholtzOp->setLambda(1.0);
     DiagPreconOp->configure(HelmholtzOp);
@@ -45,16 +50,19 @@ BOOST_FIXTURE_TEST_CASE(diagpreconcuda_tri_quad, Helmholtz2D_Tri_Quad)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
 BOOST_FIXTURE_TEST_CASE(diagpreconcuda_hex, Helmholtz3D_Hex)
 {
     Configure();
-    SetTestCase(fixtcuda_in->GetBlocks(), fixtcuda_in->GetStorage().GetCPUPtr());
-    auto HelmholtzOp  = Helmholtz<>::create(fixt_explist, "StdMat");
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    auto HelmholtzOp  = Helmholtz<>::create(fixt_explist, "CUDA");
     auto DiagPreconOp = DiagPrecon<>::create(fixt_explist, "CUDA");
     HelmholtzOp->setLambda(1.0);
     DiagPreconOp->configure(HelmholtzOp);
@@ -64,16 +72,19 @@ BOOST_FIXTURE_TEST_CASE(diagpreconcuda_hex, Helmholtz3D_Hex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-10));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
     }
 }
 
 BOOST_FIXTURE_TEST_CASE(diagpreconcuda_prism, Helmholtz3D_Prism)
 {
     Configure();
-    SetTestCase(fixtcuda_in->GetBlocks(), fixtcuda_in->GetStorage().GetCPUPtr());
-    auto HelmholtzOp  = Helmholtz<>::create(fixt_explist, "StdMat");
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    auto HelmholtzOp  = Helmholtz<>::create(fixt_explist, "CUDA");
     auto DiagPreconOp = DiagPrecon<>::create(fixt_explist, "CUDA");
     HelmholtzOp->setLambda(1.0);
     DiagPreconOp->configure(HelmholtzOp);
@@ -83,16 +94,19 @@ BOOST_FIXTURE_TEST_CASE(diagpreconcuda_prism, Helmholtz3D_Prism)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-10));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
     }
 }
 
 BOOST_FIXTURE_TEST_CASE(diagpreconcuda_pyr, Helmholtz3D_Pyr)
 {
     Configure();
-    SetTestCase(fixtcuda_in->GetBlocks(), fixtcuda_in->GetStorage().GetCPUPtr());
-    auto HelmholtzOp  = Helmholtz<>::create(fixt_explist, "StdMat");
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    auto HelmholtzOp  = Helmholtz<>::create(fixt_explist, "CUDA");
     auto DiagPreconOp = DiagPrecon<>::create(fixt_explist, "CUDA");
     HelmholtzOp->setLambda(1.0);
     DiagPreconOp->configure(HelmholtzOp);
@@ -102,16 +116,19 @@ BOOST_FIXTURE_TEST_CASE(diagpreconcuda_pyr, Helmholtz3D_Pyr)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-10));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
     }
 }
 
-BOOST_FIXTURE_TEST_CASE(diagpreconcuda_tet, Helmholtz3D_Tet)
+/*BOOST_FIXTURE_TEST_CASE(diagpreconcuda_tet, Helmholtz3D_Tet)
 {
     Configure();
-    SetTestCase(fixtcuda_in->GetBlocks(), fixtcuda_in->GetStorage().GetCPUPtr());
-    auto HelmholtzOp  = Helmholtz<>::create(fixt_explist, "StdMat");
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    auto HelmholtzOp  = Helmholtz<>::create(fixt_explist, "CUDA");
     auto DiagPreconOp = DiagPrecon<>::create(fixt_explist, "CUDA");
     HelmholtzOp->setLambda(1.0);
     DiagPreconOp->configure(HelmholtzOp);
@@ -121,9 +138,10 @@ BOOST_FIXTURE_TEST_CASE(diagpreconcuda_tet, Helmholtz3D_Tet)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-10));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
     }
-}
+}*/
 
 BOOST_AUTO_TEST_SUITE_END()
diff --git a/tests/test_dirichletcuda.cpp b/tests/test_dirichletcuda.cpp
index 9778d3ac4b13373b1e4907a5b726abaade3d293d..cc036970b42954f20741b03c815eec5743601246 100644
--- a/tests/test_dirichletcuda.cpp
+++ b/tests/test_dirichletcuda.cpp
@@ -9,7 +9,7 @@
 #include "Operators/OperatorDirBndCond.hpp"
 #include "init_dirichletfields.hpp"
 
-BOOST_AUTO_TEST_SUITE(TestDirichlet)
+BOOST_AUTO_TEST_SUITE(TestDirichletCUDA)
 
 BOOST_FIXTURE_TEST_CASE(dirichlet1dcuda_seg, Helmholtz1D_Seg)
 {
@@ -21,8 +21,9 @@ BOOST_FIXTURE_TEST_CASE(dirichlet1dcuda_seg, Helmholtz1D_Seg)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -36,8 +37,9 @@ BOOST_FIXTURE_TEST_CASE(dirichlet2dcuda_tri_quad, Helmholtz2D_Tri_Quad)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -51,8 +53,9 @@ BOOST_FIXTURE_TEST_CASE(dirichlet3dcuda_hex, Helmholtz3D_Hex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -66,8 +69,9 @@ BOOST_FIXTURE_TEST_CASE(dirichlet3dcuda_prism, Helmholtz3D_Prism)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -81,8 +85,9 @@ BOOST_FIXTURE_TEST_CASE(dirichlet3dcuda_pyr, Helmholtz3D_Pyr)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -96,8 +101,9 @@ BOOST_FIXTURE_TEST_CASE(dirichlet3dcuda_tet, Helmholtz3D_Tet)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
diff --git a/tests/test_fwdtrans.cpp b/tests/test_fwdtrans.cpp
index 2d2af43d221ce99e2f2c783e211d8cae95882867..d0f7b4237212e2164063b39302a6413defc88c33 100644
--- a/tests/test_fwdtrans.cpp
+++ b/tests/test_fwdtrans.cpp
@@ -12,11 +12,6 @@
 
 BOOST_AUTO_TEST_SUITE(TestFwdTrans)
 
-using namespace std;
-using namespace Nektar::Operators;
-using namespace Nektar::LibUtilities;
-using namespace Nektar;
-
 BOOST_FIXTURE_TEST_CASE(fwdtrans_seg, Helmholtz1D_Seg)
 {
     Configure();
diff --git a/tests/test_helmholtzcuda.cpp b/tests/test_helmholtzcuda.cpp
index 197b8b7e88ce1a0e3ef33d06d84c062c7afea8b1..afa0d1253e2071bcf16a06fdda9997c0f1cac53d 100644
--- a/tests/test_helmholtzcuda.cpp
+++ b/tests/test_helmholtzcuda.cpp
@@ -25,8 +25,9 @@ BOOST_FIXTURE_TEST_CASE(helmholtzcuda_seg, Seg)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -44,8 +45,9 @@ BOOST_FIXTURE_TEST_CASE(helmholtzcuda_quad, Quad)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -63,8 +65,9 @@ BOOST_FIXTURE_TEST_CASE(helmholtzcuda_tri, Tri)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -82,8 +85,9 @@ BOOST_FIXTURE_TEST_CASE(helmholtzcuda_square_all_elements, SquareAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -101,8 +105,9 @@ BOOST_FIXTURE_TEST_CASE(helmholtzcuda_hex, Hex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -120,8 +125,9 @@ BOOST_FIXTURE_TEST_CASE(helmholtzcuda_prism, Prism)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -139,8 +145,9 @@ BOOST_FIXTURE_TEST_CASE(helmholtzcuda_pyr, Pyr)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -158,8 +165,9 @@ BOOST_FIXTURE_TEST_CASE(helmholtzcuda_tet, Tet)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -177,8 +185,9 @@ BOOST_FIXTURE_TEST_CASE(helmholtzcuda_cube_prism_hex, CubePrismHex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -196,8 +205,9 @@ BOOST_FIXTURE_TEST_CASE(helmholtzcuda_cube_all_elements, CubeAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
diff --git a/tests/test_helmsolve.cpp b/tests/test_helmsolve.cpp
index ea4880e838c0b40f41aa3ed16722f915bb517592..d5f07ef5a137892e2e946daad4728ef2b5e91caa 100644
--- a/tests/test_helmsolve.cpp
+++ b/tests/test_helmsolve.cpp
@@ -12,11 +12,6 @@
 
 BOOST_AUTO_TEST_SUITE(TestHelmSolve)
 
-using namespace std;
-using namespace Nektar::Operators;
-using namespace Nektar::LibUtilities;
-using namespace Nektar;
-
 BOOST_FIXTURE_TEST_CASE(helmsolve_seg, Helmholtz1D_Seg)
 {
     Configure();
diff --git a/tests/test_identitycoeffcuda.cpp b/tests/test_identitycoeffcuda.cpp
index 8a0d31184dd7bb9b50a9fea02cf0219bc89bf15d..d8bd9630200ec874a7865e8423d42e8450b42d6f 100644
--- a/tests/test_identitycoeffcuda.cpp
+++ b/tests/test_identitycoeffcuda.cpp
@@ -24,8 +24,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_seg, Seg)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -42,8 +43,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_quad, Quad)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -60,8 +62,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_tri, Tri)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -78,8 +81,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_square_all_elements, SquareAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -96,8 +100,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_hex, Hex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -114,8 +119,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_prism, Prism)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -132,8 +138,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_pyr, Pyr)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -150,8 +157,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_tet, Tet)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -168,8 +176,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_cube_prism_hex, CubePrismHex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -186,8 +195,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_cube_all_elements, CubeAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
diff --git a/tests/test_identityphyscuda.cpp b/tests/test_identityphyscuda.cpp
index 1780fcbe1e98994f6541de04a5729e486c8166d7..ed6823c46b42e4564cd9e512d189ddc843fb4e7b 100644
--- a/tests/test_identityphyscuda.cpp
+++ b/tests/test_identityphyscuda.cpp
@@ -24,8 +24,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_seg, Seg)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -42,8 +43,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_quad, Quad)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -60,8 +62,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_tri, Tri)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -78,8 +81,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_square_all_elements, SquareAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -96,8 +100,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_hex, Hex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -114,8 +119,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_prism, Prism)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -132,8 +138,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_pyr, Pyr)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -150,8 +157,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_tet, Tet)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -168,8 +176,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_cube_prism_hex, CubePrismHex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -186,8 +195,9 @@ BOOST_FIXTURE_TEST_CASE(identitycuda_cube_all_elements, CubeAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
diff --git a/tests/test_ipwrtbasecuda.cpp b/tests/test_ipwrtbasecuda.cpp
index e315b3a3cd831a10c907573d5bd5eb1b6338bf22..2d1858677cdd49538cef3ec2012ca9641323102f 100644
--- a/tests/test_ipwrtbasecuda.cpp
+++ b/tests/test_ipwrtbasecuda.cpp
@@ -25,8 +25,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtbasecuda_seg, Seg)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -44,8 +45,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtbasecuda_quad, Quad)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -63,8 +65,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtbasecuda_tri, Tri)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -82,8 +85,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtbasecuda_square_all_elements, SquareAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -101,8 +105,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtbasecuda_hex, Hex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -120,8 +125,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtbasecuda_prism, Prism)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -139,8 +145,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtbasecuda_pyr, Pyr)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -158,8 +165,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtbasecuda_tet, Tet)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -177,8 +185,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtbasecuda_cube_prism_hex, CubePrismHex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -196,8 +205,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtbasecuda_cube_all_elements, CubeAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
diff --git a/tests/test_ipwrtderivbasecuda.cpp b/tests/test_ipwrtderivbasecuda.cpp
index 55261b9328145150fdf8a9ba2ef7560afaa9f577..efe3fa636dec457b56ba8b06ce57660672fdfa86 100644
--- a/tests/test_ipwrtderivbasecuda.cpp
+++ b/tests/test_ipwrtderivbasecuda.cpp
@@ -25,8 +25,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtderivbasecuda_seg, Seg)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -44,8 +45,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtderivbasecuda_quad, Quad)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -63,8 +65,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtderivbasecuda_tri, Tri)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -83,8 +86,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtderivbasecuda_square_all_elements,
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -102,8 +106,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtderivbasecuda_hex, Hex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -121,8 +126,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtderivbasecuda_prism, Prism)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -140,8 +146,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtderivbasecuda_pyr, Pyr)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -159,8 +166,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtderivbasecuda_tet, Tet)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -178,8 +186,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtderivbasecuda_cube_prism_hex, CubePrismHex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -197,8 +206,9 @@ BOOST_FIXTURE_TEST_CASE(ipwrtderivbasecuda_cube_all_elements, CubeAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
diff --git a/tests/test_masscuda.cpp b/tests/test_masscuda.cpp
index 4182592152a3155e0d2c29e8d8073f00b9140a5c..a48828c240667d1d4af8415a7c27ddcc53788fd8 100644
--- a/tests/test_masscuda.cpp
+++ b/tests/test_masscuda.cpp
@@ -14,150 +14,180 @@ BOOST_AUTO_TEST_SUITE(TestMassCUDA)
 BOOST_FIXTURE_TEST_CASE(masscuda_seg, Seg)
 {
     Configure();
-    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
-    Mass<>::create(fixt_explist, "CUDA")->apply(*fixt_in, *fixt_out);
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Mass<>::create(fixt_explist, "CUDA")->apply(*fixtcuda_in, *fixtcuda_out);
     ExpectedSolution(fixt_expected->GetBlocks(),
                      fixt_expected->GetStorage().GetCPUPtr());
-    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
 BOOST_FIXTURE_TEST_CASE(masscuda_quad, Quad)
 {
     Configure();
-    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
-    Mass<>::create(fixt_explist, "CUDA")->apply(*fixt_in, *fixt_out);
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Mass<>::create(fixt_explist, "CUDA")->apply(*fixtcuda_in, *fixtcuda_out);
     ExpectedSolution(fixt_expected->GetBlocks(),
                      fixt_expected->GetStorage().GetCPUPtr());
-    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
 BOOST_FIXTURE_TEST_CASE(masscuda_tri, Tri)
 {
     Configure();
-    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
-    Mass<>::create(fixt_explist, "CUDA")->apply(*fixt_in, *fixt_out);
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Mass<>::create(fixt_explist, "CUDA")->apply(*fixtcuda_in, *fixtcuda_out);
     ExpectedSolution(fixt_expected->GetBlocks(),
                      fixt_expected->GetStorage().GetCPUPtr());
-    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
 BOOST_FIXTURE_TEST_CASE(masscuda_square_all_elements, SquareAllElements)
 {
     Configure();
-    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
-    Mass<>::create(fixt_explist, "CUDA")->apply(*fixt_in, *fixt_out);
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Mass<>::create(fixt_explist, "CUDA")->apply(*fixtcuda_in, *fixtcuda_out);
     ExpectedSolution(fixt_expected->GetBlocks(),
                      fixt_expected->GetStorage().GetCPUPtr());
-    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
 BOOST_FIXTURE_TEST_CASE(masscuda_hex, Hex)
 {
     Configure();
-    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
-    Mass<>::create(fixt_explist, "CUDA")->apply(*fixt_in, *fixt_out);
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Mass<>::create(fixt_explist, "CUDA")->apply(*fixtcuda_in, *fixtcuda_out);
     ExpectedSolution(fixt_expected->GetBlocks(),
                      fixt_expected->GetStorage().GetCPUPtr());
-    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
 BOOST_FIXTURE_TEST_CASE(masscuda_prism, Prism)
 {
     Configure();
-    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
-    Mass<>::create(fixt_explist, "CUDA")->apply(*fixt_in, *fixt_out);
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Mass<>::create(fixt_explist, "CUDA")->apply(*fixtcuda_in, *fixtcuda_out);
     ExpectedSolution(fixt_expected->GetBlocks(),
                      fixt_expected->GetStorage().GetCPUPtr());
-    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
 BOOST_FIXTURE_TEST_CASE(masscuda_pyr, Pyr)
 {
     Configure();
-    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
-    Mass<>::create(fixt_explist, "CUDA")->apply(*fixt_in, *fixt_out);
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Mass<>::create(fixt_explist, "CUDA")->apply(*fixtcuda_in, *fixtcuda_out);
     ExpectedSolution(fixt_expected->GetBlocks(),
                      fixt_expected->GetStorage().GetCPUPtr());
-    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
 BOOST_FIXTURE_TEST_CASE(masscuda_tet, Tet)
 {
     Configure();
-    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
-    Mass<>::create(fixt_explist, "CUDA")->apply(*fixt_in, *fixt_out);
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Mass<>::create(fixt_explist, "CUDA")->apply(*fixtcuda_in, *fixtcuda_out);
     ExpectedSolution(fixt_expected->GetBlocks(),
                      fixt_expected->GetStorage().GetCPUPtr());
-    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
 BOOST_FIXTURE_TEST_CASE(masscuda_cube_prism_hex, CubePrismHex)
 {
     Configure();
-    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
-    Mass<>::create(fixt_explist, "CUDA")->apply(*fixt_in, *fixt_out);
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Mass<>::create(fixt_explist, "CUDA")->apply(*fixtcuda_in, *fixtcuda_out);
     ExpectedSolution(fixt_expected->GetBlocks(),
                      fixt_expected->GetStorage().GetCPUPtr());
-    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
 BOOST_FIXTURE_TEST_CASE(masscuda_cube_all_elements, CubeAllElements)
 {
     Configure();
-    SetTestCase(fixt_in->GetBlocks(), fixt_in->GetStorage().GetCPUPtr());
-    Mass<>::create(fixt_explist, "CUDA")->apply(*fixt_in, *fixt_out);
+    SetTestCase(
+        fixtcuda_in->GetBlocks(),
+        fixtcuda_in->template GetStorage<MemoryRegionCUDA>().GetCPUPtr());
+    Mass<>::create(fixt_explist, "CUDA")->apply(*fixtcuda_in, *fixtcuda_out);
     ExpectedSolution(fixt_expected->GetBlocks(),
                      fixt_expected->GetStorage().GetCPUPtr());
-    BOOST_TEST(fixt_out->compare(*fixt_expected, 1.0E-12));
+    BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixt_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
diff --git a/tests/test_neumanncuda.cpp b/tests/test_neumanncuda.cpp
index 33051462aa0b38d574869e347293031789b38671..1bde94249f718475101729530e0c81c56336b407 100644
--- a/tests/test_neumanncuda.cpp
+++ b/tests/test_neumanncuda.cpp
@@ -19,8 +19,9 @@ BOOST_FIXTURE_TEST_CASE(dirichlet1dcuda_seg, Helmholtz1D_Seg)
     BOOST_TEST(fixtcuda_out->compare(*fixt_out, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_out->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_out->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -32,8 +33,9 @@ BOOST_FIXTURE_TEST_CASE(dirichlet2dcuda_tri_quad, Helmholtz2D_Tri_Quad)
     BOOST_TEST(fixtcuda_out->compare(*fixt_out, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_out->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_out->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -45,8 +47,9 @@ BOOST_FIXTURE_TEST_CASE(dirichlet3dcuda_hex, Helmholtz3D_Hex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_out, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_out->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_out->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -58,8 +61,9 @@ BOOST_FIXTURE_TEST_CASE(dirichlet3dcuda_prism, Helmholtz3D_Prism)
     BOOST_TEST(fixtcuda_out->compare(*fixt_out, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_out->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_out->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -71,8 +75,9 @@ BOOST_FIXTURE_TEST_CASE(dirichlet3dcuda_pyr, Helmholtz3D_Pyr)
     BOOST_TEST(fixtcuda_out->compare(*fixt_out, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_out->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_out->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -84,8 +89,9 @@ BOOST_FIXTURE_TEST_CASE(dirichlet3dcuda_tet, Helmholtz3D_Tet)
     BOOST_TEST(fixtcuda_out->compare(*fixt_out, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_out->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_out->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
diff --git a/tests/test_nullpreconcuda.cpp b/tests/test_nullpreconcuda.cpp
index 945ad6b6d568ff5c4f51dbde17838473b9652fca..8af17c41e18177199d1de1e2e26c4fac60caf86b 100644
--- a/tests/test_nullpreconcuda.cpp
+++ b/tests/test_nullpreconcuda.cpp
@@ -11,7 +11,6 @@
 
 BOOST_AUTO_TEST_SUITE(TestNullPreconCUDA)
 
-
 BOOST_FIXTURE_TEST_CASE(nullpreconcuda_seg, Helmholtz1D_Seg)
 {
     Configure();
@@ -25,8 +24,9 @@ BOOST_FIXTURE_TEST_CASE(nullpreconcuda_seg, Helmholtz1D_Seg)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-15));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-15);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-15);
     }
 }
 
@@ -43,8 +43,9 @@ BOOST_FIXTURE_TEST_CASE(nullpreconcuda_tri_quad, Helmholtz2D_Tri_Quad)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-15));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-15);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-15);
     }
 }
 
@@ -61,8 +62,9 @@ BOOST_FIXTURE_TEST_CASE(nullpreconcuda_hex, Helmholtz3D_Hex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-15));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-15);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-15);
     }
 }
 
@@ -79,8 +81,9 @@ BOOST_FIXTURE_TEST_CASE(nullpreconcuda_prism, Helmholtz3D_Prism)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-15));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-15);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-15);
     }
 }
 
@@ -97,8 +100,9 @@ BOOST_FIXTURE_TEST_CASE(nullpreconcuda_pyr, Helmholtz3D_Pyr)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-15));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-15);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-15);
     }
 }
 
@@ -115,8 +119,9 @@ BOOST_FIXTURE_TEST_CASE(nullpreconcuda_tet, Helmholtz3D_Tet)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-15));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-15);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-15);
     }
 }
 
diff --git a/tests/test_physderivcuda.cpp b/tests/test_physderivcuda.cpp
index 5cb37c6f821b03eaae812ece60b208728ac83b45..54a160488c6f9f4be34019e03f1b574a9ffda1df 100644
--- a/tests/test_physderivcuda.cpp
+++ b/tests/test_physderivcuda.cpp
@@ -24,8 +24,9 @@ BOOST_FIXTURE_TEST_CASE(physderivcuda_seg, Seg)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -42,8 +43,9 @@ BOOST_FIXTURE_TEST_CASE(physderivcuda_quad, Quad)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -60,8 +62,9 @@ BOOST_FIXTURE_TEST_CASE(physderivcuda_tri, Tri)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-12);
     }
 }
 
@@ -78,8 +81,9 @@ BOOST_FIXTURE_TEST_CASE(physderivcuda_square_all_elements, SquareAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-12));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
     }
 }
 
@@ -96,8 +100,9 @@ BOOST_FIXTURE_TEST_CASE(physderivcuda_hex, Hex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-10));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
     }
 }
 
@@ -114,8 +119,9 @@ BOOST_FIXTURE_TEST_CASE(physderivcuda_prism, Prism)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-10));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
     }
 }
 
@@ -132,8 +138,9 @@ BOOST_FIXTURE_TEST_CASE(physderivcuda_pyr, Pyr)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-10));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
     }
 }
 
@@ -150,8 +157,9 @@ BOOST_FIXTURE_TEST_CASE(physderivcuda_tet, Tet)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-10));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
     }
 }
 
@@ -168,8 +176,9 @@ BOOST_FIXTURE_TEST_CASE(physderivcuda_cube_prism_hex, CubePrismHex)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-10));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
     }
 }
 
@@ -186,8 +195,9 @@ BOOST_FIXTURE_TEST_CASE(physderivcuda_cube_all_elements, CubeAllElements)
     BOOST_TEST(fixtcuda_out->compare(*fixt_expected, 1.0E-10));
     boost::test_tools::output_test_stream output;
     {
-        OutputIfNotMatch(fixtcuda_out->GetStorage().GetCPUPtr(),
-                         fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
+        OutputIfNotMatch(
+            fixtcuda_out->template GetStorage<MemoryRegionCUDA>().GetCPUPtr(),
+            fixt_expected->GetStorage().GetCPUPtr(), 1.0E-10);
     }
 }