diff --git a/.github/workflows/linux_test.yml b/.github/workflows/linux_test.yml
index fe31fcfa..b6fd7108 100644
--- a/.github/workflows/linux_test.yml
+++ b/.github/workflows/linux_test.yml
@@ -85,11 +85,15 @@ jobs:
           cmake -B ./build \
             -D ViennaLS_DIR=${{github.workspace}}/viennals/install/lib/cmake/ViennaLS/ \
             -D ViennaRay_DIR=${{github.workspace}}/viennaray/install/lib/cmake/ViennaRay/ \
-            -D VIENNAPS_BUILD_EXAMPLES=ON
+            -D VIENNAPS_BUILD_EXAMPLES=ON \
+            -D VIENNAPS_BUILD_APPLICATION=ON
 
       - name: Build Examples
         run: cmake --build ./build --target buildExamples
 
+      - name: Build application
+        run: cmake --build ./build --target buildApplication
+
       - name: Run
         working-directory: ${{github.workspace}}/build/Examples/InterpolationDemo
         run: ./InterpolationDemo
diff --git a/.github/workflows/macos_test.yml b/.github/workflows/macos_test.yml
index e94b426d..179de99e 100644
--- a/.github/workflows/macos_test.yml
+++ b/.github/workflows/macos_test.yml
@@ -101,11 +101,15 @@ jobs:
             -D OpenMP_C_LIB_NAMES="omp" \
             -D OpenMP_CXX_LIB_NAMES="omp" \
             -D OpenMP_omp_LIBRARY="$(brew --prefix libomp)/lib/libomp.a" \
-            -D VIENNAPS_BUILD_EXAMPLES=ON
+            -D VIENNAPS_BUILD_EXAMPLES=ON \
+            -D VIENNAPS_BUILD_APPLICATION=ON
 
       - name: Build Examples
         run: cmake --build ./build --target buildExamples
 
+      - name: Build application
+        run: cmake --build ./build --target buildApplication
+
       - name: Run
         working-directory: ${{github.workspace}}/build/Examples/InterpolationDemo
         run: ./InterpolationDemo
diff --git a/.github/workflows/windows_test.yml b/.github/workflows/windows_test.yml
index 8d2ac88f..98c5c206 100644
--- a/.github/workflows/windows_test.yml
+++ b/.github/workflows/windows_test.yml
@@ -83,11 +83,14 @@ jobs:
       # ViennaPS
       - name: Configure ViennaPS
         run: |
-          cmake -B ${{github.workspace}}/build -DViennaLS_DIR=${{github.workspace}}/viennals/install/lib/cmake/ViennaLS/ -DViennaRay_DIR=${{github.workspace}}/viennaray/install/lib/cmake/ViennaRay/ -D CMAKE_CXX_FLAGS="-openmp:llvm" -D VIENNAPS_BUILD_EXAMPLES=ON
+          cmake -B ${{github.workspace}}/build -DViennaLS_DIR=${{github.workspace}}/viennals/install/lib/cmake/ViennaLS/ -DViennaRay_DIR=${{github.workspace}}/viennaray/install/lib/cmake/ViennaRay/ -D CMAKE_CXX_FLAGS="-openmp:llvm" -D VIENNAPS_BUILD_EXAMPLES=ON -D VIENNAPS_BUILD_APPLICATION=ON
 
       - name: Build ViennaPS Examples
         run: cmake --build ${{github.workspace}}/build/Examples --config ${{env.BUILD_TYPE}}
 
+      - name: Build ViennaPS Application
+        run: cmake --build ${{github.workspace}}/build/app --config ${{env.BUILD_TYPE}}
+
       - name: Run ViennaPS Example
         working-directory: ${{github.workspace}}/build/Examples/InterpolationDemo/${{env.BUILD_TYPE}}
         shell: bash
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 018cd252..1c59f5e0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -77,6 +77,14 @@ if(VIENNAPS_BUILD_EXAMPLES)
   add_subdirectory(Examples)
 endif(VIENNAPS_BUILD_EXAMPLES)
 
+# ##################################################################################################
+# BUILD VIENNAPS APPLICATION
+# ##################################################################################################
+option(VIENNAPS_BUILD_APPLICATION "Build ViennaPS application." OFF)
+if(VIENNAPS_BUILD_APPLICATION)
+  add_subdirectory(app)
+endif(VIENNAPS_BUILD_APPLICATION)
+
 # ##################################################################################################
 # INSTALL
 # ##################################################################################################
diff --git a/Examples/GDSReader/GDSReader.cpp b/Examples/GDSReader/GDSReader.cpp
index 25d1c59c..dc98d4c8 100644
--- a/Examples/GDSReader/GDSReader.cpp
+++ b/Examples/GDSReader/GDSReader.cpp
@@ -22,7 +22,7 @@ int main(int argc, char **argv) {
   psGDSReader<NumericType, D>(mask, "mask.gds").apply();
 
   // geometry setup
-  NumericType *bounds = mask->getBounds();
+  double *bounds = mask->getBounds();
   auto geometry = psSmartPointer<psDomain<NumericType, D>>::New();
 
   // substrate plane
diff --git a/Examples/HoleEtching/HoleEtching.cpp b/Examples/HoleEtching/HoleEtching.cpp
index 3a6c5c84..1fb44e9b 100644
--- a/Examples/HoleEtching/HoleEtching.cpp
+++ b/Examples/HoleEtching/HoleEtching.cpp
@@ -28,7 +28,8 @@ int main(int argc, char *argv[]) {
       geometry, params.gridDelta /* grid delta */, params.xExtent /*x extent*/,
       params.yExtent /*y extent*/, params.holeRadius /*hole radius*/,
       params.maskHeight /* mask height*/,
-      params.taperAngle /* tapering angle in degrees */, true /*create mask*/)
+      params.taperAngle /* tapering angle in degrees */, 0 /* base height */,
+      false /* periodic boundary */, true /*create mask*/)
       .apply();
 
   SF6O2Etching<NumericType, D> model(params.totalIonFlux /*ion flux*/,
diff --git a/Examples/TrenchDeposition/TrenchDeposition.cpp b/Examples/TrenchDeposition/TrenchDeposition.cpp
index ba3126cb..25cbdf33 100644
--- a/Examples/TrenchDeposition/TrenchDeposition.cpp
+++ b/Examples/TrenchDeposition/TrenchDeposition.cpp
@@ -28,7 +28,7 @@ int main(int argc, char *argv[]) {
       geometry, params.gridDelta /* grid delta */, params.xExtent /*x extent*/,
       params.yExtent /*y extent*/, params.trenchWidth /*trench width*/,
       params.trenchHeight /*trench height*/,
-      params.taperAngle /* tapering angle */, false /*create mask*/)
+      params.taperAngle /* tapering angle */)
       .apply();
 
   // copy top layer to capture deposition
diff --git a/Examples/TrenchDepositionGeometric/TrenchDepositionGeometric.cpp b/Examples/TrenchDepositionGeometric/TrenchDepositionGeometric.cpp
index 38787275..ad069766 100644
--- a/Examples/TrenchDepositionGeometric/TrenchDepositionGeometric.cpp
+++ b/Examples/TrenchDepositionGeometric/TrenchDepositionGeometric.cpp
@@ -26,7 +26,7 @@ int main(int argc, char *argv[]) {
   auto geometry = psSmartPointer<psDomain<NumericType, D>>::New();
   psMakeTrench<NumericType, D>(geometry, params.gridDelta, params.xExtent,
                                params.yExtent, params.trenchWidth,
-                               params.trenchHeight, false /*create mask*/)
+                               params.trenchHeight)
       .apply();
 
   // copy top layer to capture deposition
diff --git a/Examples/VolumeModel/VolumeModel.cpp b/Examples/VolumeModel/VolumeModel.cpp
index 8ca38a9b..6612114f 100644
--- a/Examples/VolumeModel/VolumeModel.cpp
+++ b/Examples/VolumeModel/VolumeModel.cpp
@@ -24,6 +24,7 @@ int main(int argc, char *argv[]) {
   auto geometry = psSmartPointer<psDomain<NumericType, D>>::New();
   psMakeFin<NumericType, D>(geometry, params.gridDelta, params.xExtent,
                             params.yExtent, params.finWidth, params.finHeight,
+                            0. /* base height*/, false /*periodic boundary*/,
                             false /*create mask*/)
       .apply();
 
diff --git a/README.md b/README.md
index ee081ffc..7ee31a4e 100644
--- a/README.md
+++ b/README.md
@@ -68,6 +68,7 @@ The examples can be built using CMake (make sure all dependencies are installed/
 ```bash
 mkdir build && cd build
 cmake .. -DVIENNAPS_BUILD_EXAMPLES=ON
+make buildExamples
 ```
 
 The examples can then be executed in their respective folders with the config files
@@ -96,13 +97,26 @@ By changing the dimension of the hole etching example (_D = 2_), we can easily s
   <img src="https://raw.githubusercontent.com/ViennaTools/ViennaPS/master/data/images/sidewall_tapering.svg" width=700 style="background-color:white;">
 </div>
 
+## Application
+
+It is also possible to build an application which can parse a custom configuration file and execute pre-defined processes. The application can be built using CMake (make sure all dependencies are installed/ have been built previously):
+```bash
+mkdir build && cd build
+cmake .. -DVIENNAPS_BUILD_APPLICATION=ON
+make buildApplication
+```
+This creates 2 executables `ViennaPS2D` and `ViennaPS3D` which run processes in 2 or 3 dimensions respectively. Every configuration file can be run in 2D or 3D mode.
+
+The configuration file must obey a certain structure in order to be parsed correctly. An example for a configuration file can be seen in _SampleConfig.txt_. The configuration file is parsed line by line and each succesfully parsed line is executed immediately. A detailed documentation for the configuration file can be found in **app/README.md**.
+
+
 ## Contributing
 
 If you want to contribute to ViennaPS, make sure to follow the [LLVM Coding guidelines](https://llvm.org/docs/CodingStandards.html). Before creating a pull request, make sure ALL files have been formatted by clang-format, which can be done using the format-project.sh script in the root directory.
 
 ## Authors
 
-Current contributors: Tobias Reiter, Josip Bobinac, Xaver Klemenschits, Julius Piso
+Current contributors: Tobias Reiter, Julius Piso, Josip Bobinac, Xaver Klemenschits
 
 Contact us via: viennatools@iue.tuwien.ac.at
 
diff --git a/app/Application.hpp b/app/Application.hpp
new file mode 100644
index 00000000..ec09d170
--- /dev/null
+++ b/app/Application.hpp
@@ -0,0 +1,368 @@
+#pragma once
+
+#include <fstream>
+#include <sstream>
+
+#include <lsReader.hpp>
+
+#include <psDomain.hpp>
+#include <psGDSReader.hpp>
+#include <psMakeHole.hpp>
+#include <psMakePlane.hpp>
+#include <psMakeTrench.hpp>
+#include <psPlanarize.hpp>
+#include <psProcess.hpp>
+#include <psUtils.hpp>
+
+#include <ApplicationParameters.hpp>
+#include <ApplicationParser.hpp>
+
+#include <DirectionalEtching.hpp>
+#include <GeometricUniformDeposition.hpp>
+#include <SF6O2Etching.hpp>
+#include <SimpleDeposition.hpp>
+
+template <int D> class Application {
+  psSmartPointer<psDomain<NumericType, D>> geometry = nullptr;
+  psSmartPointer<ApplicationParameters> params = nullptr;
+  ApplicationParser parser;
+  int clArgC = 0;
+  char **clArgV;
+
+public:
+  Application(int argc, char **argv) : clArgC(argc), clArgV(argv) {}
+
+  void run() {
+    if (clArgC < 2) {
+      std::cout << "No input file specified. " << std::endl;
+      return;
+    }
+
+    std::fstream inputFile;
+    inputFile.open(clArgV[1], std::fstream::in);
+
+    if (!inputFile.is_open()) {
+      std::cout << "Could not open input file." << std::endl;
+      return;
+    }
+
+    params = psSmartPointer<ApplicationParameters>::New();
+    params->defaultParameters(true);
+    parser.setParameters(params);
+
+    int lineNumber = 0;
+    std::string line;
+    while (std::getline(inputFile, line)) {
+      if (line.empty() || line[0] == '#') // skipping empty line and comments
+        continue;
+      std::istringstream lineStream(line);
+      switch (parser.parseCommand(lineStream, lineNumber)) {
+      case CommandType::INIT:
+        runInit();
+        break;
+
+      case CommandType::GEOMETRY:
+        createGeometry();
+        break;
+
+      case CommandType::PROCESS:
+        runProcess();
+        break;
+
+      case CommandType::PLANARIZE:
+        planarizeGeometry();
+        break;
+
+      case CommandType::OUTPUT:
+        writeVTP();
+        break;
+
+      case CommandType::NONE:
+        break;
+
+      default:
+        assert(false);
+      }
+
+      lineNumber++;
+      params->defaultParameters();
+    }
+  }
+
+protected:
+  virtual void
+  runSimpleDeposition(psSmartPointer<psDomain<NumericType, D>> processGeometry,
+                      psSmartPointer<ApplicationParameters> processParams) {
+    SimpleDeposition<NumericType, D> model(processParams->sticking,
+                                           processParams->cosinePower);
+
+    psProcess<NumericType, D> process;
+    process.setDomain(processGeometry);
+    process.setProcessModel(model.getProcessModel());
+    process.setNumberOfRaysPerPoint(processParams->raysPerPoint);
+    process.setProcessDuration(processParams->processTime *
+                               processParams->rate / processParams->sticking);
+    process.setPrintIntdermediate(params->printIntermediate);
+    process.apply();
+  }
+
+  virtual void
+  runSF6O2Etching(psSmartPointer<psDomain<NumericType, D>> processGeometry,
+                  psSmartPointer<ApplicationParameters> processParams) {
+    SF6O2Etching<NumericType, D> model(
+        processParams->totalIonFlux, processParams->totalEtchantFlux,
+        processParams->totalOxygenFlux, processParams->ionEnergy,
+        processParams->A_O, processParams->maskId);
+
+    psProcess<NumericType, D> process;
+    process.setDomain(processGeometry);
+    process.setProcessModel(model.getProcessModel());
+    process.setMaxCoverageInitIterations(10);
+    process.setNumberOfRaysPerPoint(processParams->raysPerPoint);
+    process.setProcessDuration(processParams->processTime);
+    process.setPrintIntdermediate(params->printIntermediate);
+    process.apply();
+  }
+
+  virtual void runGeometriyUniformDeposition(
+      psSmartPointer<psDomain<NumericType, D>> processGeometry,
+      psSmartPointer<ApplicationParameters> processParams) {
+    GeometricUniformDeposition<NumericType, D> model(
+        processParams->processTime * processParams->rate);
+
+    psProcess<NumericType, D> process;
+    process.setDomain(processGeometry);
+    process.setProcessModel(model.getProcessModel());
+    process.setPrintIntdermediate(params->printIntermediate);
+    process.apply();
+  }
+
+  virtual void runDirectionalEtching(
+      psSmartPointer<psDomain<NumericType, D>> processGeometry,
+      psSmartPointer<ApplicationParameters> processParams) {
+
+    DirectionalEtching<NumericType, D> model(
+        getDirection(processParams->direction), processParams->directionalRate,
+        processParams->isotropicRate);
+
+    psProcess<NumericType, D> process;
+    process.setDomain(processGeometry);
+    process.setProcessModel(model.getProcessModel());
+    process.setProcessDuration(params->processTime);
+    process.setPrintIntdermediate(params->printIntermediate);
+    process.apply();
+  }
+
+private:
+  void runInit() {
+    std::cout << "\tx-Extent: " << params->xExtent
+              << "\n\ty-Extent: " << params->yExtent
+              << "\n\tResolution: " << params->gridDelta
+              << "\n\tPrint intermediate: "
+              << boolString(params->printIntermediate)
+              << "\n\tPeriodic boundary: "
+              << boolString(params->periodicBoundary) << "\n\n";
+
+    geometry = psSmartPointer<psDomain<NumericType, D>>::New();
+  }
+
+  void createGeometry() {
+    std::cout << "\tGeometry type: ";
+    switch (params->geometryType) {
+    case GeometryType::TRENCH:
+      std::cout << "Trench\n\tWidth: " << params->trenchWidth
+                << "\n\tHeight: " << params->trenchHeight
+                << "\n\tzPos: " << params->maskZPos
+                << "\n\tTapering angle: " << params->taperAngle
+                << "\n\tMask: " << boolString(params->mask) << "\n\n";
+      psMakeTrench<NumericType, D>(
+          geometry, params->gridDelta, params->xExtent, params->yExtent,
+          params->trenchWidth, params->trenchHeight, params->taperAngle,
+          params->maskZPos, params->periodicBoundary, params->mask)
+          .apply();
+      break;
+
+    case GeometryType::HOLE:
+      std::cout << "Hole\n\tRadius: " << params->holeRadius
+                << "\n\tDepth: " << params->holeDepth
+                << "\n\tzPos: " << params->maskZPos
+                << "\n\tTapering angle: " << params->taperAngle
+                << "\n\tMask: " << boolString(params->mask) << "\n\n";
+      psMakeHole<NumericType, D>(
+          geometry, params->gridDelta, params->xExtent, params->yExtent,
+          params->holeRadius, params->holeDepth, params->taperAngle,
+          params->maskZPos, params->periodicBoundary, params->mask)
+          .apply();
+      break;
+
+    case GeometryType::PLANE:
+      std::cout << "Plane"
+                << "\n\tzPos: " << params->maskZPos << "\n\n";
+      psMakePlane<NumericType, D>(geometry, params->gridDelta, params->xExtent,
+                                  params->yExtent, params->maskZPos,
+                                  params->periodicBoundary)
+          .apply();
+      break;
+    case GeometryType::GDS: {
+      std::cout << "GDS file import\n\tFile name: " << params->fileName
+                << "\n\tLayer: " << params->layers
+                << "\n\tMask height: " << params->maskHeight
+                << "\n\tzPos: " << params->maskZPos << "\n\n";
+
+      typename lsDomain<NumericType, D>::BoundaryType boundaryCons[D];
+      for (int i = 0; i < D - 1; i++) {
+        if (params->periodicBoundary) {
+          boundaryCons[i] =
+              lsDomain<NumericType, D>::BoundaryType::PERIODIC_BOUNDARY;
+        } else {
+          boundaryCons[i] =
+              lsDomain<NumericType, D>::BoundaryType::REFLECTIVE_BOUNDARY;
+        }
+      }
+      boundaryCons[D - 1] =
+          lsDomain<NumericType, D>::BoundaryType::INFINITE_BOUNDARY;
+      auto mask =
+          psSmartPointer<psGDSGeometry<NumericType, D>>::New(params->gridDelta);
+      mask->setBoundaryConditions(boundaryCons);
+      psGDSReader<NumericType, D>(mask, params->fileName).apply();
+      auto layer = mask->layerToLevelSet(params->layers, params->maskZPos,
+                                         params->maskHeight);
+      geometry->insertNextLevelSet(layer);
+      break;
+    }
+
+    case GeometryType::IMPORT: {
+      std::cout << "ViennaLS file import\n\tFile name: " << params->fileName
+                << "\n\tNumber of layers: " << params->layers << "\n";
+      for (int i = 0; i < params->layers; i++) {
+        std::string layerFileName =
+            params->fileName + "_layer" + std::to_string(i) + ".lvst";
+        std::cout << "\tReading " << layerFileName << std::endl;
+        auto layer = psSmartPointer<lsDomain<NumericType, D>>::New();
+        lsReader<NumericType, D>(layer, layerFileName).apply();
+        if (!geometry->getLevelSets()->empty() &&
+            layer->getGrid().getGridDelta() !=
+                geometry->getLevelSets()->back()->getGrid().getGridDelta()) {
+          std::cout << std::setprecision(8);
+          std::cout << "Import geometry grid does not match. Grid resolution: "
+                    << params->gridDelta << ", Import grid resolution: "
+                    << layer->getGrid().getGridDelta()
+                    << "\nCannot add geometry." << std::endl;
+          continue;
+        }
+        geometry->insertNextLevelSet(layer, false);
+      }
+      break;
+    }
+
+    default:
+      break;
+    }
+  }
+
+  void runProcess() {
+    if (geometry->getLevelSets()->empty()) {
+      std::cout << "Cannot run process on empty geometry." << std::endl;
+      return;
+    }
+
+    std::cout << "\tModel: ";
+    switch (params->processType) {
+    case ProcessType::DEPOSITION: {
+      std::cout << "Single particle deposition\n\tRate: " << params->rate
+                << "\n\tTime: " << params->processTime
+                << "\n\tSticking probability: " << params->sticking
+                << "\n\tCosine exponent: " << params->cosinePower
+                << "\n\tUsing " << params->raysPerPoint
+                << " rays per source grid point\n\n";
+
+      // copy top layer to capture deposition
+      auto topLayerCopy = psSmartPointer<lsDomain<NumericType, D>>::New(
+          geometry->getLevelSets()->back());
+      geometry->insertNextLevelSet(topLayerCopy);
+      runSimpleDeposition(geometry, params);
+      break;
+    }
+
+    case ProcessType::SF6O2ETCHING: {
+      std::cout << "SF6O2 etching\n\tTime: " << params->processTime
+                << "\n\tEtchant flux: " << params->totalEtchantFlux
+                << "\n\tOxygen flux: " << params->totalOxygenFlux
+                << "\n\tIon flux: " << params->totalIonFlux
+                << "\n\tIon energy: " << params->ionEnergy
+                << "\n\tA_O: " << params->A_O << "\n\tUsing "
+                << params->raysPerPoint << " rays per source grid point\n\n";
+      runSF6O2Etching(geometry, params);
+      break;
+    }
+
+    case ProcessType::GEOMETRICUNIFORMDEPOSITION: {
+      std::cout << "Geometric uniform deposition\n\tTime: "
+                << params->processTime << "\n\tRate: " << params->rate
+                << "\n\n";
+      runGeometriyUniformDeposition(geometry, params);
+      break;
+    }
+
+    case ProcessType::DIRECTIONALETCHING: {
+      std::cout << "Directional etching\n\tTime: " << params->processTime
+                << "\n\tDirectional rate: " << params->directionalRate
+                << "\n\tIsotropic rate: " << params->isotropicRate << "\n\n";
+      runDirectionalEtching(geometry, params);
+      break;
+    }
+
+    case ProcessType::NONE:
+      std::cout << "Process model could not be parsed. Skipping line."
+                << std::endl;
+      break;
+
+    default:
+      assert(false);
+      break;
+    }
+  }
+
+  void planarizeGeometry() {
+    psPlanarize<NumericType, D>(geometry, params->maskZPos).apply();
+  }
+
+  void writeVTP() {
+    if (geometry->getLevelSets()->empty()) {
+      std::cout << "Cannot write empty geometry." << std::endl;
+      return;
+    }
+    std::cout << "\tOut file name: " << params->fileName << ".vtp\n\n";
+    geometry->printSurface(params->fileName + ".vtp");
+  }
+
+  std::array<NumericType, 3> getDirection(const std::string &directionString) {
+    std::array<NumericType, 3> direction = {0};
+
+    if (directionString == "negZ") {
+      int i = 2;
+      if constexpr (D == 2)
+        i--;
+      direction[i] = -1.;
+    } else if (directionString == "posZ") {
+      int i = 2;
+      if constexpr (D == 2)
+        i--;
+      direction[i] = 1.;
+    } else if (directionString == "negY") {
+      direction[1] = -1.;
+    } else if (directionString == "posY") {
+      direction[1] = 1.;
+    } else if (directionString == "negX") {
+      direction[0] = -1.;
+    } else if (directionString == "posX") {
+      direction[0] = 1.;
+    } else {
+      std::cout << "Invalid direction: " << directionString << std::endl;
+    }
+
+    return direction;
+  }
+
+  std::string boolString(const int in) { return in == 0 ? "false" : "true"; }
+};
\ No newline at end of file
diff --git a/app/Application2D.cpp b/app/Application2D.cpp
new file mode 100644
index 00000000..e4eb56f5
--- /dev/null
+++ b/app/Application2D.cpp
@@ -0,0 +1,9 @@
+#include "Application.hpp"
+
+int main(int argc, char **argv) {
+
+  Application<2> app(argc, argv);
+  app.run();
+
+  return 0;
+}
\ No newline at end of file
diff --git a/app/Application3D.cpp b/app/Application3D.cpp
new file mode 100644
index 00000000..ac757a30
--- /dev/null
+++ b/app/Application3D.cpp
@@ -0,0 +1,9 @@
+#include "Application.hpp"
+
+int main(int argc, char **argv) {
+
+  Application<3> app(argc, argv);
+  app.run();
+
+  return 0;
+}
\ No newline at end of file
diff --git a/app/ApplicationParameters.hpp b/app/ApplicationParameters.hpp
new file mode 100644
index 00000000..dac5acc1
--- /dev/null
+++ b/app/ApplicationParameters.hpp
@@ -0,0 +1,109 @@
+#pragma once
+
+#include <string>
+#include <unordered_map>
+
+#include <psUtils.hpp>
+
+enum class CommandType { NONE, INIT, GEOMETRY, PROCESS, OUTPUT, PLANARIZE };
+
+enum class GeometryType { NONE, TRENCH, HOLE, PLANE, GDS, IMPORT };
+
+enum class ProcessType {
+  NONE,
+  SF6O2ETCHING,
+  DEPOSITION,
+  GEOMETRICUNIFORMDEPOSITION,
+  DIRECTIONALETCHING
+};
+
+#ifdef VIENNAPS_USE_DOUBLE
+using NumericType = double;
+#else
+using NumericType = float;
+#endif
+
+struct ApplicationParameters {
+  int printIntermediate = 0;
+  GeometryType geometryType = GeometryType::NONE;
+  ProcessType processType = ProcessType::NONE;
+
+  // Domain
+  NumericType gridDelta = 0.02;
+  NumericType xExtent = 1.0;
+  NumericType yExtent = 1.0;
+  int periodicBoundary = 0;
+
+  // Geometry
+  int mask = 0;
+  int maskId = 0;
+  NumericType taperAngle = 0.;
+  NumericType maskZPos = 0.;
+  // MAKE Trench
+  NumericType trenchWidth = 0.2;
+  NumericType trenchHeight = 0.5;
+  // MAKE hole
+  NumericType holeDepth = 0.2;
+  NumericType holeRadius = 0.2;
+  // GDS / IMPORT
+  int layers = 0;
+  std::string fileName = "";
+  NumericType maskHeight = 0.1;
+
+  // Process
+  NumericType processTime = 1;
+  int raysPerPoint = 3000;
+  // SF6O2Etching
+  NumericType totalEtchantFlux = 4.5e16;
+  NumericType totalOxygenFlux = 1e18;
+  NumericType totalIonFlux = 2e16;
+  NumericType ionEnergy = 100;
+  NumericType A_O = 3.;
+  // Deposition
+  NumericType rate = 1.;
+  NumericType sticking = 1.;
+  NumericType cosinePower = 1.;
+  // Directional etching
+  NumericType directionalRate = 1.;
+  NumericType isotropicRate = 0.;
+  std::string direction = "negZ";
+
+  ApplicationParameters() {}
+
+  void defaultParameters(bool all = false) {
+    mask = 0;
+    maskId = 0;
+    taperAngle = 0.;
+    maskZPos = 0.;
+    trenchWidth = 0.2;
+    trenchHeight = 0.5;
+    holeDepth = 0.2;
+    holeRadius = 0.2;
+    layers = 0;
+    fileName = "";
+    maskHeight = 0.1;
+    processTime = 1;
+    raysPerPoint = 3000;
+    totalEtchantFlux = 4.5e16;
+    totalOxygenFlux = 1e18;
+    totalIonFlux = 2e16;
+    ionEnergy = 100;
+    A_O = 3.;
+    rate = 1.;
+    sticking = 1.;
+    cosinePower = 1.;
+    directionalRate = 1.;
+    isotropicRate = 0.;
+    direction = "negZ";
+
+    if (all) {
+      printIntermediate = 0;
+      geometryType = GeometryType::NONE;
+      processType = ProcessType::NONE;
+      gridDelta = 0.02;
+      xExtent = 1.0;
+      yExtent = 1.0;
+      periodicBoundary = 0;
+    }
+  }
+};
diff --git a/app/ApplicationParser.hpp b/app/ApplicationParser.hpp
new file mode 100644
index 00000000..b97d079c
--- /dev/null
+++ b/app/ApplicationParser.hpp
@@ -0,0 +1,183 @@
+#pragma once
+
+#include <sstream>
+
+#include <psSmartPointer.hpp>
+
+#include <ApplicationParameters.hpp>
+
+class ApplicationParser {
+private:
+  psSmartPointer<ApplicationParameters> params;
+
+public:
+  ApplicationParser() {}
+
+  ApplicationParser(
+      const psSmartPointer<ApplicationParameters> passedParameters)
+      : params(passedParameters) {}
+
+  void
+  setParameters(const psSmartPointer<ApplicationParameters> passedParameters) {
+    params = passedParameters;
+  }
+
+  CommandType parseCommand(std::istringstream &stream, const int lineNumber) {
+    std::string command;
+    stream >> command;
+    if (command == "INIT") {
+      std::cout << "Initializing ...\n";
+      parseInit(stream);
+      return CommandType::INIT;
+    } else if (command == "GEOMETRY") {
+      std::cout << "Creating geometry ..." << std::endl;
+      parseGeometry(stream);
+      return CommandType::GEOMETRY;
+    } else if (command == "PROCESS") {
+      std::cout << "Starting process ..." << std::endl;
+      parseProcess(stream);
+      return CommandType::PROCESS;
+    } else if (command == "OUTPUT") {
+      std::cout << "Writing geometry to file ..." << std::endl;
+      parseOutput(stream);
+      return CommandType::OUTPUT;
+    } else if (command == "PLANARIZE") {
+      std::cout << "Planarizing ..." << std::endl;
+      parsePlanarize(stream);
+      return CommandType::PLANARIZE;
+    } else {
+      std::cout << "Unknown command in config file. Skipping line "
+                << lineNumber << std::endl;
+      return CommandType::NONE;
+    }
+  }
+
+private:
+  void parseInit(std::istringstream &stream) {
+    auto config = parseLineStream(stream);
+    psUtils::AssignItems(
+        config, psUtils::Item{"xExtent", params->xExtent},
+        psUtils::Item{"yExtent", params->yExtent},
+        psUtils::Item{"resolution", params->gridDelta},
+        psUtils::Item{"printIntermediate", params->printIntermediate},
+        psUtils::Item{"periodic", params->periodicBoundary});
+  }
+
+  void parseGeometry(std::istringstream &stream) {
+    std::string type;
+    stream >> type;
+    auto config = parseLineStream(stream);
+    if (type == "Trench") {
+      params->geometryType = GeometryType::TRENCH;
+      psUtils::AssignItems(config, psUtils::Item{"width", params->trenchWidth},
+                           psUtils::Item{"depth", params->trenchHeight},
+                           psUtils::Item{"zPos", params->maskZPos},
+                           psUtils::Item{"tapering", params->taperAngle},
+                           psUtils::Item{"mask", params->mask});
+    } else if (type == "Hole") {
+      params->geometryType = GeometryType::HOLE;
+      psUtils::AssignItems(config, psUtils::Item{"radius", params->holeRadius},
+                           psUtils::Item{"depth", params->holeDepth},
+                           psUtils::Item{"zPos", params->maskZPos},
+                           psUtils::Item{"tapering", params->taperAngle},
+                           psUtils::Item{"mask", params->mask});
+    } else if (type == "Plane") {
+      params->geometryType = GeometryType::PLANE;
+      psUtils::AssignItems(config, psUtils::Item{"zPos", params->maskZPos});
+    } else if (type == "GDS") {
+      params->geometryType = GeometryType::GDS;
+      psUtils::AssignItems(config, psUtils::Item{"file", params->fileName},
+                           psUtils::Item{"layer", params->layers},
+                           psUtils::Item{"zPos", params->maskZPos},
+                           psUtils::Item{"maskHeight", params->maskHeight});
+    } else if (type == "Import") {
+      params->geometryType = GeometryType::IMPORT;
+      psUtils::AssignItems(config, psUtils::Item{"file", params->fileName},
+                           psUtils::Item{"layers", params->layers});
+    } else {
+      params->geometryType = GeometryType::NONE;
+      std::cout << "Invalid mask type." << std::endl;
+      exit(1);
+    }
+  }
+
+  void parseProcess(std::istringstream &stream) {
+    std::string model;
+    stream >> model;
+    auto config = parseLineStream(stream);
+    if (model == "Deposition") {
+      params->processType = ProcessType::DEPOSITION;
+      psUtils::AssignItems(config, psUtils::Item{"rate", params->rate},
+                           psUtils::Item{"time", params->processTime},
+                           psUtils::Item{"sticking", params->sticking},
+                           psUtils::Item{"cosineExponent", params->cosinePower},
+                           psUtils::Item{"raysPerPoint", params->raysPerPoint});
+    } else if (model == "SF6O2Etching") {
+      params->processType = ProcessType::SF6O2ETCHING;
+      psUtils::AssignItems(
+          config, psUtils::Item{"time", params->processTime},
+          psUtils::Item{"ionFlux", params->totalIonFlux},
+          psUtils::Item{"ionEnergy", params->ionEnergy},
+          psUtils::Item{"etchantFlux", params->totalEtchantFlux},
+          psUtils::Item{"oxygenFlux", params->totalOxygenFlux},
+          psUtils::Item{"A_O", params->A_O},
+          psUtils::Item{"raysPerPoint", params->raysPerPoint});
+    } else if (model == "GeometricUniformDeposition") {
+      params->processType = ProcessType::GEOMETRICUNIFORMDEPOSITION;
+      psUtils::AssignItems(config, psUtils::Item{"rate", params->rate},
+                           psUtils::Item{"time", params->processTime});
+    } else if (model == "DirectionalEtching") {
+      params->processType = ProcessType::DIRECTIONALETCHING;
+      psUtils::AssignItems(
+          config, psUtils::Item{"direction", params->direction},
+          psUtils::Item{"directionalRate", params->directionalRate},
+          psUtils::Item{"isotropicRate", params->isotropicRate},
+          psUtils::Item{"time", params->processTime});
+    } else {
+      params->processType = ProcessType::NONE;
+      std::cout << "Invalid process model: " << model << std::endl;
+    }
+  }
+
+  void parsePlanarize(std::istringstream &stream) {
+    auto config = parseLineStream(stream);
+    psUtils::AssignItems(config, psUtils::Item{"height", params->maskZPos});
+  }
+
+  void parseOutput(std::istringstream &stream) { stream >> params->fileName; }
+
+  std::unordered_map<std::string, std::string>
+  parseLineStream(std::istringstream &input) {
+    // Regex to find trailing and leading whitespaces
+    const auto wsRegex = std::regex("^ +| +$|( ) +");
+
+    // Regular expression for extracting key and value separated by '=' as two
+    // separate capture groups
+    const auto keyValueRegex = std::regex(
+        R"rgx([ \t]*([0-9a-zA-Z_\-\.+]+)[ \t]*=[ \t]*([0-9a-zA-Z_\-\.+]+).*$)rgx");
+
+    // Reads a simple config file containing a single <key>=<value> pair per
+    // line and returns the content as an unordered map
+    std::unordered_map<std::string, std::string> paramMap;
+    std::string expression;
+    while (input >> expression) {
+      // Remove trailing and leading whitespaces
+      expression = std::regex_replace(expression, wsRegex, "$1");
+      // Skip this expression if it is marked as a comment
+      if (expression.rfind('#') == 0 || expression.empty())
+        continue;
+
+      // Extract key and value
+      std::smatch smatch;
+      if (std::regex_search(expression, smatch, keyValueRegex)) {
+        if (smatch.size() < 3) {
+          std::cerr << "Malformed expression:\n " << expression;
+          continue;
+        }
+
+        paramMap.insert({smatch[1], smatch[2]});
+      }
+    }
+    return paramMap;
+  }
+};
diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt
new file mode 100644
index 00000000..43f1adf7
--- /dev/null
+++ b/app/CMakeLists.txt
@@ -0,0 +1,48 @@
+cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
+
+project(ViennaPSApplication)
+
+find_package(ViennaPS CONFIG PATHS ${ViennaPS_BINARY_DIR})
+if(NOT ViennaPS_FOUND)
+  message(
+    FATAL_ERROR
+      "One of the dependencies hasn't been built yet and thus can't be used by the buildExamples target.\n\
+  Try the following steps to fix the problem:\n\
+  1. Turn off the buildExamples target by disabling the VIENNAPS_BUILD_EXAMPLES option when configuring the project\n\
+  (e.g. with the `-D VIENNAPS_BUILD_EXAMPLES=OFF` flag)\n\
+  2. Build the dependencies target: e.g. `make buildDependencies` when using Makefiles\n\
+  3. Turn on the buildExamples target again by re-enabling the VIENNAPS_BUILD_EXAMPLES option\n\
+  (e.g. with the `-D VIENNAPS_BUILD_EXAMPLES=ON` flag)\n\
+  4. Now the configure step should pick up the built dependencies and proceed as intended.")
+endif()
+
+list(PREPEND VIENNAPS_INCLUDE_DIRS ${VIENNAPS_BUILD_INCLUDE_DIRS})
+
+# Add subfolders to the include directories, since during installation the directory structure is
+# flattened, but in the source tree - which we are using here - it is not.
+list(PREPEND VIENNAPS_INCLUDE_DIRS ${VIENNAPS_BUILD_INCLUDE_DIRS}/CellSet
+     ${VIENNAPS_BUILD_INCLUDE_DIRS}/Geometries ${VIENNAPS_BUILD_INCLUDE_DIRS}/Models
+     ${VIENNAPS_BUILD_INCLUDE_DIRS}/Compact)
+
+set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
+
+option(VIENNAPS_USE_DOUBLE "Use double precision floating point numbers in application." OFF)
+if(VIENNAPS_USE_DOUBLE)
+  add_compile_definitions(VIENNAPS_USE_DOUBLE)
+endif(VIENNAPS_USE_DOUBLE)
+
+set_property(DIRECTORY PROPERTY EXCLUDE_FROM_ALL TRUE)
+
+add_custom_target(buildApplication)
+
+add_executable(ViennaPS2D ${PROJECT_SOURCE_DIR}/Application2D.cpp)
+target_include_directories(ViennaPS2D PUBLIC ${VIENNAPS_INCLUDE_DIRS} ${PROJECT_SOURCE_DIR})
+target_link_libraries(ViennaPS2D ${VIENNAPS_LIBRARIES})
+add_dependencies(buildApplication ViennaPS2D)
+
+add_executable(ViennaPS3D ${PROJECT_SOURCE_DIR}/Application3D.cpp)
+target_include_directories(ViennaPS3D PUBLIC ${VIENNAPS_INCLUDE_DIRS} ${PROJECT_SOURCE_DIR})
+target_link_libraries(ViennaPS3D ${VIENNAPS_LIBRARIES})
+add_dependencies(buildApplication ViennaPS3D)
+
+configure_file(SampleConfig.txt ${CMAKE_BINARY_DIR}/SampleConfig.txt COPYONLY)
diff --git a/app/README.md b/app/README.md
new file mode 100644
index 00000000..5194c936
--- /dev/null
+++ b/app/README.md
@@ -0,0 +1,156 @@
+# Configuration File Documentation
+
+The configuration file must obey a certain structure in order to be parsed correctly. An example for a configuration file can be seen in **SampleConfig.txt**. The configuration file is parsed line by line and each succesfully parsed line is executed immediately.
+
+## Commands
+Every line has to start with a command statement, followed by a list of parameters for the command. Possible commands are:
+
+- **INIT**: Initialize the simulation domain.
+- **GEOMETRY** _\<GeometryType>_: Create or import a geometry. Possible options for creating geometries are: _Trench_, _Hole_ and _Plane_. It is also possible to import geometries from _.lvst_ files by specifying _Import_ or to import layers from GDSII file format by secifying _GDS_ (only possible in 3D mode). Parameters for the geometries are described below.
+- **PROCESS** _\<ProcessType>_: Run a process. Possible process types are: _Deposition_, _GeometricUniformDeposition_, _SF6O2Etching_ and _DirectionalEtching_. Parameters for the processes are described below.
+- **PLANARIZE**: Planarize the geometry at a given height.
+- **OUTPUT** _\<fileName>_: Print the surface of the geometry in _fileName.vtp_ file format.
+
+## Parameters
+
+Every parameter is specified by _paramName=paramValue_ (there must not be any blank space before or after the _=_ sign) after the command statement.
+All parameters which are parsed additional to a command are described below. For any parameter which is not specified, a default value is assigned.
+
+---
+**INIT**
+<dl>
+  <dt>xExtent</dt>
+  <dd>width of the domain in x-direction (numeric value, default: 1)</dd>
+  <dt>yExtent</dt>
+  <dd>width of domain in y-direction (this is ignored in 2D mode, numeric value, default: 1)</dd>
+  <dt>resolution</dt>
+  <dd>distance between grid points in the domain (numeric value, default: 0.02)</dd>
+  <dt>printIntermediate</dt>
+  <dd>output intermediate disk meshes for each process step (boolean, default: 0)</dd>
+  <dt>periodic</dt>
+  <dd>use periodic boundary conditions (boolean, default: 0)</dd>
+</dl>
+
+---
+**GEOMETRY Trench**
+<dl>
+  <dt>width</dt>
+  <dd> width of the trench (numeric value, default: 0.2)</dd>
+  <dt>depth</dt>
+  <dd>depth of the trench (numeric value, default: 0.5)</dd>
+  <dt>zPos</dt>
+  <dd>offset of the trench in z-direction (y-direction in 2D mode, numeric value, default: 0)</dd>
+  <dt>tapering</dt>
+  <dd>tapering angle of the trench in degrees (numeric value, default: 0)</dd>
+  <dt>mask</dt>
+  <dd>creates a mask from the trench, such that the bottom represents a different material and a selective etching or deposition process can be used (boolean, default: 0)
+</dl>
+
+---
+**GEOMETRY Hole**
+<dl>
+  <dt>radius</dt>
+  <dd> radius of the hole (numeric value, default: 0.2)</dd>
+  <dt>depth</dt>
+  <dd>depth of the hole (numeric value, default: 0.2)</dd>
+  <dt>zPos</dt>
+  <dd>offset of the hole in z-direction (y-direction in 2D mode, numeric value, default: 0)</dd>
+  <dt>tapering</dt>
+  <dd>tapering angle of the hole in degrees (numeric value, default: 0)</dd>
+  <dt>mask</dt>
+  <dd>creates a mask from the hole, such that the bottom represents a different material and a selective etching or deposition process can be used (boolean, default: 0)
+</dl>
+
+---
+**GEOMETRY Plane**
+<dl>
+  <dt>zPos</dt>
+  <dd>offset of the plane in z-direction (y-direction in 2D mode, numeric value, default: 0)</dd>
+</dl>
+
+---
+**GEOMETRY Import**
+<dl>
+  <dt>file</dt>
+  <dd>file name of ViennaLS geometry files. The file name is assumed to be appended with "_layer(i).lvst" (string, no default value) </dd>
+  <dt>layers</dt>
+  <dd>number of layers to read (integer value, default: 0)</dd>
+</dl>
+
+---
+**GEOMETRY GDS**
+<dl>
+  <dt>file</dt>
+  <dd>file name of GDSII file (string, no default value) </dd>
+  <dt>layer</dt>
+  <dd>specify which layer of the file should be read (integer value, default: 0)</dd>
+  <dt>maskHeight</dt>
+  <dd>height of the layer in z-direction (numeric value, default: 0.1)</dd>
+  <dt>zPos</dt>
+  <dd>offset of the layer in z-direction (numeric value, default: 0)</dd>
+</dl>
+
+---
+**PROCESS Deposition**
+<dl>
+  <dt>time</dt>
+  <dd>process time (numeric value, default: 1)</dd>
+  <dt>rate</dt>
+  <dd>deposition rate (numeric value, default: 1)</dd>
+  <dt>sticking</dt>
+  <dd>sticking coefficient (numeric value, default: 1)</dd>
+  <dt>cosineExponent</dt>
+  <dd>Exponent for cosine distribution of initial ray directions (numeric value, default: 1)</dd>
+  <dt>raysPerPoint</dt>
+  <dd>number of rays traced per grid point in the surface geometry (integer value, default: 3000)</dd>
+</dl>
+
+---
+**PROCESS GeometricUniformDeposition**
+<dl>
+  <dt>time</dt>
+  <dd>process time (numeric value, default: 1)</dd>
+  <dt>rate</dt>
+  <dd>deposition rate (numeric value, default: 1)</dd>
+</dl>
+
+---
+**PROCESS SF6O2Etching**
+<dl>
+  <dt>time</dt>
+  <dd>process time (numeric value, default: 1)</dd>
+  <dt>ionFlux</dt>
+  <dd>total flux of ions in plasma (numeric value, default: 2e16)</dd>
+  <dt>ionEnergy</dt>
+  <dd>mean ion energy (numeric value, default: 100)</dd>
+  <dt>etchantFlux</dt>
+  <dd>total flux of etchant species in plasma (numeric value, default: 4.5e16)</dd>
+  <dt>oxygenFlux</dt>
+  <dd>total flux of oxygen in plasma (numeric value, default: 1e18)</dd>
+  <dt>A_O</dt>
+  <dd>factor for ion etching yield on oxygen (numeric value, default: 3)</dd>
+  <dt>raysPerPoint</dt>
+  <dd>number of rays traced per grid point in the surface geometry (integer value, default: 3000)</dd>
+</dl>
+
+---
+**PROCESS DirectionalEtching**
+<dl>
+  <dt>time</dt>
+  <dd>process time (numeric value, default: 1)</dd>
+  <dt>direction</dt>
+  <dd>primal direction of etching (string, default: negZ, negY in 2D mode) 
+  <dt>directionalRate</dt>
+  <dd>etching rate in primal direction (numeric value, default: 1)</dd>
+  <dt>isotropicRate</dt>
+  <dd>isotropic etching rate (numeric value, default: 0)</dd>
+</dl>
+
+---
+**PLANARIZE**
+<dl>
+  <dt>height</dt>
+  <dd>height of planarization on the z-axis (y-axis in 2D, numeric value, default: 0)</dd>
+</dl>
+
+---
diff --git a/app/SampleConfig.txt b/app/SampleConfig.txt
new file mode 100644
index 00000000..d4be7b89
--- /dev/null
+++ b/app/SampleConfig.txt
@@ -0,0 +1,5 @@
+INIT xExtent=1 yExtent=1 resolution=0.01 printIntermediate=0 periodic=0
+GEOMETRY Trench width=0.4 depth=0.5 zPos=0.5 tapering=10
+OUTPUT initial
+PROCESS GeometricUniformDeposition rate=0.1 time=1
+OUTPUT result
\ No newline at end of file
diff --git a/include/Geometries/psMakeFin.hpp b/include/Geometries/psMakeFin.hpp
index c8947b73..096db440 100644
--- a/include/Geometries/psMakeFin.hpp
+++ b/include/Geometries/psMakeFin.hpp
@@ -15,9 +15,11 @@ template <class NumericType, int D> class psMakeFin {
   NumericType gridDelta = .25;
   NumericType xExtent = 10;
   NumericType yExtent = 7;
+  NumericType baseHeight = 0.;
 
   NumericType finWidth = 7;
   NumericType finHeight = 17.5;
+  bool periodicBoundary = false;
   bool makeMask = true;
 
   psMakeFin(PSPtrType passedDomain) : domain(passedDomain) {}
@@ -25,10 +27,12 @@ template <class NumericType, int D> class psMakeFin {
   psMakeFin(PSPtrType passedDomain, const NumericType passedGridDelta,
             const NumericType passedXExtent, const NumericType passedYExtent,
             const NumericType passedFinWidth, const NumericType passedFinHeight,
-            const bool passedMakeMask = true)
+            const NumericType passedBaseHeight = 0.,
+            const bool passedPeriodic = false, const bool passedMakeMask = true)
       : domain(passedDomain), gridDelta(passedGridDelta),
         xExtent(passedXExtent), yExtent(passedYExtent),
         finWidth(passedFinWidth), finHeight(passedFinHeight),
+        baseHeight(passedBaseHeight), periodicBoundary(passedPeriodic),
         makeMask(passedMakeMask) {}
 
   void apply() {
diff --git a/include/Geometries/psMakeHole.hpp b/include/Geometries/psMakeHole.hpp
index f5591575..aa328b01 100644
--- a/include/Geometries/psMakeHole.hpp
+++ b/include/Geometries/psMakeHole.hpp
@@ -21,11 +21,13 @@ template <class NumericType, int D> class psMakeHole {
   NumericType gridDelta = .25;
   NumericType xExtent = 10;
   NumericType yExtent = 10;
+  NumericType baseHeight = 0.;
 
   NumericType holeRadius = 3;
   NumericType taperAngle = 0; // tapering angle in degrees
   NumericType holeDepth = 6;
   bool makeMask = true;
+  bool periodicBoundary = false;
 
   psMakeHole(PSPtrType passedDomain) : domain(passedDomain) {}
 
@@ -33,22 +35,15 @@ template <class NumericType, int D> class psMakeHole {
              const NumericType passedXExtent, const NumericType passedYExtent,
              const NumericType passedHoleRadius,
              const NumericType passedHoleDepth,
+             const NumericType passedTaperAngle = 0.,
+             const NumericType passedBaseHeight = 0.,
+             const bool passedPeriodicBoundary = false,
              const bool passedMakeMask = true)
       : domain(passedDomain), gridDelta(passedGridDelta),
         xExtent(passedXExtent), yExtent(passedYExtent),
         holeRadius(passedHoleRadius), holeDepth(passedHoleDepth),
-        makeMask(passedMakeMask) {}
-
-  psMakeHole(PSPtrType passedDomain, const NumericType passedGridDelta,
-             const NumericType passedXExtent, const NumericType passedYExtent,
-             const NumericType passedHoleRadius,
-             const NumericType passedHoleDepth,
-             const NumericType passedTaperAngle,
-             const bool passedMakeMask = true)
-      : domain(passedDomain), gridDelta(passedGridDelta),
-        xExtent(passedXExtent), yExtent(passedYExtent),
-        holeRadius(passedHoleRadius), taperAngle(passedTaperAngle),
-        holeDepth(passedHoleDepth), makeMask(passedMakeMask) {}
+        taperAngle(passedTaperAngle), baseHeight(passedBaseHeight),
+        periodicBoundary(passedPeriodicBoundary), makeMask(passedMakeMask) {}
 
   void apply() {
     if (D != 3) {
@@ -71,18 +66,24 @@ template <class NumericType, int D> class psMakeHole {
     if constexpr (D == 3) {
       bounds[2] = -yExtent / 2.;
       bounds[3] = yExtent / 2.;
-      bounds[4] = -gridDelta;
-      bounds[5] = holeDepth + gridDelta;
+      bounds[4] = baseHeight - gridDelta;
+      bounds[5] = baseHeight + holeDepth + gridDelta;
     } else {
-      bounds[2] = -gridDelta;
-      bounds[3] = holeDepth + gridDelta;
+      bounds[2] = baseHeight - gridDelta;
+      bounds[3] = baseHeight + holeDepth + gridDelta;
     }
 
     typename lsDomain<NumericType, D>::BoundaryType boundaryCons[D];
 
-    for (int i = 0; i < D - 1; i++)
-      boundaryCons[i] =
-          lsDomain<NumericType, D>::BoundaryType::REFLECTIVE_BOUNDARY;
+    for (int i = 0; i < D - 1; i++) {
+      if (periodicBoundary) {
+        boundaryCons[i] =
+            lsDomain<NumericType, D>::BoundaryType::PERIODIC_BOUNDARY;
+      } else {
+        boundaryCons[i] =
+            lsDomain<NumericType, D>::BoundaryType::REFLECTIVE_BOUNDARY;
+      }
+    }
     boundaryCons[D - 1] =
         lsDomain<NumericType, D>::BoundaryType::INFINITE_BOUNDARY;
 
@@ -91,19 +92,20 @@ template <class NumericType, int D> class psMakeHole {
     NumericType normal[D] = {0.};
     NumericType origin[D] = {0.};
     normal[D - 1] = 1.;
+    origin[D - 1] = baseHeight;
     lsMakeGeometry<NumericType, D>(
         substrate, lsSmartPointer<lsPlane<NumericType, D>>::New(origin, normal))
         .apply();
 
     // mask layer
     auto mask = LSPtrType::New(bounds, boundaryCons, gridDelta);
-    origin[D - 1] = holeDepth;
+    origin[D - 1] = holeDepth + baseHeight;
     lsMakeGeometry<NumericType, D>(
         mask, lsSmartPointer<lsPlane<NumericType, D>>::New(origin, normal))
         .apply();
 
     auto maskAdd = LSPtrType::New(bounds, boundaryCons, gridDelta);
-    origin[D - 1] = 0.;
+    origin[D - 1] = baseHeight;
     normal[D - 1] = -1.;
     lsMakeGeometry<NumericType, D>(
         maskAdd, lsSmartPointer<lsPlane<NumericType, D>>::New(origin, normal))
@@ -115,7 +117,7 @@ template <class NumericType, int D> class psMakeHole {
 
     // cylinder cutout
     normal[D - 1] = 1.;
-    origin[D - 1] = 0.;
+    origin[D - 1] = baseHeight;
 
     NumericType topRadius = holeRadius;
     if (taperAngle) {
diff --git a/include/Geometries/psMakePlane.hpp b/include/Geometries/psMakePlane.hpp
new file mode 100644
index 00000000..40994ac6
--- /dev/null
+++ b/include/Geometries/psMakePlane.hpp
@@ -0,0 +1,74 @@
+#pragma once
+
+#include <lsBooleanOperation.hpp>
+#include <lsDomain.hpp>
+#include <lsFromSurfaceMesh.hpp>
+#include <lsMakeGeometry.hpp>
+
+#include <psDomain.hpp>
+
+/**
+  Creates a plane in z(3D)/y(2D) direction.
+*/
+template <class NumericType, int D> class psMakePlane {
+  using LSPtrType = psSmartPointer<lsDomain<NumericType, D>>;
+  using PSPtrType = psSmartPointer<psDomain<NumericType, D>>;
+
+public:
+  PSPtrType domain = nullptr;
+
+  NumericType gridDelta = .25;
+  NumericType xExtent = 10;
+  NumericType yExtent = 10;
+  NumericType height = 0.;
+  bool periodicBoundary = false;
+
+  psMakePlane(PSPtrType passedDomain) : domain(passedDomain) {}
+
+  psMakePlane(PSPtrType passedDomain, const NumericType passedGridDelta,
+              const NumericType passedXExtent, const NumericType passedYExtent,
+              const NumericType passedHeight, const bool passedPeriodic = false)
+      : domain(passedDomain), gridDelta(passedGridDelta),
+        xExtent(passedXExtent), yExtent(passedYExtent), height(passedHeight),
+        periodicBoundary(passedPeriodic) {}
+
+  void apply() {
+    domain->clear();
+    double bounds[2 * D];
+    bounds[0] = -xExtent / 2.;
+    bounds[1] = xExtent / 2.;
+
+    if constexpr (D == 3) {
+      bounds[2] = -yExtent / 2.;
+      bounds[3] = yExtent / 2.;
+      bounds[4] = -gridDelta;
+      bounds[5] = gridDelta;
+    } else {
+      bounds[2] = -gridDelta;
+      bounds[3] = gridDelta;
+    }
+
+    typename lsDomain<NumericType, D>::BoundaryType boundaryCons[D];
+
+    for (int i = 0; i < D - 1; i++) {
+      if (periodicBoundary) {
+        boundaryCons[i] =
+            lsDomain<NumericType, D>::BoundaryType::PERIODIC_BOUNDARY;
+      } else {
+        boundaryCons[i] =
+            lsDomain<NumericType, D>::BoundaryType::REFLECTIVE_BOUNDARY;
+      }
+    }
+    boundaryCons[D - 1] =
+        lsDomain<NumericType, D>::BoundaryType::INFINITE_BOUNDARY;
+
+    auto substrate = LSPtrType::New(bounds, boundaryCons, gridDelta);
+    NumericType normal[D] = {0.};
+    NumericType origin[D] = {0.};
+    normal[D - 1] = 1.;
+    origin[D - 1] = height;
+    lsMakeGeometry<NumericType, D>(
+        substrate, lsSmartPointer<lsPlane<NumericType, D>>::New(origin, normal))
+        .apply();
+  }
+};
\ No newline at end of file
diff --git a/include/Geometries/psMakeTrench.hpp b/include/Geometries/psMakeTrench.hpp
index 3730d346..f1933133 100644
--- a/include/Geometries/psMakeTrench.hpp
+++ b/include/Geometries/psMakeTrench.hpp
@@ -20,34 +20,27 @@ template <class NumericType, int D> class psMakeTrench {
   NumericType gridDelta = .25;
   NumericType xExtent = 20;
   NumericType yExtent = 14;
+  NumericType baseHeight = 0.;
 
   NumericType trenchWidth = 7;
   NumericType taperingAngle = 0;
   NumericType trenchDepth = 17.5;
-  bool makeMask = true;
-
-  psMakeTrench(PSPtrType passedDomain) : domain(passedDomain) {}
+  bool periodicBoundary = false;
+  bool makeMask = false;
 
   psMakeTrench(PSPtrType passedDomain, const NumericType passedGridDelta,
                const NumericType passedXExtent, const NumericType passedYExtent,
                const NumericType passedTrenchWidth,
                const NumericType passedTrenchHeight,
-               const bool passedMakeMask = true)
+               const NumericType passedTaperingAngle = 0.,
+               const NumericType passedBaseHeight = 0.,
+               const bool passedPeriodicBoundary = false,
+               const bool passedMakeMask = false)
       : domain(passedDomain), gridDelta(passedGridDelta),
         xExtent(passedXExtent), yExtent(passedYExtent),
         trenchWidth(passedTrenchWidth), trenchDepth(passedTrenchHeight),
-        makeMask(passedMakeMask) {}
-
-  psMakeTrench(PSPtrType passedDomain, const NumericType passedGridDelta,
-               const NumericType passedXExtent, const NumericType passedYExtent,
-               const NumericType passedTrenchWidth,
-               const NumericType passedTrenchHeight,
-               const NumericType passedTaperingAngle,
-               const bool passedMakeMask = true)
-      : domain(passedDomain), gridDelta(passedGridDelta),
-        xExtent(passedXExtent), yExtent(passedYExtent),
-        trenchWidth(passedTrenchWidth), taperingAngle(passedTaperingAngle),
-        trenchDepth(passedTrenchHeight), makeMask(passedMakeMask) {}
+        taperingAngle(passedTaperingAngle), baseHeight(passedBaseHeight),
+        periodicBoundary(passedPeriodicBoundary), makeMask(passedMakeMask) {}
 
   void apply() {
     domain->clear();
@@ -67,9 +60,15 @@ template <class NumericType, int D> class psMakeTrench {
 
     typename lsDomain<NumericType, D>::BoundaryType boundaryCons[D];
 
-    for (int i = 0; i < D - 1; i++)
-      boundaryCons[i] =
-          lsDomain<NumericType, D>::BoundaryType::REFLECTIVE_BOUNDARY;
+    for (int i = 0; i < D - 1; i++) {
+      if (periodicBoundary) {
+        boundaryCons[i] =
+            lsDomain<NumericType, D>::BoundaryType::PERIODIC_BOUNDARY;
+      } else {
+        boundaryCons[i] =
+            lsDomain<NumericType, D>::BoundaryType::REFLECTIVE_BOUNDARY;
+      }
+    }
     boundaryCons[D - 1] =
         lsDomain<NumericType, D>::BoundaryType::INFINITE_BOUNDARY;
 
@@ -77,18 +76,19 @@ template <class NumericType, int D> class psMakeTrench {
     NumericType normal[D] = {0.};
     NumericType origin[D] = {0.};
     normal[D - 1] = 1.;
+    origin[D - 1] = baseHeight;
     lsMakeGeometry<NumericType, D>(
         substrate, lsSmartPointer<lsPlane<NumericType, D>>::New(origin, normal))
         .apply();
 
     auto mask = LSPtrType::New(bounds, boundaryCons, gridDelta);
-    origin[D - 1] = trenchDepth;
+    origin[D - 1] = trenchDepth + baseHeight;
     lsMakeGeometry<NumericType, D>(
         mask, lsSmartPointer<lsPlane<NumericType, D>>::New(origin, normal))
         .apply();
 
     auto maskAdd = LSPtrType::New(bounds, boundaryCons, gridDelta);
-    origin[D - 1] = 0.;
+    origin[D - 1] = baseHeight;
     normal[D - 1] = -1.;
     lsMakeGeometry<NumericType, D>(
         maskAdd, lsSmartPointer<lsPlane<NumericType, D>>::New(origin, normal))
@@ -112,9 +112,12 @@ template <class NumericType, int D> class psMakeTrench {
         mesh->nodes[0][0] = -trenchWidth / 2.;
         mesh->nodes[1][0] = trenchWidth / 2.;
         mesh->nodes[2][0] = trenchWidth / 2. + offset;
-        mesh->nodes[2][1] = trenchDepth;
         mesh->nodes[3][0] = -trenchWidth / 2. - offset;
-        mesh->nodes[3][1] = trenchDepth;
+
+        mesh->nodes[0][1] = baseHeight;
+        mesh->nodes[1][1] = baseHeight;
+        mesh->nodes[2][1] = trenchDepth + baseHeight;
+        mesh->nodes[3][1] = trenchDepth + baseHeight;
 
         mesh->insertNextLine(std::array<unsigned, 2>{0, 3});
         mesh->insertNextLine(std::array<unsigned, 2>{3, 2});
@@ -127,32 +130,36 @@ template <class NumericType, int D> class psMakeTrench {
           mesh->insertNextNode(node);
         }
         mesh->nodes[0][0] = -trenchWidth / 2.;
-        mesh->nodes[0][1] = -yExtent / 2.;
+        mesh->nodes[0][1] = -yExtent / 2. - gridDelta;
+        mesh->nodes[0][2] = baseHeight;
 
         mesh->nodes[1][0] = trenchWidth / 2.;
-        mesh->nodes[1][1] = -yExtent / 2.;
+        mesh->nodes[1][1] = -yExtent / 2. - gridDelta;
+        mesh->nodes[1][2] = baseHeight;
 
         mesh->nodes[2][0] = trenchWidth / 2.;
-        mesh->nodes[2][1] = yExtent / 2.;
+        mesh->nodes[2][1] = yExtent / 2. + gridDelta;
+        mesh->nodes[2][2] = baseHeight;
 
         mesh->nodes[3][0] = -trenchWidth / 2.;
-        mesh->nodes[3][1] = yExtent / 2.;
+        mesh->nodes[3][1] = yExtent / 2. + gridDelta;
+        mesh->nodes[3][2] = baseHeight;
 
         mesh->nodes[4][0] = -trenchWidth / 2. - offset;
-        mesh->nodes[4][1] = -yExtent / 2.;
-        mesh->nodes[4][2] = trenchDepth;
+        mesh->nodes[4][1] = -yExtent / 2. - gridDelta;
+        mesh->nodes[4][2] = trenchDepth + baseHeight;
 
         mesh->nodes[5][0] = trenchWidth / 2. + offset;
-        mesh->nodes[5][1] = -yExtent / 2.;
-        mesh->nodes[5][2] = trenchDepth;
+        mesh->nodes[5][1] = -yExtent / 2. - gridDelta;
+        mesh->nodes[5][2] = trenchDepth + baseHeight;
 
         mesh->nodes[6][0] = trenchWidth / 2. + offset;
-        mesh->nodes[6][1] = yExtent / 2.;
-        mesh->nodes[6][2] = trenchDepth;
+        mesh->nodes[6][1] = yExtent / 2. + gridDelta;
+        mesh->nodes[6][2] = trenchDepth + baseHeight;
 
         mesh->nodes[7][0] = -trenchWidth / 2. - offset;
-        mesh->nodes[7][1] = yExtent / 2.;
-        mesh->nodes[7][2] = trenchDepth;
+        mesh->nodes[7][1] = yExtent / 2. + gridDelta;
+        mesh->nodes[7][2] = trenchDepth + baseHeight;
 
         mesh->insertNextTriangle(std::array<unsigned, 3>{0, 3, 1});
         mesh->insertNextTriangle(std::array<unsigned, 3>{1, 3, 2});
@@ -184,11 +191,11 @@ template <class NumericType, int D> class psMakeTrench {
       if constexpr (D == 3) {
         minPoint[1] = -yExtent / 2.;
         maxPoint[1] = yExtent / 2.;
-        minPoint[2] = 0.;
-        maxPoint[2] = trenchDepth;
+        minPoint[2] = baseHeight;
+        maxPoint[2] = trenchDepth + baseHeight;
       } else {
-        minPoint[1] = 0.;
-        maxPoint[1] = trenchDepth;
+        minPoint[1] = baseHeight;
+        maxPoint[1] = trenchDepth + baseHeight;
       }
       lsMakeGeometry<NumericType, D>(
           cutout,
diff --git a/include/Models/DirectionalEtching.hpp b/include/Models/DirectionalEtching.hpp
new file mode 100644
index 00000000..2618670e
--- /dev/null
+++ b/include/Models/DirectionalEtching.hpp
@@ -0,0 +1,83 @@
+#pragma once
+
+#include <psProcessModel.hpp>
+#include <psSurfaceModel.hpp>
+#include <psVelocityField.hpp>
+
+// Directional etch for one material
+template <class NumericType, int D>
+class DirectionalEtchVelocityField : public psVelocityField<NumericType> {
+  const std::array<NumericType, 3> direction;
+  const NumericType dirVel = 1.;
+  const NumericType isoVel = 0.;
+  const int maskId;
+
+public:
+  DirectionalEtchVelocityField(std::array<NumericType, 3> dir,
+                               const NumericType dVel, const NumericType iVel,
+                               const int mask = 0)
+      : direction(dir), dirVel(dVel), isoVel(iVel), maskId(mask) {}
+
+  std::array<NumericType, 3>
+  getVectorVelocity(const std::array<NumericType, 3> &coordinate, int material,
+                    const std::array<NumericType, 3> &normalVector,
+                    unsigned long) override {
+    if (material != maskId) {
+      std::array<NumericType, 3> dir(direction);
+      for (unsigned i = 0; i < D; ++i) {
+        if (dir[i] == 0.) {
+          dir[i] -= isoVel * (normalVector[i] < 0 ? -1 : 1);
+        } else {
+          dir[i] *= dirVel;
+        }
+      }
+      return dir;
+    } else {
+      return {0};
+    }
+  }
+
+  NumericType getDissipationAlpha(
+      int /*direction*/, int /*material*/,
+      const std::array<NumericType, 3> & /*centralDifferences*/) {
+    return -1;
+  }
+};
+
+template <typename NumericType, int D>
+class DirectionalEtchingSurfaceModel : public psSurfaceModel<NumericType> {
+public:
+  psSmartPointer<std::vector<NumericType>>
+  calculateVelocities(psSmartPointer<psPointData<NumericType>> Rates,
+                      const std::vector<NumericType> &materialIds) override {
+    return nullptr;
+  }
+};
+
+template <typename NumericType, int D> class DirectionalEtching {
+  psSmartPointer<psProcessModel<NumericType, D>> processModel = nullptr;
+
+public:
+  DirectionalEtching(const std::array<NumericType, 3> direction,
+                     const NumericType directionalVelocity = 1.,
+                     const NumericType isotropicVelocity = 0.) {
+    processModel = psSmartPointer<psProcessModel<NumericType, D>>::New();
+
+    // surface model
+    auto surfModel =
+        psSmartPointer<DirectionalEtchingSurfaceModel<NumericType, D>>::New();
+
+    // velocity field
+    auto velField =
+        psSmartPointer<DirectionalEtchVelocityField<NumericType, D>>::New(
+            direction, directionalVelocity, isotropicVelocity, 0);
+
+    processModel->setSurfaceModel(surfModel);
+    processModel->setVelocityField(velField);
+    processModel->setProcessName("DirectionalEtching");
+  }
+
+  psSmartPointer<psProcessModel<NumericType, D>> getProcessModel() {
+    return processModel;
+  }
+};
\ No newline at end of file
diff --git a/include/Models/GeometricUniformDeposition.hpp b/include/Models/GeometricUniformDeposition.hpp
new file mode 100644
index 00000000..bec0e84a
--- /dev/null
+++ b/include/Models/GeometricUniformDeposition.hpp
@@ -0,0 +1,48 @@
+#pragma once
+
+#include <lsGeometricAdvect.hpp>
+#include <lsGeometricAdvectDistributions.hpp>
+
+#include <psGeometricModel.hpp>
+#include <psProcessModel.hpp>
+#include <psSmartPointer.hpp>
+
+// Simple geometric model that implements a
+template <typename NumericType, int D>
+class GeometricUniformDepositionModel
+    : public psGeometricModel<NumericType, D> {
+  NumericType layerThickness;
+
+  using psGeometricModel<NumericType, D>::domain;
+
+public:
+  GeometricUniformDepositionModel(NumericType passedLayerThickness)
+      : layerThickness(passedLayerThickness) {}
+
+  void apply() {
+    auto dist = psSmartPointer<lsSphereDistribution<double, D>>::New(
+        layerThickness, domain->getGrid().getGridDelta());
+    lsGeometricAdvect<NumericType, D>(domain->getLevelSets()->back(), dist)
+        .apply();
+  }
+};
+
+template <typename NumericType, int D> class GeometricUniformDeposition {
+  psSmartPointer<psProcessModel<NumericType, D>> processModel = nullptr;
+
+public:
+  GeometricUniformDeposition(const NumericType layerThickness = 1.) {
+    processModel = psSmartPointer<psProcessModel<NumericType, D>>::New();
+
+    auto geomModel =
+        psSmartPointer<GeometricUniformDepositionModel<NumericType, D>>::New(
+            layerThickness);
+
+    processModel->setGeometricModel(geomModel);
+    processModel->setProcessName("GeometricUniformDeposition");
+  }
+
+  psSmartPointer<psProcessModel<NumericType, D>> getProcessModel() {
+    return processModel;
+  }
+};
diff --git a/include/Models/SF6O2Etching.hpp b/include/Models/SF6O2Etching.hpp
index 7bdc2a59..7db0c3f2 100644
--- a/include/Models/SF6O2Etching.hpp
+++ b/include/Models/SF6O2Etching.hpp
@@ -162,7 +162,9 @@ class SF6O2Ion : public rayParticle<SF6O2Ion<NumericType, D>, NumericType> {
     assert(cosTheta >= 0 && "Hit backside of disc");
     assert(cosTheta <= 1 + 1e6 && "Error in calculating cos theta");
 
-    const double angle = std::acos(std::max(std::min(cosTheta, 1.), 0.));
+    const double angle =
+        std::acos(std::max(std::min(cosTheta, static_cast<NumericType>(1.)),
+                           static_cast<NumericType>(0.)));
 
     if (cosTheta > 0.5) {
       f_Si_theta = 1.;
@@ -209,7 +211,8 @@ class SF6O2Ion : public rayParticle<SF6O2Ion<NumericType, D>, NumericType> {
     assert(cosTheta <= 1 + 1e-6 && "Error in calculating cos theta");
 
     const NumericType incAngle =
-        std::acos(std::max(std::min(cosTheta, 1.), 0.));
+        std::acos(std::max(std::min(cosTheta, static_cast<NumericType>(1.)),
+                           static_cast<NumericType>(0.)));
 
     NumericType Eref_peak = 0;
 
diff --git a/include/psGDSGeometry.hpp b/include/psGDSGeometry.hpp
index e9912dff..e6aa09c2 100644
--- a/include/psGDSGeometry.hpp
+++ b/include/psGDSGeometry.hpp
@@ -22,7 +22,7 @@ template <class NumericType, int D = 3> class psGDSGeometry {
   std::array<NumericType, 2> minBounds;
   std::array<NumericType, 2> maxBounds;
 
-  double bounds[2 * D];
+  double bounds[6];
   NumericType gridDelta;
   typename lsDomain<NumericType, D>::BoundaryType boundaryCons[3] = {
       lsDomain<NumericType, D>::BoundaryType::REFLECTIVE_BOUNDARY,
@@ -30,10 +30,24 @@ template <class NumericType, int D = 3> class psGDSGeometry {
       lsDomain<NumericType, D>::BoundaryType::INFINITE_BOUNDARY};
 
 public:
-  psGDSGeometry() {}
+  psGDSGeometry() {
+    if constexpr (D == 2) {
+      lsMessage::getInstance()
+          .addWarning("Cannot import 2D geometry from GDS file.")
+          .print();
+      return;
+    }
+  }
 
   psGDSGeometry(const NumericType passedGridDelta)
-      : gridDelta(passedGridDelta) {}
+      : gridDelta(passedGridDelta) {
+    if constexpr (D == 2) {
+      lsMessage::getInstance()
+          .addWarning("Cannot import 2D geometry from GDS file.")
+          .print();
+      return;
+    }
+  }
 
   void setGridDelta(const NumericType passedGridDelta) {
     gridDelta = passedGridDelta;
@@ -48,8 +62,8 @@ template <class NumericType, int D = 3> class psGDSGeometry {
   }
 
   void setBoundaryConditions(
-      typename lsDomain<NumericType, D>::BoundaryType passedBoundaryCons[D]) {
-    for (int i = 0; i < D; i++)
+      typename lsDomain<NumericType, D>::BoundaryType passedBoundaryCons[3]) {
+    for (int i = 0; i < 3; i++)
       boundaryCons[i] = passedBoundaryCons[i];
   }
 
@@ -115,16 +129,15 @@ template <class NumericType, int D = 3> class psGDSGeometry {
             if (sref.angle > 0.) {
               lsTransformMesh<NumericType>(
                   preBuiltStrMesh, lsTransformEnum::ROTATION,
-                  hrleVectorType<NumericType, 3>{0., 0., 1.},
-                  deg2rad(sref.angle))
+                  hrleVectorType<double, 3>{0., 0., 1.}, deg2rad(sref.angle))
                   .apply();
             }
 
             if (sref.magnification > 0.) {
               lsTransformMesh<NumericType>(
                   preBuiltStrMesh, lsTransformEnum::SCALE,
-                  hrleVectorType<NumericType, 3>{sref.magnification,
-                                                 sref.magnification, 1.})
+                  hrleVectorType<double, 3>{sref.magnification,
+                                            sref.magnification, 1.})
                   .apply();
             }
 
@@ -137,8 +150,8 @@ template <class NumericType, int D = 3> class psGDSGeometry {
 
             lsTransformMesh<NumericType>(
                 preBuiltStrMesh, lsTransformEnum::TRANSLATION,
-                hrleVectorType<NumericType, 3>{sref.refPoint[0],
-                                               sref.refPoint[1], 0.})
+                hrleVectorType<double, 3>{sref.refPoint[0], sref.refPoint[1],
+                                          0.})
                 .apply();
 
             strMesh->append(*preBuiltStrMesh);
@@ -158,8 +171,8 @@ template <class NumericType, int D = 3> class psGDSGeometry {
     if (mask) {
       auto topPlane = psSmartPointer<lsDomain<NumericType, D>>::New(
           bounds, boundaryCons, gridDelta);
-      NumericType normal[D] = {0., 0., 1.};
-      NumericType origin[D] = {0., 0., baseHeight + height};
+      NumericType normal[3] = {0., 0., 1.};
+      NumericType origin[3] = {0., 0., baseHeight + height};
       lsMakeGeometry<NumericType, D>(
           topPlane,
           lsSmartPointer<lsPlane<NumericType, D>>::New(origin, normal))
@@ -316,7 +329,7 @@ template <class NumericType, int D = 3> class psGDSGeometry {
     return {minBounds, maxBounds};
   }
 
-  NumericType *getBounds() { return bounds; }
+  double *getBounds() { return bounds; }
 
   void addBox(psSmartPointer<lsDomain<NumericType, D>> levelSet,
               psGDSElement<NumericType> &element, const NumericType baseHeight,
@@ -396,7 +409,7 @@ template <class NumericType, int D = 3> class psGDSGeometry {
 
     // sidewalls
     for (unsigned i = 0; i < numPointsFlat; i++) {
-      std::array<NumericType, D> offsetPoint = element.pointCloud[i];
+      std::array<NumericType, 3> offsetPoint = element.pointCloud[i];
       offsetPoint[0] += xOffset;
       offsetPoint[1] += yOffset;
       offsetPoint[2] = baseHeight;
@@ -408,7 +421,7 @@ template <class NumericType, int D = 3> class psGDSGeometry {
 
     for (unsigned i = 0; i < numPointsFlat; i++) {
       unsigned upPoint = i + numPointsFlat;
-      std::array<NumericType, D> offsetPoint = element.pointCloud[i];
+      std::array<NumericType, 3> offsetPoint = element.pointCloud[i];
       offsetPoint[0] += xOffset;
       offsetPoint[1] += yOffset;
       offsetPoint[2] = baseHeight + height;
diff --git a/include/psGDSReader.hpp b/include/psGDSReader.hpp
index 7e8eee65..1110b9a9 100644
--- a/include/psGDSReader.hpp
+++ b/include/psGDSReader.hpp
@@ -17,7 +17,7 @@
 #define endian_swap_short(w) (((w & 0xff) << 8) | ((w & 0xff00) >> 8))
 #endif
 
-template <typename NumericType, int D> class psGDSReader {
+template <typename NumericType, int D = 3> class psGDSReader {
   using PSPtrType = psSmartPointer<psGDSGeometry<NumericType, D>>;
 
   FILE *filePtr = nullptr;
@@ -34,6 +34,13 @@ template <typename NumericType, int D> class psGDSReader {
   void setFileName(std::string passedFileName) { fileName = passedFileName; }
 
   void apply() {
+    if constexpr (D == 2) {
+      lsMessage::getInstance()
+          .addWarning("Cannot import 2D geometry from GDS file.")
+          .print();
+      return;
+    }
+
     parseFile();
     geometry->checkReferences();
     geometry->calculateBoundingBoxes();
diff --git a/include/psPlanarize.hpp b/include/psPlanarize.hpp
new file mode 100644
index 00000000..c1a800d4
--- /dev/null
+++ b/include/psPlanarize.hpp
@@ -0,0 +1,34 @@
+#pragma once
+
+#include <psDomain.hpp>
+
+#include <lsBooleanOperation.hpp>
+#include <lsMakeGeometry.hpp>
+
+template <class NumericType, int D> class psPlanarize {
+  psSmartPointer<psDomain<NumericType, D>> domain;
+  NumericType cutoffPosition = 0.;
+
+public:
+  psPlanarize() {}
+  psPlanarize(psSmartPointer<psDomain<NumericType, D>> passedDomain,
+              const NumericType passedCutoff)
+      : domain(passedDomain), cutoffPosition(passedCutoff) {}
+
+  void apply() {
+    NumericType origin[D] = {0.};
+    origin[D - 1] = cutoffPosition;
+    NumericType normal[D] = {0.};
+    normal[D - 1] = -1.;
+    auto plane = lsSmartPointer<lsDomain<NumericType, D>>::New(
+        domain->getLevelSets()->back()->getGrid());
+    lsMakeGeometry<NumericType, D>(
+        plane, lsSmartPointer<lsPlane<NumericType, D>>::New(origin, normal))
+        .apply();
+    for (auto &ls : *domain->getLevelSets()) {
+      lsBooleanOperation<NumericType, D>(
+          ls, plane, lsBooleanOperationEnum::RELATIVE_COMPLEMENT)
+          .apply();
+    }
+  }
+};
\ No newline at end of file
diff --git a/include/psProcess.hpp b/include/psProcess.hpp
index 19674356..5cb8d359 100644
--- a/include/psProcess.hpp
+++ b/include/psProcess.hpp
@@ -52,6 +52,10 @@ template <typename NumericType, int D> class psProcess {
     integrationScheme = passedIntegrationScheme;
   }
 
+  void setPrintIntdermediate(const bool passedPrint) {
+    printIntermediate = passedPrint;
+  }
+
   void apply() {
     /* ---------- Process Setup --------- */
     auto name = model->getProcessName();
@@ -142,8 +146,7 @@ template <typename NumericType, int D> class psProcess {
     }
 
     // Determine whether there are process parameters used in ray tracing
-    if (model->getSurfaceModel())
-      model->getSurfaceModel()->initializeProcessParameters();
+    model->getSurfaceModel()->initializeProcessParameters();
     const bool useProcessParams =
         model->getSurfaceModel()->getProcessParameters() != nullptr;
 
@@ -232,8 +235,9 @@ template <typename NumericType, int D> class psProcess {
             diskMesh->getCellData().insertNextScalarData(
                 *Rates->getScalarData(idx), label);
           }
-          printDiskMesh(diskMesh, name + "_covIinit_" +
-                                      std::to_string(iterations) + ".vtp");
+          if (printIntermediate)
+            printDiskMesh(diskMesh, name + "_covIinit_" +
+                                        std::to_string(iterations) + ".vtp");
           std::cerr << "\r"
                     << "Iteration: " << iterations + 1 << " / "
                     << maxIterations;
@@ -324,7 +328,9 @@ template <typename NumericType, int D> class psProcess {
         diskMesh->getCellData().insertNextScalarData(*Rates->getScalarData(idx),
                                                      label);
       }
-      printDiskMesh(diskMesh, name + "_" + std::to_string(counter++) + ".vtp");
+      if (printIntermediate)
+        printDiskMesh(diskMesh,
+                      name + "_" + std::to_string(counter++) + ".vtp");
 #endif
       // apply volume model
       if (useVolumeModel) {
@@ -474,6 +480,7 @@ template <typename NumericType, int D> class psProcess {
   bool useRandomSeeds = true;
   size_t maxIterations = 20;
   bool coveragesInitialized = false;
+  bool printIntermediate = true;
 };
 
 #endif
\ No newline at end of file
diff --git a/include/psVelocityField.hpp b/include/psVelocityField.hpp
index e563002d..e5196037 100644
--- a/include/psVelocityField.hpp
+++ b/include/psVelocityField.hpp
@@ -5,9 +5,6 @@
 #include <vector>
 
 template <typename NumericType> class psVelocityField {
-private:
-  psSmartPointer<std::vector<NumericType>> velocities = nullptr;
-
 public:
   psVelocityField() {}