diff --git a/MeshLib/MeshGenerators/LayeredVolume.cpp b/MeshLib/MeshGenerators/LayeredVolume.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bd4825fa1b6d390777a5c5c376a6e1072b9f82d8
--- /dev/null
+++ b/MeshLib/MeshGenerators/LayeredVolume.cpp
@@ -0,0 +1,219 @@
+/**
+ * \file LayeredVolume.cpp
+ * \author Karsten Rink
+ * \date 2014-04-11
+ * \brief Implementation of the LayeredVolume class.
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "LayeredVolume.h"
+
+#include <fstream>
+#include <numeric>
+
+#include "Vector3.h"
+
+#include "GEOObjects.h"
+#include "PointVec.h"
+#include "Mesh.h"
+#include "convertMeshToGeo.h"
+#include "Elements/Element.h"
+#include "Elements/Tri.h"
+#include "Elements/Quad.h"
+#include "MeshGenerators/MeshLayerMapper.h"
+#include "MeshQuality/MeshValidation.h"
+#include "MeshEditing/ElementExtraction.h"
+
+LayeredVolume::LayeredVolume()
+: _invalid_value(-9999), _mesh(nullptr)
+{
+}
+
+bool LayeredVolume::createGeoVolumes(const MeshLib::Mesh &mesh, const std::vector<std::string> &raster_paths, double noDataReplacementValue)
+{
+ if (mesh.getDimension() != 2 || !allRastersExist(raster_paths))
+ return false;
+
+ std::vector<GeoLib::Raster const*const> rasters;
+ rasters.reserve(raster_paths.size());
+ for (auto path = raster_paths.begin(); path != raster_paths.end(); ++path)
+ rasters.push_back(GeoLib::Raster::getRasterFromASCFile(*path));
+
+ const bool result = this->createGeoVolumes(mesh, rasters, noDataReplacementValue);
+ std::for_each(rasters.begin(), rasters.end(), [](GeoLib::Raster const*const raster){ delete raster; });
+ return result;
+}
+
+
+bool LayeredVolume::createGeoVolumes(const MeshLib::Mesh &mesh, const std::vector<GeoLib::Raster const*const> &rasters, double noDataReplacementValue)
+{
+ if (mesh.getDimension() != 2)
+ return false;
+
+ // remove line elements, only tri + quad remain
+ MeshLib::ElementExtraction ex(mesh);
+ ex.searchByElementType(MeshElemType::LINE);
+ MeshLib::Mesh* mesh_layer (ex.removeMeshElements("MeshLayer"));
+ if (mesh_layer==nullptr)
+ mesh_layer = new MeshLib::Mesh(mesh);
+
+ // map each layer and attach to subsurface mesh
+ const std::size_t nRasters (rasters.size());
+ std::vector<GeoLib::Point> in_region_points(nRasters-1, GeoLib::Point(0,0,0));
+ for (size_t i=0; i<nRasters; ++i)
+ {
+ const double replacement_value = (i==0) ? noDataReplacementValue : _invalid_value;
+ if (!MeshLayerMapper::LayerMapping(*mesh_layer, *rasters[i], 0, 0, _invalid_value))
+ {
+ this->cleanUpOnError();
+ return false;
+ }
+ this->addLayerToMesh(*mesh_layer, i);
+ }
+ // close boundaries between layers
+ this->addLayerBoundaries(*mesh_layer, nRasters);
+ delete mesh_layer;
+
+ this->removeCongruentElements(nRasters, mesh.getNElements());
+ _mesh = new MeshLib::Mesh("BoundaryMesh", _nodes, _elements);
+ MeshLib::MeshValidation::removeUnusedMeshNodes(*_mesh);
+ return true;
+}
+
+void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned layer_id)
+{
+ const std::vector<MeshLib::Node*> &layer_nodes (mesh_layer.getNodes());
+ const std::size_t nNodes (layer_nodes.size());
+ const std::size_t node_id_offset (_nodes.size());
+ const std::size_t last_layer_offset (node_id_offset-nNodes);
+
+ for (std::size_t i=0; i<nNodes; ++i)
+ {
+ if (layer_id > 0 &&
+ ((*layer_nodes[i])[2] == _invalid_value ||
+ (*_nodes[last_layer_offset+i])[2]-(*layer_nodes[i])[2] < 1))
+ _nodes.push_back(new MeshLib::Node(*_nodes[last_layer_offset+i]));
+ else
+ _nodes.push_back(new MeshLib::Node(layer_nodes[i]->getCoords(), _nodes.size()));
+ }
+
+ const std::vector<MeshLib::Element*> &layer_elements (mesh_layer.getElements());
+ for (MeshLib::Element* elem : layer_elements)
+ {
+ if (elem->getGeomType() == MeshElemType::TRIANGLE)
+ {
+ std::array<MeshLib::Node*,3> tri_nodes = { _nodes[node_id_offset+elem->getNode(0)->getID()],
+ _nodes[node_id_offset+elem->getNode(1)->getID()],
+ _nodes[node_id_offset+elem->getNode(2)->getID()] };
+ _elements.push_back(new MeshLib::Tri(tri_nodes, layer_id));
+ }
+ else if (elem->getGeomType() == MeshElemType::QUAD)
+ {
+ std::array<MeshLib::Node*,4> quad_nodes = { _nodes[node_id_offset+elem->getNode(0)->getID()],
+ _nodes[node_id_offset+elem->getNode(1)->getID()],
+ _nodes[node_id_offset+elem->getNode(2)->getID()],
+ _nodes[node_id_offset+elem->getNode(3)->getID()] };
+ _elements.push_back(new MeshLib::Quad(quad_nodes, layer_id));
+ }
+ }
+}
+
+void LayeredVolume::addLayerBoundaries(const MeshLib::Mesh &layer, std::size_t nLayers)
+{
+ const unsigned nLayerBoundaries (nLayers-1);
+ const std::size_t nNodes (layer.getNNodes());
+ const std::vector<MeshLib::Element*> &layer_elements (layer.getElements());
+ for (MeshLib::Element* elem : layer_elements)
+ {
+ const std::size_t nElemNodes (elem->getNNodes());
+ for (unsigned i=0; i<nElemNodes; ++i)
+ if (elem->getNeighbor(i) == nullptr)
+ for (unsigned j=0; j<nLayerBoundaries; ++j)
+ {
+ const std::size_t offset (j*nNodes);
+ MeshLib::Node* n0 = _nodes[offset + elem->getNode(i)->getID()];
+ MeshLib::Node* n1 = _nodes[offset + elem->getNode((i+1)%nElemNodes)->getID()];
+ MeshLib::Node* n2 = _nodes[offset + nNodes + elem->getNode((i+1)%nElemNodes)->getID()];
+ MeshLib::Node* n3 = _nodes[offset + nNodes + elem->getNode(i)->getID()];
+
+ if (MathLib::Vector3(*n1, *n2).getLength() > std::numeric_limits<double>::epsilon())
+ {
+ const std::array<MeshLib::Node*,3> tri_nodes = { n0, n1, n2 };
+ _elements.push_back(new MeshLib::Tri(tri_nodes, nLayers+j));
+ }
+ if (MathLib::Vector3(*n0, *n3).getLength() > std::numeric_limits<double>::epsilon())
+ {
+ const std::array<MeshLib::Node*,3> tri_nodes = { n0, n2, n3 };
+ _elements.push_back(new MeshLib::Tri(tri_nodes, nLayers+j));
+ }
+ }
+ }
+}
+
+void LayeredVolume::removeCongruentElements(std::size_t nLayers, std::size_t nElementsPerLayer)
+{
+ for (std::size_t i=1; i<nLayers; ++i)
+ {
+ const std::size_t upper_offset ((i-1) * nElementsPerLayer);
+ const std::size_t lower_offset ( i * nElementsPerLayer);
+ for (std::size_t j=0; j<nElementsPerLayer; ++j)
+ {
+ MeshLib::Element const*const high (_elements[upper_offset+j]);
+ MeshLib::Element *low (_elements[lower_offset+j]);
+
+ unsigned count(0);
+ const std::size_t nElemNodes (low->getNNodes());
+ for (std::size_t k=0; k<nElemNodes; ++k)
+ if (high->getNode(k)->getID() == low->getNode(k)->getID())
+ {
+ low->setNode(k, _nodes[high->getNodeIndex(k)]);
+ count++;
+ }
+
+ if (count == nElemNodes)
+ {
+ delete _elements[upper_offset+j];
+ _elements[upper_offset+j] = nullptr;
+ }
+ else
+ {
+ MeshLib::Node attr = high->getCenterOfGravity();
+ _attribute_points.push_back(MeshLib::Node(attr[0], attr[1], (attr[2] + low->getCenterOfGravity()[2])/2.0, low->getValue()));
+ }
+ }
+ }
+ auto elem_vec_end = std::remove(_elements.begin(), _elements.end(), nullptr);
+ _elements.erase(elem_vec_end, _elements.end());
+}
+
+bool LayeredVolume::exportToGeometry(GeoLib::GEOObjects &geo_objects) const
+{
+ if (_mesh == nullptr)
+ return false;
+ MeshLib::convertMeshToGeo(*_mesh, geo_objects, std::numeric_limits<double>::min());
+ return true;
+}
+
+bool LayeredVolume::allRastersExist(const std::vector<std::string> &raster_paths) const
+{
+ for (auto raster = raster_paths.begin(); raster != raster_paths.end(); ++raster)
+ {
+ std::ifstream file_stream (*raster, std::ifstream::in);
+ if (!file_stream.good())
+ return false;
+ file_stream.close();
+ }
+ return true;
+}
+
+void LayeredVolume::cleanUpOnError()
+{
+ std::for_each(_nodes.begin(), _nodes.end(), [](MeshLib::Node *node) { delete node; });
+ std::for_each(_elements.begin(), _elements.end(), [](MeshLib::Element *elem) { delete elem; });
+}
diff --git a/MeshLib/MeshGenerators/LayeredVolume.h b/MeshLib/MeshGenerators/LayeredVolume.h
new file mode 100644
index 0000000000000000000000000000000000000000..3ccd09527dd7a0f9f29b0c341349413472c79f39
--- /dev/null
+++ b/MeshLib/MeshGenerators/LayeredVolume.h
@@ -0,0 +1,93 @@
+/**
+ * \file LayeredVolume.h
+ * \author Karsten Rink
+ * \date 2014-04-11
+ * \brief Definition of the LayeredVolume class
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef LAYEREDVOLUME_H
+#define LAYEREDVOLUME_H
+
+#include <string>
+#include <vector>
+
+#include "Raster.h"
+#include "Node.h"
+
+namespace GeoLib {
+ class GEOObjects;
+ class Surface;
+}
+
+namespace MeshLib {
+ class Mesh;
+ class Element;
+}
+
+/**
+ * \brief Creates a volume geometry from 2D mesh layers based on raster data.
+ */
+class LayeredVolume
+{
+public:
+ LayeredVolume();
+ ~LayeredVolume() {}
+
+ /**
+ * Constructs a subsurface representation of a mesh using only 2D elements (i.e. layer boundaries are represented by surfaces)
+ * @param mesh The 2D surface mesh that is used as a basis for the subsurface mesh
+ * @param rasters Containing all the raster-data for the subsurface layers from top to bottom (starting with the DEM)
+ * @param noDataReplacementValue Default z-coordinate if there are mesh nodes not located on the DEM raster (i.e. raster_paths[0])
+ * @result true if the subsurface representation has been created, false if there was an error
+ */
+ bool createGeoVolumes(const MeshLib::Mesh &mesh, const std::vector<GeoLib::Raster const*const> &rasters, double noDataReplacementValue = 0.0);
+
+ /**
+ * Constructs a subsurface representation of a mesh using only 2D elements (i.e. layer boundaries are represented by surfaces)
+ * @param mesh The 2D surface mesh that is used as a basis for the subsurface mesh
+ * @param raster_paths Containing all the raster-file-names for the subsurface layers from top to bottom (starting with the DEM)
+ * @param noDataReplacementValue Default z-coordinate if there are mesh nodes not located on the DEM raster (i.e. raster_paths[0])
+ * @result true if the subsurface representation has been created, false if there was an error
+ */
+ bool createGeoVolumes(const MeshLib::Mesh &mesh, const std::vector<std::string> &raster_paths, double noDataReplacementValue = 0.0);
+
+ /// Returns the subsurface representation that can be used as input for the TetGenInterface
+ MeshLib::Mesh* getMesh() const { return _mesh; }
+
+ /// Returns the region attribute vector necessary for assigning region attributes via TetGen
+ std::vector<MeshLib::Node> getAttributePoints() { return _attribute_points; }
+
+ /// Converts the subsurface representation to geometry objects and adds them to the geo storage object
+ bool exportToGeometry(GeoLib::GEOObjects &geo_objects) const;
+
+private:
+ /// Adds another layer to the subsurface mesh
+ void addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned layer_id);
+
+ /// Creates boundary surfaces between the mapped layers to make the volumes watertight
+ void addLayerBoundaries(const MeshLib::Mesh &layer, std::size_t nLayers);
+
+ /// Removes duplicate 2D elements (possible due to outcroppings)
+ void removeCongruentElements(std::size_t nLayers, std::size_t nElementsPerLayer);
+
+ /// Checks if all raster files actually exist
+ bool allRastersExist(const std::vector<std::string> &raster_paths) const;
+
+ /// Cleans up the already created object in case of an error
+ void cleanUpOnError();
+
+ const double _invalid_value;
+ std::vector<MeshLib::Node*> _nodes;
+ std::vector<MeshLib::Element*> _elements;
+ std::vector<MeshLib::Node> _attribute_points;
+ MeshLib::Mesh* _mesh;
+};
+
+#endif //LAYEREDVOLUME_H