diff --git a/Gui/DataView/MeshLayerEditDialog.cpp b/Gui/DataView/MeshLayerEditDialog.cpp
index 5f7b9461f188714e1450c8ece9f8f82198830efd..3e1d0bb8f40ab1c017dc3e035de8150ef68e3d12 100644
--- a/Gui/DataView/MeshLayerEditDialog.cpp
+++ b/Gui/DataView/MeshLayerEditDialog.cpp
@@ -188,16 +188,19 @@ MeshLib::Mesh* MeshLayerEditDialog::createPrismMesh()
layer_thickness.push_back(thickness);
}
- MeshLayerMapper const mapper;
- MeshLib::Mesh* new_mesh = mapper.createStaticLayers(*_msh, layer_thickness);
+ MeshLayerMapper mapper;
+ MeshLib::Mesh* new_mesh (nullptr);
if (_use_rasters)
{
std::vector<std::string> raster_paths;
for (int i=nLayers; i>=0; --i)
raster_paths.push_back(this->_edits[i]->text().toStdString());
- new_mesh = mapper.createRasterLayers(*_msh, raster_paths);
+ if (mapper.createLayers(*_msh, raster_paths))
+ new_mesh= mapper.getMesh("SubsurfaceMesh");
}
+ else
+ new_mesh = mapper.createStaticLayers(*_msh, layer_thickness);
return new_mesh;
}
@@ -218,9 +221,9 @@ MeshLib::Mesh* MeshLayerEditDialog::createTetMesh()
for (int i=nLayers; i>=0; --i)
raster_paths[i] = this->_edits[i]->text().toStdString();
LayeredVolume lv;
- lv.createGeoVolumes(*_msh, raster_paths);
+ lv.createLayers(*_msh, raster_paths);
- tg_mesh = lv.getMesh();
+ tg_mesh = lv.getMesh("SubsurfaceMesh");
QString file_path("");
if (tg_mesh)
diff --git a/MeshLib/MeshGenerators/LayeredVolume.cpp b/MeshLib/MeshGenerators/LayeredVolume.cpp
index 2f8033302b3f90319ee9a98f95b275653a995aa1..509d314c23f2c9a2baa1eebd21ef55acbbcac944 100644
--- a/MeshLib/MeshGenerators/LayeredVolume.cpp
+++ b/MeshLib/MeshGenerators/LayeredVolume.cpp
@@ -14,49 +14,17 @@
#include "LayeredVolume.h"
-#include <fstream>
-#include <numeric>
-
#include "Vector3.h"
-#include "FileIO/AsciiRasterInterface.h"
+#include "Raster.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"
-const double LayeredVolume::_invalid_value = -9999;
-
-LayeredVolume::LayeredVolume()
-: _elevation_epsilon(0.0001), _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*> rasters;
- rasters.reserve(raster_paths.size());
- for (auto path = raster_paths.begin(); path != raster_paths.end(); ++path)
- rasters.push_back(FileIO::AsciiRasterInterface::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*> &rasters, double noDataReplacementValue)
+bool LayeredVolume::createRasterLayers(const MeshLib::Mesh &mesh, const std::vector<GeoLib::Raster const*> &rasters, double noDataReplacementValue)
{
if (mesh.getDimension() != 2)
return false;
@@ -73,38 +41,36 @@ bool LayeredVolume::createGeoVolumes(const MeshLib::Mesh &mesh, const std::vecto
mesh_layer = new MeshLib::Mesh(mesh);
// map each layer and attach to subsurface mesh
- MeshLayerMapper const mapper;
const std::size_t nRasters (rasters.size());
for (size_t i=0; i<nRasters; ++i)
{
- const double replacement_value = (i==(nRasters-1)) ? noDataReplacementValue : _invalid_value;
- if (!mapper.layerMapping(*mesh_layer, *rasters[i], replacement_value))
+ const double replacement_value = (i==(nRasters-1)) ? noDataReplacementValue : rasters[i]->getNoDataValue();
+ if (!MeshLayerMapper::layerMapping(*mesh_layer, *rasters[i], replacement_value))
{
this->cleanUpOnError();
return false;
}
- this->addLayerToMesh(*mesh_layer, i);
+ this->addLayerToMesh(*mesh_layer, i, *rasters[i]);
}
// close boundaries between layers
this->addLayerBoundaries(*mesh_layer, nRasters);
this->removeCongruentElements(nRasters, mesh_layer->getNElements());
delete mesh_layer;
- _mesh = new MeshLib::Mesh("BoundaryMesh", _nodes, _elements);
- MeshLib::MeshValidation::removeUnusedMeshNodes(*_mesh);
return true;
}
-void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned layer_id)
+void LayeredVolume::addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned layer_id, GeoLib::Raster const& raster)
{
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);
+ const double no_data_value (raster.getNoDataValue());
for (std::size_t i=0; i<nNodes; ++i)
{
if (layer_id > 0 &&
- ((*layer_nodes[i])[2] == _invalid_value ||
+ ((*layer_nodes[i])[2] == no_data_value ||
(*_nodes[last_layer_offset+i])[2]-(*layer_nodes[i])[2] < _elevation_epsilon))
_nodes.push_back(new MeshLib::Node(*_nodes[last_layer_offset+i]));
else
@@ -199,36 +165,3 @@ void LayeredVolume::removeCongruentElements(std::size_t nLayers, std::size_t nEl
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;
-}
-
-double LayeredVolume::calcEpsilon(const GeoLib::Raster &high, const GeoLib::Raster &low)
-{
- const double max (*std::max_element(high.begin(), high.end()));
- const double min (*std::min_element( low.begin(), low.end()));
- return ((max-min)*1e-07);
-}
-
-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
index 1801f81c90f2bb4925bfe99c1d66e9c401c5a5da..172503a64bf3f7bf510061fd38cf91182e69a167 100644
--- a/MeshLib/MeshGenerators/LayeredVolume.h
+++ b/MeshLib/MeshGenerators/LayeredVolume.h
@@ -18,26 +18,21 @@
#include <string>
#include <vector>
-#include "Raster.h"
#include "Node.h"
+#include "SubsurfaceMapper.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
+class LayeredVolume : public SubsurfaceMapper
{
public:
- LayeredVolume();
+ LayeredVolume() {}
~LayeredVolume() {}
/**
@@ -47,29 +42,14 @@ public:
* @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*> &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 bottom to top (starting with the bottom of the oldest layer and ending 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; }
-
+ bool createRasterLayers(const MeshLib::Mesh &mesh, const std::vector<GeoLib::Raster const*> &rasters, double noDataReplacementValue = 0.0);
+
/// 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);
+ void addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned layer_id, GeoLib::Raster const& raster);
/// Creates boundary surfaces between the mapped layers to make the volumes watertight
void addLayerBoundaries(const MeshLib::Mesh &layer, std::size_t nLayers);
@@ -77,21 +57,7 @@ private:
/// Removes duplicate 2D elements (possible due to outcroppings)
void removeCongruentElements(std::size_t nLayers, std::size_t nElementsPerLayer);
- /// Calculates a data-dependent epsilon value
- double calcEpsilon(const GeoLib::Raster &high, const GeoLib::Raster &low);
-
- /// 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();
-
- static const double _invalid_value;
- double _elevation_epsilon;
- std::vector<MeshLib::Node*> _nodes;
- std::vector<MeshLib::Element*> _elements;
std::vector<MeshLib::Node> _attribute_points;
- MeshLib::Mesh* _mesh;
};
#endif //LAYEREDVOLUME_H
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.cpp b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
index 20eb5c939241d81a88c93a603305af7ed6264a62..2264191178d387fcb7f22d1fadb08d1b4ebbb828 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.cpp
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
@@ -12,28 +12,21 @@
*
*/
-// stl
+#include "MeshLayerMapper.h"
+
#include <algorithm>
-// ThirdParty/logog
#include "logog/include/logog.hpp"
-#include "MeshLayerMapper.h"
-
#include "FileIO/AsciiRasterInterface.h"
-// GeoLib
#include "Raster.h"
-#include "Mesh.h"
-#include "Node.h"
-#include "Elements/Element.h"
#include "Elements/Tet.h"
#include "Elements/Hex.h"
#include "Elements/Pyramid.h"
#include "Elements/Prism.h"
#include "MeshSurfaceExtraction.h"
-#include "MeshQuality/MeshValidation.h"
#include "MathTools.h"
@@ -44,7 +37,6 @@ const unsigned MeshLayerMapper::_pyramid_base[3][4] =
{0, 3, 4, 1}, // Point 6 missing
};
-
MeshLib::Mesh* MeshLayerMapper::createStaticLayers(MeshLib::Mesh const& mesh, std::vector<float> const& layer_thickness_vector, std::string const& mesh_name) const
{
std::vector<float> thickness;
@@ -114,27 +106,12 @@ MeshLib::Mesh* MeshLayerMapper::createStaticLayers(MeshLib::Mesh const& mesh, st
return new MeshLib::Mesh(mesh_name, new_nodes, new_elems);
}
-MeshLib::Mesh* MeshLayerMapper::createRasterLayers(MeshLib::Mesh const& mesh, std::vector<std::string> const& raster_paths, std::string const& mesh_name) const
-{
- if (mesh.getDimension() != 2 || !allRastersExist(raster_paths))
- return false;
-
- std::vector<GeoLib::Raster const*> rasters;
- rasters.reserve(raster_paths.size());
- for (auto path = raster_paths.begin(); path != raster_paths.end(); ++path)
- rasters.push_back(FileIO::AsciiRasterInterface::getRasterFromASCFile(*path));
-
- MeshLib::Mesh* result = this->createRasterLayers(mesh, rasters, mesh_name);
- std::for_each(rasters.begin(), rasters.end(), [](GeoLib::Raster const*const raster){ delete raster; });
- return result;
-}
-
-MeshLib::Mesh* MeshLayerMapper::createRasterLayers(MeshLib::Mesh const& mesh, std::vector<GeoLib::Raster const*> const& rasters, std::string const& mesh_name) const
+bool MeshLayerMapper::createRasterLayers(MeshLib::Mesh const& mesh, std::vector<GeoLib::Raster const*> const& rasters, double noDataReplacementValue)
{
const std::size_t nLayers(rasters.size());
if (nLayers < 1 || mesh.getDimension() != 2)
{
- ERR("MeshLayerMapper::createStaticLayers(): A 2D mesh with nLayers > 0 is required as input.");
+ ERR("MeshLayerMapper::createRasterLayers(): A 2D mesh and at least two rasters required as input.");
return nullptr;
}
@@ -142,26 +119,22 @@ MeshLib::Mesh* MeshLayerMapper::createRasterLayers(MeshLib::Mesh const& mesh, st
if (layerMapping(*dem_mesh, *rasters.back(), 0))
{
std::size_t const nNodes = mesh.getNNodes();
- std::vector<MeshLib::Node*> new_nodes;
- new_nodes.reserve(nLayers * nNodes);
+ _nodes.reserve(nLayers * nNodes);
// number of triangles in the original mesh
std::size_t const nElems (std::count_if(mesh.getElements().begin(), mesh.getElements().end(),
[](MeshLib::Element const* elem) { return (elem->getGeomType() == MeshElemType::TRIANGLE);}));
- std::vector<MeshLib::Element*> new_elems;
- new_elems.reserve(nElems * (nLayers-1));
+ _elements.reserve(nElems * (nLayers-1));
for (std::size_t i=0; i<nLayers; ++i)
- addLayerToMesh(*dem_mesh, i, *rasters[i], new_nodes, new_elems);
-
- MeshLib::Mesh* result = new MeshLib::Mesh(mesh_name, new_nodes, new_elems);
- MeshLib::MeshValidation::removeUnusedMeshNodes(*result);
- return result;
+ addLayerToMesh(*dem_mesh, i, *rasters[i]);
+ return true;
}
+ return false;
}
-void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer_id, GeoLib::Raster const& raster, std::vector<MeshLib::Node*> &new_nodes, std::vector<MeshLib::Element*> &new_elems) const
+void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned layer_id, GeoLib::Raster const& raster)
{
std::size_t const nNodes = dem_mesh.getNNodes();
std::vector<MeshLib::Node*> const& nodes = dem_mesh.getNodes();
@@ -177,10 +150,10 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay
if ((layer_id > 0) &&
((std::abs(elevation - no_data_value) < std::numeric_limits<double>::epsilon()) ||
- (elevation - (*new_nodes[last_layer_node_offset + i])[2] < std::numeric_limits<double>::epsilon())))
- new_nodes.push_back(new MeshLib::Node((*nodes[i])[0], (*nodes[i])[1], elevation, new_nodes[last_layer_node_offset +i]->getID()));
+ (elevation - (*_nodes[last_layer_node_offset + i])[2] < std::numeric_limits<double>::epsilon())))
+ _nodes.push_back(new MeshLib::Node((*nodes[i])[0], (*nodes[i])[1], elevation, _nodes[last_layer_node_offset +i]->getID()));
else
- new_nodes.push_back(new MeshLib::Node((*nodes[i])[0], (*nodes[i])[1], elevation, (layer_id * nNodes) + i));
+ _nodes.push_back(new MeshLib::Node((*nodes[i])[0], (*nodes[i])[1], elevation, (layer_id * nNodes) + i));
}
if (layer_id == 0)
@@ -198,8 +171,8 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay
std::array<MeshLib::Node*, 6> new_elem_nodes;
for (unsigned j=0; j<3; ++j)
{
- new_elem_nodes[j] = new_nodes[last_layer_node_offset + elem->getNodeIndex(j)];
- new_elem_nodes[node_counter] = (new_nodes[last_layer_node_offset + elem->getNodeIndex(j) + nNodes]);
+ new_elem_nodes[j] = _nodes[last_layer_node_offset + elem->getNodeIndex(j)];
+ new_elem_nodes[node_counter] = (_nodes[last_layer_node_offset + elem->getNodeIndex(j) + nNodes]);
if (new_elem_nodes[j]->getID() != new_elem_nodes[node_counter]->getID())
node_counter++;
else
@@ -209,7 +182,7 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay
switch (node_counter)
{
case 6:
- new_elems.push_back(new MeshLib::Prism(new_elem_nodes, layer_id));
+ _elements.push_back(new MeshLib::Prism(new_elem_nodes, layer_id));
break;
case 5:
std::array<MeshLib::Node*, 5> pyramid_nodes;
@@ -218,12 +191,12 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay
pyramid_nodes[2] = new_elem_nodes[_pyramid_base[missing_idx][2]];
pyramid_nodes[3] = new_elem_nodes[_pyramid_base[missing_idx][3]];
pyramid_nodes[4] = new_elem_nodes[missing_idx];
- new_elems.push_back(new MeshLib::Pyramid(pyramid_nodes, layer_id));
+ _elements.push_back(new MeshLib::Pyramid(pyramid_nodes, layer_id));
break;
case 4:
std::array<MeshLib::Node*, 4> tet_nodes;
std::copy(new_elem_nodes.begin(), new_elem_nodes.begin() + node_counter, tet_nodes.begin());
- new_elems.push_back(new MeshLib::Tet(tet_nodes, layer_id));
+ _elements.push_back(new MeshLib::Tet(tet_nodes, layer_id));
break;
default:
continue;
@@ -231,7 +204,7 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay
}
}
-bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, std::string const& rasterfile, double noDataReplacementValue = 0.0) const
+bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, std::string const& rasterfile, double noDataReplacementValue = 0.0)
{
const GeoLib::Raster *raster(FileIO::AsciiRasterInterface::getRasterFromASCFile(rasterfile));
if (! raster) {
@@ -243,7 +216,7 @@ bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, std::string const& r
return result;
}
-bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, GeoLib::Raster const& raster, double noDataReplacementValue = 0.0) const
+bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, GeoLib::Raster const& raster, double noDataReplacementValue = 0.0)
{
if (new_mesh.getDimension() != 2)
{
@@ -281,16 +254,3 @@ bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, GeoLib::Raster const
return true;
}
-
-bool MeshLayerMapper::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;
-}
-
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.h b/MeshLib/MeshGenerators/MeshLayerMapper.h
index 6eeb1bc1202d7497b7ee6570acdf2eecb0c3c05c..fa164ac36b87a1474a3ccb0861d5cf8dc63ec3ed 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.h
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.h
@@ -15,25 +15,16 @@
#ifndef MESHLAYERMAPPER_H
#define MESHLAYERMAPPER_H
-#include <string>
-#include "GeoLib/Raster.h"
-
-class QImage;
-
-namespace MeshLib {
- class Mesh;
- class Node;
- class Element;
-}
+#include "SubsurfaceMapper.h"
/**
-* \brief Manipulating and adding layers to an existing mesh
-*/
-class MeshLayerMapper
+ * \brief Manipulating and adding prism element layers to an existing 2D mesh
+ */
+class MeshLayerMapper : public SubsurfaceMapper
{
public:
- MeshLayerMapper() {};
- ~MeshLayerMapper() {};
+ MeshLayerMapper() {}
+ ~MeshLayerMapper() {}
/**
* Based on a 2D triangle-or quad mesh this method creates a 3D mesh with with a given number of prism- or hex-layers
@@ -47,36 +38,30 @@ public:
/**
* Based on a 2D triangle mesh this method creates a 3D mesh with with a given number of prism-layers.
- * Note: This method would theoretically also work with quad meshes. However, this is discouraged as quad elements will most likely not
- * be coplanar after the mapping process and result in invaled elements.
+ * Note: While this method would technically also work with quad meshes, this is discouraged as quad elements will most likely not
+ * be coplanar after the mapping process which result in invaled mesh elements.
* \param mesh The 2D triangle mesh that is the basis for the new 3D prism mesh
- * \param nLayers The number of layers of prism/hex elements that will be extruded from the triangle/quad elements of the original mesh
- * \param thickness The thickness of each of these newly added layers
- * \param mesh_name The name of the newly created mesh
- * \return A mesh with the requested number of layers of prism/hex elements
+ * \param rasters Containing all the raster-data for the subsurface layers from bottom to top (starting with the bottom of the oldest layer and ending with the DEM)
+ * \param noDataReplacementValue Default z-coordinate if there are mesh nodes not located on the DEM raster (i.e. raster_paths[0])
+ * \return A mesh with the requested number of layers of prism elements (also including Tet- & Pyramid-elements in case of degenerated prisms)
*/
- MeshLib::Mesh* createRasterLayers(MeshLib::Mesh const& mesh, std::vector<std::string> const& raster_paths, std::string const& mesh_name = "SubsurfaceMesh") const;
-
- MeshLib::Mesh* createRasterLayers(MeshLib::Mesh const& mesh, std::vector<GeoLib::Raster const*> const& rasters, std::string const& mesh_name = "SubsurfaceMesh") const;
+ bool createRasterLayers(MeshLib::Mesh const& mesh, std::vector<GeoLib::Raster const*> const& rasters, double noDataReplacementValue = 0.0);
/**
* Maps the elevation of nodes of a given 2D mesh according to the raster specified by the file path.
* At locations wher no information is given, node elevation is set to noDataReplacementValue.
*/
- bool layerMapping(MeshLib::Mesh &mesh, const std::string &rasterfile, double noDataReplacementValue) const;
+ static bool layerMapping(MeshLib::Mesh &mesh, const std::string &rasterfile, double noDataReplacementValue);
/**
* Maps the elevation of nodes of a given 2D mesh according to the raster. At locations wher no
* information is given, node elevation is set to noDataReplacementValue.
*/
- bool layerMapping(MeshLib::Mesh &mesh, const GeoLib::Raster &raster, double noDataReplacementValue) const;
+ static bool layerMapping(MeshLib::Mesh &mesh, const GeoLib::Raster &raster, double noDataReplacementValue);
private:
/// Adds another layer to a subsurface mesh
- void addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned layer_id, GeoLib::Raster const& raster, std::vector<MeshLib::Node*> &new_nodes, std::vector<MeshLib::Element*> &new_elems) const;
-
- /// Checks if all raster files actually exist
- bool allRastersExist(const std::vector<std::string> &raster_paths) const;
+ void addLayerToMesh(const MeshLib::Mesh &mesh_layer, unsigned layer_id, GeoLib::Raster const& raster);
static const unsigned _pyramid_base[3][4];
};
diff --git a/MeshLib/MeshGenerators/SubsurfaceMapper.cpp b/MeshLib/MeshGenerators/SubsurfaceMapper.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2d5a223fc12f87b90b39d040b05d9f7691e76aeb
--- /dev/null
+++ b/MeshLib/MeshGenerators/SubsurfaceMapper.cpp
@@ -0,0 +1,82 @@
+/**
+ * \file SubsurfaceMapper.cpp
+ * \author Karsten Rink
+ * \date 2014-09-18
+ * \brief Implementation of the SubsurfaceMapper class.
+ *
+ * \copyright
+ * Copyright (c) 2012-2014, 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 "SubsurfaceMapper.h"
+
+#include <vector>
+#include <fstream>
+
+#include "FileIO/AsciiRasterInterface.h"
+
+#include "Raster.h"
+
+#include "Mesh.h"
+#include "Node.h"
+#include "Elements/Element.h"
+#include "MeshQuality/MeshValidation.h"
+
+SubsurfaceMapper::SubsurfaceMapper()
+: _elevation_epsilon(0.0001)
+{
+}
+
+bool SubsurfaceMapper::createLayers(MeshLib::Mesh const& mesh, std::vector<std::string> const& raster_paths, double noDataReplacementValue)
+{
+ if (mesh.getDimension() != 2 || !allRastersExist(raster_paths))
+ return false;
+
+ std::vector<GeoLib::Raster const*> rasters;
+ rasters.reserve(raster_paths.size());
+ for (auto path = raster_paths.begin(); path != raster_paths.end(); ++path)
+ rasters.push_back(FileIO::AsciiRasterInterface::getRasterFromASCFile(*path));
+
+ bool result = createRasterLayers(mesh, rasters, noDataReplacementValue);
+ std::for_each(rasters.begin(), rasters.end(), [](GeoLib::Raster const*const raster){ delete raster; });
+ return result;
+}
+
+MeshLib::Mesh* SubsurfaceMapper::getMesh(std::string const& mesh_name) const
+{
+ if (_nodes.empty() || _elements.empty())
+ return nullptr;
+
+ MeshLib::Mesh* result (new MeshLib::Mesh(mesh_name, _nodes, _elements));
+ MeshLib::MeshValidation::removeUnusedMeshNodes(*result);
+ return result;
+}
+
+double SubsurfaceMapper::calcEpsilon(GeoLib::Raster const& high, GeoLib::Raster const& low)
+{
+ const double max (*std::max_element(high.begin(), high.end()));
+ const double min (*std::min_element( low.begin(), low.end()));
+ return ((max-min)*1e-07);
+}
+
+bool SubsurfaceMapper::allRastersExist(std::vector<std::string> const& 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 SubsurfaceMapper::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/SubsurfaceMapper.h b/MeshLib/MeshGenerators/SubsurfaceMapper.h
new file mode 100644
index 0000000000000000000000000000000000000000..ea29d23afc160b0fcb043027b1550231f0558e12
--- /dev/null
+++ b/MeshLib/MeshGenerators/SubsurfaceMapper.h
@@ -0,0 +1,79 @@
+/**
+ * \file SubsurfaceMapper.h
+ * \author Karsten Rink
+ * \date 2014-09-18
+ * \brief Definition of the SubsurfaceMapper class
+ *
+ * \copyright
+ * Copyright (c) 2012-2014, 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 SUBSURFACEMAPPER_H
+#define SUBSURFACEMAPPER_H
+
+#include <string>
+#include <vector>
+
+namespace GeoLib {
+ class Raster;
+}
+
+namespace MeshLib {
+ class Mesh;
+ class Node;
+ class Element;
+}
+
+/**
+ * \brief Base class for creation of 3D subsurface meshes based on raster data
+ */
+class SubsurfaceMapper
+{
+public:
+ /**
+ * Returns a subsurface representation of a region represented by a 2D by reading raster files and calling the appropriate construction method.
+ * @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 bottom to top (starting with the bottom of the oldest layer and ending 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
+ */
+ virtual bool createLayers(MeshLib::Mesh const& mesh, std::vector<std::string> const& raster_paths, double noDataReplacementValue = 0.0) final;
+
+ /**
+ * Constructs a subsurface representation based on a 2D mesh and a number of rasters representing subsurface layer boundaries.
+ * @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 bottom to top (starting with the bottom of the oldest layer and ending 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
+ */
+ virtual bool createRasterLayers(MeshLib::Mesh const& mesh, std::vector<GeoLib::Raster const*> const& rasters, double noDataReplacementValue) = 0;
+
+ /// Returns a mesh of the subsurface representation
+ MeshLib::Mesh* getMesh(std::string const& mesh_name) const;
+
+protected:
+ SubsurfaceMapper();
+ ~SubsurfaceMapper() {}
+
+ /// Adds another layer to the subsurface mesh
+ virtual void addLayerToMesh(MeshLib::Mesh const& mesh_layer, unsigned layer_id, GeoLib::Raster const& raster) = 0;
+
+ /// Calculates a data-dependent epsilon value
+ double calcEpsilon(GeoLib::Raster const& high, GeoLib::Raster const& low);
+
+ /// Checks if all raster files actually exist
+ bool allRastersExist(std::vector<std::string> const& raster_paths) const;
+
+ /// Cleans up the already created objects in case of an error
+ void cleanUpOnError();
+
+ double _elevation_epsilon;
+ std::vector<MeshLib::Node*> _nodes;
+ std::vector<MeshLib::Element*> _elements;
+};
+
+#endif //SUBSURFACEMAPPER_H
diff --git a/MeshLib/Node.h b/MeshLib/Node.h
index a59a637036c6235fcc4fa59eadd0997bdc7fefb5..527066b0d64044804a77353aedc7769f3483d50b 100644
--- a/MeshLib/Node.h
+++ b/MeshLib/Node.h
@@ -36,7 +36,7 @@ class Element;
class Node : public GeoLib::PointWithID
{
/* friend functions: */
- friend bool MeshLayerMapper::layerMapping(MeshLib::Mesh &mesh, const GeoLib::Raster &raster, double noDataReplacementValue) const;
+ friend static bool MeshLayerMapper::layerMapping(MeshLib::Mesh &mesh, const GeoLib::Raster &raster, double noDataReplacementValue);
friend MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, bool keep3dMeshIds);
/* friend classes: */