Skip to content
Snippets Groups Projects
Commit d53c27b3 authored by Karsten Rink's avatar Karsten Rink
Browse files

simplyfied layerMapping method to only handle single layered meshes

parent aa07dea5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
MeshLib/MeshGenerators/LayeredVolume.cpp
+ 1
− 1
View file @ d53c27b3
......@@ -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;
......
MeshLib/MeshGenerators/MeshLayerMapper.cpp
+ 17
− 167
View file @ d53c27b3
......@@ -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)
......
MeshLib/MeshGenerators/MeshLayerMapper.h
+ 18
− 18
View file @ d53c27b3
......@@ -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
......
MeshLib/Node.h
+ 1
− 3
View file @ d53c27b3
......@@ -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: */
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment