diff --git a/.gitignore b/.gitignore
index 1fd2b13..130eeec 100644
--- a/.gitignore
+++ b/.gitignore
@@ -14,6 +14,7 @@ src/*.json
 
 .vs/
 .vscode/settings.json
+.vscode/settings.dev.json
 .vscode/launch.json
 .vscode/tasks.json
 
diff --git a/.vscode/settings.dev.json b/.vscode/settings.dev.json
index 5303634..9d190ee 100644
--- a/.vscode/settings.dev.json
+++ b/.vscode/settings.dev.json
@@ -8,12 +8,8 @@
   "C_Cpp.intelliSenseEngine": "disabled",
   "clangd.arguments": ["-log=verbose",
                           "-pretty",
-                          "--background-index",
+                          "--background-index", 
                           "--query-driver=/usr/bin/g++"
   ],
-  "clangd.fallbackFlags": ["-std=c++17"],
-	"lldb.showDisassembly": "never",
-	"lldb.launch.preRunCommands": [
-		"settings set target.import-std-module true"
-	]
+  "clangd.fallbackFlags": ["-std=c++17"]
 }
diff --git a/src/app/launcher.cpp b/src/app/launcher.cpp
index 6ce1de5..c9bb54c 100644
--- a/src/app/launcher.cpp
+++ b/src/app/launcher.cpp
@@ -150,8 +150,9 @@ meshlib::meshers::StaircaseMesherOptions readStaircaseMesherOptions(const nlohma
     }
 
     meshlib::meshers::StaircaseMesherOptions res;
-    
-    res.isVolume = isVolume;
+    if (isVolume){
+        res.volumeGroups.insert(0);
+    }
     if (mesherConfig.contains("options") && 
         mesherConfig["options"].contains("compress")) {
         res.compress = mesherConfig["options"]["compress"];
@@ -160,7 +161,7 @@ meshlib::meshers::StaircaseMesherOptions readStaircaseMesherOptions(const nlohma
     return res;
 }
 
-meshlib::meshers::ConformalMesherOptions readConformalMesherOptions(const nlohmann::json& fileData, const std::optional<nlohmann::json>& override)
+meshlib::meshers::ConformalMesherOptions readConformalMesherOptions(const nlohmann::json& fileData, bool isVolume, const std::optional<nlohmann::json>& override)
 {    
     nlohmann::json mesherConfig;
     if (override.has_value()) {
@@ -170,6 +171,9 @@ meshlib::meshers::ConformalMesherOptions readConformalMesherOptions(const nlohma
     }
     
     meshlib::meshers::ConformalMesherOptions res;
+    if (isVolume){
+        res.volumeGroups.insert(0);
+    }
     if (mesherConfig.contains("options")) {
         res.snapperOptions.edgePoints = mesherConfig["options"]["edgePoints"];
         res.snapperOptions.forbiddenLength = mesherConfig["options"]["forbiddenLength"];
@@ -204,7 +208,7 @@ std::unique_ptr<meshlib::meshers::MesherBase> buildMesher(const Mesh& in, const
             readStaircaseMesherOptions(fileData, objDef.isVolume, objDef.mesherOverride)
         });
     } else if (mesherType == meshlib::app::conformal_mesher) {
-        return std::make_unique<meshlib::meshers::ConformalMesher>(meshlib::meshers::ConformalMesher{in, readConformalMesherOptions(fileData, objDef.mesherOverride)});
+        return std::make_unique<meshlib::meshers::ConformalMesher>(meshlib::meshers::ConformalMesher{in, readConformalMesherOptions(fileData, objDef.isVolume, objDef.mesherOverride)});
     } else {
         throw std::runtime_error("Unsupported mesher type");
     }
diff --git a/src/app/vtkIO.cpp b/src/app/vtkIO.cpp
index d55a561..70a105c 100644
--- a/src/app/vtkIO.cpp
+++ b/src/app/vtkIO.cpp
@@ -3,6 +3,7 @@
 #include <vtkCellData.h>
 #include <vtkCellType.h>
 #include <vtkTriangle.h>
+#include <vtkTetra.h>
 #include <vtkQuad.h>
 #include <vtkLine.h>
 #include <vtkVertex.h>
@@ -76,6 +77,7 @@ Element vtkCellToElement(vtkCell* cell)
     vtkVertex* vertex = nullptr;
     vtkLine* line = nullptr;
     vtkTriangle* triangle = nullptr;
+    vtkTetra* tetra = nullptr;
 
     switch (cell->GetCellType()) {
     case VTK_VERTEX:
@@ -102,6 +104,17 @@ Element vtkCellToElement(vtkCell* cell)
         };
         elem.type = meshlib::Element::Type::Surface;
         break;
+
+    case VTK_TETRA:
+        tetra = vtkTetra::SafeDownCast(cell);
+        elem.vertices = {
+            CoordinateId(tetra->GetPointIds()->GetId(0)),
+            CoordinateId(tetra->GetPointIds()->GetId(1)),
+            CoordinateId(tetra->GetPointIds()->GetId(2)),
+            CoordinateId(tetra->GetPointIds()->GetId(3))
+        };
+        elem.type = meshlib::Element::Type::Volume;
+        break;
     }
     
     return elem;
@@ -131,6 +144,7 @@ Mesh vtuToMesh(vtkUnstructuredGrid* vtu)
         }
     } else {
         mesh.groups.resize(1);
+        auto k = vtu->GetNumberOfCells();
         mesh.groups[0].elements.reserve(vtu->GetNumberOfCells());
         for (vtkIdType i = 0; i < vtu->GetNumberOfCells(); i++) {
             mesh.groups[0].elements.push_back(
@@ -169,8 +183,13 @@ vtkSmartPointer<vtkStringArray> toVTKGroupNamesArray(const Mesh& mesh)
     groupNamesArray->SetName("groupNames");
     groupNamesArray->SetNumberOfComponents(1);
     
-    for (const auto& group : mesh.groups) {
-        groupNamesArray->InsertNextValue(group.name.c_str());
+    // for (const auto& group : mesh.groups) {
+    //     groupNamesArray->InsertNextValue(group.name.c_str());
+    // }
+    for (auto g = 0; g < mesh.groups.size(); g++) {
+        for (auto e = 0; e < mesh.groups[g].elements.size(); e++) {
+            groupNamesArray->InsertNextValue( mesh.groups[g].name.c_str() );
+        }
     }
     
     return groupNamesArray;
diff --git a/src/cgal/filler/Filler.cpp b/src/cgal/filler/Filler.cpp
index ae72d0a..60c23bf 100644
--- a/src/cgal/filler/Filler.cpp
+++ b/src/cgal/filler/Filler.cpp
@@ -253,6 +253,29 @@ void sliceAlignedByGrid(
 	);
 }
 
+void sliceAlignedByGridAndRemove(
+	Filler::GridSlices& slices,
+	const Polyhedron& m,
+	const Grid& g,
+	const Priority& priority)
+{
+	
+	auto polygons{ buildGridPlanesPolygons(makeFacesCCWOriented(m), g)};
+	const std::array<Axis, 3> axis{ X, Y, Z };
+
+	std::for_each(
+#ifdef TESSELLATOR_EXECUTION_POLICIES
+		std::execution::par,
+#endif
+		axis.begin(), axis.end(),
+		[&](const auto& x) {
+			for (const auto& [i, polygon] : polygons[x]) {
+				slices[x][i].remove(polygon, priority);
+			}
+		}
+	);
+}
+
 Priority Filler::getGroupPriority(const GroupId& gId) const
 {
 	if (gId < groupPriorities_.size()) {
@@ -508,7 +531,8 @@ FillerPolyhedrons buildFillerPolyhedrons(
 Filler::Filler(
 	const Mesh& volumeMesh,
 	const Mesh& surfaceMesh,
-	const std::vector<Priority>& groupPriorities)
+	const std::vector<Priority>& groupPriorities,
+	const FillerMode& fillerMode)
 {
 	utils::meshTools::checkNoNullAreasExist(volumeMesh);
 	utils::meshTools::checkNoNullAreasExist(surfaceMesh);
@@ -547,7 +571,11 @@ Filler::Filler(
 		log("Slicing surfaces", 2);
 		sliceNonAlignedByGrid(slices_, fP.surfaces, grid_, pr, SlicingMode::Surface);
 		log("Slicing aligned", 2);
-		sliceAlignedByGrid(slices_, fP.aligned, grid_, pr);
+		if (fillerMode == FillerMode::insideAndOutside) {
+			sliceAlignedByGrid(slices_, fP.aligned, grid_, pr);
+		} else if (fillerMode == FillerMode::onlyInside){
+			sliceAlignedByGridAndRemove(slices_, fP.aligned, grid_, pr);
+		}
 		log("Building segments arrays", 2);
 		buildSegmentsArray(segmentsArray_, fP.aligned, grid_, pr);
 		buildSegmentsArray(segmentsArray_, fP.volumes, grid_, pr);
diff --git a/src/cgal/filler/Filler.h b/src/cgal/filler/Filler.h
index 2b2d5af..ed87558 100644
--- a/src/cgal/filler/Filler.h
+++ b/src/cgal/filler/Filler.h
@@ -9,17 +9,26 @@
 
 namespace meshlib::cgal::filler {
 
+enum class FillerMode{
+	insideAndOutside,
+	onlyInside
+};
+
+
 class Filler {
-public:
+	public:
 	using Slices = std::map<SliceNumber, Slice>;
 	using GridSlices = std::array<Slices, 3>;
 	using SegmentsArray = std::map<ArrayIndex, Segments>;
 	using GridSegmentsArray = std::array<SegmentsArray, 3>;
-
+	
+	FillerMode mode = FillerMode::insideAndOutside;
+	
 	Filler(
 		const Mesh& volumeMesh, 
 		const Mesh& surfaceMesh = Mesh(),
-		const std::vector<Priority>& groupPriorities = std::vector<Priority>());
+		const std::vector<Priority>& groupPriorities = std::vector<Priority>(),
+		const FillerMode& mode = FillerMode::insideAndOutside);
 	Filler(const Filler&) = delete;
 	Filler(Filler&&) = default;
 	Filler& operator=(const Filler&) = delete;
@@ -32,7 +41,7 @@ class Filler {
 	FillingState getFillingState(const CellIndex&) const;
 	
 	Mesh getMeshFilling() const;
-
+	GridSlices getSlices() const {return slices_;};
 private:
 	GridSlices slices_;
 	GridSegmentsArray segmentsArray_;
@@ -42,6 +51,7 @@ class Filler {
 
 	void mergeGroupsWithSamePriority(Groups& vGroups, Groups& sGroups);
 
+
 };
 
 
diff --git a/src/cgal/filler/Slice.cpp b/src/cgal/filler/Slice.cpp
index 9a674d7..5133ba6 100644
--- a/src/cgal/filler/Slice.cpp
+++ b/src/cgal/filler/Slice.cpp
@@ -183,6 +183,21 @@ void Slice::add(const Polylines2& polylines, const Priority& pr)
 	}
 }
 
+void Slice::remove(const Polylines2& polylines, const Priority& pr)
+{
+	SliceData& sd = data_[pr];
+	
+	for (const auto& p : polylines) {
+		if (p.size() == 1) {
+			continue;
+		}
+		auto r{ removeSegmentsContainedInAnyAxis(p) };
+		for (const auto& rr : r){
+			sd.lines.erase(std::find(sd.lines.begin(), sd.lines.end(), rr));
+		}
+	}
+}
+
 FillingState::FillingState(const FillingType& t) :
 	type{ t },
 	priority_{ 0 }
@@ -261,6 +276,18 @@ void Slice::add(const HPolygonSet& polygons, const Priority& pr)
 	removeInSuperiorPriorities(pr);
 }
 
+void Slice::remove(const HPolygonSet& polygons, const Priority& pr)
+{
+	if (polygons.isEmpty()) {
+		return;
+	}
+
+	SliceData& sd = data_[pr];
+	sd.surfaces.difference(polygons);
+
+	removeInSuperiorPriorities(pr);
+}
+
 void Slice::mergeLines(const Slice& lhs)
 {
 	for (const auto& [pr, sd] : lhs.data_) {
diff --git a/src/cgal/filler/Slice.h b/src/cgal/filler/Slice.h
index b61883c..b04bf76 100644
--- a/src/cgal/filler/Slice.h
+++ b/src/cgal/filler/Slice.h
@@ -85,8 +85,10 @@ class Slice {
 	FillingState getFillingState(const ArrayIndex&) const;
 
 	void add(const Polylines2&, const Priority&);
+	void remove(const Polylines2&, const Priority&);
 	void addAsPolygon(const Polylines2&, const Priority&);
 	void add(const HPolygonSet&, const Priority&);
+	void remove(const HPolygonSet&, const Priority&);
 	void mergeLines(const Slice& lhs);
 	void buildSearchMap();
 	void buildTriangulations();
diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt
index faab34a..205c3a6 100644
--- a/src/core/CMakeLists.txt
+++ b/src/core/CMakeLists.txt
@@ -10,4 +10,7 @@ add_library(tessellator-core
     "Staircaser.cpp"
 )
 
-target_link_libraries(tessellator-core tessellator-utils)
\ No newline at end of file
+target_link_libraries(tessellator-core 
+    tessellator-utils
+    tessellator-cgal
+    CGAL::CGAL)
\ No newline at end of file
diff --git a/src/core/Slicer.cpp b/src/core/Slicer.cpp
index fb5d7a2..8246fff 100644
--- a/src/core/Slicer.cpp
+++ b/src/core/Slicer.cpp
@@ -38,7 +38,6 @@ void orient(const Coordinates& coords,
     }
 }
 
-
 Slicer::Slicer(const Mesh& input, const std::vector<Element::Type>& dimensionPolicy, const SlicerOptions& opts) :
     GridTools(input.grid),
     opts_(opts)
diff --git a/src/core/Slicer.h b/src/core/Slicer.h
index e0bae73..b51ae15 100644
--- a/src/core/Slicer.h
+++ b/src/core/Slicer.h
@@ -24,8 +24,6 @@ class Slicer : public utils::GridTools {
 
     Slicer(const Mesh&, const std::vector<Element::Type>& dimensionPolicy = {}, const SlicerOptions& opts = SlicerOptions());
     Mesh getMesh() const { return mesh_; };
-
-
     static Elements buildTrianglesFromPath(const std::vector<Coordinate>&, const std::vector<CoordinateId>&);
 
 private:
@@ -53,6 +51,11 @@ class Slicer : public utils::GridTools {
         const Cell&, const Coordinate&,
         const Cell&, const Coordinate&) const;
 
+    // IdSet buildGroupIntersectionsWithGridPlanes(
+    //     Coordinates& sCoords,
+    //     const std::vector<Element>& elements);
+
+
 };
 
 }
diff --git a/src/meshers/MesherBaseOptions.h b/src/meshers/MesherBaseOptions.h
index 7403619..55138c1 100644
--- a/src/meshers/MesherBaseOptions.h
+++ b/src/meshers/MesherBaseOptions.h
@@ -7,7 +7,6 @@ namespace meshlib::meshers {
 
 class MesherBaseOptions {
 public:
-    bool isVolume = false;
     std::set<GroupId> volumeGroups{};
 };
 
diff --git a/src/meshers/StaircaseMesher.cpp b/src/meshers/StaircaseMesher.cpp
index 6e9034d..d2d9006 100644
--- a/src/meshers/StaircaseMesher.cpp
+++ b/src/meshers/StaircaseMesher.cpp
@@ -9,6 +9,7 @@
 #include "core/Compressor.h"
 
 #include "cgal/filler/Filler.h"
+// #include "cgal/Manifolder.h"
 
 #include "utils/RedundancyCleaner.h"
 #include "utils/MeshTools.h"
@@ -26,14 +27,39 @@ StaircaseMesher::StaircaseMesher(const Mesh& inputMesh, int decimalPlacesInColla
     opts_(opts)
 {
     log("Preparing surfaces.");
-    surfaceMesh_ = buildMeshFilteringElements(inputMesh, isNotTetrahedron);
-
+    surfaceMesh_ = MesherBase::buildSurfaceMesh(inputMesh, opts_.volumeGroups);
     log("Processing surface mesh.");
     process(surfaceMesh_);
+
+    log("Preparing volumes");
+    volumeMesh_ = MesherBase::buildVolumeMesh(inputMesh, opts_.volumeGroups);
+    fillMesh(volumeMesh_);
+    log("Processing volume mesh.");
+    process(volumeMesh_);
+
+    mergeMesh(surfaceMesh_, volumeMesh_);
     
-    log("Surface mesh built succesfully.", 1);
+    log("Mesh built succesfully.", 1);
 }
 
+void StaircaseMesher::fillMesh(Mesh& m){
+    if (m.countElems() == 0) return;
+    meshlib::cgal::filler::FillerMode mode = meshlib::cgal::filler::FillerMode::onlyInside;
+    if (m.groups[0].elements[0].isTetrahedron()) {
+        mode = meshlib::cgal::filler::FillerMode::insideAndOutside;
+    }
+    auto filling = m;
+    utils::meshTools::convertToRelativeCoordinates(filling);
+    filling = meshlib::cgal::filler::Filler(filling, Mesh(), std::vector<Priority>(), mode).getMeshFilling();
+    utils::meshTools::convertToAbsoluteCoordinates(filling);
+    if (mode == meshlib::cgal::filler::FillerMode::insideAndOutside){
+        m = filling;
+    } else if (mode == meshlib::cgal::filler::FillerMode::onlyInside){
+        mergeMesh(m, filling);
+    }
+}
+
+
 Mesh StaircaseMesher::buildSurfaceMesh(const Mesh& inputMesh, const Mesh & volumeSurface)
 {
     auto resultMesh = buildMeshFilteringElements(inputMesh, isNotTetrahedron);
@@ -41,28 +67,40 @@ Mesh StaircaseMesher::buildSurfaceMesh(const Mesh& inputMesh, const Mesh & volum
     return resultMesh;
 }
 
+std::vector<std::string> getGroupNames(const Groups& groups){
+    std::vector<std::string> names;
+    names.reserve(groups.size());
+    for (auto gId{0}; gId < groups.size(); ++gId) {        
+        names.push_back(groups[gId].name);
+    }
+    return names;
+}
+
+void copyGroupNames(Mesh& m, const std::vector<std::string>& names){
+    for (auto gId{0}; gId < m.groups.size(); ++gId) {        
+        m.groups[gId].name = names[gId];
+    }
+}
+
+static Mesh toAbsolute(const Mesh& m) 
+{
+    auto r{ m };
+    r.coordinates = 
+        utils::GridTools{ m.grid }.relativeToAbsolute(m.coordinates);
+    return r;
+}
+
 void StaircaseMesher::process(Mesh& mesh) const
 {
-    
+    const auto groupNames = getGroupNames(mesh.groups);
     const auto slicingGrid{ buildSlicingGrid(originalGrid_, enlargedGrid_) };
-    
     if (mesh.countElems() == 0) {
-        mesh.grid = slicingGrid;
+        // mesh.grid = slicingGrid;
         return;
     }
 
     auto dimensions = getHighestDimensionByGroup(mesh);
 
-    if (opts_.isVolume){
-        if (meshTools::isAClosedTopology(mesh.groups[0].elements)){
-            meshlib::cgal::filler::Filler f{ mesh };
-            auto filling = f.getMeshFilling();
-            mergeMesh(mesh, filling);
-        } else {
-            throw std::runtime_error("Input object marked to be meshed as a volume, but surface is not closed");
-        }
-    }
-
     log("Slicing.", 1);
     mesh.grid = slicingGrid;
     mesh = Slicer{ mesh, dimensions }.getMesh();
@@ -113,6 +151,8 @@ void StaircaseMesher::process(Mesh& mesh) const
     logNumberOfQuads(countMeshElementsIf(mesh, isQuad));
     logNumberOfLines(countMeshElementsIf(mesh, isLine));
 
+    copyGroupNames(mesh, groupNames);
+
 }
 
 
diff --git a/src/meshers/StaircaseMesher.h b/src/meshers/StaircaseMesher.h
index 3cabd4b..c6be343 100644
--- a/src/meshers/StaircaseMesher.h
+++ b/src/meshers/StaircaseMesher.h
@@ -17,9 +17,11 @@ class StaircaseMesher : public MesherBase {
 	int decimalPlacesInCollapser_;
 
 	Mesh surfaceMesh_;
+	Mesh volumeMesh_;
 	StaircaseMesherOptions opts_;
 
 	virtual Mesh buildSurfaceMesh(const Mesh& inputMesh, const Mesh& volumeSurface);
+	static void fillMesh(Mesh& inputMesh);
 	void process(Mesh&) const;
 
 };
diff --git a/src/utils/Geometry.cpp b/src/utils/Geometry.cpp
index fe2fc4e..2015763 100644
--- a/src/utils/Geometry.cpp
+++ b/src/utils/Geometry.cpp
@@ -70,6 +70,17 @@ std::vector<ElementsView> Geometry::buildDisjointSmoothSets(
 }
 
 
+QuaV Geometry::asQuaV(const Element& el, const std::vector<Coordinate>& co) {
+    if (el.vertices.size() != 4) {
+        throw std::logic_error("Invalid conversion from element to QuaV");
+    }
+    QuaV res;
+    for (std::size_t i = 0; i < el.vertices.size(); i++) {
+        res[i] = co[el.vertices[i]];
+    }
+    return res;
+}
+
 TriV Geometry::asTriV(const Element& el, const std::vector<Coordinate>& co) {
     if (el.vertices.size() != 3) {
         throw std::logic_error("Invalid conversion from element to TriV");
@@ -229,6 +240,18 @@ double Geometry::area(const TriV& tri) {
     return ((tri[0] - tri[1]) ^ (tri[1] - tri[2])).norm() / 2.0;
 }
 
+double Geometry::area(const QuaV& qua) {
+    const Coordinates cs{ qua.begin(), qua.end() };
+
+    VecD crossSum;
+    for (std::size_t i = 0; i < qua.size(); ++i) {
+        const auto& p = qua[i];
+        const auto& q = qua[(i + 1) % qua.size()];
+        crossSum += (p ^ q);
+    }
+    return 0.5*crossSum.norm();
+}
+
 
 }
 }
diff --git a/src/utils/Geometry.h b/src/utils/Geometry.h
index 113593b..e185f0c 100644
--- a/src/utils/Geometry.h
+++ b/src/utils/Geometry.h
@@ -28,6 +28,7 @@ class Geometry {
     static bool areAdjacentLines(const Element&, const Element&);
     
 
+    static QuaV asQuaV(const Element&, const Coordinates&);
     static TriV asTriV(const Element&, const Coordinates&);
     static LinV asLinV(const Element&, const Coordinates&);
     
@@ -38,6 +39,7 @@ class Geometry {
     static VecD getCentroid(const Element&, const std::vector<Coordinate>&);
     static VecD getCentroid(const TriV&);
     static double area(const TriV& tri);
+    static double area(const QuaV& qua);
     static bool isDegenerate(const TriV& tri, const double& areaTolerance = NORM_TOLERANCE);
     static bool areCollinear(const Coordinates&);
     template <std::size_t N>
diff --git a/src/utils/MeshTools.cpp b/src/utils/MeshTools.cpp
index 02bbe57..00e0f1b 100644
--- a/src/utils/MeshTools.cpp
+++ b/src/utils/MeshTools.cpp
@@ -273,12 +273,22 @@ void checkNoNullAreasExist(const Mesh& m)
                     msg << info(e, m) << std::endl;
                 }
             }
-            else if (Geometry::area(Geometry::asTriV(e, m.coordinates)) == 0.0) {
-                nullAreas = true;
-                msg << std::endl;
-                msg << "Group: " << &g - &m.groups.front()
-                    << ", Element: " << &e - &g.elements.front() << std::endl;
-                msg << info(e, m) << std::endl;
+            else if (e.isTriangle()){
+                if (Geometry::area(Geometry::asTriV(e, m.coordinates)) == 0.0) {
+                    nullAreas = true;
+                    msg << std::endl;
+                    msg << "Group: " << &g - &m.groups.front()
+                        << ", Element: " << &e - &g.elements.front() << std::endl;
+                    msg << info(e, m) << std::endl;
+                }
+            } else if (e.isQuad()){
+                if (Geometry::area(Geometry::asQuaV(e, m.coordinates)) == 0.0) {
+                    nullAreas = true;
+                    msg << std::endl;
+                    msg << "Group: " << &g - &m.groups.front()
+                        << ", Element: " << &e - &g.elements.front() << std::endl;
+                    msg << info(e, m) << std::endl;
+                }
             }
         }
     }
@@ -299,6 +309,17 @@ void convertToAbsoluteCoordinates(Mesh& m)
     );
 }
 
+void convertToRelativeCoordinates(Mesh& m)
+{
+    GridTools gT{ m.grid };
+
+    std::transform(
+        m.coordinates.begin(), m.coordinates.end(),
+        m.coordinates.begin(),
+        [&](const auto& v) { return gT.getRelative(v); }
+    );
+}
+
 void checkSlicedMeshInvariants(const Mesh& m)
 {
     checkNoCellsAreCrossed(m);
@@ -319,6 +340,7 @@ Mesh buildMeshFilteringElements(
             inElems.begin(), inElems.end(), 
             std::back_inserter(r.groups[gId].elements), 
             filter);
+        r.groups[gId].name = in.groups[gId].name;
     }
     return r;
 }
@@ -378,7 +400,7 @@ void mergeMesh(Mesh& lMesh, const Mesh& iMesh)
     assert(lMesh.groups.size() == iMesh.groups.size());
 
     auto coordCount{ lMesh.coordinates.size() };
-    
+    if (iMesh.countElems() == 0) return;
     lMesh.coordinates.insert(lMesh.coordinates.end(),
         iMesh.coordinates.begin(), iMesh.coordinates.end());
 
diff --git a/src/utils/MeshTools.h b/src/utils/MeshTools.h
index 5fc4546..9bfed81 100644
--- a/src/utils/MeshTools.h
+++ b/src/utils/MeshTools.h
@@ -35,6 +35,7 @@ void reduceGrid(Mesh&, const Grid&);
 Mesh reduceGrid(const Mesh& m, const Grid& g);
 
 void convertToAbsoluteCoordinates(Mesh&);
+void convertToRelativeCoordinates(Mesh&);
 	
 void checkSlicedMeshInvariants(Mesh& m);
 	
diff --git a/test/MeshFixtures.h b/test/MeshFixtures.h
index 65c77ad..7eda7c0 100644
--- a/test/MeshFixtures.h
+++ b/test/MeshFixtures.h
@@ -38,6 +38,7 @@ static Mesh buildNonManifoldPatchMesh(double stepSize)
     return m;
 }
 
+
 static Mesh buildTetAndTriMesh(double stepSize) 
 {
     Mesh m;
diff --git a/test/app/launcherTest.cpp b/test/app/launcherTest.cpp
index 5b3f3e4..056f5ee 100644
--- a/test/app/launcherTest.cpp
+++ b/test/app/launcherTest.cpp
@@ -71,7 +71,7 @@ TEST_F(LauncherTest, builds_staircased_mesher_default)
 
     EXPECT_NO_THROW(auto staircaseMesher = dynamic_cast<meshlib::meshers::StaircaseMesher&>(*mesher));
     const auto & options = dynamic_cast<meshlib::meshers::StaircaseMesher&>(*mesher).getOptions();
-    EXPECT_EQ(options.isVolume, false);
+    EXPECT_EQ(options.volumeGroups.size(), 0);
     EXPECT_EQ(options.compress, false);
 }
 
@@ -95,7 +95,7 @@ TEST_F(LauncherTest, builds_staircased_mesher_without_compression)
 
     EXPECT_NO_THROW(auto staircaseMesher = dynamic_cast<meshlib::meshers::StaircaseMesher&>(*mesher));
     const auto & options = dynamic_cast<meshlib::meshers::StaircaseMesher&>(*mesher).getOptions();
-    EXPECT_EQ(options.isVolume, false);
+    EXPECT_EQ(options.volumeGroups.size(), 0);
     EXPECT_EQ(options.compress, false);
 }
 
@@ -119,7 +119,7 @@ TEST_F(LauncherTest, builds_staircased_mesher_with_compression)
     
     EXPECT_NO_THROW(auto staircaseMesher = dynamic_cast<meshlib::meshers::StaircaseMesher&>(*mesher));
     const auto & options = dynamic_cast<meshlib::meshers::StaircaseMesher&>(*mesher).getOptions();
-    EXPECT_EQ(options.isVolume, false);
+    EXPECT_EQ(options.volumeGroups.size(), 0);
     EXPECT_EQ(options.compress, true);
 }
 
@@ -315,7 +315,7 @@ TEST_F(LauncherTest, builds_staircased_mesher_with_override)
     
     EXPECT_NO_THROW(auto staircaseMesher = dynamic_cast<meshlib::meshers::StaircaseMesher&>(*mesher));
     const auto & options = dynamic_cast<meshlib::meshers::StaircaseMesher&>(*mesher).getOptions();
-    EXPECT_EQ(options.isVolume, false);
+    EXPECT_EQ(options.volumeGroups.size(), 0);
     EXPECT_EQ(options.compress, true);
 }
 
diff --git a/test/cgal/ManifolderTest.cpp b/test/cgal/ManifolderTest.cpp
index cedd1c3..c05f1c5 100644
--- a/test/cgal/ManifolderTest.cpp
+++ b/test/cgal/ManifolderTest.cpp
@@ -87,6 +87,7 @@ TEST_F(ManifolderTest, volume_and_surface)
 	ASSERT_EQ(1, r.countElems());
 }
 
+
 TEST_F(ManifolderTest, closed_surface)
 {
 	Mesh m = buildCubeSurfaceMesh(1.0);
diff --git a/test/cgal/filler/FillerTest.cpp b/test/cgal/filler/FillerTest.cpp
index b2230e6..e8fbfa7 100644
--- a/test/cgal/filler/FillerTest.cpp
+++ b/test/cgal/filler/FillerTest.cpp
@@ -155,6 +155,14 @@ class FillerTest : public ::testing::Test {
         return r;
     }
 
+    static Mesh toAbsolute(const Mesh& m) 
+    {
+        auto r{ m };
+        r.coordinates = 
+            utils::GridTools{ m.grid }.relativeToAbsolute(m.coordinates);
+        return r;
+    }
+
     static bool allAreSimple(const FaceFilling& ff)
     {
         for (const auto [pr, ss] : ff.tris) {
@@ -337,6 +345,59 @@ TEST_F(FillerTest, parallelogram_as_surface)
     EXPECT_EQ(0, countPWHs(f.getFaceFilling({ Cell({0, 0, 1}), Z })));
 }
 
+TEST_F(FillerTest, fill_cube1x1x1_size1_grid)
+{
+    Mesh m = buildCubeSurfaceMesh(01.0);
+        
+    Mesh out;
+    ASSERT_NO_THROW(out = Slicer{m}.getMesh());
+    EXPECT_EQ(12, countMeshElementsIf(out, isTriangle));
+
+    Mesh filled = Filler{out}.getMeshFilling();
+    EXPECT_EQ(12, countMeshElementsIf(filled, isTriangle));
+}
+
+TEST_F(FillerTest, fill_cube1x1x1_size05_grid)
+{
+    //filling with unstruc. triangles
+    Mesh m = buildCubeSurfaceMesh(0.5);
+    Mesh filled = Filler{Slicer{buildCubeSurfaceMesh(0.5) }.getMesh()}.getMeshFilling();
+    Mesh filled_no_slicing = Filler{toRelative(buildCubeSurfaceMesh(0.5))}.getMeshFilling();
+    EXPECT_EQ(18, countMeshElementsIf(filled, isTriangle));
+    EXPECT_EQ(18, countMeshElementsIf(filled_no_slicing, isTriangle));
+        
+    //slicing hull
+    Mesh out_1;
+    ASSERT_NO_THROW(out_1 = Slicer{buildCubeSurfaceMesh(0.5) }.getMesh());
+    EXPECT_EQ(48, countMeshElementsIf(out_1, isTriangle));
+
+    //filling and slicing
+    Mesh out_2;
+    ASSERT_NO_THROW(out_2 = Slicer{toAbsolute(filled) }.getMesh());
+    EXPECT_EQ(72, countMeshElementsIf(out_2, isTriangle));
+
+    ASSERT_NO_THROW(out_2 = Slicer{toAbsolute(filled_no_slicing)}.getMesh());
+    EXPECT_EQ(72, countMeshElementsIf(out_2, isTriangle));
+
+}
+
+TEST_F(FillerTest, fill_cube1x1x1_size025_grid)
+{
+    Mesh m = buildCubeSurfaceMesh(0.25);
+    Mesh filled = Filler{Slicer{buildCubeSurfaceMesh(0.25) }.getMesh()}.getMeshFilling();
+    EXPECT_EQ(30, countMeshElementsIf(filled, isTriangle));
+        
+    Mesh out_1;
+    ASSERT_NO_THROW(out_1 = Slicer{buildCubeSurfaceMesh(0.25) }.getMesh());
+    EXPECT_EQ(192, countMeshElementsIf(out_1, isTriangle));
+
+    Mesh out_2;
+    ASSERT_NO_THROW(out_2 = Slicer{toAbsolute(filled) }.getMesh());
+    EXPECT_EQ(480, countMeshElementsIf(out_2, isTriangle));
+
+}
+
+
 TEST_F(FillerTest, planeXY_mesh_filling)
 {
     Filler f{ Slicer{ buildPlaneXYMesh(1.0) }.getMesh() };
diff --git a/test/core/SlicerTest.cpp b/test/core/SlicerTest.cpp
index 2367b93..70e18b4 100644
--- a/test/core/SlicerTest.cpp
+++ b/test/core/SlicerTest.cpp
@@ -1090,4 +1090,20 @@ TEST_F(SlicerTest, sphere_case_patch_contour_check_2)
     #endif
 }
 
+// TEST_F(SlicerTest, cube1x1x1_size05_grid_fill)
+// {
+//     Mesh m = buildCubeSurfaceMesh(0.5);
+        
+//     Mesh out;
+//     ASSERT_NO_THROW(out = Slicer{m}.getMesh());
+//     EXPECT_EQ(48, countMeshElementsIf(out, isTriangle));
+//     EXPECT_FALSE(containsDegenerateTriangles(out));
+//     EXPECT_EQ(countContours(m), countContours(out));
+    
+//     Mesh filled;
+//     ASSERT_NO_THROW(filled = Slicer::fill(out));
+//     EXPECT_EQ(72, countMeshElementsIf(out, isTriangle));
+// }
+
+
 }
diff --git a/test/meshers/StaircaseMesherTest.cpp b/test/meshers/StaircaseMesherTest.cpp
index 6d1dfdb..334cacd 100644
--- a/test/meshers/StaircaseMesherTest.cpp
+++ b/test/meshers/StaircaseMesherTest.cpp
@@ -369,6 +369,47 @@ TEST_F(StaircaseMesherTest, testStaircaseWithCompression)
     EXPECT_EQ(0, countMeshElementsIf(compressedMesh, isNode));
 }
 
+TEST_F(StaircaseMesherTest, mesh_tetrahedron_volume_2x2){
+
+	Mesh m = buildCubeVolumeMesh(0.5);
+    meshlib::meshers::StaircaseMesherOptions opts;
+    opts.volumeGroups.insert(0);
+    // opts.isVolume = true;
+    auto staircasedMesh = StaircaseMesher{m, 4, opts }.mesh();
+
+    EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTriangle));
+    EXPECT_EQ(36, countMeshElementsIf(staircasedMesh, isQuad));
+    EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTetrahedron));
+
+}
+
+TEST_F(StaircaseMesherTest, mesh_surface_volume_2x2){
+
+	Mesh m = buildCubeSurfaceMesh(0.5);
+    meshlib::meshers::StaircaseMesherOptions opts;
+    opts.volumeGroups.insert(0);
+    // opts.isVolume = true;
+    auto staircasedMesh = StaircaseMesher{m, 4, opts }.mesh();
+
+    EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTriangle));
+    EXPECT_EQ(36, countMeshElementsIf(staircasedMesh, isQuad));
+    EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTetrahedron));
+
+}
+
+TEST_F(StaircaseMesherTest, mesh_surface_not_volume_2x2){
+
+	Mesh m = buildCubeSurfaceMesh(0.5);
+    meshlib::meshers::StaircaseMesherOptions opts;
+    // opts.isVolume = true;
+    auto staircasedMesh = StaircaseMesher{m, 4, opts }.mesh();
+
+    EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTriangle));
+    EXPECT_EQ(24, countMeshElementsIf(staircasedMesh, isQuad));
+    EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTetrahedron));
+
+}
+
 #if APP_LOADED
 
 TEST_F(StaircaseMesherTest, fills_closed_volume_with_quads)
@@ -380,10 +421,10 @@ TEST_F(StaircaseMesherTest, fills_closed_volume_with_quads)
     mesh.grid[Z] = utils::GridTools::linspace(-100.0, 100.0, 51);
 
     meshlib::meshers::StaircaseMesherOptions opts;
-    opts.isVolume = false;
+    // opts.isVolume = false;
     auto staircasedMesh = StaircaseMesher{mesh, 4, opts }.mesh();
     
-    opts.isVolume = true;
+    opts.volumeGroups.insert(0);
     auto staircasedMeshVolume = StaircaseMesher{mesh, 4, opts }.mesh();
     
     EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTriangle));