diff --git a/MeshLib/MeshGenerators/LayeredVolume.cpp b/MeshLib/MeshGenerators/LayeredVolume.cpp
index 4866b8222a067f8f6fa200b08a475d51f7c34913..2f8033302b3f90319ee9a98f95b275653a995aa1 100644
--- a/MeshLib/MeshGenerators/LayeredVolume.cpp
+++ b/MeshLib/MeshGenerators/LayeredVolume.cpp
@@ -78,7 +78,7 @@ bool LayeredVolume::createGeoVolumes(const MeshLib::Mesh &mesh, const std::vecto
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], 0, 0, replacement_value))
+ if (!mapper.layerMapping(*mesh_layer, *rasters[i], replacement_value))
{
this->cleanUpOnError();
return false;
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.cpp b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
index cb16ee54144a2421691f18e6380ac60dec8736a6..20eb5c939241d81a88c93a603305af7ed6264a62 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.cpp
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
@@ -139,7 +139,7 @@ MeshLib::Mesh* MeshLayerMapper::createRasterLayers(MeshLib::Mesh const& mesh, st
}
MeshLib::Mesh* dem_mesh (new MeshLib::Mesh(mesh));
- if (layerMapping(*dem_mesh, *rasters.back(), 0, 0, 0))
+ if (layerMapping(*dem_mesh, *rasters.back(), 0))
{
std::size_t const nNodes = mesh.getNNodes();
std::vector<MeshLib::Node*> new_nodes;
@@ -167,6 +167,7 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay
std::vector<MeshLib::Node*> const& nodes = dem_mesh.getNodes();
int const last_layer_node_offset = (layer_id-1) * nNodes;
unsigned const layer_node_offset = layer_id * nNodes;
+ double const no_data_value (raster.getNoDataValue());
// add nodes for new layer
for (std::size_t i=0; i<nNodes; ++i)
@@ -174,11 +175,12 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay
// min of dem elevation and layer elevation
double const elevation = std::min(raster.interpolateValueAtPoint(*nodes[i]), (*nodes[i])[2]);
- if ((layer_id == 0) ||
- (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, (layer_id * nNodes) + i));
- else
+ 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()));
+ else
+ new_nodes.push_back(new MeshLib::Node((*nodes[i])[0], (*nodes[i])[1], elevation, (layer_id * nNodes) + i));
}
if (layer_id == 0)
@@ -229,27 +231,24 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay
}
}
-
-bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, const std::string &rasterfile,
- const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue = 0.0) const
+bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, std::string const& rasterfile, double noDataReplacementValue = 0.0) const
{
const GeoLib::Raster *raster(FileIO::AsciiRasterInterface::getRasterFromASCFile(rasterfile));
if (! raster) {
- ERR("MshLayerMapper::LayerMapping - could not read raster file %s", rasterfile.c_str());
+ ERR("MshLayerMapper::layerMapping - could not read raster file %s", rasterfile.c_str());
return false;
}
- const bool result = layerMapping(new_mesh, *raster, nLayers, layer_id, noDataReplacementValue);
+ const bool result = layerMapping(new_mesh, *raster, noDataReplacementValue);
delete raster;
return result;
}
-bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, const GeoLib::Raster &raster,
- const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue = 0.0) const
+bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, GeoLib::Raster const& raster, double noDataReplacementValue = 0.0) const
{
- if (nLayers < layer_id)
- {
- ERR("MshLayerMapper::LayerMapping() - Mesh has only %d Layers, cannot assign layer %d.", nLayers, layer_id);
- return false;
+ if (new_mesh.getDimension() != 2)
+ {
+ ERR("MshLayerMapper::layerMapping - requires 2D mesh");
+ return false;
}
const double x0(raster.getOrigin()[0]);
@@ -262,14 +261,10 @@ bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, const GeoLib::Raster
const std::pair<double, double> yDim(y0, y0 + raster.getNRows() * delta); // extension in y-dimension
const size_t nNodes (new_mesh.getNNodes());
- const size_t nNodesPerLayer (nNodes / (nLayers+1));
-
- const size_t firstNode (layer_id * nNodesPerLayer);
- const size_t lastNode (firstNode + nNodesPerLayer);
-
+
const double half_delta = 0.5*delta;
const std::vector<MeshLib::Node*> &nodes = new_mesh.getNodes();
- for (unsigned i = firstNode; i < lastNode; ++i)
+ for (unsigned i = 0; i < nNodes; ++i)
{
if (!raster.isPntOnRaster(*nodes[i]))
{
@@ -287,151 +282,6 @@ bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, const GeoLib::Raster
return true;
}
-MeshLib::Mesh* MeshLayerMapper::blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster) const
-{
- // construct surface mesh from DEM
- const MathLib::Vector3 dir(0,0,1);
- MeshLib::Mesh* dem = MeshLib::MeshSurfaceExtraction::getMeshSurface(mesh, dir);
- this->layerMapping(*dem, dem_raster, 0, 0);
- const std::vector<MeshLib::Node*> &dem_nodes (dem->getNodes());
-
- const std::vector<MeshLib::Node*> &mdl_nodes (mesh.getNodes());
- const size_t nNodes (mesh.getNNodes());
- const size_t nNodesPerLayer = nNodes / (nLayers+1);
- std::vector<bool> is_surface_node(nNodes, false);
- std::vector<bool> nodes_below_surface(nNodes, false);
-
- // check if bottom layer nodes are below DEM
- const unsigned bottom_firstNode = nLayers * nNodesPerLayer;
- const unsigned bottom_lastNode = bottom_firstNode + nNodesPerLayer;
- for (unsigned i = bottom_firstNode; i < bottom_lastNode; ++i)
- {
- nodes_below_surface[i]=true;
- const double* coords = mdl_nodes[i]->getCoords();
- const double* dem_coords = dem_nodes[i-bottom_firstNode]->getCoords();
- if (coords[2] >= dem_coords[2])
- {
- WARN("Node %d (in bottom-layer) is above surface node %d. (%f, %f)", i, (i-bottom_firstNode), coords[2], dem_coords[2]);
- is_surface_node[i] = true;
- }
- }
-
- // for all other layers:
- // if node < dem-node: do nothing
- // if node > dem-node:
- // if first node above surface: map to dem and mark as surface node
- // else remove node (i.e. don't copy it)
- for (int layer_id=nLayers-1; layer_id>=0; --layer_id)
- {
- const size_t firstNode = layer_id * nNodesPerLayer;
- const size_t lastNode = firstNode + nNodesPerLayer;
-
- for(unsigned i = firstNode; i < lastNode; ++i)
- {
- if (is_surface_node[i+nNodesPerLayer])
- is_surface_node[i]=true;
- else
- {
- nodes_below_surface[i]=true;
- MeshLib::Node* node (mdl_nodes[i]);
- const double* coords = node->getCoords();
- const double* dem_coords = dem_nodes[i-firstNode]->getCoords();
- if (coords[2] > dem_coords[2])
- {
- node->updateCoordinates(dem_coords[0], dem_coords[1], dem_coords[2]);
- is_surface_node[i] = true;
- }
- }
- }
- }
- delete dem; // no longer needed
-
- // copy valid nodes to new node vector
- std::vector<MeshLib::Node*> new_nodes;
- std::vector<int> node_index_map(nNodes, -1);
- size_t node_count(0);
- for (unsigned j=0; j<nNodes; ++j)
- if (nodes_below_surface[j])
- {
- new_nodes.push_back(new MeshLib::Node(mdl_nodes[j]->getCoords(), mdl_nodes[j]->getID()));
- node_index_map[j]=node_count++;
- }
-
- // copy elements (old elements need to have at least 4 nodes remaining and form a 3d element
- const std::vector<MeshLib::Element*> &mdl_elements (mesh.getElements());
- const size_t nElems (mesh.getNElements());
- std::vector<MeshLib::Element*> new_elements;
- for (unsigned j=0; j<nElems; ++j)
- {
- const MeshLib::Element* elem = mdl_elements[j];
-
- size_t count(0);
- for (unsigned i=0; i<6; ++i) // check top surface of prism
- if (nodes_below_surface[elem->getNode(i)->getID()]) ++count;
-
- if (count==6) // copy prism elements if all six nodes are valid
- {
- MeshLib::Node** e_nodes = new MeshLib::Node*[count];
- for (unsigned i=0; i<6; ++i)
- e_nodes[i] = new_nodes[node_index_map[elem->getNode(i)->getID()]];
-
- MeshLib::Element* prism (new MeshLib::Prism(e_nodes, elem->getValue()));
- new_elements.push_back(prism);
- }
- else if (count==5) // change the current element to two tetrahedra if only five nodes are valid
- {
- MeshLib::Node** e_nodes = new MeshLib::Node*[count];
- unsigned top_idx(6);
- for (unsigned i=3; i<6; ++i) // find node that has been cut
- if (!nodes_below_surface[elem->getNode(i)->getID()])
- top_idx = i-3;
-
- // construct pyramid element based on missing node
- unsigned idx1 ((top_idx+1)%3);
- unsigned idx2 ((top_idx+2)%3);
- e_nodes[0] = new_nodes[node_index_map[elem->getNode(idx1)->getID()]];
- e_nodes[1] = new_nodes[node_index_map[elem->getNode(idx1+3)->getID()]];
- e_nodes[2] = new_nodes[node_index_map[elem->getNode(idx2+3)->getID()]];
- e_nodes[3] = new_nodes[node_index_map[elem->getNode(idx2)->getID()]];
- e_nodes[4] = new_nodes[node_index_map[elem->getNode(top_idx)->getID()]];
-
- MeshLib::Element* pyr (new MeshLib::Pyramid(e_nodes, elem->getValue()));
- new_elements.push_back(pyr);
- }
- else if (count==4) // change the current element to a tetrahedron if only four nodes are valid
- {
- MeshLib::Node** e_nodes = new MeshLib::Node*[count];
- for (unsigned i=0; i<3; ++i) // first three nodes are the bottom-face
- {
- unsigned idx (elem->getNode(i)->getID());
- if (nodes_below_surface[idx])
- e_nodes[i] = new_nodes[node_index_map[idx]];
- else
- e_nodes[i] = NULL;
- }
-
- if (e_nodes[0] && e_nodes[1] && e_nodes[2]) //make sure that the 4 remaining nodes don't form a quad
- {
- for (unsigned i=3; i<6; ++i) // last node
- {
- unsigned idx (elem->getNode(i)->getID());
- if (nodes_below_surface[idx])
- {
- e_nodes[3] = new_nodes[node_index_map[idx]];
- break;
- }
- }
-
- MeshLib::Element* tet (new MeshLib::Tet(e_nodes, elem->getValue()));
- new_elements.push_back(tet);
- }
- else delete e_nodes;
- }
- // else remove element, if less than four nodes are valid
- }
- return new MeshLib::Mesh("SubsurfaceMesh", new_nodes, new_elements);
-}
-
bool MeshLayerMapper::allRastersExist(const std::vector<std::string> &raster_paths) const
{
for (auto raster = raster_paths.begin(); raster != raster_paths.end(); ++raster)
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.h b/MeshLib/MeshGenerators/MeshLayerMapper.h
index 370e035061afcc7b137da31c704cb50be776d7b7..6eeb1bc1202d7497b7ee6570acdf2eecb0c3c05c 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.h
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.h
@@ -36,7 +36,7 @@ public:
~MeshLayerMapper() {};
/**
- * Based on a triangle-or quad mesh this method creates a 3D mesh with with a given number of prism- or hex-layers
+ * Based on a 2D triangle-or quad mesh this method creates a 3D mesh with with a given number of prism- or hex-layers
* \param mesh The triangle/quad mesh that is the basis for the new prism/hex 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
@@ -45,34 +45,34 @@ public:
*/
MeshLib::Mesh* createStaticLayers(MeshLib::Mesh const& mesh, std::vector<float> const& layer_thickness_vector, std::string const& mesh_name = "SubsurfaceMesh") const;
+ /**
+ * 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.
+ * \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
+ */
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;
/**
- * Maps the z-values of nodes in the designated layer of the given mesh according to the given raster.
- * Note: This only results in a valid mesh if the layers don't intersect each other.
- */
- bool layerMapping(MeshLib::Mesh &mesh, const std::string &rasterfile,
- const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue) const;
-
- /**
- * Maps the z-values of nodes in the designated layer of the given mesh according to the given raster.
- * Note: This only results in a valid mesh if the layers don't intersect each other.
+ * 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 GeoLib::Raster &raster,
- const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue) const;
+ bool layerMapping(MeshLib::Mesh &mesh, const std::string &rasterfile, double noDataReplacementValue) const;
/**
- * Blends a mesh with the surface given by dem_raster. Nodes and elements above the surface are either removed or adapted to fit the surface.
- * Note: It is unlikely but possible that the new nodes vector contains (very few) nodes that are not part of any element. This problem is
- * remedied at the end of method upon creating the actual mesh from the new node- and element-vector as the mesh-constructor checks for such
- * nodes and removes them. This note is just to call this issue to attention in case this methods is changed.
+ * 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.
*/
- MeshLib::Mesh* blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster) const;
+ bool layerMapping(MeshLib::Mesh &mesh, const GeoLib::Raster &raster, double noDataReplacementValue) const;
private:
- /// Adds another layer to the subsurface mesh
+ /// 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
diff --git a/MeshLib/Node.h b/MeshLib/Node.h
index ce4e783e122fd861561ea030d0d75b762148108d..a59a637036c6235fcc4fa59eadd0997bdc7fefb5 100644
--- a/MeshLib/Node.h
+++ b/MeshLib/Node.h
@@ -36,9 +36,7 @@ class Element;
class Node : public GeoLib::PointWithID
{
/* friend functions: */
- friend bool MeshLayerMapper::layerMapping(MeshLib::Mesh &mesh, const GeoLib::Raster &raster, const unsigned nLayers,
- const unsigned layer_id, double noDataReplacementValue) const;
- friend MeshLib::Mesh* MeshLayerMapper::blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster) const;
+ friend bool MeshLayerMapper::layerMapping(MeshLib::Mesh &mesh, const GeoLib::Raster &raster, double noDataReplacementValue) const;
friend MeshLib::Mesh* MeshSurfaceExtraction::getMeshSurface(const MeshLib::Mesh &mesh, const MathLib::Vector3 &dir, bool keep3dMeshIds);
/* friend classes: */