diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp
index 3f751709..79e51ef9 100644
--- a/libsrc/meshing/boundarylayer.cpp
+++ b/libsrc/meshing/boundarylayer.cpp
@@ -1422,10 +1422,19 @@ void BoundaryLayerTool ::FixVolumeElements() {
 
       if(params.project_boundaries.has_value())
         {
-          regex pattern(*params.project_boundaries);
-          for(int i = 1; i<=mesh.GetNFD(); i++)
-            if(regex_match(mesh.GetFaceDescriptor(i).GetBCName(), pattern))
-              par_project_boundaries.Append(i);
+          auto proj_bnd = *params.project_boundaries;
+          if(string* s = get_if<string>(&proj_bnd); s)
+            {
+              regex pattern(*s);
+              for(int i = 1; i<=mesh.GetNFD(); i++)
+                if(regex_match(mesh.GetFaceDescriptor(i).GetBCName(), pattern))
+                  par_project_boundaries.Append(i);
+            }
+          else
+           {
+              for(auto id : *get_if<std::vector<int>>(&proj_bnd))
+                par_project_boundaries.Append(id);
+           }
         }
 
       if(double* height = get_if<double>(&params.thickness); height)
@@ -1449,10 +1458,14 @@ void BoundaryLayerTool ::FixVolumeElements() {
             if(regex_match(mesh.GetMaterial(i), pattern))
               domains.SetBit(i);
         }
+      else if(int *idomain = get_if<int>(&params.domain); idomain)
+        {
+          domains.SetBit(*idomain);
+        }
       else
         {
-          auto idomain = *get_if<int>(&params.domain);
-          domains.SetBit(idomain);
+          for (auto i : *get_if<std::vector<int>>(&params.domain))
+            domains.SetBit(i);
         }
   }
 
diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp
index a14d9ca5..7f49d334 100644
--- a/libsrc/meshing/meshtype.cpp
+++ b/libsrc/meshing/meshtype.cpp
@@ -2907,42 +2907,30 @@ namespace netgen
 
   ostream & operator<< (ostream & ost, const BoundaryLayerParameters & mp)
   {
+    auto print_vec = [&ost](auto & v) {
+        ost << "[";
+        for (const auto & val : v)
+          ost << val << " ";
+        ost << "]";
+    };
     ost << "BoundaryLayerParameters" << endl;
     ost << "  boundary: ";
     switch(mp.boundary.index())
       {
-        case 0:
-        ost << std::get<0>(mp.boundary);
-        break;
-        case 1:
-        ost << std::get<1>(mp.boundary);
-        break;
-        case 2:
-        ost << "[";
-        for (auto val : std::get<2>(mp.boundary))
-          ost << val << " ";
-        ost << "]";
-        break;
+        case 0: ost << std::get<0>(mp.boundary); break;
+        case 1: ost << std::get<1>(mp.boundary); break;
+        case 2: print_vec(std::get<2>(mp.boundary)); break;
       }
     ost << "\n  thickness: ";
     switch(mp.thickness.index())
       {
-        case 0:
-        ost << std::get<0>(mp.thickness);
-        break;
-        case 1:
-        ost << "[";
-        for (auto val : std::get<1>(mp.thickness))
-          ost << val << " ";
-        ost << "]";
-        break;
+        case 0: ost << std::get<0>(mp.thickness); break;
+        case 1: print_vec(std::get<1>(mp.thickness)); break;
       }
     ost <<"\n  new_material: ";
     switch(mp.new_material.index())
       {
-        case 0:
-        ost << std::get<0>(mp.new_material);
-        break;
+        case 0: ost << std::get<0>(mp.new_material); break;
         case 1:
         for (const auto & [key, value] : std::get<1>(mp.new_material))
           ost << key << " -> " << value << ", ";
@@ -2951,15 +2939,22 @@ namespace netgen
     ost << "\n  domain: ";
     switch(mp.domain.index())
       {
-        case 0:
-        ost << std::get<0>(mp.domain);
-        break;
-        case 1:
-        ost << std::get<1>(mp.domain);
-        break;
+        case 0: ost << std::get<0>(mp.domain); break;
+        case 1: ost << std::get<1>(mp.domain); break;
+        case 2: print_vec(std::get<2>(mp.domain)); break;
       }
     ost << "\n  outside: " << mp.outside;
-    ost << "\n  project_boundaries: " << (mp.project_boundaries ? ToString(*mp.project_boundaries) : "nullopt");
+    ost << "\n  project_boundaries: ";
+    if(mp.project_boundaries) {
+      auto & proj_bnd = *mp.project_boundaries;
+      switch(proj_bnd.index())
+        {
+          case 0: ost << std::get<0>(proj_bnd); break;
+          case 1: print_vec(std::get<1>(proj_bnd)); break;
+        }
+    }
+    else
+      ost << "nullopt";
     ost << "\n  grow_edges: " << mp.grow_edges;
     ost << "\n  limit_growth_vectors: " << mp.limit_growth_vectors;
     ost << "\n  sides_keep_surfaceindex: " << mp.sides_keep_surfaceindex;
diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp
index ec189ad6..178dd471 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -1279,9 +1279,9 @@ namespace netgen
     std::variant<string, int, std::vector<int>> boundary;
     std::variant<double, std::vector<double>> thickness;
     std::variant<string, std::map<string, string>> new_material;
-    std::variant<string, int> domain;
+    std::variant<string, int, std::vector<int>> domain;
     bool outside;
-    std::optional<string> project_boundaries;
+    std::optional<std::variant<string, std::vector<int>>> project_boundaries;
     bool grow_edges;
     bool limit_growth_vectors;
     bool sides_keep_surfaceindex;
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index 6937c905..0500ac5c 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -1472,8 +1472,8 @@ py::arg("point_tolerance") = -1.)
     .def ("BoundaryLayer", [](Mesh & self, variant<string, int, std::vector<int>> boundary,
                               variant<double, std::vector<double>> thickness,
                               variant<string, map<string, string>> material,
-                              variant<string, int> domain, bool outside,
-                              optional<string> project_boundaries,
+                              variant<string, int, std::vector<int>> domain, bool outside,
+                              optional<variant<string, std::vector<int>>> project_boundaries,
                               bool grow_edges, bool limit_growth_vectors,
                               bool sides_keep_surfaceindex,
                               bool keep_surfaceindex)
@@ -1725,9 +1725,9 @@ py::arg("point_tolerance") = -1.)
         std::variant<string, int, std::vector<int>> boundary,
         std::variant<double, std::vector<double>> thickness,
         std::variant<string, std::map<string, string>> new_material,
-        std::variant<string, int> domain,
+        std::variant<string, int, std::vector<int>> domain,
         bool outside,
-        std::optional<string> project_boundaries,
+        std::optional<std::variant<string, std::vector<int>>> project_boundaries,
         bool grow_edges,
         bool limit_growth_vectors,
         bool sides_keep_surfaceindex,