From bde82bbbb3bb340bd2369f1ca8ea0a6e13a993c3 Mon Sep 17 00:00:00 2001
From: Karsten Rink <karsten.rink@ufz.de>
Date: Thu, 18 Sep 2014 14:52:26 +0200
Subject: [PATCH] reworked prims-layer-generation to work from top to bottom
and prefer older (lower) layers in cases of conflicting information
---
MeshLib/MeshGenerators/MeshLayerMapper.cpp | 183 ++++++++++++++++++---
MeshLib/MeshGenerators/MeshLayerMapper.h | 11 ++
2 files changed, 170 insertions(+), 24 deletions(-)
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.cpp b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
index 38ea9700502..cb16ee54144 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.cpp
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
@@ -33,9 +33,18 @@
#include "Elements/Pyramid.h"
#include "Elements/Prism.h"
#include "MeshSurfaceExtraction.h"
+#include "MeshQuality/MeshValidation.h"
#include "MathTools.h"
+const unsigned MeshLayerMapper::_pyramid_base[3][4] =
+{
+ {1, 3, 4, 2}, // Point 4 missing
+ {2, 4, 3, 0}, // Point 5 missing
+ {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;
@@ -52,7 +61,7 @@ MeshLib::Mesh* MeshLayerMapper::createStaticLayers(MeshLib::Mesh const& mesh, st
return nullptr;
}
- const size_t nNodes = mesh.getNNodes();
+ const std::size_t nNodes = mesh.getNNodes();
// count number of 2d elements in the original mesh
const std::size_t nElems (std::count_if(mesh.getElements().begin(), mesh.getElements().end(),
[](MeshLib::Element const* elem) { return (elem->getDimension() == 2);}));
@@ -75,36 +84,152 @@ MeshLib::Mesh* MeshLayerMapper::createStaticLayers(MeshLib::Mesh const& mesh, st
[&z_offset](MeshLib::Node* node){ return new MeshLib::Node((*node)[0], (*node)[1], (*node)[2]-z_offset); });
// starting with 2nd layer create prism or hex elements connecting the last layer with the current one
- if (layer_id > 0)
- {
- node_offset -= nNodes;
- const unsigned mat_id (nLayers - layer_id);
+ if (layer_id == 0)
+ continue;
- for (unsigned i = 0; i < nOrgElems; ++i)
- {
- const MeshLib::Element* sfc_elem( elems[i] );
- if (sfc_elem->getDimension() < 2) // ignore line-elements
- continue;
+ node_offset -= nNodes;
+ const unsigned mat_id (nLayers - layer_id);
+
+ for (unsigned i = 0; i < nOrgElems; ++i)
+ {
+ const MeshLib::Element* sfc_elem( elems[i] );
+ if (sfc_elem->getDimension() < 2) // ignore line-elements
+ continue;
- const unsigned nElemNodes(sfc_elem->getNNodes());
- MeshLib::Node** e_nodes = new MeshLib::Node*[2*nElemNodes];
+ const unsigned nElemNodes(sfc_elem->getNNodes());
+ MeshLib::Node** e_nodes = new MeshLib::Node*[2*nElemNodes];
- for (unsigned j=0; j<nElemNodes; ++j)
- {
- const unsigned node_id = sfc_elem->getNode(j)->getID() + node_offset;
- e_nodes[j] = new_nodes[node_id+nNodes];
- e_nodes[j+nElemNodes] = new_nodes[node_id];
- }
- if (sfc_elem->getGeomType() == MeshElemType::TRIANGLE) // extrude triangles to prism
- new_elems.push_back (new MeshLib::Prism(e_nodes, mat_id));
- else if (sfc_elem->getGeomType() == MeshElemType::QUAD) // extrude quads to hexes
- new_elems.push_back (new MeshLib::Hex(e_nodes, mat_id));
+ for (unsigned j=0; j<nElemNodes; ++j)
+ {
+ const unsigned node_id = sfc_elem->getNode(j)->getID() + node_offset;
+ e_nodes[j] = new_nodes[node_id+nNodes];
+ e_nodes[j+nElemNodes] = new_nodes[node_id];
}
+ if (sfc_elem->getGeomType() == MeshElemType::TRIANGLE) // extrude triangles to prism
+ new_elems.push_back (new MeshLib::Prism(e_nodes, mat_id));
+ else if (sfc_elem->getGeomType() == MeshElemType::QUAD) // extrude quads to hexes
+ new_elems.push_back (new MeshLib::Hex(e_nodes, mat_id));
}
}
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
+{
+ 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.");
+ return nullptr;
+ }
+
+ MeshLib::Mesh* dem_mesh (new MeshLib::Mesh(mesh));
+ if (layerMapping(*dem_mesh, *rasters.back(), 0, 0, 0))
+ {
+ std::size_t const nNodes = mesh.getNNodes();
+ std::vector<MeshLib::Node*> new_nodes;
+ new_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));
+
+ 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;
+
+ }
+}
+
+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
+{
+ std::size_t const nNodes = dem_mesh.getNNodes();
+ 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;
+
+ // add nodes for new layer
+ for (std::size_t i=0; i<nNodes; ++i)
+ {
+ // 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
+ new_nodes.push_back(new MeshLib::Node((*nodes[i])[0], (*nodes[i])[1], elevation, new_nodes[last_layer_node_offset +i]->getID()));
+ }
+
+ if (layer_id == 0)
+ return;
+
+ std::vector<MeshLib::Element*> const& elems = dem_mesh.getElements();
+ std::size_t const nElems (dem_mesh.getNElements());
+
+ for (std::size_t i=0; i<nElems; ++i)
+ {
+ MeshLib::Element* elem (elems[i]);
+ if (elem->getGeomType() != MeshElemType::TRIANGLE)
+ continue;
+ unsigned node_counter(3), missing_idx(0);
+ 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]);
+ if (new_elem_nodes[j]->getID() != new_elem_nodes[node_counter]->getID())
+ node_counter++;
+ else
+ missing_idx = j;
+ }
+
+ switch (node_counter)
+ {
+ case 6:
+ new_elems.push_back(new MeshLib::Prism(new_elem_nodes, layer_id));
+ break;
+ case 5:
+ std::array<MeshLib::Node*, 5> pyramid_nodes;
+ pyramid_nodes[0] = new_elem_nodes[_pyramid_base[missing_idx][0]];
+ pyramid_nodes[1] = new_elem_nodes[_pyramid_base[missing_idx][1]];
+ 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));
+ 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));
+ break;
+ default:
+ continue;
+ }
+ }
+}
+
+
bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, const std::string &rasterfile,
const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue = 0.0) const
{
@@ -154,7 +279,7 @@ bool MeshLayerMapper::layerMapping(MeshLib::Mesh &new_mesh, const GeoLib::Raster
}
double elevation (raster.interpolateValueAtPoint(*nodes[i]));
- if (elevation - no_data < std::numeric_limits<double>::epsilon())
+ if (std::abs(elevation - no_data) < std::numeric_limits<double>::epsilon())
elevation = noDataReplacementValue;
nodes[i]->updateCoordinates((*nodes[i])[0], (*nodes[i])[1], elevation);
}
@@ -307,5 +432,15 @@ MeshLib::Mesh* MeshLayerMapper::blendLayersWithSurface(MeshLib::Mesh &mesh, cons
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)
+ {
+ 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 aab59cb115c..370e035061a 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.h
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.h
@@ -23,6 +23,7 @@ class QImage;
namespace MeshLib {
class Mesh;
class Node;
+ class Element;
}
/**
@@ -44,6 +45,10 @@ public:
*/
MeshLib::Mesh* createStaticLayers(MeshLib::Mesh const& mesh, std::vector<float> const& layer_thickness_vector, std::string const& mesh_name = "SubsurfaceMesh") const;
+ 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.
@@ -67,7 +72,13 @@ public:
MeshLib::Mesh* blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster) const;
private:
+ /// Adds another layer to the 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;
+ static const unsigned _pyramid_base[3][4];
};
#endif //MESHLAYERMAPPER_H
--
GitLab