diff --git a/gpu/benchmark/Benchmark.hpp b/gpu/benchmark/Benchmark.hpp
index c9a7bbf6..31e67990 100644
--- a/gpu/benchmark/Benchmark.hpp
+++ b/gpu/benchmark/Benchmark.hpp
@@ -6,16 +6,36 @@
 #include <psDomain.hpp>
 #include <rayParticle.hpp>
 
-#define MAKE_GEO Trench
+#define MAKE_GEO Hole
 #define DEFAULT_GRID_DELTA 0.1
 #define DEFAULT_STICKING 0.1
-#define DIM 2
+#define DIM 3
 #define FIXED_RAYS false
 
-constexpr int particleType = 1;
 using TranslatorType = std::unordered_map<unsigned long, unsigned long>;
 using namespace viennaps;
 
+template <class NumericType, int N>
+consteval std::array<NumericType, N> linspace(NumericType start,
+                                              NumericType end) {
+  std::array<NumericType, N> arr{};
+  NumericType step = (end - start) / static_cast<NumericType>(N - 1);
+  for (int i = 0; i < N; ++i) {
+    arr[i] = start + i * step;
+  }
+  return arr;
+}
+
+constexpr int particleType = 0;
+constexpr bool runDisk = false;
+constexpr bool runLine = false;
+constexpr bool runTriangle = true;
+
+constexpr auto gridDeltaValues = linspace<float, 8>(0.08f, 0.4f);
+constexpr int numRuns = 10;
+constexpr int raysPerPoint = 1000;
+constexpr int numRays = int(1.4e8);
+
 template <class NumericType>
 auto Trench(NumericType gridDelta = DEFAULT_GRID_DELTA) {
   NumericType xExtent = 20.;
@@ -42,16 +62,6 @@ auto Hole(NumericType gridDelta = DEFAULT_GRID_DELTA) {
   return domain;
 }
 
-template <class NumericType, int N>
-std::array<NumericType, N> linspace(NumericType start, NumericType end) {
-  std::array<NumericType, N> arr{};
-  NumericType step = (end - start) / static_cast<NumericType>(N - 1);
-  for (int i = 0; i < N; ++i) {
-    arr[i] = start + i * step;
-  }
-  return arr;
-}
-
 template <typename NumericType> auto getIBEParameters() {
   viennaps::IBEParameters<NumericType> params;
   params.redepositionRate = 1.0;
@@ -199,8 +209,9 @@ void setupTriangleGeometry(
     elementKdTree_->setPoints(triangleCenters);
     elementKdTree_->build();
   } else {
-    triangleMesh = CreateTriangleMesh(
-        static_cast<float>(domain->getGridDelta()), surfaceMesh_);
+    CopyTriangleMesh(static_cast<float>(domain->getGridDelta()), surfaceMesh_,
+                     triangleMesh);
+    // preserves surfaceMesh_
   }
 
   rayTracer_.setGeometry(triangleMesh);
diff --git a/gpu/benchmark/CPU_Benchmark.cpp b/gpu/benchmark/CPU_Benchmark.cpp
index 9c5fa51c..1a253617 100644
--- a/gpu/benchmark/CPU_Benchmark.cpp
+++ b/gpu/benchmark/CPU_Benchmark.cpp
@@ -7,27 +7,18 @@
 
 #include "Benchmark.hpp"
 
-using namespace viennaps;
-
 int main() {
   omp_set_num_threads(16);
   using NumericType = float;
   constexpr int D = DIM;
 
-  // const std::array<NumericType, 4> gridDeltaValues = {0.05f, 0.1f, 0.2f,
-  // 0.4f};
-  auto gridDeltaValues = linspace<NumericType, 8>(0.01f, 0.4f);
-  const int numRuns = 20;
-  const int raysPerPoint = 1000;
-  const int numRays = int(1.4e8);
-
   auto particle = makeCPUParticle<NumericType, D>();
   std::vector<std::unique_ptr<viennaray::AbstractParticle<NumericType>>>
       particles;
   particles.push_back(particle->clone());
   std::string fluxLabel = particleType == 0 ? "flux" : "ionFlux";
 
-  { // Disk
+  if constexpr (runDisk) { // Disk
     std::ofstream file("CPU_Benchmark_Disk.txt");
     file << "Meshing;Tracing;Postprocessing;GridDelta\n";
 
@@ -112,7 +103,7 @@ int main() {
     file.close();
   }
 
-  { // Triangle
+  if constexpr (runTriangle) { // Triangle
     std::ofstream file("CPU_Benchmark_Triangle.txt");
     file << "Meshing;Tracing;Postprocessing;GridDelta\n";
 
@@ -125,7 +116,7 @@ int main() {
     tracer.setUseRandomSeeds(false);
     tracer.setParticleType(particle);
 
-    std::vector<std::string> dataLabels = particle->getLocalDataLabels();
+    const auto &dataLabels = particle->getLocalDataLabels();
 
     for (int i = 0; i < gridDeltaValues.size(); i++) {
       std::cout << "  Grid Delta: " << gridDeltaValues[i] << "\n";
diff --git a/gpu/benchmark/GPU_Benchmark.cpp b/gpu/benchmark/GPU_Benchmark.cpp
index 1df0f15b..9379c11d 100644
--- a/gpu/benchmark/GPU_Benchmark.cpp
+++ b/gpu/benchmark/GPU_Benchmark.cpp
@@ -13,21 +13,10 @@
 
 #include "Benchmark.hpp"
 
-using namespace viennaps;
-
 int main() {
   omp_set_num_threads(16);
   using NumericType = float;
   constexpr int D = DIM;
-
-  const NumericType cSticking = DEFAULT_STICKING;
-  // const std::array<NumericType, 4> gridDeltaValues = {0.05f, 0.1f, 0.2f,
-  // 0.4f};
-  auto gridDeltaValues = linspace<NumericType, 8>(0.01f, 0.4f);
-  const int numRuns = 20;
-  const int raysPerPoint = 1000;
-  const int numRays = int(1.4e8);
-
   auto context = DeviceContext::createContext();
 
   CudaBuffer deviceParamsBuffer;
@@ -36,7 +25,7 @@ int main() {
     deviceParamsBuffer.allocUploadSingle(deviceParams);
   }
 
-  { // Triangle
+  if constexpr (runTriangle) { // Triangle
     std::ofstream file("GPU_Benchmark_Triangle.txt");
     file << "Meshing;Tracing;Postprocessing;GridDelta\n";
 
@@ -55,8 +44,7 @@ int main() {
     }
     tracer.prepareParticlePrograms();
 
-    std::vector<std::string> dataLabels =
-        std::get<0>(particleConfig).dataLabels;
+    const auto &dataLabels = std::get<0>(particleConfig).dataLabels;
 
     std::cout << "Starting Triangle Benchmark\n";
 
@@ -138,7 +126,7 @@ int main() {
     file.close();
   }
 
-  { // Disk
+  if constexpr (runDisk) { // Disk
     std::ofstream file("GPU_Benchmark_Disk.txt");
     file << "Meshing;Tracing;Postprocessing;GridDelta\n";
 
@@ -186,7 +174,6 @@ int main() {
       for (int j = 0; j < numRuns; j++) {
         std::cout << "    Process Step: " << j + 1 << "\n";
         advectionKernel.prepareLS();
-        file << cSticking << ";";
 
         Timer timer;
 
@@ -235,7 +222,7 @@ int main() {
     file.close();
   }
 
-  if constexpr (D == 2) { // Line
+  if constexpr (D == 2 && runLine) { // Line
     std::ofstream file("GPU_Benchmark_Line.txt");
     file << "Meshing;Tracing;Postprocessing;GridDelta\n";
 
diff --git a/gpu/benchmark/holeCompare.cpp b/gpu/benchmark/holeCompare.cpp
index 94a0fe11..e2e73e8d 100644
--- a/gpu/benchmark/holeCompare.cpp
+++ b/gpu/benchmark/holeCompare.cpp
@@ -23,7 +23,7 @@ void runBenchmark(SmartPointer<Domain<NumericType, D>> geometry,
   process.apply();
   timer.finish();
 
-  std::cout << "Flux Engine: " << to_string(fluxEngine)
+  std::cout << "Flux Engine: " << util::toString(fluxEngine)
             << ", Time taken: " << timer.currentDuration / 1e9 << " s"
             << std::endl;
   copy->saveSurfaceMesh(outputFilename);
@@ -61,7 +61,7 @@ int main(int argc, char *argv[]) {
 
     for (const auto &fluxEngine : fluxEngines) {
       runBenchmark<NumericType, D>(geometry, model, fluxEngine,
-                                   "result_PE_" + to_string(fluxEngine) +
+                                   "result_PE_" + util::toString(fluxEngine) +
                                        ".vtp");
     }
   }
@@ -72,7 +72,7 @@ int main(int argc, char *argv[]) {
 
     for (const auto &fluxEngine : fluxEngines) {
       runBenchmark<NumericType, D>(geometry, model, fluxEngine,
-                                   "result_SP_" + to_string(fluxEngine) +
+                                   "result_SP_" + util::toString(fluxEngine) +
                                        ".vtp");
     }
   }
@@ -90,7 +90,7 @@ int main(int argc, char *argv[]) {
 
     for (const auto &fluxEngine : fluxEngines) {
       runBenchmark<NumericType, D>(geometry, model, fluxEngine,
-                                   "result_IBE_" + to_string(fluxEngine) +
+                                   "result_IBE_" + util::toString(fluxEngine) +
                                        ".vtp");
     }
   }
@@ -114,7 +114,7 @@ int main(int argc, char *argv[]) {
 
     for (const auto &fluxEngine : fluxEngines) {
       runBenchmark<NumericType, D>(geometry, model, fluxEngine,
-                                   "result_MP_" + to_string(fluxEngine) +
+                                   "result_MP_" + util::toString(fluxEngine) +
                                        ".vtp");
     }
   }
diff --git a/gpu/benchmark/holeCompareGrid.cpp b/gpu/benchmark/holeCompareGrid.cpp
index b77ecdc3..be4052db 100644
--- a/gpu/benchmark/holeCompareGrid.cpp
+++ b/gpu/benchmark/holeCompareGrid.cpp
@@ -28,7 +28,7 @@ void runBenchmark(SmartPointer<Domain<NumericType, D>> geometry,
   process.apply();
   timer.finish();
 
-  std::cout << "Flux Engine: " << to_string(fluxEngine)
+  std::cout << "Flux Engine: " << util::toString(fluxEngine)
             << ", Time taken: " << timer.currentDuration / 1e9 << " s"
             << std::endl;
   copy->saveSurfaceMesh(outputFilename);
@@ -65,7 +65,7 @@ int main(int argc, char *argv[]) {
 
       for (const auto &fluxEngine : fluxEngines) {
         runBenchmark<NumericType, D>(geometry, model, fluxEngine,
-                                     "result_PE_" + to_string(fluxEngine) +
+                                     "result_PE_" + util::toString(fluxEngine) +
                                          "_" + std::to_string(i) + ".vtp",
                                      0.2);
       }
@@ -92,8 +92,9 @@ int main(int argc, char *argv[]) {
 
       for (const auto &fluxEngine : fluxEngines) {
         runBenchmark<NumericType, D>(geometry, model, fluxEngine,
-                                     "result_IBE_" + to_string(fluxEngine) +
-                                         "_" + std::to_string(i) + ".vtp");
+                                     "result_IBE_" +
+                                         util::toString(fluxEngine) + "_" +
+                                         std::to_string(i) + ".vtp");
       }
     }
   }
diff --git a/include/viennaps/process/psCPUTriangleEngine.hpp b/include/viennaps/process/psCPUTriangleEngine.hpp
index 7a5bc47d..f51f220d 100644
--- a/include/viennaps/process/psCPUTriangleEngine.hpp
+++ b/include/viennaps/process/psCPUTriangleEngine.hpp
@@ -114,8 +114,8 @@ class CPUTriangleEngine final : public FluxEngine<NumericType, D> {
       elementKdTree_->setPoints(triangleCenters);
       elementKdTree_->build();
     } else {
-      triangleMesh = CreateTriangleMesh(
-          static_cast<float>(context.domain->getGridDelta()), surfaceMesh_);
+      CopyTriangleMesh(static_cast<float>(context.domain->getGridDelta()),
+                       surfaceMesh_, triangleMesh);
     }
 
     rayTracer_.setGeometry(triangleMesh);
diff --git a/include/viennaps/process/psGPUTriangleEngine.hpp b/include/viennaps/process/psGPUTriangleEngine.hpp
index 8fc095a6..95bbc66f 100644
--- a/include/viennaps/process/psGPUTriangleEngine.hpp
+++ b/include/viennaps/process/psGPUTriangleEngine.hpp
@@ -136,8 +136,9 @@ class GPUTriangleEngine final : public FluxEngine<NumericType, D> {
       elementKdTree_->setPoints(triangleCenters);
       elementKdTree_->build();
     } else {
-      triangleMesh = CreateTriangleMesh(
-          static_cast<float>(context.domain->getGridDelta()), surfaceMesh_);
+      CopyTriangleMesh(static_cast<float>(context.domain->getGridDelta()),
+                       surfaceMesh_, triangleMesh);
+      // preserves surfaceMesh_
     }
 
     rayTracer_.setGeometry(triangleMesh);
@@ -228,7 +229,8 @@ class GPUTriangleEngine final : public FluxEngine<NumericType, D> {
         downloadCoverages(d_coverages, surfaceMesh_->getCellData(), coverages,
                           surfaceMesh_->getElements<3>().size());
       }
-      saveResultsToPointData(surfaceMesh_->getCellData());
+      saveResultsToPointData(
+          surfaceMesh_->getCellData()); // save fluxes to elements
       viennals::VTKWriter<float>(
           surfaceMesh_, context.intermediateOutputPath +
                             context.getProcessName() + "_flux_" +
@@ -269,22 +271,20 @@ class GPUTriangleEngine final : public FluxEngine<NumericType, D> {
     assert(numPoints == surfaceMesh_->getElements<3>().size());
     auto const &results = rayTracer_.getResults();
     auto particles = rayTracer_.getParticles();
+    const auto &dataLabels = model_->getParticleDataLabels();
 
-    int offset = 0;
-    for (int pIdx = 0; pIdx < particles.size(); pIdx++) {
-      for (int dIdx = 0; dIdx < particles[pIdx].dataLabels.size(); dIdx++) {
-        auto name = particles[pIdx].dataLabels[dIdx];
-        assert(offset + dIdx < results.size());
-        const auto &data = results[offset + dIdx];
+    assert(dataLabels.size() == results.size());
 
-        std::vector<float> values(numPoints);
-        for (unsigned i = 0; i < numPoints; ++i) {
-          values[i] = static_cast<float>(data[i]);
-        }
+    for (int dIdx = 0; dIdx < dataLabels.size(); dIdx++) {
+      const auto &name = dataLabels[dIdx];
+      const auto &data = results[dIdx];
 
-        pointData.insertReplaceScalarData(std::move(values), name);
+      std::vector<float> values(numPoints);
+      for (unsigned i = 0; i < numPoints; ++i) {
+        values[i] = static_cast<float>(data[i]);
       }
-      offset += particles[pIdx].dataLabels.size();
+
+      pointData.insertReplaceScalarData(std::move(values), name);
     }
   }
 
diff --git a/include/viennaps/process/psProcess.hpp b/include/viennaps/process/psProcess.hpp
index 154e3886..360c42c3 100644
--- a/include/viennaps/process/psProcess.hpp
+++ b/include/viennaps/process/psProcess.hpp
@@ -99,7 +99,7 @@ template <typename NumericType, int D> class Process {
     VIENNACORE_LOG_DEBUG("Using strategy: " + std::string(strategy->name()));
 
     if (strategy->requiresFluxEngine()) {
-      VIENNACORE_LOG_DEBUG("Setting up " + to_string(fluxEngineType_) +
+      VIENNACORE_LOG_DEBUG("Setting up " + util::toString(fluxEngineType_) +
                            " flux engine for strategy.");
       strategy->setFluxEngine(createFluxEngine());
     }
@@ -191,7 +191,7 @@ template <typename NumericType, int D> class Process {
     assert(fluxEngineType_ != FluxEngineType::AUTO &&
            "Flux engine type must be specified before creation.");
     VIENNACORE_LOG_DEBUG("Creating flux engine of type: " +
-                         to_string(fluxEngineType_));
+                         util::toString(fluxEngineType_));
     // Create CPU engine
     if (fluxEngineType_ == FluxEngineType::CPU_DISK) {
       return std::make_unique<CPUDiskEngine<NumericType, D>>();
@@ -295,7 +295,7 @@ template <typename NumericType, int D> class Process {
         fluxEngineType_ = FluxEngineType::CPU_DISK;
       }
       VIENNACORE_LOG_DEBUG("Auto-selected flux engine type: " +
-                           to_string(fluxEngineType_));
+                           util::toString(fluxEngineType_));
     }
 
 #ifdef VIENNACORE_COMPILE_GPU
diff --git a/include/viennaps/process/psProcessModel.hpp b/include/viennaps/process/psProcessModel.hpp
index 2ed6a471..b82274b7 100644
--- a/include/viennaps/process/psProcessModel.hpp
+++ b/include/viennaps/process/psProcessModel.hpp
@@ -174,8 +174,6 @@ class ProcessModelCPU : public ProcessModelBase<NumericType, D> {
 #ifdef VIENNACORE_COMPILE_GPU
 namespace viennaps::gpu {
 
-using namespace viennacore;
-
 template <class NumericType, int D>
 class ProcessModelGPU : public ProcessModelBase<NumericType, D> {
 private:
@@ -187,7 +185,7 @@ class ProcessModelGPU : public ProcessModelBase<NumericType, D> {
   std::vector<viennaray::gpu::CallableConfig> callableMap_;
 
 public:
-  CudaBuffer processData;
+  viennacore::CudaBuffer processData;
   auto &getParticleTypes() { return particles; }
   auto getProcessDataDPtr() const { return processData.dPointer(); }
   bool useMaterialIds() const { return materialIds; }
diff --git a/include/viennaps/psCreateSurfaceMesh.hpp b/include/viennaps/psCreateSurfaceMesh.hpp
index fa59d06c..e850ed37 100644
--- a/include/viennaps/psCreateSurfaceMesh.hpp
+++ b/include/viennaps/psCreateSurfaceMesh.hpp
@@ -24,7 +24,18 @@ CreateTriangleMesh(const float gridDelta,
   triangleMesh.maximumExtent = mesh->maximumExtent;
   triangleMesh.normals = *mesh->getCellData().getVectorData("Normals");
 
-  return triangleMesh;
+  return std::move(triangleMesh);
+}
+
+inline void CopyTriangleMesh(const float gridDelta,
+                             const SmartPointer<viennals::Mesh<float>> &mesh,
+                             viennaray::TriangleMesh &triangleMesh) {
+  triangleMesh.gridDelta = gridDelta;
+  triangleMesh.triangles = mesh->triangles;
+  triangleMesh.nodes = mesh->nodes;
+  triangleMesh.minimumExtent = mesh->minimumExtent;
+  triangleMesh.maximumExtent = mesh->maximumExtent;
+  triangleMesh.normals = *mesh->getCellData().getVectorData("Normals");
 }
 
 template <class LsNT, class MeshNT = LsNT, int D = 3> class CreateSurfaceMesh {
diff --git a/include/viennaps/psUtil.hpp b/include/viennaps/psUtil.hpp
index 2049d42a..65a5b92d 100644
--- a/include/viennaps/psUtil.hpp
+++ b/include/viennaps/psUtil.hpp
@@ -191,11 +191,9 @@ template <typename T> [[nodiscard]] std::string toString(const T &value) {
   if constexpr (std::is_same_v<T, bool>)
     return value ? "true" : "false";
   else if constexpr (std::is_same_v<T, viennals::SpatialSchemeEnum>)
-    return convertSpatialSchemeToString(
-        static_cast<viennals::SpatialSchemeEnum>(value));
+    return convertSpatialSchemeToString(value);
   else if constexpr (std::is_same_v<T, viennals::TemporalSchemeEnum>)
-    return convertTemporalSchemeToString(
-        static_cast<viennals::TemporalSchemeEnum>(value));
+    return convertTemporalSchemeToString(value);
   else if constexpr (std::is_same_v<T, viennaps::Material>) {
     std::string mat = viennaps::to_string_view(value);
     return mat;
@@ -228,9 +226,3 @@ hexToRGBArray(const uint32_t hexColor) {
   return rgb;
 }
 }; // namespace viennacore::util
-
-namespace viennaps {
-[[nodiscard]] inline std::string to_string(FluxEngineType type) {
-  return viennacore::util::convertFluxEngineTypeToString(type);
-}
-} // namespace viennaps
diff --git a/tests/fluxEngines/fluxEngines.cpp b/tests/fluxEngines/fluxEngines.cpp
index c7b3a434..69d5dc06 100644
--- a/tests/fluxEngines/fluxEngines.cpp
+++ b/tests/fluxEngines/fluxEngines.cpp
@@ -40,7 +40,7 @@ template <class T, int D> void RunTest() {
           1.0, 1.0, 1.0);
 
   for (const auto &engineType : engineTypes) {
-    std::cout << "Testing flux engine: " << viennaps::to_string(engineType)
+    std::cout << "Testing flux engine: " << viennaps::util::toString(engineType)
               << " in " << D << "D" << std::endl;
     auto domain = viennaps::Domain<T, D>::New(1.0, 10.0, 10.0);
 
@@ -58,7 +58,7 @@ template <class T, int D> void RunTest() {
 
     domain->saveSurfaceMesh("fluxEngineTest_" +
                             std::to_string(static_cast<int>(D)) + "D_" +
-                            viennaps::to_string(engineType) + ".vtp");
+                            viennaps::util::toString(engineType) + ".vtp");
 
     checkSurfaceHeight<T, D>(domain, 1.0);
   }
diff --git a/tests/rngSeed/rngSeed.cpp b/tests/rngSeed/rngSeed.cpp
index 2fbf0008..2c634acc 100644
--- a/tests/rngSeed/rngSeed.cpp
+++ b/tests/rngSeed/rngSeed.cpp
@@ -33,8 +33,8 @@ template <typename T, int D> void testRngSeed() {
                                     // here, also when rngSeed is set
 
   for (const auto &engineType : engineTypes) {
-    std::cout << "Testing rng seed: " << to_string(engineType) << " in " << D
-              << "D" << std::endl;
+    std::cout << "Testing rng seed: " << util::toString(engineType) << " in "
+              << D << "D" << std::endl;
 
     std::vector<T> flux1, flux2;
     {