diff --git a/examples/SiGeSelectiveEtching/SiGeEtching.py b/examples/SiGeSelectiveEtching/SiGeEtching.py
index 1c5a1029..5caeb64c 100644
--- a/examples/SiGeSelectiveEtching/SiGeEtching.py
+++ b/examples/SiGeSelectiveEtching/SiGeEtching.py
@@ -40,36 +40,29 @@
 modelParams.Passivation.A_C_ie = params["A_C"]
 
 # Use Material enum
-modelParams.gamma_F = {
-    ps.Material.Mask: 0.0,
-    ps.Material.Si: 0.1,
-    ps.Material.SiGe: 0.1,
-}
-modelParams.gamma_F_oxidized = {
-    ps.Material.Mask: 0.0,
-    ps.Material.Si: 0.1,
-    ps.Material.SiGe: 0.1,
-}
-modelParams.gamma_O = {
-    ps.Material.Mask: 0.0,
-    ps.Material.Si: 0.7,
-    ps.Material.SiGe: 0.7,
-}
-modelParams.gamma_O_passivated = {
-    ps.Material.Mask: 0.0,
-    ps.Material.Si: 0.7,
-    ps.Material.SiGe: 0.7,
-}
-modelParams.gamma_C = {
-    ps.Material.Mask: 0.0,
-    ps.Material.Si: 0.7,
-    ps.Material.SiGe: 0.7,
-}
-modelParams.gamma_C_oxidized = {
-    ps.Material.Mask: 0.0,
-    ps.Material.Si: 0.7,
-    ps.Material.SiGe: 0.7,
-}
+modelParams.gamma_F.set(ps.Material.Mask, 0.0)
+modelParams.gamma_F.set(ps.Material.Si, 0.1)
+modelParams.gamma_F.set(ps.Material.SiGe, 0.1)
+
+modelParams.gamma_F_oxidized.set(ps.Material.Mask, 0.0)
+modelParams.gamma_F_oxidized.set(ps.Material.Si, 0.1)
+modelParams.gamma_F_oxidized.set(ps.Material.SiGe, 0.1)
+
+modelParams.gamma_O.set(ps.Material.Mask, 0.0)
+modelParams.gamma_O.set(ps.Material.Si, 0.7)
+modelParams.gamma_O.set(ps.Material.SiGe, 0.7)
+
+modelParams.gamma_O_passivated.set(ps.Material.Mask, 0.0)
+modelParams.gamma_O_passivated.set(ps.Material.Si, 0.7)
+modelParams.gamma_O_passivated.set(ps.Material.SiGe, 0.7)
+
+modelParams.gamma_C.set(ps.Material.Mask, 0.0)
+modelParams.gamma_C.set(ps.Material.Si, 0.7)
+modelParams.gamma_C.set(ps.Material.SiGe, 0.7)
+
+modelParams.gamma_C_oxidized.set(ps.Material.Mask, 0.0)
+modelParams.gamma_C_oxidized.set(ps.Material.Si, 0.7)
+modelParams.gamma_C_oxidized.set(ps.Material.SiGe, 0.7)
 
 model = ps.CF4O2Etching(modelParams)
 parameters = model.getParameters()
diff --git a/examples/atomicLayerDeposition/atomicLayerDeposition.py b/examples/atomicLayerDeposition/atomicLayerDeposition.py
index 9fec9359..1626ecc9 100644
--- a/examples/atomicLayerDeposition/atomicLayerDeposition.py
+++ b/examples/atomicLayerDeposition/atomicLayerDeposition.py
@@ -1,127 +1,86 @@
 import viennaps as ps
 import viennals as ls
-import matplotlib.pyplot as plt
-import numpy as np
 
 ps.setDimension(2)
 ls.setDimension(2)
 
-ps.Logger.setLogLevel(ps.LogLevel.DEBUG)
-
-
-def run_simulation(gd, filename):
-    params = ps.readConfigFile("config.txt")
-    geometry = ps.Domain()
-
-    # Create the geometry
-    boundaryCons = [
-        ps.BoundaryType.REFLECTIVE_BOUNDARY,
-        ps.BoundaryType.INFINITE_BOUNDARY,
-    ]
-    gridDelta = gd
-    bounds = [
-        0.0,
-        params["openingWidth"] / 2.0 + params["xPad"] + params["gapLength"],
-        -gridDelta,
-        params["openingDepth"] + params["gapHeight"] + gridDelta,
-    ]
-
-    substrate = ls.Domain(bounds, boundaryCons, gridDelta)
-    normal = [0.0, 1.0]
-    origin = [0.0, params["openingDepth"] + params["gapHeight"]]
-    ls.MakeGeometry(substrate, ls.Plane(origin, normal)).apply()
-
-    geometry.insertNextLevelSetAsMaterial(substrate, ps.Material.Si)
-
-    vertBox = ls.Domain(bounds, boundaryCons, gridDelta)
-    minPoint = [-gridDelta, 0.0]
-    maxPoint = [
-        params["openingWidth"] / 2.0,
-        params["gapHeight"] + params["openingDepth"] + gridDelta,
-    ]
-    ls.MakeGeometry(vertBox, ls.Box(minPoint, maxPoint)).apply()
-
-    geometry.applyBooleanOperation(vertBox, ls.BooleanOperationEnum.RELATIVE_COMPLEMENT)
-
-    horiBox = ls.Domain(bounds, boundaryCons, gridDelta)
-    minPoint = [params["openingWidth"] / 2.0 - gridDelta, 0.0]
-    maxPoint = [params["openingWidth"] / 2.0 + params["gapLength"], params["gapHeight"]]
-    ls.MakeGeometry(horiBox, ls.Box(minPoint, maxPoint)).apply()
-    geometry.applyBooleanOperation(horiBox, ls.BooleanOperationEnum.RELATIVE_COMPLEMENT)
-
-    geometry.saveSurfaceMesh("SingleParticleALD_initial")
-
-    model = ps.SingleParticleProcess(rate=1.0, stickingProbability=1e-4)
-
-    geometry.duplicateTopLevelSet(ps.Material.Al2O3)
-
-    gasMFP = ps.constants.gasMeanFreePath(
-        params["pressure"], params["temperature"], params["diameter"]
-    )
-    print("Mean free path: ", gasMFP, " um")
-
-    model = ps.SingleParticleALD(
-        stickingProbability=params["stickingProbability"],
-        numCycles=int(params["numCycles"]),
-        growthPerCycle=params["growthPerCycle"],
-        totalCycles=int(params["totalCycles"]),
-        coverageTimeStep=params["coverageTimeStep"],
-        evFlux=params["evFlux"],
-        inFlux=params["inFlux"],
-        s0=params["s0"],
-        gasMFP=gasMFP,
-    )
-    model.setProcessName(filename)
-
-    alpParams = ps.AtomicLayerProcessParameters()
-    alpParams.pulseTime = params["pulseTime"]
-    alpParams.coverageTimeStep = params["coverageTimeStep"]
-    alpParams.numCycles = 1
-
-    p = ps.Process(geometry, model, 1.0)
-    p.setParameters(alpParams)
-    p.setFluxEngineType(ps.FluxEngineType.GPU_TRIANGLE)
-    p.apply()
-
-    # flux = p.calculateFlux()
-
-    # points = np.array(flux.getNodes())
-    # fluxValues = np.array(flux.getCellData().getScalarData("particleFlux", True))
-
-    # cuts = points[:, 1] < 0.25
-    # points = points[cuts]
-    # fluxValues = fluxValues[cuts]
-
-    # ps.ls.VTKWriter(flux, filename).apply()
-
-    mesh = ps.ls.Mesh()
-    ps.ToDiskMesh(geometry, mesh).apply()
-    points = np.array(mesh.getNodes())
-    cuts = points[:, 1] < 0.25
-    points = points[cuts]
-    return points
-    # return points, fluxValues
-
-
-p = run_simulation(0.1, "flux_0p1")
-plt.plot(p[:, 0], p[:, 1], "--")
-p = run_simulation(0.01, "flux_0p01")
-plt.plot(p[:, 0], p[:, 1], ".-")
-p = run_simulation(0.05, "flux_0p05")
-plt.plot(p[:, 0], p[:, 1], "-")
-plt.show()
-
-
-# rayParams = ps.RayTracingParameters()
-# rayParams.raysPerPoint = int(params["numRaysPerPoint"])
-
-# ALP = ps.Process(geometry, model)
-# ALP.setParameters(rayParams)
-# ALP.setParameters(alpParams)
-# ALP.apply()
+params = ps.readConfigFile("config.txt")
+geometry = ps.Domain()
+
+# Create the geometry
+boundaryCons = [
+    ps.BoundaryType.REFLECTIVE_BOUNDARY,
+    ps.BoundaryType.INFINITE_BOUNDARY,
+]
+gridDelta = params["gridDelta"]
+bounds = [
+    0.0,
+    params["openingWidth"] / 2.0 + params["xPad"] + params["gapLength"],
+    -gridDelta,
+    params["openingDepth"] + params["gapHeight"] + gridDelta,
+]
+
+substrate = ls.Domain(bounds, boundaryCons, gridDelta)
+normal = [0.0, 1.0]
+origin = [0.0, params["openingDepth"] + params["gapHeight"]]
+ls.MakeGeometry(substrate, ls.Plane(origin, normal)).apply()
+
+geometry.insertNextLevelSetAsMaterial(substrate, ps.Material.Si)
+
+vertBox = ls.Domain(bounds, boundaryCons, gridDelta)
+minPoint = [-gridDelta, 0.0]
+maxPoint = [
+    params["openingWidth"] / 2.0,
+    params["gapHeight"] + params["openingDepth"] + gridDelta,
+]
+ls.MakeGeometry(vertBox, ls.Box(minPoint, maxPoint)).apply()
+
+geometry.applyBooleanOperation(vertBox, ls.BooleanOperationEnum.RELATIVE_COMPLEMENT)
+
+horiBox = ls.Domain(bounds, boundaryCons, gridDelta)
+minPoint = [params["openingWidth"] / 2.0 - gridDelta, 0.0]
+maxPoint = [params["openingWidth"] / 2.0 + params["gapLength"], params["gapHeight"]]
+ls.MakeGeometry(horiBox, ls.Box(minPoint, maxPoint)).apply()
+geometry.applyBooleanOperation(horiBox, ls.BooleanOperationEnum.RELATIVE_COMPLEMENT)
+
+geometry.saveSurfaceMesh("SingleParticleALD_initial")
+
+model = ps.SingleParticleProcess(rate=1.0, stickingProbability=1e-4)
+
+geometry.duplicateTopLevelSet(ps.Material.Al2O3)
+
+gasMFP = ps.constants.gasMeanFreePath(
+    params["pressure"], params["temperature"], params["diameter"]
+)
+print("Mean free path: ", gasMFP, " um")
+
+model = ps.SingleParticleALD(
+    stickingProbability=params["stickingProbability"],
+    numCycles=int(params["numCycles"]),
+    growthPerCycle=params["growthPerCycle"],
+    totalCycles=int(params["totalCycles"]),
+    coverageTimeStep=params["coverageTimeStep"],
+    evFlux=params["evFlux"],
+    inFlux=params["inFlux"],
+    s0=params["s0"],
+    gasMFP=gasMFP,
+)
+
+alpParams = ps.AtomicLayerProcessParameters()
+alpParams.pulseTime = params["pulseTime"]
+alpParams.coverageTimeStep = params["coverageTimeStep"]
+alpParams.numCycles = int(params["numCycles"])
+
+rayParams = ps.RayTracingParameters()
+rayParams.raysPerPoint = int(params["numRaysPerPoint"])
+
+ALP = ps.Process(geometry, model)
+ALP.setParameters(rayParams)
+ALP.setParameters(alpParams)
+ALP.apply()
 
 # ## TODO: Implement MeasureProfile in Python
 # #   MeasureProfile<NumericType, D>(domain, params.get("gapHeight") / 2.)
 # #       .save(params.get<std::string>("outputFile"));
 
-# geometry.saveVolumeMesh("SingleParticleALD_final")
+geometry.saveSurfaceMesh("SingleParticleALD_final")
diff --git a/examples/blazedGratingsEtching/blazedGratingsEtching.py b/examples/blazedGratingsEtching/blazedGratingsEtching.py
index fa09aff8..4f88ff91 100644
--- a/examples/blazedGratingsEtching/blazedGratingsEtching.py
+++ b/examples/blazedGratingsEtching/blazedGratingsEtching.py
@@ -64,7 +64,8 @@
 ibeParams = ps.IBEParameters()
 ibeParams.exponent = params["exponent"]
 ibeParams.meanEnergy = params["meanEnergy"]
-ibeParams.materialPlaneWaferRate = {ps.Material.SiO2: 1, ps.Material.Mask: 1 / 11}
+ibeParams.materialPlaneWaferRate.set(ps.Material.SiO2, 1.0)
+ibeParams.materialPlaneWaferRate.set(ps.Material.Mask, 1.0 / 11.0)
 ibeParams.cos4Yield.isDefined = True
 ibeParams.cos4Yield.a1 = params["yieldFactor"]
 ibeParams.cos4Yield.a2 = -1.55
diff --git a/examples/emulation/FinFET.cpp b/examples/emulation/FinFET.cpp
index e1798149..b088bbb9 100644
--- a/examples/emulation/FinFET.cpp
+++ b/examples/emulation/FinFET.cpp
@@ -1,7 +1,5 @@
 #include <viennaps.hpp>
 
-#include <lsReader.hpp>
-
 using namespace viennaps;
 
 constexpr int D = 3;
@@ -11,19 +9,9 @@ using DomainType = SmartPointer<Domain<NumericType, D>>;
 using IsotropicProcessGeometric = SmartPointer<SphereDistribution<double, D>>;
 constexpr bool volumeOutput = false;
 
-void writeVolume(DomainType domain) {
-  if (!volumeOutput)
-    return;
-  static int volumeNum = 0;
-  std::cout << "Writing volume mesh ..." << std::flush;
-  domain->saveVolumeMesh("FinFET_" + std::to_string(volumeNum), 0.05);
-  std::cout << " done" << std::endl;
-  ++volumeNum;
-}
-
 void writeSurface(DomainType domain) {
   static int outputNum = 0;
-  domain->saveSurfaceMesh("FinFET_" + std::to_string(outputNum) + ".vtp", true);
+  domain->saveSurfaceMesh("FinFET_" + std::to_string(outputNum) + ".vtp");
   ++outputNum;
 }
 
@@ -60,7 +48,7 @@ int main() {
   { // DP-Depo
     std::cout << "DP-Depo ..." << std::flush;
     const NumericType thickness = 4; // nm
-    domain->duplicateTopLevelSet(Material::Metal);
+    domain->duplicateTopLevelSet(Material::SiO2_HM);
     auto dist = IsotropicProcessGeometric::New(thickness);
     Process<NumericType, D>(domain, dist).apply();
     std::cout << " done" << std::endl;
@@ -81,19 +69,21 @@ int main() {
   // Remove mask with boolean operation
   domain->removeMaterial(Material::Mask);
   writeSurface(domain);
-  writeVolume(domain);
 
   // pattern si
   {
     std::cout << "Si-Patterning ..." << std::flush;
     const NumericType etchDepth = 90.; // nm
-    Vec3D<NumericType> direction = {0, 0, 1};
-    auto model = SmartPointer<DirectionalProcess<NumericType, D>>::New(
-        direction, 1.1, 0.1, Material::Metal, false);
+    Vec3D<NumericType> direction = {0, 0, -1};
+    DirectionalProcess<NumericType, D>::RateMap rates;
+    rates.set(Material::Si,
+              {1., 0.1});       // directional and isotropic rates for silicon
+    rates.setDefault({0., 0.}); // default rates for non-silicon materials
+    auto model =
+        SmartPointer<DirectionalProcess<NumericType, D>>::New(direction, rates);
     Process<NumericType, D>(domain, model, etchDepth).apply();
     std::cout << " done" << std::endl;
   }
-  writeVolume(domain);
   writeSurface(domain);
 
   // Remove DP mask (metal)
@@ -113,7 +103,7 @@ int main() {
 
   // CMP at 80
   Planarize<NumericType, D>(domain, 80.0).apply();
-  writeVolume(domain);
+  writeSurface(domain);
 
   // pattern STI material
   {
@@ -124,7 +114,6 @@ int main() {
     std::cout << " done" << std::endl;
   }
   writeSurface(domain);
-  writeVolume(domain);
 
   // deposit gate material
   {
@@ -179,7 +168,6 @@ int main() {
   // Remove mask
   domain->removeTopLevelSet();
   writeSurface(domain);
-  writeVolume(domain);
 
   // Spacer Deposition and Etch
   { // spacer depo
@@ -194,22 +182,20 @@ int main() {
 
   { // spacer etch
     std::cout << "Spacer Etch ..." << std::flush;
-    auto ls = domain->getLevelSets()[domain->getLevelSets().size() - 2];
+    auto mask = domain->getLevelSets()[domain->getLevelSets().size() - 2];
     auto dist = SmartPointer<BoxDistribution<double, D>>::New(
-        std::array<NumericType, D>{-gridDelta, -gridDelta, -50}, ls);
+        std::array<NumericType, D>{-gridDelta, -gridDelta, -50}, mask);
     Process<NumericType, D>(domain, dist).apply();
     std::cout << " done" << std::endl;
   }
   writeSurface(domain);
-  writeVolume(domain);
 
   // isotropic etch (fin-release)
   {
     std::cout << "Fin-Release ..." << std::flush;
-    std::vector<Material> masks = {Material::PolySi, Material::SiO2,
-                                   Material::Si3N4};
-    auto model =
-        SmartPointer<IsotropicProcess<NumericType, D>>::New(-1., masks);
+    auto model = SmartPointer<IsotropicProcess<NumericType, D>>::New(0.0);
+    model->setMaterialRate(Material::Si, -1.);
+
     AdvectionParameters advectionParams;
     advectionParams.spatialScheme =
         viennals::SpatialSchemeEnum::LAX_FRIEDRICHS_2ND_ORDER;
@@ -219,22 +205,21 @@ int main() {
     std::cout << " done" << std::endl;
   }
   writeSurface(domain);
-  writeVolume(domain);
 
   // source/drain epitaxy
   {
     std::cout << "S/D Epitaxy ..." << std::flush;
     domain->duplicateTopLevelSet(Material::SiGe);
-    Logger::setLogLevel(LogLevel::INFO);
     AdvectionParameters advectionParams;
     advectionParams.spatialScheme =
         viennals::SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER;
     lsInternal::StencilLocalLaxFriedrichsScalar<NumericType, D,
                                                 1>::setMaxDissipation(100);
 
-    std::vector<std::pair<Material, NumericType>> material = {
-        {Material::Si, 1.}, {Material::SiGe, 1.}};
-    auto model = SmartPointer<SelectiveEpitaxy<NumericType, D>>::New(material);
+    auto model = SmartPointer<SelectiveEpitaxy<NumericType, D>>::New();
+    model->setMaterialRate(Material::Si, 2.);
+    model->setMaterialRate(Material::SiGe, 2.);
+
     Process<NumericType, D> proc(domain, model, 14.);
     proc.setParameters(advectionParams);
     proc.apply();
@@ -242,7 +227,6 @@ int main() {
     std::cout << " done" << std::endl;
   }
   writeSurface(domain);
-  writeVolume(domain);
 
   // deposit dielectric
   {
@@ -257,12 +241,11 @@ int main() {
 
   // CMP at 90
   Planarize<NumericType, D>(domain, 90.0).apply();
-  writeVolume(domain);
+  writeSurface(domain);
 
-  // now remove gate and add new gate materials
+  // now remove dummy gate and add new gate materials
   domain->removeMaterial(Material::PolySi);
   writeSurface(domain);
-  writeVolume(domain);
 
   // now deposit TiN and PolySi as replacement gate
   {
@@ -283,10 +266,11 @@ int main() {
     Process<NumericType, D>(domain, dist).apply();
     std::cout << " done" << std::endl;
   }
+  writeSurface(domain);
 
   // CMP at 90
   Planarize<NumericType, D>(domain, 90.0).apply();
-  writeVolume(domain);
+  writeSurface(domain);
 
   domain->saveVolumeMesh("FinFET_Final", 0.05);
 }
diff --git a/examples/emulation/FinFET.py b/examples/emulation/FinFET.py
index 13c18761..fc4c000d 100644
--- a/examples/emulation/FinFET.py
+++ b/examples/emulation/FinFET.py
@@ -1,21 +1,9 @@
 import viennaps as ps
 
 ps.setDimension(3)
-volumeOutput = False
-volumeNum = 0
 surfaceNum = 0
 
 
-def writeVolume(domain):
-    if not volumeOutput:
-        return
-    global volumeNum
-    print("Writing volume mesh ...", end="", flush=True)
-    domain.saveVolumeMesh("FinFET_" + str(volumeNum))
-    print(" done")
-    volumeNum += 1
-
-
 def writeSurface(domain):
     global surfaceNum
     domain.saveSurfaceMesh("FinFET_" + str(surfaceNum) + ".vtp", True)
@@ -38,7 +26,6 @@ def writeSurface(domain):
 ps.MakePlane(domain, 70.0, ps.Material.Si).apply()
 writeSurface(domain)
 
-
 # Add double patterning mask
 box = ps.ls.Domain(domain.getGrid())
 minPoints = [30.0, -10.0, 69.9]
@@ -60,7 +47,7 @@ def writeSurface(domain):
 print("DP-Patterning ...", end="", flush=True)
 etchDepth = 6.0  # nm
 dist = ps.BoxDistribution(
-    [-gridDelta, -gridDelta, -etchDepth], domain.getLevelSets()[0]
+    halfAxes=[-gridDelta, -gridDelta, -etchDepth], mask=domain.getLevelSets()[0]
 )
 ps.Process(domain, dist, 0).apply()
 print(" done")
@@ -69,22 +56,16 @@ def writeSurface(domain):
 # Remove mask with boolean operation
 domain.removeMaterial(ps.Material.Mask)
 writeSurface(domain)
-writeVolume(domain)
 
 # pattern si
 print("Si-Patterning ...", end="", flush=True)
 etchDepth = 90.0  # nm
-direction = [0.0, 0.0, 1.0]
 model = ps.DirectionalProcess(
-    direction=direction,
-    directionalVelocity=1.1,
-    isotropicVelocity=0.1,
-    maskMaterial=ps.Material.Metal,
-    calculateVisibility=False,
+    direction=[0.0, 0.0, -1.0],
+    materialRates={ps.Material.Si: [1.0, 0.1]},
 )
 ps.Process(domain, model, etchDepth).apply()
 print(" done")
-writeVolume(domain)
 writeSurface(domain)
 
 # Remove DP mask (metal)
@@ -102,7 +83,7 @@ def writeSurface(domain):
 
 # CMP at 80
 ps.Planarize(domain, 80.0).apply()
-writeVolume(domain)
+writeSurface(domain)
 
 # pattern STI material
 print("STI Patterning ...", end="", flush=True)
@@ -110,7 +91,6 @@ def writeSurface(domain):
 ps.Process(domain, dist, 0).apply()
 print(" done")
 writeSurface(domain)
-writeVolume(domain)
 
 # deposit gate material
 print("Gate Deposition HfO2 ...", end="", flush=True)
@@ -131,6 +111,7 @@ def writeSurface(domain):
 
 # CMP at 150
 ps.Planarize(domain, 150.0).apply()
+writeSurface(domain)
 
 # dummy gate mask addition
 box = ps.ls.Domain(domain.getGrid())
@@ -158,7 +139,6 @@ def writeSurface(domain):
 # Remove mask
 domain.removeTopLevelSet()
 writeSurface(domain)
-writeVolume(domain)
 
 # Spacer Deposition and Etch
 print("Spacer Deposition and Etch ...", end="", flush=True)
@@ -176,12 +156,11 @@ def writeSurface(domain):
 ps.Process(domain, dist, 0).apply()
 print(" done")
 writeSurface(domain)
-writeVolume(domain)
 
 # isotropic etch (fin-release)
 print("Fin-Release ...", end="", flush=True)
-masks = [ps.Material.PolySi, ps.Material.SiO2, ps.Material.Si3N4]
-model = ps.IsotropicProcess(rate=-1.0, maskMaterial=masks)
+model = ps.IsotropicProcess(rate=0.0)
+model.setMaterialRate(ps.Material.Si, -1.0)
 advParams = ps.AdvectionParameters()
 advParams.spatialScheme = ps.SpatialScheme.LAX_FRIEDRICHS_2ND_ORDER
 process = ps.Process(domain, model, 5.0)
@@ -189,25 +168,21 @@ def writeSurface(domain):
 process.apply()
 print(" done")
 writeSurface(domain)
-writeVolume(domain)
 
 # source/drain epitaxy
 print("S/D Epitaxy ...", end="", flush=True)
 domain.duplicateTopLevelSet(ps.Material.SiGe)
 advectionParams = ps.AdvectionParameters()
 advectionParams.spatialScheme = ps.SpatialScheme.STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER
-ps.StencilLocalLaxFriedrichsScalar.setMaxDissipation(1000)
-material = [
-    (ps.Material.Si, 1.0),
-    (ps.Material.SiGe, 1.0),
-]
-model = ps.SelectiveEpitaxy(materialRates=material)
+ps.StencilLocalLaxFriedrichsScalar.setMaxDissipation(100)
+model = ps.SelectiveEpitaxy()
+model.setMaterialRate(ps.Material.SiGe, 2.0)
+model.setMaterialRate(ps.Material.Si, 2.0)
 process = ps.Process(domain, model, 14.0)
 process.setParameters(advectionParams)
 process.apply()
 print(" done")
 writeSurface(domain)
-writeVolume(domain)
 
 # deposit dielectric
 print("Dielectric Deposition ...", end="", flush=True)
@@ -220,12 +195,11 @@ def writeSurface(domain):
 
 # CMP at 90
 ps.Planarize(domain, 90.0).apply()
-writeVolume(domain)
+writeSurface(domain)
 
 # now remove gate and add new gate materials
 domain.removeMaterial(ps.Material.PolySi)
 writeSurface(domain)
-writeVolume(domain)
 
 # now deposit TiN and PolySi as replacement gate
 print("Gate Deposition TiN ...", end="", flush=True)
@@ -242,9 +216,10 @@ def writeSurface(domain):
 dist = ps.SphereDistribution(radius=thickness)
 ps.Process(domain, dist, 0).apply()
 print(" done")
+writeSurface(domain)
 
 # CMP at 90
 ps.Planarize(domain, 90.0).apply()
-writeVolume(domain)
+writeSurface(domain)
 
 domain.saveVolumeMesh("FinFET_Final", 0.05)
diff --git a/examples/emulation/stackedNanowire.py b/examples/emulation/stackedNanowire.py
index 8dedc603..5dfedeba 100644
--- a/examples/emulation/stackedNanowire.py
+++ b/examples/emulation/stackedNanowire.py
@@ -36,7 +36,7 @@
 geo.apply()
 
 domain.insertNextLevelSetAsMaterial(mask, ps.Material.Mask)
-domain.saveVolumeMesh("StackedNanowire_" + str(n))
+domain.saveSurfaceMesh("StackedNanowire_" + str(n))
 n += 1
 
 # Double patterning
@@ -53,7 +53,7 @@
 ps.Process(domain, directional, 20.0).apply()
 
 domain.removeMaterial(ps.Material.Mask)
-domain.saveVolumeMesh("StackedNanowire_" + str(n))
+domain.saveSurfaceMesh("StackedNanowire_" + str(n))
 n += 1
 
 # Pattern Si/SiGe/Si
@@ -68,7 +68,7 @@
 ps.Process(domain, directional, 30.0).apply()
 
 domain.removeMaterial(ps.Material.Metal)
-domain.saveVolumeMesh("StackedNanowire_" + str(n))
+domain.saveSurfaceMesh("StackedNanowire_" + str(n))
 n += 1
 
 # Deposit dummy gate material
@@ -79,7 +79,7 @@
 # CMP at 80nm height
 ps.Planarize(domain, 80.0).apply()
 
-domain.saveVolumeMesh("StackedNanowire_" + str(n))
+domain.saveSurfaceMesh("StackedNanowire_" + str(n))
 n += 1
 
 # Dummy gate mask addition
@@ -101,70 +101,52 @@
 ps.Process(tmpDomain, geometricEtch, 1.0).apply()
 
 domain.removeMaterial(ps.Material.Mask)
-domain.saveVolumeMesh("StackedNanowire_" + str(n))
+domain.saveSurfaceMesh("StackedNanowire_" + str(n))
 n += 1
 
 # Spacer Deposition
 print("Depositing spacer material ...")
 domain.duplicateTopLevelSet(ps.Material.Si3N4)
 ps.Process(domain, growth, 12.0).apply()
+domain.saveSurfaceMesh("StackedNanowire_" + str(n))
+n += 1
 
 # Spacer patterning
 print("Patterning spacer ...")
 directional = ps.DirectionalProcess(
-    direction=[0.0, 0.0, 1.0],
-    directionalVelocity=1.0,
-    isotropicVelocity=-0.05,
-    maskMaterial=[
-        ps.Material.PolySi,
-        ps.Material.Si,
-        ps.Material.SiGe,
-        ps.Material.SiO2,
-    ],
-    calculateVisibility=False,
+    direction=[0.0, 0.0, -1.0], materialRates={ps.Material.Si3N4: [1.0, 0.0]}
 )
 ps.Process(domain, directional, 40.0).apply()
 
-domain.saveVolumeMesh("StackedNanowire_" + str(n))
+domain.saveSurfaceMesh("StackedNanowire_" + str(n))
 n += 1
 
 # Fin patterning
 print("Patterning Fin ...")
-rateSet1 = ps.RateSet(
-    direction=[0.0, 0.0, 1.0],
-    directionalVelocity=1.0,
-    isotropicVelocity=0.0,
-    maskMaterials=[ps.Material.PolySi, ps.Material.SiO2, ps.Material.Si3N4],
-    calculateVisibility=False,
-)
-rateSet2 = ps.RateSet(
-    direction=[0.0, 0.0, 1.0],
-    directionalVelocity=0.1,
-    isotropicVelocity=0.0,
-    maskMaterials=[
-        ps.Material.SiO2,
-        ps.Material.Si3N4,
-        ps.Material.Si,
-        ps.Material.SiGe,
-    ],
-    calculateVisibility=False,
+directional = ps.DirectionalProcess(
+    direction=[0.0, 0.0, -1.0],
+    materialRates={
+        ps.Material.PolySi: [0.1, 0.0],
+        ps.Material.Si: [1.0, 0.0],
+        ps.Material.SiGe: [1.0, 0.0],
+    },
+    defaultIsotropicRate=-0.005,
 )
-directional = ps.DirectionalProcess([rateSet1, rateSet2])
 ps.Process(domain, directional, 21.0).apply()
 
-domain.saveVolumeMesh("StackedNanowire_" + str(n))
+domain.saveSurfaceMesh("StackedNanowire_" + str(n))
 n += 1
 
 # SD Epitaxy
 print("SD Epitaxy ...")
-domain.duplicateTopLevelSet(ps.Material.Metal)
-epitaxy = ps.IsotropicProcess(
-    rate=1.0,
-    maskMaterial=[ps.Material.PolySi, ps.Material.SiO2, ps.Material.Si3N4],
-)
+epiMaterial = ps.MaterialMap.fromString("EpiSi")
+domain.duplicateTopLevelSet(epiMaterial)
+epitaxy = ps.IsotropicProcess(0.0)
+epitaxy.setMaterialRate(ps.Material.Si, 1.0)
+epitaxy.setMaterialRate(epiMaterial, 1.0)
 ps.Process(domain, epitaxy, 9.0).apply()
 
-domain.saveVolumeMesh("StackedNanowire_" + str(n))
+domain.saveSurfaceMesh("StackedNanowire_" + str(n))
 n += 1
 
 # Dielectric deposition
@@ -174,30 +156,22 @@
 
 ps.Planarize(domain, 60).apply()
 
-domain.saveVolumeMesh("StackedNanowire_" + str(n))
+domain.saveSurfaceMesh("StackedNanowire_" + str(n))
 n += 1
 
 # Remove dummy gate
 print("Removing dummy gate ...")
 domain.removeMaterial(ps.Material.PolySi)
-domain.saveVolumeMesh("StackedNanowire_" + str(n))
+domain.saveSurfaceMesh("StackedNanowire_" + str(n))
 n += 1
 
 # remove SiGe interlayer
 print("Removing SiGe interlayer ...")
-etch = ps.IsotropicProcess(
-    rate=-1.0,
-    maskMaterial=[
-        ps.Material.Si,
-        ps.Material.SiO2,
-        ps.Material.Si3N4,
-        ps.Material.Dielectric,
-        ps.Material.Metal,
-    ],
-)
+etch = ps.IsotropicProcess(rate=0.0)  # default rate is 0.0
+etch.setMaterialRate(ps.Material.SiGe, -1.0)
 ps.Process(domain, etch, 10.0).apply()
 
-domain.saveVolumeMesh("StackedNanowire_" + str(n))
+domain.saveSurfaceMesh("StackedNanowire_" + str(n))
 n += 1
 
 # gate dielectric, gate metal, gate contact deposition
@@ -213,4 +187,5 @@
 
 ps.Planarize(domain, 47.5).apply()
 
-domain.saveVolumeMesh("StackedNanowire_" + str(n))
+domain.saveSurfaceMesh("StackedNanowire_" + str(n))
+domain.saveVolumeMesh("StackedNanowire_Final")
diff --git a/examples/exampleProcess/velocityField.hpp b/examples/exampleProcess/velocityField.hpp
index cf2207aa..7b2ecef4 100644
--- a/examples/exampleProcess/velocityField.hpp
+++ b/examples/exampleProcess/velocityField.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include <materials/psMaterials.hpp>
+#include <materials/psMaterialMap.hpp>
 #include <process/psVelocityField.hpp>
 #include <vector>
 
diff --git a/examples/holeEtching/holeEtching.py b/examples/holeEtching/holeEtching.py
index 31ad612b..30c60b24 100644
--- a/examples/holeEtching/holeEtching.py
+++ b/examples/holeEtching/holeEtching.py
@@ -89,7 +89,8 @@ def run_simulation(intermediate_velocities, suffix):
         output_file += suffix
 
     geometry.saveSurfaceMesh(
-        filename=output_file, addInterfaces=True, boolMaterials=True
+        filename=output_file,
+        addInterfaces=True,
     )
 
 
diff --git a/examples/ionBeamEtching/ionBeamEtching.py b/examples/ionBeamEtching/ionBeamEtching.py
index 60e5685c..6ba8250c 100644
--- a/examples/ionBeamEtching/ionBeamEtching.py
+++ b/examples/ionBeamEtching/ionBeamEtching.py
@@ -35,7 +35,7 @@
 # copy top layer to capture deposition
 geometry.duplicateTopLevelSet(ps.Material.Polymer)
 
-ibeParams = ps.IonBeamEtching.defaultParamters()
+ibeParams = ps.IonBeamEtching.defaultParameters()
 ibeParams.tiltAngle = params["angle"]
 ibeParams.exponent = params["exponent"]
 ibeParams.thetaRMin = 0.0
diff --git a/examples/selectiveEpitaxy/config.txt b/examples/selectiveEpitaxy/config.txt
index d7f121f4..73d9fdf0 100644
--- a/examples/selectiveEpitaxy/config.txt
+++ b/examples/selectiveEpitaxy/config.txt
@@ -12,6 +12,6 @@ etchTime = 4.0
 
 # Process
 processTime=5.
-epitaxyRate=1.
+epitaxyRate=2.
 R111=0.4
 R100=1.0
\ No newline at end of file
diff --git a/gpu/models/PlasmaEtching.cuh b/gpu/models/PlasmaEtching.cuh
index 79cc431b..f7cae010 100644
--- a/gpu/models/PlasmaEtching.cuh
+++ b/gpu/models/PlasmaEtching.cuh
@@ -6,7 +6,7 @@
 #include "raygLaunchParams.hpp"
 #include "raygReflection.hpp"
 
-#include "materials/psMaterials.hpp"
+#include "materials/psMaterialMap.hpp"
 #include "models/psPipelineParameters.hpp"
 #include "models/psPlasmaEtchingParameters.hpp"
 
diff --git a/include/viennaps/materials/psMaterials.hpp b/include/viennaps/materials/psMaterialMap.hpp
similarity index 100%
rename from include/viennaps/materials/psMaterials.hpp
rename to include/viennaps/materials/psMaterialMap.hpp
diff --git a/include/viennaps/materials/psMaterialValueMap.hpp b/include/viennaps/materials/psMaterialValueMap.hpp
index 2699724d..84fc0848 100644
--- a/include/viennaps/materials/psMaterialValueMap.hpp
+++ b/include/viennaps/materials/psMaterialValueMap.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include "psMaterials.hpp"
+#include "psMaterial.hpp"
 
 namespace viennaps {
 
@@ -11,12 +11,20 @@ template <class T> class MaterialValueMap {
   MaterialValueMap() = default;
 
   // generic constructor (works with map, vector<pair>, initializer_list, etc.)
-  template <class MapLike> explicit MaterialValueMap(const MapLike &mapLike) {
+  template <class MapLike>
+  explicit MaterialValueMap(const MapLike &mapLike, T defaultValue = T{})
+      : default_(defaultValue) {
     for (const auto &[key, value] : mapLike) {
       set(key, value);
     }
   }
 
+  static MaterialValueMap fromDefault(T defaultValue) {
+    MaterialValueMap map;
+    map.setDefault(defaultValue);
+    return map;
+  }
+
   void set(Material material, const T &value) {
     if (material.isBuiltIn()) {
       auto idx = toIndex(material.builtIn());
@@ -79,6 +87,14 @@ template <class T> class MaterialValueMap {
     return getBuiltInValue(idx);
   }
 
+  const T getEntryByIndex(std::size_t idx) const {
+    if (auto it = begin().goToIndex(idx); it != end()) {
+      return (*it).value;
+    } else {
+      throw std::out_of_range("Index out of range in MaterialValueMap.");
+    }
+  }
+
   void setDefault(const T &v) { default_ = v; }
 
   [[nodiscard]] const T &getDefault() const { return default_; }
@@ -93,6 +109,15 @@ template <class T> class MaterialValueMap {
 
   [[nodiscard]] bool has(BuiltInMaterial m) const { return has(Material(m)); }
 
+  bool empty() const {
+    for (const auto &b : isSet_) {
+      if (b) {
+        return false;
+      }
+    }
+    return customValues_.empty();
+  }
+
   // remove value -> fallback to default
   void clear(Material material) {
     if (material.isBuiltIn()) {
@@ -113,55 +138,103 @@ template <class T> class MaterialValueMap {
 
   // ================= ITERATOR =================
   struct Entry {
-    BuiltInMaterial material;
-    const T *value;
-    bool set;
-
-    [[nodiscard]] bool isSet() const { return set; }
-    [[nodiscard]] const T &getValue() const { return *value; }
-    [[nodiscard]] BuiltInMaterial getMaterial() const { return material; }
+    Material material;
+    const T &value;
   };
 
   class Iterator {
   public:
-    Iterator(const MaterialValueMap *map, std::size_t idx)
-        : map_(map), idx_(idx) {
+    using CustomIterator =
+        typename std::unordered_map<Material::ValueType, T>::const_iterator;
+
+    Iterator(const MaterialValueMap *map, bool isEnd)
+        : map_(map), phase_(isEnd ? Phase::End : Phase::BuiltIn), idx_(0) {
+      if (!isEnd) {
+        advanceToValid();
+      }
+    }
+
+    Iterator goToIndex(std::size_t targetIdx) {
+      idx_ = 0;
+      phase_ = Phase::BuiltIn;
       advanceToValid();
+
+      size_t idx = 0;
+      while (phase_ != Phase::End && idx != targetIdx) {
+        ++(*this);
+        ++idx;
+      }
+      return *this;
     }
 
     Iterator &operator++() {
-      ++idx_;
-      advanceToValid();
+      if (phase_ == Phase::BuiltIn) {
+        ++idx_;
+        advanceToValid();
+      } else if (phase_ == Phase::Custom) {
+        ++customIt_;
+        if (customIt_ == map_->customValues_.cend()) {
+          phase_ = Phase::End;
+        }
+      }
       return *this;
     }
 
-    bool operator!=(const Iterator &other) const { return idx_ != other.idx_; }
+    bool operator!=(const Iterator &other) const {
+      if (phase_ != other.phase_) {
+        return true;
+      }
+
+      switch (phase_) {
+      case Phase::BuiltIn:
+        return idx_ != other.idx_;
+      case Phase::Custom:
+        return customIt_ != other.customIt_;
+      case Phase::End:
+        return false;
+      }
+      return false;
+    }
 
     Entry operator*() const {
-      return Entry{static_cast<BuiltInMaterial>(idx_), &map_->values_[idx_],
-                   map_->isSet_[idx_]};
+      if (phase_ == Phase::BuiltIn) {
+        return Entry{Material(static_cast<BuiltInMaterial>(idx_)),
+                     map_->values_[idx_]};
+      }
+
+      return Entry{Material::custom(customIt_->first), customIt_->second};
     }
 
   private:
+    enum class Phase { BuiltIn, Custom, End };
+
     void advanceToValid() {
       while (idx_ < map_->values_.size() && !map_->isSet_[idx_]) {
         ++idx_;
       }
+
+      if (idx_ >= map_->values_.size()) {
+        phase_ = Phase::Custom;
+        customIt_ = map_->customValues_.cbegin();
+        if (customIt_ == map_->customValues_.cend()) {
+          phase_ = Phase::End;
+        }
+      }
     }
 
     const MaterialValueMap *map_;
+    Phase phase_;
     std::size_t idx_;
+    CustomIterator customIt_{};
   };
 
-  // iterators over built-in materials with user-provided values
-  Iterator begin() const { return Iterator(this, 0); }
+  // iterators over set built-in materials and custom materials
+  Iterator begin() const { return Iterator(this, false); }
 
-  Iterator end() const { return Iterator(this, values_.size()); }
+  Iterator end() const { return Iterator(this, true); }
 
   const auto &getCustomMaterialValues() const { return customValues_; }
 
-  auto &customMaterialRates() { return customValues_; }
-
 private:
   static constexpr std::size_t toIndex(BuiltInMaterial m) {
     return static_cast<std::uint16_t>(m);
diff --git a/include/viennaps/models/psCF4O2Etching.hpp b/include/viennaps/models/psCF4O2Etching.hpp
index 5dc4cc9a..b62ec9df 100644
--- a/include/viennaps/models/psCF4O2Etching.hpp
+++ b/include/viennaps/models/psCF4O2Etching.hpp
@@ -370,8 +370,10 @@ class CF4O2Etchant final
       const auto &phi_O = globalData->getVectorData(1)[primID];
       // F surface coverage on oxidized SiGe
       const auto &phi_C = globalData->getVectorData(2)[primID];
-      NumericType gamma_F = sticking(materialId);
-      NumericType gamma_FO = stickingOxidized(materialId);
+      NumericType gamma_F =
+          params.gamma_F.get(Material::fromLegacyId(materialId));
+      NumericType gamma_FO =
+          params.gamma_F_oxidized.get(Material::fromLegacyId(materialId));
       S_eff =
           gamma_F * std::max(1. - phi_F - phi_O - phi_C, 0.) + gamma_FO * phi_O;
     }
@@ -392,8 +394,10 @@ class CF4O2Etchant final
     // F surface coverage on oxidized SiGe
     const auto &phi_C = globalData->getVectorData(2)[primID];
     // Obtain the sticking probability
-    NumericType gamma_F = sticking(materialId);
-    NumericType gamma_FO = stickingOxidized(materialId);
+    NumericType gamma_F =
+        params.gamma_F.get(Material::fromLegacyId(materialId));
+    NumericType gamma_FO =
+        params.gamma_F_oxidized.get(Material::fromLegacyId(materialId));
     NumericType S_eff =
         gamma_F * std::max(1. - phi_F - phi_O - phi_C, 0.) + gamma_FO * phi_O;
 
@@ -405,27 +409,6 @@ class CF4O2Etchant final
   [[nodiscard]] std::vector<std::string> getLocalDataLabels() const override {
     return {"etchantFlux"};
   }
-
-private:
-  NumericType sticking(const int materialId) const {
-    auto material = static_cast<Material>(materialId);
-    auto gamma = params.gamma_F.find(material);
-    if (gamma != params.gamma_F.end())
-      return gamma->second;
-
-    // default value
-    return 1.0;
-  }
-
-  // Sticking probability for oxidized SiGe or Si
-  NumericType stickingOxidized(const int materialId) const {
-    auto material = static_cast<Material>(materialId);
-    auto gamma = params.gamma_F_oxidized.find(material);
-    if (gamma != params.gamma_F_oxidized.end())
-      return gamma->second;
-
-    return 1.0; // Default
-  }
 };
 
 template <typename NumericType, int D>
@@ -448,8 +431,10 @@ class CF4O2Oxygen final
       const auto &phi_F = globalData->getVectorData(0)[primID];
       const auto &phi_O = globalData->getVectorData(1)[primID];
       const auto &phi_C = globalData->getVectorData(2)[primID];
-      NumericType gamma_O = sticking(materialId);
-      NumericType gamma_OC = stickingPassivated(materialId);
+      NumericType gamma_O =
+          params.gamma_O.get(Material::fromLegacyId(materialId));
+      NumericType gamma_OC =
+          params.gamma_O_passivated.get(Material::fromLegacyId(materialId));
       S_eff =
           gamma_O * std::max(1. - phi_O - phi_F - phi_C, 0.) + gamma_OC * phi_C;
     }
@@ -466,8 +451,10 @@ class CF4O2Oxygen final
     const auto &phi_F = globalData->getVectorData(0)[primID];
     const auto &phi_O = globalData->getVectorData(1)[primID];
     const auto &phi_C = globalData->getVectorData(2)[primID];
-    NumericType gamma_O = sticking(materialId);
-    NumericType gamma_OC = stickingPassivated(materialId);
+    NumericType gamma_O =
+        params.gamma_O.get(Material::fromLegacyId(materialId));
+    NumericType gamma_OC =
+        params.gamma_O_passivated.get(Material::fromLegacyId(materialId));
     NumericType S_eff =
         gamma_O * std::max(1. - phi_O - phi_F - phi_C, 0.) + gamma_OC * phi_C;
 
@@ -479,27 +466,6 @@ class CF4O2Oxygen final
   [[nodiscard]] std::vector<std::string> getLocalDataLabels() const override {
     return {"oxygenFlux"};
   }
-
-private:
-  NumericType sticking(const int materialId) const {
-    auto material = static_cast<Material>(materialId);
-    auto gamma = params.gamma_O.find(material);
-    if (gamma != params.gamma_O.end())
-      return gamma->second;
-
-    // default value
-    return 1.0;
-  }
-
-  NumericType stickingPassivated(const int materialId) const {
-    auto material = static_cast<Material>(materialId);
-    auto gamma = params.gamma_O_passivated.find(material);
-    if (gamma != params.gamma_O_passivated.end())
-      return gamma->second;
-
-    // default value
-    return 1.0;
-  }
 };
 
 template <typename NumericType, int D>
@@ -522,8 +488,10 @@ class CF4O2Polymer final
       const auto &phi_F = globalData->getVectorData(0)[primID];
       const auto &phi_O = globalData->getVectorData(1)[primID];
       const auto &phi_C = globalData->getVectorData(2)[primID];
-      NumericType gamma_C = sticking(materialId);
-      NumericType gamma_CO = stickingOxidized(materialId);
+      NumericType gamma_C =
+          params.gamma_C.get(Material::fromLegacyId(materialId));
+      NumericType gamma_CO =
+          params.gamma_C_oxidized.get(Material::fromLegacyId(materialId));
       S_eff =
           gamma_C * std::max(1. - phi_O - phi_F - phi_C, 0.) + gamma_CO * phi_O;
     }
@@ -540,8 +508,10 @@ class CF4O2Polymer final
     const auto &phi_F = globalData->getVectorData(0)[primID];
     const auto &phi_O = globalData->getVectorData(1)[primID];
     const auto &phi_C = globalData->getVectorData(2)[primID];
-    NumericType gamma_C = sticking(materialId);
-    NumericType gamma_CO = stickingOxidized(materialId);
+    NumericType gamma_C =
+        params.gamma_C.get(Material::fromLegacyId(materialId));
+    NumericType gamma_CO =
+        params.gamma_C_oxidized.get(Material::fromLegacyId(materialId));
     NumericType S_eff =
         gamma_C * std::max(1. - phi_O - phi_F - phi_C, 0.) + gamma_CO * phi_O;
 
@@ -553,27 +523,6 @@ class CF4O2Polymer final
   [[nodiscard]] std::vector<std::string> getLocalDataLabels() const override {
     return {"polymerFlux"};
   }
-
-private:
-  NumericType sticking(const int materialId) const {
-    auto material = static_cast<Material>(materialId);
-    auto gamma = params.gamma_C.find(material);
-    if (gamma != params.gamma_C.end())
-      return gamma->second;
-
-    // default value
-    return 1.0;
-  }
-
-  NumericType stickingOxidized(const int materialId) const {
-    auto material = static_cast<Material>(materialId);
-    auto gamma = params.gamma_C_oxidized.find(material);
-    if (gamma != params.gamma_C_oxidized.end())
-      return gamma->second;
-
-    // default value
-    return 1.0;
-  }
 };
 } // namespace impl
 
diff --git a/include/viennaps/models/psCF4O2Parameters.hpp b/include/viennaps/models/psCF4O2Parameters.hpp
index e021c0aa..d3e5a09a 100644
--- a/include/viennaps/models/psCF4O2Parameters.hpp
+++ b/include/viennaps/models/psCF4O2Parameters.hpp
@@ -1,6 +1,7 @@
 #pragma once
 
-#include "../materials/psMaterials.hpp"
+#include "../materials/psMaterialMap.hpp"
+#include "../materials/psMaterialValueMap.hpp"
 #include "../psConstants.hpp"
 #include <unordered_map>
 
@@ -14,18 +15,44 @@ template <typename NumericType> struct CF4O2Parameters {
   NumericType polymerFlux = 1.0e2;
 
   // sticking probabilities
-  std::unordered_map<Material, NumericType> gamma_F = {
-      {Material::Mask, 0.7}, {Material::Si, 0.7}, {Material::SiGe, 0.7}};
-  std::unordered_map<Material, NumericType> gamma_F_oxidized = {
-      {Material::Mask, 0.3}, {Material::Si, 0.3}, {Material::SiGe, 0.3}};
-  std::unordered_map<Material, NumericType> gamma_O = {
-      {Material::Mask, 1.0}, {Material::Si, 1.0}, {Material::SiGe, 1.0}};
-  std::unordered_map<Material, NumericType> gamma_O_passivated = {
-      {Material::Mask, 0.3}, {Material::Si, 0.3}, {Material::SiGe, 0.3}};
-  std::unordered_map<Material, NumericType> gamma_C = {
-      {Material::Mask, 1.0}, {Material::Si, 1.0}, {Material::SiGe, 1.0}};
-  std::unordered_map<Material, NumericType> gamma_C_oxidized = {
-      {Material::Mask, 0.3}, {Material::Si, 0.3}, {Material::SiGe, 0.3}};
+  MaterialValueMap<NumericType> gamma_F = MaterialValueMap<NumericType>(
+      std::initializer_list<std::pair<Material, NumericType>>{
+          {Material::Mask, 0.7}, {Material::Si, 0.7}, {Material::SiGe, 0.7}},
+      1.0);
+
+  MaterialValueMap<NumericType> gamma_F_oxidized =
+      MaterialValueMap<NumericType>(
+          std::initializer_list<std::pair<Material, NumericType>>{
+              {Material::Mask, 0.3},
+              {Material::Si, 0.3},
+              {Material::SiGe, 0.3}},
+          1.0);
+
+  MaterialValueMap<NumericType> gamma_O = MaterialValueMap<NumericType>(
+      std::initializer_list<std::pair<Material, NumericType>>{
+          {Material::Mask, 1.0}, {Material::Si, 1.0}, {Material::SiGe, 1.0}},
+      1.0);
+
+  MaterialValueMap<NumericType> gamma_O_passivated =
+      MaterialValueMap<NumericType>(
+          std::initializer_list<std::pair<Material, NumericType>>{
+              {Material::Mask, 0.3},
+              {Material::Si, 0.3},
+              {Material::SiGe, 0.3}},
+          1.0);
+
+  MaterialValueMap<NumericType> gamma_C = MaterialValueMap<NumericType>(
+      std::initializer_list<std::pair<Material, NumericType>>{
+          {Material::Mask, 1.0}, {Material::Si, 1.0}, {Material::SiGe, 1.0}},
+      1.0);
+
+  MaterialValueMap<NumericType> gamma_C_oxidized =
+      MaterialValueMap<NumericType>(
+          std::initializer_list<std::pair<Material, NumericType>>{
+              {Material::Mask, 0.3},
+              {Material::Si, 0.3},
+              {Material::SiGe, 0.3}},
+          1.0);
 
   NumericType etchStopDepth = std::numeric_limits<NumericType>::lowest();
   bool fluxIncludeSticking = false;
@@ -136,28 +163,28 @@ template <typename NumericType> struct CF4O2Parameters {
 
     // sticking probabilities
     for (const auto &gamma : gamma_F) {
-      metaData["gamma_F_" + MaterialMap::toString(gamma.first)] = {
-          gamma.second};
+      metaData["gamma_F_" + MaterialMap::toString(gamma.material)] = {
+          gamma.value};
     }
     for (const auto &gamma : gamma_F_oxidized) {
-      metaData["gamma_F_oxidized_" + MaterialMap::toString(gamma.first)] = {
-          gamma.second};
+      metaData["gamma_F_oxidized_" + MaterialMap::toString(gamma.material)] = {
+          gamma.value};
     }
     for (const auto &gamma : gamma_O) {
-      metaData["gamma_O_" + MaterialMap::toString(gamma.first)] = {
-          gamma.second};
+      metaData["gamma_O_" + MaterialMap::toString(gamma.material)] = {
+          gamma.value};
     }
     for (const auto &gamma : gamma_O_passivated) {
-      metaData["gamma_O_passivated_" + MaterialMap::toString(gamma.first)] = {
-          gamma.second};
+      metaData["gamma_O_passivated_" + MaterialMap::toString(gamma.material)] =
+          {gamma.value};
     }
     for (const auto &gamma : gamma_C) {
-      metaData["gamma_C_" + MaterialMap::toString(gamma.first)] = {
-          gamma.second};
+      metaData["gamma_C_" + MaterialMap::toString(gamma.material)] = {
+          gamma.value};
     }
     for (const auto &gamma : gamma_C_oxidized) {
-      metaData["gamma_C_oxidized_" + MaterialMap::toString(gamma.first)] = {
-          gamma.second};
+      metaData["gamma_C_oxidized_" + MaterialMap::toString(gamma.material)] = {
+          gamma.value};
     }
 
     return metaData;
diff --git a/include/viennaps/models/psCSVFileProcess.hpp b/include/viennaps/models/psCSVFileProcess.hpp
index 790295af..9f2ae779 100644
--- a/include/viennaps/models/psCSVFileProcess.hpp
+++ b/include/viennaps/models/psCSVFileProcess.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include "../materials/psMaterials.hpp"
+#include "../materials/psMaterialMap.hpp"
 #include "../process/psProcessModel.hpp"
 #include "../process/psSurfaceModel.hpp"
 #include "../psDomain.hpp"
diff --git a/include/viennaps/models/psDirectionalProcess.hpp b/include/viennaps/models/psDirectionalProcess.hpp
index 7894f1f3..c35cefc7 100644
--- a/include/viennaps/models/psDirectionalProcess.hpp
+++ b/include/viennaps/models/psDirectionalProcess.hpp
@@ -173,6 +173,11 @@ class DirectionalProcess : public ProcessModelCPU<NumericType, D> {
   using RateSet = impl::RateSet<NumericType>;
   using RateMap = MaterialValueMap<std::pair<NumericType, NumericType>>;
 
+  DirectionalProcess(const Vec3D<NumericType> &direction,
+                     RateMap materialRates) {
+    initialize(std::move(materialRates), direction);
+  }
+
   DirectionalProcess(
       const Vec3D<NumericType> &direction,
       std::unordered_map<Material, std::pair<NumericType, NumericType>>
diff --git a/include/viennaps/models/psFaradayCageEtching.hpp b/include/viennaps/models/psFaradayCageEtching.hpp
index 2dff5220..a4b03bcb 100644
--- a/include/viennaps/models/psFaradayCageEtching.hpp
+++ b/include/viennaps/models/psFaradayCageEtching.hpp
@@ -3,7 +3,7 @@
 #include "psIonBeamEtching.hpp"
 #include "psPipelineParameters.hpp"
 
-#include "../materials/psMaterials.hpp"
+#include "../materials/psMaterialMap.hpp"
 #include "../process/psProcessModel.hpp"
 #include "../psUtil.hpp"
 
@@ -257,23 +257,18 @@ class FaradayCageEtching final : public ProcessModelGPU<NumericType, D> {
 
     // surface model
     // adjust rate here since we trace two particles per point
-    auto surfaceModelParams = params_.ibeParams;
-    surfaceModelParams.planeWaferRate *= 0.5;
-    surfaceModelParams.materialPlaneWaferRate.setDefault(
-        surfaceModelParams.planeWaferRate);
-    for (auto entry : surfaceModelParams.materialPlaneWaferRate) {
-      surfaceModelParams.materialPlaneWaferRate.set(entry.getMaterial(),
-                                                    entry.getValue() * 0.5);
+    surfaceModelParams_ = params_.ibeParams;
+    surfaceModelParams_.planeWaferRate *= 0.5;
+    surfaceModelParams_.materialPlaneWaferRate.setDefault(
+        surfaceModelParams_.planeWaferRate);
+    for (auto entry : params_.ibeParams.materialPlaneWaferRate) {
+      surfaceModelParams_.materialPlaneWaferRate.set(entry.material,
+                                                     entry.value * 0.5);
     }
-    auto &customMaterialRates =
-        surfaceModelParams.materialPlaneWaferRate.customMaterialRates();
-    for (auto &mat : customMaterialRates) {
-      mat.second *= 0.5;
-    }
-    surfaceModelParams.redepositionRate *= 0.5;
+    surfaceModelParams_.redepositionRate *= 0.5;
     auto surfModel =
         SmartPointer<::viennaps::impl::IBESurfaceModel<NumericType>>::New(
-            surfaceModelParams);
+            surfaceModelParams_);
 
     // velocity field
     auto velField = SmartPointer<DefaultVelocityField<NumericType, D>>::New();
@@ -289,6 +284,7 @@ class FaradayCageEtching final : public ProcessModelGPU<NumericType, D> {
 
 private:
   FaradayCageParameters<NumericType> params_;
+  IBEParameters<NumericType> surfaceModelParams_;
 };
 } // namespace gpu
 #endif
diff --git a/include/viennaps/models/psFluorocarbonEtching.hpp b/include/viennaps/models/psFluorocarbonEtching.hpp
index 9f3a5b73..6b83d82e 100644
--- a/include/viennaps/models/psFluorocarbonEtching.hpp
+++ b/include/viennaps/models/psFluorocarbonEtching.hpp
@@ -2,7 +2,7 @@
 
 #include <cmath>
 
-#include "../materials/psMaterials.hpp"
+#include "../materials/psMaterialMap.hpp"
 #include "../process/psProcessModel.hpp"
 #include "../psConstants.hpp"
 #include "../psUnits.hpp"
diff --git a/include/viennaps/models/psHBrO2Etching.hpp b/include/viennaps/models/psHBrO2Etching.hpp
index a410845e..9edcc752 100644
--- a/include/viennaps/models/psHBrO2Etching.hpp
+++ b/include/viennaps/models/psHBrO2Etching.hpp
@@ -44,11 +44,7 @@ class HBrO2Etching final : public ProcessModelGPU<NumericType, D> {
     etchant.dataLabels.push_back("etchantFlux");
     etchant.cosineExponent = 1.f;
     for (auto entry : params.beta_E) {
-      etchant.materialSticking[static_cast<int>(entry.getMaterial())] =
-          entry.getValue();
-    }
-    for (auto &entry : params.beta_E.getCustomMaterialValues()) {
-      etchant.materialSticking[entry.first] = entry.second;
+      etchant.materialSticking[static_cast<int>(entry.material)] = entry.value;
     }
 
     viennaray::gpu::Particle<NumericType> oxygen;
@@ -56,11 +52,7 @@ class HBrO2Etching final : public ProcessModelGPU<NumericType, D> {
     oxygen.dataLabels.push_back("passivationFlux");
     oxygen.cosineExponent = 1.f;
     for (auto entry : params.beta_P) {
-      oxygen.materialSticking[static_cast<int>(entry.getMaterial())] =
-          entry.getValue();
-    }
-    for (auto &entry : params.beta_P.getCustomMaterialValues()) {
-      oxygen.materialSticking[entry.first] = entry.second;
+      oxygen.materialSticking[static_cast<int>(entry.material)] = entry.value;
     }
 
     // surface model
@@ -190,13 +182,9 @@ class HBrO2Etching : public ProcessModelCPU<NumericType, D> {
     // sticking probabilities
     defParams.beta_E.set(Material::Si, 0.1);
     defParams.beta_E.set(Material::Mask, 0.1);
-    defParams.beta_E.setDefault(1.0);
 
     defParams.beta_P.set(Material::Si, 1.);
     defParams.beta_P.set(Material::Mask, 1.);
-    defParams.beta_P.setDefault(1.0);
-
-    defParams.rateFactors.setDefault(1.0);
 
     defParams.etchStopDepth = std::numeric_limits<NumericType>::lowest();
 
diff --git a/include/viennaps/models/psIonBeamEtching.hpp b/include/viennaps/models/psIonBeamEtching.hpp
index f42ab96d..9738a9cc 100644
--- a/include/viennaps/models/psIonBeamEtching.hpp
+++ b/include/viennaps/models/psIonBeamEtching.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include "../materials/psMaterials.hpp"
+#include "../materials/psMaterialMap.hpp"
 #include "../process/psProcessModel.hpp"
 #include "../psConstants.hpp"
 #include "psIonBeamParameters.hpp"
diff --git a/include/viennaps/models/psIonBeamParameters.hpp b/include/viennaps/models/psIonBeamParameters.hpp
index 85e7420c..f4eea13a 100644
--- a/include/viennaps/models/psIonBeamParameters.hpp
+++ b/include/viennaps/models/psIonBeamParameters.hpp
@@ -1,7 +1,7 @@
 #pragma once
 
+#include "../materials/psMaterialMap.hpp"
 #include "../materials/psMaterialValueMap.hpp"
-#include "../materials/psMaterials.hpp"
 
 #include <functional>
 #include <unordered_map>
@@ -51,8 +51,8 @@ template <typename NumericType> struct IBEParameters {
 
     processData["Default PWR"] = {planeWaferRate};
     for (const auto &entry : materialPlaneWaferRate) {
-      processData[MaterialMap::toString(entry.getMaterial()) + " PWR"] = {
-          entry.getValue()};
+      processData[MaterialMap::toString(entry.material) + " PWR"] = {
+          entry.value};
     }
 
     processData["Mean Energy"] = {meanEnergy};
diff --git a/include/viennaps/models/psIsotropicProcess.hpp b/include/viennaps/models/psIsotropicProcess.hpp
index cd77fc90..b750bc1b 100644
--- a/include/viennaps/models/psIsotropicProcess.hpp
+++ b/include/viennaps/models/psIsotropicProcess.hpp
@@ -73,9 +73,9 @@ class IsotropicProcess : public ProcessModelCPU<NumericType, D> {
 
     // store process data (before moving materialRates)
     addMetaData("Isotropic Rate", materialRates.getDefault());
-    for (const auto &materialRate : materialRates) {
-      addMetaData("Rate " + MaterialMap::toString(materialRate.getMaterial()),
-                  materialRate.getValue());
+    for (auto materialRate : materialRates) {
+      addMetaData("Rate " + MaterialMap::toString(materialRate.material),
+                  materialRate.value);
     }
 
     // velocity field
diff --git a/include/viennaps/models/psMultiParticleProcess.hpp b/include/viennaps/models/psMultiParticleProcess.hpp
index 298c812d..192bdf76 100644
--- a/include/viennaps/models/psMultiParticleProcess.hpp
+++ b/include/viennaps/models/psMultiParticleProcess.hpp
@@ -1,6 +1,7 @@
 #pragma once
 
-#include "../materials/psMaterials.hpp"
+#include "../materials/psMaterialMap.hpp"
+#include "../materials/psMaterialValueMap.hpp"
 #include "../process/psProcessModel.hpp"
 #include "../psConstants.hpp"
 #include "../psUtil.hpp"
@@ -170,14 +171,12 @@ template <typename NumericType, int D>
 class DiffuseParticle
     : public viennaray::Particle<DiffuseParticle<NumericType, D>, NumericType> {
 public:
-  DiffuseParticle(NumericType stickingProbability, const std::string &dataLabel)
-      : stickingProbability_(stickingProbability), dataLabel_(dataLabel) {}
-
-  DiffuseParticle(NumericType stickingProbability,
-                  std::unordered_map<Material, NumericType> materialSticking,
-                  const std::string &dataLabel)
-      : materialSticking_(materialSticking),
-        stickingProbability_(stickingProbability), dataLabel_(dataLabel) {}
+  DiffuseParticle(
+      NumericType defaultStickingProbability,
+      const std::unordered_map<Material, NumericType> &materialSticking,
+      const std::string &dataLabel)
+      : materialSticking_(materialSticking, defaultStickingProbability),
+        dataLabel_(dataLabel) {}
 
   void surfaceCollision(NumericType rayWeight, const Vec3D<NumericType> &rayDir,
                         const Vec3D<NumericType> &geomNormal,
@@ -193,13 +192,7 @@ class DiffuseParticle
                     const int materialId,
                     const viennaray::TracingData<NumericType> *,
                     RNG &rngState) override final {
-    NumericType sticking = stickingProbability_;
-    if (auto mat =
-            materialSticking_.find(MaterialMap::mapToMaterial(materialId));
-        mat != materialSticking_.end()) {
-      sticking = mat->second;
-    }
-
+    auto sticking = materialSticking_.get(Material::fromLegacyId(materialId));
     auto direction =
         viennaray::ReflectionDiffuse<NumericType, D>(geomNormal, rngState);
     return std::pair<NumericType, Vec3D<NumericType>>{sticking, direction};
@@ -209,8 +202,7 @@ class DiffuseParticle
   }
 
 public:
-  const std::unordered_map<Material, NumericType> materialSticking_;
-  const NumericType stickingProbability_;
+  const MaterialValueMap<NumericType> materialSticking_;
   const std::string dataLabel_;
 };
 } // namespace impl
@@ -310,11 +302,11 @@ class MultiParticleProcess final : public ProcessModelGPU<NumericType, D> {
         if (!neutralAdded[i] && label == neutrals[i].dataLabel_) {
           const auto &neutral = neutrals[i];
           if (neutral.materialSticking_.empty()) {
-            addNeutralParticle(neutral.stickingProbability_,
+            addNeutralParticle(neutral.materialSticking_.getDefault(),
                                neutral.dataLabel_);
           } else {
             addNeutralParticle(neutral.materialSticking_,
-                               neutral.stickingProbability_,
+                               neutral.materialSticking_.getDefault(),
                                neutral.dataLabel_);
           }
           neutralAdded[i] = true;
@@ -345,10 +337,12 @@ class MultiParticleProcess final : public ProcessModelGPU<NumericType, D> {
           "conversion. Appending neutral particle with label: " +
           neutral.dataLabel_);
       if (neutral.materialSticking_.empty()) {
-        addNeutralParticle(neutral.stickingProbability_, neutral.dataLabel_);
+        addNeutralParticle(neutral.materialSticking_.getDefault(),
+                           neutral.dataLabel_);
       } else {
         addNeutralParticle(neutral.materialSticking_,
-                           neutral.stickingProbability_, neutral.dataLabel_);
+                           neutral.materialSticking_.getDefault(),
+                           neutral.dataLabel_);
       }
     }
 
@@ -378,10 +372,9 @@ class MultiParticleProcess final : public ProcessModelGPU<NumericType, D> {
                          std::to_string(stickingProbability));
   }
 
-  void addNeutralParticle(
-      const std::unordered_map<Material, NumericType> &materialSticking,
-      NumericType defaultStickingProbability = 1.,
-      const std::string &label = "neutralFlux") {
+  void addNeutralParticle(const MaterialValueMap<NumericType> &materialSticking,
+                          NumericType defaultStickingProbability = 1.,
+                          const std::string &label = "neutralFlux") {
     std::string dataLabel = label + std::to_string(fluxDataLabels_.size());
     fluxDataLabels_.push_back(dataLabel);
 
@@ -389,8 +382,8 @@ class MultiParticleProcess final : public ProcessModelGPU<NumericType, D> {
     particle.name = "Neutral";
     particle.sticking = defaultStickingProbability;
     particle.dataLabels.push_back(dataLabel);
-    for (auto &mat : materialSticking) {
-      particle.materialSticking[static_cast<int>(mat.first)] = mat.second;
+    for (auto entry : materialSticking) {
+      particle.materialSticking[static_cast<int>(entry.material)] = entry.value;
     }
 
     this->insertNextParticleType(particle);
@@ -535,12 +528,16 @@ class MultiParticleProcess : public ProcessModelCPU<NumericType, D> {
                           const std::string &label = "neutralFlux") {
     std::string dataLabel = label + std::to_string(fluxDataLabels_.size());
     fluxDataLabels_.push_back(dataLabel);
+    std::unordered_map<Material, NumericType>
+        materialSticking; // empty map will use default sticking for all
+                          // materials
     auto particle = std::make_unique<impl::DiffuseParticle<NumericType, D>>(
-        stickingProbability, dataLabel);
+        stickingProbability, materialSticking, dataLabel);
     this->insertNextParticleType(particle);
 
     addStickingData(stickingProbability);
-    neutralParticles_.emplace_back(stickingProbability, dataLabel);
+    neutralParticles_.emplace_back(stickingProbability, materialSticking,
+                                   dataLabel);
   }
 
   void addNeutralParticle(
diff --git a/include/viennaps/models/psPlasmaEtchingParameters.hpp b/include/viennaps/models/psPlasmaEtchingParameters.hpp
index 20634f95..1ef043d5 100644
--- a/include/viennaps/models/psPlasmaEtchingParameters.hpp
+++ b/include/viennaps/models/psPlasmaEtchingParameters.hpp
@@ -14,10 +14,14 @@ template <typename NumericType> struct PlasmaEtchingParameters {
   NumericType passivationFlux = 1.0e2;
 
   // sticking probabilities
-  MaterialValueMap<NumericType> beta_E;
-  MaterialValueMap<NumericType> beta_P;
-
-  MaterialValueMap<NumericType> rateFactors;
+  MaterialValueMap<NumericType> beta_E =
+      MaterialValueMap<NumericType>::fromDefault(1.0);
+  MaterialValueMap<NumericType> beta_P =
+      MaterialValueMap<NumericType>::fromDefault(1.0);
+
+  MaterialValueMap<NumericType> rateFactors =
+      MaterialValueMap<NumericType>::fromDefault(
+          1.0); // default to 1.0 for all materials
   NumericType etchStopDepth = std::numeric_limits<NumericType>::lowest();
 
   // Mask
@@ -101,12 +105,12 @@ template <typename NumericType> struct PlasmaEtchingParameters {
     processData["Passivation Flux"] = {passivationFlux};
 
     for (auto entry : beta_E) {
-      processData["Beta_E " + MaterialMap::toString(entry.getMaterial())] =
-          std::vector<double>{entry.getValue()};
+      processData["Beta_E " + MaterialMap::toString(entry.material)] =
+          std::vector<double>{entry.value};
     }
     for (auto entry : beta_P) {
-      processData["Beta_P " + MaterialMap::toString(entry.getMaterial())] =
-          std::vector<double>{entry.getValue()};
+      processData["Beta_P " + MaterialMap::toString(entry.material)] =
+          std::vector<double>{entry.value};
     }
 
     if (etchStopDepth != std::numeric_limits<NumericType>::lowest())
diff --git a/include/viennaps/models/psSF6C4F8Etching.hpp b/include/viennaps/models/psSF6C4F8Etching.hpp
index 0272b3ec..ebd8d4ca 100644
--- a/include/viennaps/models/psSF6C4F8Etching.hpp
+++ b/include/viennaps/models/psSF6C4F8Etching.hpp
@@ -43,11 +43,7 @@ class SF6C4F8Etching final : public ProcessModelGPU<NumericType, D> {
     etchant.dataLabels.push_back("etchantFlux");
     etchant.cosineExponent = 1.f;
     for (auto entry : params.beta_E) {
-      etchant.materialSticking[static_cast<int>(entry.getMaterial())] =
-          entry.getValue();
-    }
-    for (auto &entry : params.beta_E.getCustomMaterialValues()) {
-      etchant.materialSticking[entry.first] = entry.second;
+      etchant.materialSticking[static_cast<int>(entry.material)] = entry.value;
     }
 
     // No oxygen/passivation particle for SF6C4F8
@@ -175,11 +171,8 @@ class SF6C4F8Etching : public ProcessModelCPU<NumericType, D> {
     defParams.beta_E.set(Material::Si, 0.7);
     defParams.beta_E.set(Material::Mask, 0.02);
     defParams.beta_E.set(Material::Polymer, 0.02);
-    defParams.beta_E.setDefault(1.0);
     // No beta_P needed since passivationFlux = 0
 
-    defParams.rateFactors.setDefault(1.0);
-
     defParams.etchStopDepth = std::numeric_limits<NumericType>::lowest();
 
     // Mask
diff --git a/include/viennaps/models/psSF6O2Etching.hpp b/include/viennaps/models/psSF6O2Etching.hpp
index ac8a91a8..a2ffd93a 100644
--- a/include/viennaps/models/psSF6O2Etching.hpp
+++ b/include/viennaps/models/psSF6O2Etching.hpp
@@ -43,11 +43,7 @@ class SF6O2Etching final : public ProcessModelGPU<NumericType, D> {
     etchant.dataLabels.push_back("etchantFlux");
     etchant.cosineExponent = 1.f;
     for (auto entry : params.beta_E) {
-      etchant.materialSticking[static_cast<int>(entry.getMaterial())] =
-          entry.getValue();
-    }
-    for (auto &entry : params.beta_E.getCustomMaterialValues()) {
-      etchant.materialSticking[entry.first] = entry.second;
+      etchant.materialSticking[static_cast<int>(entry.material)] = entry.value;
     }
 
     viennaray::gpu::Particle<NumericType> oxygen;
@@ -55,11 +51,7 @@ class SF6O2Etching final : public ProcessModelGPU<NumericType, D> {
     oxygen.dataLabels.push_back("passivationFlux");
     oxygen.cosineExponent = 1.f;
     for (auto entry : params.beta_P) {
-      oxygen.materialSticking[static_cast<int>(entry.getMaterial())] =
-          entry.getValue();
-    }
-    for (auto &entry : params.beta_P.getCustomMaterialValues()) {
-      oxygen.materialSticking[entry.first] = entry.second;
+      oxygen.materialSticking[static_cast<int>(entry.material)] = entry.value;
     }
 
     // surface model
@@ -190,13 +182,9 @@ class SF6O2Etching : public ProcessModelCPU<NumericType, D> {
     // sticking probabilities
     defParams.beta_E.set(Material::Si, 0.7);
     defParams.beta_E.set(Material::Mask, 0.7);
-    defParams.beta_E.setDefault(1.0);
 
     defParams.beta_P.set(Material::Si, 1.);
     defParams.beta_P.set(Material::Mask, 1.);
-    defParams.beta_P.setDefault(1.0);
-
-    defParams.rateFactors.setDefault(1.0);
 
     defParams.etchStopDepth = std::numeric_limits<NumericType>::lowest();
 
diff --git a/include/viennaps/models/psSelectiveEpitaxy.hpp b/include/viennaps/models/psSelectiveEpitaxy.hpp
index 3e617744..336bc91a 100644
--- a/include/viennaps/models/psSelectiveEpitaxy.hpp
+++ b/include/viennaps/models/psSelectiveEpitaxy.hpp
@@ -1,6 +1,7 @@
 #pragma once
 
-#include "../materials/psMaterials.hpp"
+#include "../materials/psMaterialMap.hpp"
+#include "../materials/psMaterialValueMap.hpp"
 #include "../process/psProcessModel.hpp"
 
 namespace viennaps {
@@ -17,30 +18,26 @@ class EpitaxyVelocityField : public VelocityField<NumericType, D> {
       (D > 2) ? 0.5773502691896257 : 0.7071067811865476;
   static constexpr double high = 1.0;
   const double factor;
+  const Vec3D<NumericType> nvFactors;
 
-  const std::vector<std::pair<Material, NumericType>> &materials;
+  const MaterialValueMap<NumericType> &materialRates;
 
 public:
-  EpitaxyVelocityField(
-      const std::vector<std::pair<Material, NumericType>> &passedmaterials,
-      NumericType r111, NumericType r100)
+  EpitaxyVelocityField(const MaterialValueMap<NumericType> &materials,
+                       NumericType r111, NumericType r100,
+                       const Vec3D<NumericType> &rates)
       : R111(r111), R100(r100), factor((R100 - R111) / (high - low)),
-        materials(passedmaterials) {}
-
-  NumericType getScalarVelocity(const Vec3D<NumericType> &coordinate,
-                                int material, const Vec3D<NumericType> &nv,
-                                unsigned long pointID) override {
-    for (auto const &epitaxyMaterial : materials) {
-      if (MaterialMap::isMaterial(material, epitaxyMaterial.first)) {
-        double vel = std::max(std::abs(nv[0]), std::abs(nv[D - 1]));
-        vel = (vel - low) * factor + R111;
-
-        if (std::abs(nv[0]) < std::abs(nv[D - 1])) {
-          vel *= 2.;
-        }
+        materialRates(materials), nvFactors(rates) {}
 
-        return -vel * epitaxyMaterial.second;
-      }
+  NumericType getScalarVelocity(const Vec3D<NumericType> &, int material,
+                                const Vec3D<NumericType> &nv,
+                                unsigned long) override {
+
+    auto rate = materialRates.get(Material::fromLegacyId(material));
+    if (rate > 0) {
+      double vel = MaxElement(Abs(nvFactors * nv));
+      vel = (vel - low) * factor + R111;
+      return std::min(-vel * rate, 0.0);
     }
 
     // not an epitaxy material
@@ -53,23 +50,35 @@ class EpitaxyVelocityField : public VelocityField<NumericType, D> {
 template <typename NumericType, int D>
 class SelectiveEpitaxy : public ProcessModelCPU<NumericType, D> {
 public:
+  SelectiveEpitaxy(NumericType rate111 = 0.5, NumericType rate100 = 1.)
+      : SelectiveEpitaxy(std::vector<std::pair<Material, NumericType>>{},
+                         rate111, rate100) {}
+
   // The constructor expects the materials where epitaxy is allowed including
   // the corresponding rates.
   SelectiveEpitaxy(
       const std::vector<std::pair<Material, NumericType>> pMaterials,
-      NumericType rate111 = 0.5, NumericType rate100 = 1.)
-      : materials(pMaterials) {
-    // default surface model
-    auto surfModel = SmartPointer<SurfaceModel<NumericType>>::New();
+      NumericType rate111 = 0.5, NumericType rate100 = 1.) {
+    for (const auto &[material, rate] : pMaterials) {
+      materialRates_.set(material, rate);
+    }
 
-    // velocity field
-    auto velField =
-        SmartPointer<impl::EpitaxyVelocityField<NumericType, D>>::New(
-            materials, rate111, rate100);
+    // default nvFactors
+    nvFactors_.fill(0.);
+    nvFactors_[0] = 0.5;
+    nvFactors_[D - 1] = 1.0;
 
-    this->setSurfaceModel(surfModel);
-    this->setVelocityField(velField);
-    this->setProcessName("SelectiveEpitaxy");
+    setup(rate111, rate100);
+  }
+
+  SelectiveEpitaxy(const Vec3D<NumericType> &nvFactors,
+                   NumericType rate111 = 0.5, NumericType rate100 = 1.)
+      : nvFactors_(nvFactors) {
+    setup(rate111, rate100);
+  }
+
+  void setMaterialRate(Material material, NumericType rate) {
+    materialRates_.set(material, rate);
   }
 
   void initialize(SmartPointer<Domain<NumericType, D>> domain,
@@ -107,12 +116,6 @@ class SelectiveEpitaxy : public ProcessModelCPU<NumericType, D> {
               maskLayer, lsCopy,
               viennals::BooleanOperationEnum::RELATIVE_COMPLEMENT)
               .apply();
-
-          // auto mesh = viennals::Mesh<NumericType>::New();
-          // viennals::ToSurfaceMesh(lsCopy, mesh).apply();
-          // viennals::VTKWriter<NumericType>(mesh, "SelectiveEpitaxyMask_" +
-          //                                            std::to_string(i))
-          //     .apply();
         }
       }
 
@@ -138,16 +141,27 @@ class SelectiveEpitaxy : public ProcessModelCPU<NumericType, D> {
   }
 
 private:
-  std::vector<std::pair<Material, NumericType>> materials;
+  MaterialValueMap<NumericType> materialRates_;
+  Vec3D<NumericType> nvFactors_;
   SmartPointer<Domain<NumericType, D>> domainCopy;
   bool firstInit = true;
 
+  void setup(const NumericType rate111, const NumericType rate100) {
+    // default surface model
+    auto surfModel = SmartPointer<SurfaceModel<NumericType>>::New();
+
+    // velocity field
+    auto velField =
+        SmartPointer<impl::EpitaxyVelocityField<NumericType, D>>::New(
+            materialRates_, rate111, rate100, nvFactors_);
+
+    this->setSurfaceModel(surfModel);
+    this->setVelocityField(velField);
+    this->setProcessName("SelectiveEpitaxy");
+  }
+
   bool isEpitaxyMaterial(const Material &material) const {
-    for (const auto &epitaxyMaterial : materials) {
-      if (MaterialMap::isMaterial(material, epitaxyMaterial.first))
-        return true;
-    }
-    return false;
+    return materialRates_.get(material) != NumericType(0);
   }
 };
 
diff --git a/include/viennaps/models/psSingleParticleProcess.hpp b/include/viennaps/models/psSingleParticleProcess.hpp
index ba6c34e4..be2a23e9 100644
--- a/include/viennaps/models/psSingleParticleProcess.hpp
+++ b/include/viennaps/models/psSingleParticleProcess.hpp
@@ -129,8 +129,8 @@ class SingleParticleProcess final : public ProcessModelGPU<NumericType, D> {
     this->processMetaData["Source Exponent"] =
         std::vector<double>{sourceDistributionPower_};
     for (auto rate : materialRates_) {
-      this->processMetaData[MaterialMap::toString(rate.getMaterial()) +
-                            " Rate"] = std::vector<double>{rate.getValue()};
+      this->processMetaData[MaterialMap::toString(rate.material) + " Rate"] =
+          std::vector<double>{rate.value};
     }
   }
 };
@@ -226,8 +226,8 @@ class SingleParticleProcess : public ProcessModelCPU<NumericType, D> {
     this->processMetaData["Source Exponent"] =
         std::vector<double>{sourceDistributionPower_};
     for (auto mrate : materialRates_) {
-      this->processMetaData[MaterialMap::toString(mrate.getMaterial()) +
-                            " Rate"] = std::vector<double>{mrate.getValue()};
+      this->processMetaData[MaterialMap::toString(mrate.material) + " Rate"] =
+          std::vector<double>{mrate.value};
     }
   }
 };
diff --git a/include/viennaps/models/psWetEtching.hpp b/include/viennaps/models/psWetEtching.hpp
index 33528fdc..d25b9945 100644
--- a/include/viennaps/models/psWetEtching.hpp
+++ b/include/viennaps/models/psWetEtching.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include "../materials/psMaterials.hpp"
+#include "../materials/psMaterialMap.hpp"
 #include "../process/psProcessModel.hpp"
 
 #include <vcVectorType.hpp>
diff --git a/include/viennaps/process/psTranslationField.hpp b/include/viennaps/process/psTranslationField.hpp
index a0219c81..1e777153 100644
--- a/include/viennaps/process/psTranslationField.hpp
+++ b/include/viennaps/process/psTranslationField.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include "../materials/psMaterials.hpp"
+#include "../materials/psMaterialMap.hpp"
 #include "psVelocityField.hpp"
 
 #include <lsVelocityField.hpp>
diff --git a/include/viennaps/psDomain.hpp b/include/viennaps/psDomain.hpp
index e7430bf0..d98f7e93 100644
--- a/include/viennaps/psDomain.hpp
+++ b/include/viennaps/psDomain.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include "materials/psMaterials.hpp"
+#include "materials/psMaterialMap.hpp"
 #include "psDomainSetup.hpp"
 #include "psPreCompileMacros.hpp"
 #include "psSurfacePointValuesToLevelSet.hpp"
diff --git a/include/viennaps/psUtil.hpp b/include/viennaps/psUtil.hpp
index 10472f8d..7a9a7d42 100644
--- a/include/viennaps/psUtil.hpp
+++ b/include/viennaps/psUtil.hpp
@@ -10,7 +10,7 @@
 #include <string>
 #include <unordered_map>
 
-#include "materials/psMaterials.hpp"
+#include "materials/psMaterialMap.hpp"
 
 namespace viennaps {
 enum class FluxEngineType {
diff --git a/python/pyWrapDimension.hpp b/python/pyWrapDimension.hpp
index f8e15b37..3c1ad6b5 100644
--- a/python/pyWrapDimension.hpp
+++ b/python/pyWrapDimension.hpp
@@ -1106,10 +1106,18 @@ template <int D> void bindApi(py::module &module) {
   // Selective Epitaxy Process
   py::class_<SelectiveEpitaxy<T, D>, SmartPointer<SelectiveEpitaxy<T, D>>>(
       module, "SelectiveEpitaxy", processModel)
+      .def(py::init(&SmartPointer<SelectiveEpitaxy<T, D>>::template New<T, T>),
+           py::arg("rate111") = 0.5, py::arg("rate100") = 1.0)
       .def(py::init(&SmartPointer<SelectiveEpitaxy<T, D>>::template New<
                     const std::vector<std::pair<Material, T>>, T, T>),
            py::arg("materialRates"), py::arg("rate111") = 0.5,
-           py::arg("rate100") = 1.0);
+           py::arg("rate100") = 1.0)
+      .def(py::init(&SmartPointer<SelectiveEpitaxy<T, D>>::template New<
+                    const Vec3D<T> &, T, T>),
+           py::arg("nvFactors"), py::arg("rate111") = 0.5,
+           py::arg("rate100") = 1.0)
+      .def("setMaterialRate", &SelectiveEpitaxy<T, D>::setMaterialRate,
+           py::arg("material"), py::arg("rate"));
 
   // Single Particle ALD
   py::class_<SingleParticleALD<T, D>, SmartPointer<SingleParticleALD<T, D>>>(
@@ -1541,7 +1549,7 @@ template <int D> void bindApi(py::module &module) {
                &gpu::MultiParticleProcess<T, D>::addNeutralParticle),
            py::arg("stickingProbability"), py::arg("label") = "neutralFlux")
       .def("addNeutralParticle",
-           py::overload_cast<const std::unordered_map<Material, T> &, T,
+           py::overload_cast<const MaterialValueMap<T> &, T,
                              const std::string &>(
                &gpu::MultiParticleProcess<T, D>::addNeutralParticle),
            py::arg("materialSticking"),
diff --git a/python/scripts/run_examples_from_list.py b/python/scripts/run_examples_from_list.py
new file mode 100644
index 00000000..e3f86010
--- /dev/null
+++ b/python/scripts/run_examples_from_list.py
@@ -0,0 +1,181 @@
+import argparse
+import os
+import subprocess
+import sys
+from typing import Dict, List
+
+
+# List of examples to run
+EXAMPLES: List[Dict[str, object]] = [
+    {
+        "name": "Trench_Deposition",
+        "dir": "examples/trenchDeposition",
+        "script": "trenchDeposition.py",
+        "args": ["-D", "2", "config.txt"],
+    },
+    {
+        "name": "Hole_Etching",
+        "dir": "examples/holeEtching",
+        "script": "holeEtching.py",
+        "args": ["-D", "2", "config.txt"],
+    },
+    {
+        "name": "Blazed_Gratings_Etching",
+        "dir": "examples/blazedGratingsEtching",
+        "script": "blazedGratingsEtching.py",
+        "args": ["config.txt"],
+    },
+    {
+        "name": "Bosch_simulate",
+        "dir": "examples/boschProcess",
+        "script": "boschProcessSimulate.py",
+        "args": ["-D", "2", "config.txt"],
+    },
+    {
+        "name": "Bosch_emulate",
+        "dir": "examples/boschProcess",
+        "script": "boschProcessEmulate.py",
+        "args": ["-D", "2", "config.txt"],
+    },
+    {
+        "name": "Ion_Beam_Etching",
+        "dir": "examples/ionBeamEtching",
+        "script": "ionBeamEtching.py",
+        "args": ["-D", "2", "config.txt"],
+    },
+    {
+        "name": "Faraday_Cage_Etching",
+        "dir": "examples/faradayCageEtching",
+        "script": "faradayCageEtching.py",
+        "args": ["-D", "2", "config.txt"],
+    },
+    {
+        "name": "Selective_Epitaxy",
+        "dir": "examples/selectiveEpitaxy",
+        "script": "selectiveEpitaxy.py",
+        "args": ["-D", "3", "config.txt"],
+    },
+    {
+        "name": "Simple_Etching",
+        "dir": "examples/simpleEtching",
+        "script": "simpleEtching.py",
+        "args": [],
+    },
+    {
+        "name": "Stack_Etching",
+        "dir": "examples/stackEtching",
+        "script": "stackEtching.py",
+        "args": ["config.txt"],
+    },
+    {
+        "name": "Single_TEOS",
+        "dir": "examples/TEOSTrenchDeposition",
+        "script": "singleTEOS.py",
+        "args": ["-D", "2", "singleTEOS_config.txt"],
+    },
+    {
+        "name": "Multi_TEOS",
+        "dir": "examples/TEOSTrenchDeposition",
+        "script": "multiTEOS.py",
+        "args": ["-D", "2", "multiTEOS_config.txt"],
+    },
+    {
+        "name": "Sputter_Deposition",
+        "dir": "examples/sputterDeposition",
+        "script": "sputterDeposition.py",
+        "args": ["-D", "3", "config3D.txt"],
+    },
+    {
+        "name": "SiGe_Selective_Etching",
+        "dir": "examples/SiGeSelectiveEtching",
+        "script": "SiGeEtching.py",
+        "args": [],
+    },
+    {
+        "name": "Cantilever_Wet_Etching",
+        "dir": "examples/cantileverWetEtching",
+        "script": "cantileverWetEtching.py",
+        "args": [],
+    },
+]
+
+
+def run_example(source_dir: str, example: Dict[str, object], python_exe: str) -> bool:
+    example_name = str(example["name"])
+    rel_dir = str(example["dir"])
+    script = str(example["script"])
+    args = [str(arg) for arg in example.get("args", [])]
+
+    working_dir = os.path.join(source_dir, rel_dir)
+    script_path = os.path.join(working_dir, script)
+
+    if not os.path.isdir(working_dir):
+        print(f"[FAIL] {example_name}: example directory not found: {working_dir}")
+        return False
+
+    if not os.path.isfile(script_path):
+        print(f"[FAIL] {example_name}: script not found: {script_path}")
+        return False
+
+    cmd = [python_exe, script] + args
+    print(f"[RUN ] {example_name}: {' '.join(cmd)}")
+
+    try:
+        subprocess.run(cmd, cwd=working_dir, check=True)
+        print(f"[ OK ] {example_name}")
+        return True
+    except subprocess.CalledProcessError as error:
+        print(f"[FAIL] {example_name}: exited with code {error.returncode}")
+        return False
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(
+        description="Run configured ViennaPS examples from a source folder"
+    )
+    parser.add_argument(
+        "source_dir",
+        nargs="?",
+        default=".",
+        help="Path to ViennaPS source folder containing .venv and examples/",
+    )
+    parser.add_argument(
+        "--keep-going",
+        action="store_true",
+        help="Continue with remaining examples after a failure",
+    )
+    args = parser.parse_args()
+
+    source_dir = os.path.abspath(args.source_dir)
+    python_exe = os.path.join(source_dir, ".venv", "bin", "python")
+
+    if not os.path.isfile(python_exe):
+        print(f"Python executable not found: {python_exe}")
+        return 2
+
+    print(f"Using Python: {python_exe}")
+    print(f"Source folder: {source_dir}")
+    print(f"Examples to run: {len(EXAMPLES)}")
+
+    passed = 0
+    failed = 0
+
+    for example in EXAMPLES:
+        success = run_example(source_dir, example, python_exe)
+        if success:
+            passed += 1
+        else:
+            failed += 1
+            if not args.keep_going:
+                print("Stopping after first failure. Use --keep-going to continue.")
+                break
+
+    print("\nSummary")
+    print(f"Passed: {passed}")
+    print(f"Failed: {failed}")
+
+    return 0 if failed == 0 else 1
+
+
+if __name__ == "__main__":
+    sys.exit(main())
diff --git a/python/viennaps/_core/__init__.pyi b/python/viennaps/_core/__init__.pyi
index 8e9fbc82..9ef42d98 100644
--- a/python/viennaps/_core/__init__.pyi
+++ b/python/viennaps/_core/__init__.pyi
@@ -8,8 +8,8 @@ import typing
 import viennals._core
 from viennaps import d2
 import viennaps.d2
-from viennaps import d3
 import viennaps.d3
+from viennaps import d3
 from . import constants
 from . import gpu
 from . import util
@@ -195,6 +195,12 @@ class CF4O2Parameters:
     Si: CF4O2ParametersSi
     SiGe: CF4O2ParametersSiGe
     fluxIncludeSticking: bool
+    gamma_C: MaterialValueMap
+    gamma_C_oxidized: MaterialValueMap
+    gamma_F: MaterialValueMap
+    gamma_F_oxidized: MaterialValueMap
+    gamma_O: MaterialValueMap
+    gamma_O_passivated: MaterialValueMap
     def __init__(self) -> None:
         ...
     @property
@@ -210,42 +216,6 @@ class CF4O2Parameters:
     def etchantFlux(self, arg0: typing.SupportsFloat | typing.SupportsIndex) -> None:
         ...
     @property
-    def gamma_C(self) -> dict[Material, float]:
-        ...
-    @gamma_C.setter
-    def gamma_C(self, arg0: collections.abc.Mapping[Material, typing.SupportsFloat | typing.SupportsIndex]) -> None:
-        ...
-    @property
-    def gamma_C_oxidized(self) -> dict[Material, float]:
-        ...
-    @gamma_C_oxidized.setter
-    def gamma_C_oxidized(self, arg0: collections.abc.Mapping[Material, typing.SupportsFloat | typing.SupportsIndex]) -> None:
-        ...
-    @property
-    def gamma_F(self) -> dict[Material, float]:
-        ...
-    @gamma_F.setter
-    def gamma_F(self, arg0: collections.abc.Mapping[Material, typing.SupportsFloat | typing.SupportsIndex]) -> None:
-        ...
-    @property
-    def gamma_F_oxidized(self) -> dict[Material, float]:
-        ...
-    @gamma_F_oxidized.setter
-    def gamma_F_oxidized(self, arg0: collections.abc.Mapping[Material, typing.SupportsFloat | typing.SupportsIndex]) -> None:
-        ...
-    @property
-    def gamma_O(self) -> dict[Material, float]:
-        ...
-    @gamma_O.setter
-    def gamma_O(self, arg0: collections.abc.Mapping[Material, typing.SupportsFloat | typing.SupportsIndex]) -> None:
-        ...
-    @property
-    def gamma_O_passivated(self) -> dict[Material, float]:
-        ...
-    @gamma_O_passivated.setter
-    def gamma_O_passivated(self, arg0: collections.abc.Mapping[Material, typing.SupportsFloat | typing.SupportsIndex]) -> None:
-        ...
-    @property
     def ionFlux(self) -> float:
         ...
     @ionFlux.setter
diff --git a/python/viennaps/d2/__init__.pyi b/python/viennaps/d2/__init__.pyi
index 0177d786..f464981c 100644
--- a/python/viennaps/d2/__init__.pyi
+++ b/python/viennaps/d2/__init__.pyi
@@ -890,8 +890,17 @@ class SF6O2Etching(ProcessModel):
     def setParameters(self, arg0: viennaps._core.PlasmaEtchingParameters) -> None:
         ...
 class SelectiveEpitaxy(ProcessModel):
+    @typing.overload
+    def __init__(self, rate111: typing.SupportsFloat | typing.SupportsIndex = 0.5, rate100: typing.SupportsFloat | typing.SupportsIndex = 1.0) -> None:
+        ...
+    @typing.overload
     def __init__(self, materialRates: collections.abc.Sequence[tuple[viennaps._core.Material, typing.SupportsFloat | typing.SupportsIndex]], rate111: typing.SupportsFloat | typing.SupportsIndex = 0.5, rate100: typing.SupportsFloat | typing.SupportsIndex = 1.0) -> None:
         ...
+    @typing.overload
+    def __init__(self, nvFactors: typing.Annotated[collections.abc.Sequence[typing.SupportsFloat | typing.SupportsIndex], "FixedSize(3)"], rate111: typing.SupportsFloat | typing.SupportsIndex = 0.5, rate100: typing.SupportsFloat | typing.SupportsIndex = 1.0) -> None:
+        ...
+    def setMaterialRate(self, material: viennaps._core.Material, rate: typing.SupportsFloat | typing.SupportsIndex) -> None:
+        ...
 class SingleParticleALD(ProcessModel):
     def __init__(self, stickingProbability: typing.SupportsFloat | typing.SupportsIndex, numCycles: typing.SupportsInt | typing.SupportsIndex, growthPerCycle: typing.SupportsFloat | typing.SupportsIndex, totalCycles: typing.SupportsInt | typing.SupportsIndex, coverageTimeStep: typing.SupportsFloat | typing.SupportsIndex, evFlux: typing.SupportsFloat | typing.SupportsIndex, inFlux: typing.SupportsFloat | typing.SupportsIndex, s0: typing.SupportsFloat | typing.SupportsIndex, gasMFP: typing.SupportsFloat | typing.SupportsIndex) -> None:
         ...
diff --git a/python/viennaps/d2/gpu.pyi b/python/viennaps/d2/gpu.pyi
index 8c09ae89..caa5a1c4 100644
--- a/python/viennaps/d2/gpu.pyi
+++ b/python/viennaps/d2/gpu.pyi
@@ -25,7 +25,7 @@ class MultiParticleProcess(ProcessModelGPU):
     def addNeutralParticle(self, stickingProbability: typing.SupportsFloat | typing.SupportsIndex, label: str = 'neutralFlux') -> None:
         ...
     @typing.overload
-    def addNeutralParticle(self, materialSticking: collections.abc.Mapping[viennaps._core.Material, typing.SupportsFloat | typing.SupportsIndex], defaultStickingProbability: typing.SupportsFloat | typing.SupportsIndex = 1.0, label: str = 'neutralFlux') -> None:
+    def addNeutralParticle(self, materialSticking: viennaps._core.MaterialValueMap, defaultStickingProbability: typing.SupportsFloat | typing.SupportsIndex = 1.0, label: str = 'neutralFlux') -> None:
         ...
     def setRateFunction(self, arg0: collections.abc.Callable[[collections.abc.Sequence[typing.SupportsFloat | typing.SupportsIndex], viennaps._core.Material], float]) -> None:
         ...
diff --git a/python/viennaps/d3/__init__.pyi b/python/viennaps/d3/__init__.pyi
index 476bfe78..fc196ea6 100644
--- a/python/viennaps/d3/__init__.pyi
+++ b/python/viennaps/d3/__init__.pyi
@@ -890,8 +890,17 @@ class SF6O2Etching(ProcessModel):
     def setParameters(self, arg0: viennaps._core.PlasmaEtchingParameters) -> None:
         ...
 class SelectiveEpitaxy(ProcessModel):
+    @typing.overload
+    def __init__(self, rate111: typing.SupportsFloat | typing.SupportsIndex = 0.5, rate100: typing.SupportsFloat | typing.SupportsIndex = 1.0) -> None:
+        ...
+    @typing.overload
     def __init__(self, materialRates: collections.abc.Sequence[tuple[viennaps._core.Material, typing.SupportsFloat | typing.SupportsIndex]], rate111: typing.SupportsFloat | typing.SupportsIndex = 0.5, rate100: typing.SupportsFloat | typing.SupportsIndex = 1.0) -> None:
         ...
+    @typing.overload
+    def __init__(self, nvFactors: typing.Annotated[collections.abc.Sequence[typing.SupportsFloat | typing.SupportsIndex], "FixedSize(3)"], rate111: typing.SupportsFloat | typing.SupportsIndex = 0.5, rate100: typing.SupportsFloat | typing.SupportsIndex = 1.0) -> None:
+        ...
+    def setMaterialRate(self, material: viennaps._core.Material, rate: typing.SupportsFloat | typing.SupportsIndex) -> None:
+        ...
 class SingleParticleALD(ProcessModel):
     def __init__(self, stickingProbability: typing.SupportsFloat | typing.SupportsIndex, numCycles: typing.SupportsInt | typing.SupportsIndex, growthPerCycle: typing.SupportsFloat | typing.SupportsIndex, totalCycles: typing.SupportsInt | typing.SupportsIndex, coverageTimeStep: typing.SupportsFloat | typing.SupportsIndex, evFlux: typing.SupportsFloat | typing.SupportsIndex, inFlux: typing.SupportsFloat | typing.SupportsIndex, s0: typing.SupportsFloat | typing.SupportsIndex, gasMFP: typing.SupportsFloat | typing.SupportsIndex) -> None:
         ...
diff --git a/python/viennaps/d3/gpu.pyi b/python/viennaps/d3/gpu.pyi
index 59ff81af..dce1f224 100644
--- a/python/viennaps/d3/gpu.pyi
+++ b/python/viennaps/d3/gpu.pyi
@@ -25,7 +25,7 @@ class MultiParticleProcess(ProcessModelGPU):
     def addNeutralParticle(self, stickingProbability: typing.SupportsFloat | typing.SupportsIndex, label: str = 'neutralFlux') -> None:
         ...
     @typing.overload
-    def addNeutralParticle(self, materialSticking: collections.abc.Mapping[viennaps._core.Material, typing.SupportsFloat | typing.SupportsIndex], defaultStickingProbability: typing.SupportsFloat | typing.SupportsIndex = 1.0, label: str = 'neutralFlux') -> None:
+    def addNeutralParticle(self, materialSticking: viennaps._core.MaterialValueMap, defaultStickingProbability: typing.SupportsFloat | typing.SupportsIndex = 1.0, label: str = 'neutralFlux') -> None:
         ...
     def setRateFunction(self, arg0: collections.abc.Callable[[collections.abc.Sequence[typing.SupportsFloat | typing.SupportsIndex], viennaps._core.Material], float]) -> None:
         ...
diff --git a/tests/materialMap/materialMap.cpp b/tests/materialMap/materialMap.cpp
index 3fda8ba5..6ea10d06 100644
--- a/tests/materialMap/materialMap.cpp
+++ b/tests/materialMap/materialMap.cpp
@@ -1,5 +1,5 @@
+#include "materials/psMaterialMap.hpp"
 #include "materials/psMaterialValueMap.hpp"
-#include "materials/psMaterials.hpp"
 
 #include <vcTestAsserts.hpp>
 
diff --git a/tests/materialValueMap/CMakeLists.txt b/tests/materialValueMap/CMakeLists.txt
new file mode 100644
index 00000000..21dbcfd9
--- /dev/null
+++ b/tests/materialValueMap/CMakeLists.txt
@@ -0,0 +1,7 @@
+project(materialValueMap LANGUAGES CXX)
+
+add_executable(${PROJECT_NAME} "${PROJECT_NAME}.cpp")
+target_link_libraries(${PROJECT_NAME} PRIVATE ViennaPS)
+
+add_dependencies(ViennaPS_Tests ${PROJECT_NAME})
+add_test(NAME ${PROJECT_NAME} COMMAND $<TARGET_FILE:${PROJECT_NAME}>)
diff --git a/tests/materialValueMap/materialValueMap.cpp b/tests/materialValueMap/materialValueMap.cpp
new file mode 100644
index 00000000..a1bd79a9
--- /dev/null
+++ b/tests/materialValueMap/materialValueMap.cpp
@@ -0,0 +1,69 @@
+#include <materials/psMaterialMap.hpp>
+#include <materials/psMaterialValueMap.hpp>
+#include <vcTestAsserts.hpp>
+
+#include <algorithm>
+#include <vector>
+
+using namespace viennaps;
+
+void testBuiltInMaterial() {
+  MaterialValueMap<int> map;
+  map.setDefault(-1);
+  map.set(Material::Si, 1);
+  map.set(Material::SiO2, 2);
+
+  VC_TEST_ASSERT(map.get(Material::Si) == 1);
+  VC_TEST_ASSERT(map.get(Material::SiO2) == 2);
+  VC_TEST_ASSERT(map.get(Material::Si3N4) == -1); // default
+}
+
+void testCustomMaterial() {
+  MaterialValueMap<int> map;
+  map.setDefault(0);
+  Material customMat = MaterialMap::fromString("custom");
+  map.set(customMat, 3);
+
+  VC_TEST_ASSERT(map.get(customMat) == 3);
+  VC_TEST_ASSERT(map.get(Material::Si) == 0); // default
+
+  Material customMat2 = MaterialMap::fromString("custom2");
+  VC_TEST_ASSERT(map.get(customMat2) == 0); // default
+
+  // test operator[]
+  map[customMat2] = 4;
+  VC_TEST_ASSERT(map.get(customMat2) == 4);
+}
+
+void testIteration() {
+  MaterialValueMap<int> map;
+  auto customMat = MaterialMap::fromString("custom");
+  map.set(Material::Si, 1);
+  map.set(Material::SiO2, 2);
+  map.set(customMat, 7);
+
+  std::vector<std::pair<Material, int>> expected = {
+      {Material::Si, 1},
+      {Material::SiO2, 2},
+      {customMat, 7},
+  };
+
+  size_t idx = 0;
+  for (const auto &[material, value] : map) {
+    VC_TEST_ASSERT(idx < expected.size());
+    VC_TEST_ASSERT(material == expected[idx].first);
+    VC_TEST_ASSERT(value == expected[idx].second);
+    ++idx;
+  }
+
+  VC_TEST_ASSERT(map.getEntryByIndex(0) == 1);
+  VC_TEST_ASSERT(map.getEntryByIndex(1) == 2);
+  VC_TEST_ASSERT(map.getEntryByIndex(2) == 7);
+}
+
+int main() {
+  testBuiltInMaterial();
+  testCustomMaterial();
+  testIteration();
+  return 0;
+}
\ No newline at end of file
diff --git a/tests/readAndWrite/readAndWrite.cpp b/tests/readAndWrite/readAndWrite.cpp
index c363678a..1938c8fa 100644
--- a/tests/readAndWrite/readAndWrite.cpp
+++ b/tests/readAndWrite/readAndWrite.cpp
@@ -10,7 +10,7 @@
 
 #include <geometries/psMakeHole.hpp>
 #include <geometries/psMakeTrench.hpp>
-#include <materials/psMaterials.hpp>
+#include <materials/psMaterialMap.hpp>
 #include <psDomain.hpp>
 #include <psReader.hpp>
 #include <psWriter.hpp>