diff --git a/include/viennaps/psDomain.hpp b/include/viennaps/psDomain.hpp
index 131e69e..6d49132 100644
--- a/include/viennaps/psDomain.hpp
+++ b/include/viennaps/psDomain.hpp
@@ -169,6 +169,60 @@ template <class NumericType, int D> class Domain {
     }
   }
 
+  void removeLevelSet(unsigned int idx, bool removeWrapped = true) {
+    if (idx >= levelSets_.size()) {
+      Logger::getInstance()
+          .addWarning("Trying to remove non-existing Level-Set from domain.")
+          .print();
+      return;
+    }
+
+    if (materialMap_) {
+      auto newMatMap = materialMapType::New();
+      for (std::size_t i = 0; i < levelSets_.size(); i++) {
+        if (i == idx)
+          continue;
+
+        newMatMap->insertNextMaterial(materialMap_->getMaterialAtIdx(i));
+      }
+      materialMap_ = newMatMap;
+    }
+
+    if (removeWrapped) {
+      auto remove = levelSets_.at(idx);
+
+      for (int i = idx - 1; i >= 0; i--) {
+        viennals::BooleanOperation<NumericType, D>(
+            remove, levelSets_.at(i),
+            viennals::BooleanOperationEnum::RELATIVE_COMPLEMENT)
+            .apply();
+      }
+
+      for (int i = idx + 1; i < levelSets_.size(); i++) {
+        viennals::BooleanOperation<NumericType, D>(
+            levelSets_.at(i), remove,
+            viennals::BooleanOperationEnum::RELATIVE_COMPLEMENT)
+            .apply();
+      }
+    }
+
+    levelSets_.erase(levelSets_.begin() + idx);
+    materialMapCheck();
+  }
+
+  void removeMaterial(const Material material) {
+    if (!materialMap_) {
+      return;
+    }
+
+    for (int i = 0; i < materialMap_->size(); i++) {
+      if (materialMap_->getMaterialAtIdx(i) == material) {
+        removeLevelSet(i);
+        i -= 1;
+      }
+    }
+  }
+
   // Generate the Cell-Set from the Level-Sets in the domain. The Cell-Set can
   // be used to store and track volume data.
   void generateCellSet(const NumericType position, const Material coverMaterial,
diff --git a/python/pyWrap.cpp b/python/pyWrap.cpp
index d11635d..45db51c 100644
--- a/python/pyWrap.cpp
+++ b/python/pyWrap.cpp
@@ -829,7 +829,7 @@ PYBIND11_MODULE(VIENNAPS_MODULE_NAME, module) {
            pybind11::arg("direction"),
            pybind11::arg("directionalVelocity") = 1.,
            pybind11::arg("isotropicVelocity") = 0.,
-           pybind11::arg("maskMaterial") = Material::Mask)
+           pybind11::arg("maskMaterial") = Material::None)
       .def(pybind11::init<const std::array<T, 3> &, const T, const T,
                           const std::vector<Material>>(),
            pybind11::arg("direction"), pybind11::arg("directionalVelocity"),
@@ -1144,6 +1144,9 @@ PYBIND11_MODULE(VIENNAPS_MODULE_NAME, module) {
            "by the ray tracer, is averaged over the surface point neighbors.")
       .def("disableFluxSmoothing", &Process<T, D>::disableFluxSmoothing,
            "Disable flux smoothing")
+      .def("setRayTracingDiskRadius", &Process<T, D>::setRayTracingDiskRadius,
+           "Set the radius of the disk used for ray tracing. This disk is used "
+           "for the intersection calculations at each surface point.")
       .def("enableRandomSeeds", &Process<T, D>::enableRandomSeeds,
            "Enable random seeds for the ray tracer. This will make the process "
            "results non-deterministic.")
@@ -1170,9 +1173,13 @@ PYBIND11_MODULE(VIENNAPS_MODULE_NAME, module) {
            pybind11::arg("levelSet"), pybind11::arg("material"),
            pybind11::arg("wrapLowerLevelSet") = true,
            "Insert a level set to domain as a material.")
-      .def("duplicateTopLevelSet", &Domain<T, D>::duplicateTopLevelSet)
+      .def("duplicateTopLevelSet", &Domain<T, D>::duplicateTopLevelSet,
+           "Duplicate the top level set. Should be used before a deposition "
+           "process.")
       .def("removeTopLevelSet", &Domain<T, D>::removeTopLevelSet)
       .def("applyBooleanOperation", &Domain<T, D>::applyBooleanOperation)
+      .def("removeLevelSet", &Domain<T, D>::removeLevelSet)
+      .def("removeMaterial", &Domain<T, D>::removeMaterial)
       .def("setMaterialMap", &Domain<T, D>::setMaterialMap)
       .def("getMaterialMap", &Domain<T, D>::getMaterialMap)
       .def("generateCellSet", &Domain<T, D>::generateCellSet,
diff --git a/tests/domain/domain.cpp b/tests/domain/domain.cpp
index 8317775..c4e0e47 100644
--- a/tests/domain/domain.cpp
+++ b/tests/domain/domain.cpp
@@ -92,6 +92,71 @@ template <class NumericType, int D> void RunTest() {
     VC_TEST_ASSERT(domainCopy->getMaterialMap().get() !=
                    domain->getMaterialMap().get());
   }
+
+  // remove level sets
+  {
+    // two plane geometries
+    ls::BoundaryConditionEnum<D> boundaryCondition[D];
+    double bounds[2 * D];
+
+    for (int i = 0; i < D; ++i) {
+      bounds[2 * i] = -1.;
+      bounds[2 * i + 1] = 1.;
+      boundaryCondition[i] = ls::BoundaryConditionEnum<D>::REFLECTIVE_BOUNDARY;
+    }
+    boundaryCondition[D - 1] = ls::BoundaryConditionEnum<D>::INFINITE_BOUNDARY;
+
+    NumericType origin[D] = {0.};
+    NumericType normal[D] = {0.};
+    normal[D - 1] = 1.;
+
+    auto plane1 = lsDomainType::New(bounds, boundaryCondition, 0.2);
+    ls::MakeGeometry<NumericType, D>(
+        plane1, SmartPointer<ls::Plane<NumericType, D>>::New(origin, normal))
+        .apply();
+
+    origin[D - 1] = 1.;
+    auto plane2 = lsDomainType::New(bounds, boundaryCondition, 0.2);
+    ls::MakeGeometry<NumericType, D>(
+        plane2, SmartPointer<ls::Plane<NumericType, D>>::New(origin, normal))
+        .apply();
+
+    {
+      auto domain = psDomainType::New();
+      domain->insertNextLevelSetAsMaterial(plane1, ps::Material::Si);
+      domain->insertNextLevelSetAsMaterial(plane2, ps::Material::SiO2);
+
+      domain->removeTopLevelSet();
+      VC_TEST_ASSERT(domain->getLevelSets().size() == 1);
+    }
+
+    {
+      auto domain = psDomainType::New();
+      domain->insertNextLevelSetAsMaterial(plane1, ps::Material::Si);
+      domain->insertNextLevelSetAsMaterial(plane2, ps::Material::SiO2);
+
+      domain->removeLevelSet(1);
+      VC_TEST_ASSERT(domain->getLevelSets().size() == 1);
+    }
+
+    {
+      auto domain = psDomainType::New();
+      domain->insertNextLevelSetAsMaterial(plane1, ps::Material::Si);
+      domain->insertNextLevelSetAsMaterial(plane2, ps::Material::SiO2);
+
+      domain->removeMaterial(ps::Material::Si);
+      VC_TEST_ASSERT(domain->getLevelSets().size() == 1);
+    }
+
+    {
+      auto domain = psDomainType::New();
+      domain->insertNextLevelSetAsMaterial(plane1, ps::Material::Si);
+      domain->insertNextLevelSetAsMaterial(plane2, ps::Material::Si);
+
+      domain->removeMaterial(ps::Material::Si);
+      VC_TEST_ASSERT(domain->getLevelSets().empty());
+    }
+  }
 }
 
 } // namespace viennacore