diff --git a/Gui/DataView/CMakeLists.txt b/Gui/DataView/CMakeLists.txt
index 661fc1fff5528694a73bd868a4028c2caa75b075..21a66b87756a4b3c9e3f9d0845b7a0ebbfc52dc6 100644
--- a/Gui/DataView/CMakeLists.txt
+++ b/Gui/DataView/CMakeLists.txt
@@ -28,7 +28,6 @@ set( SOURCES
ModelTreeItem.cpp
MeshLayerEditDialog.cpp
MshItem.cpp
- MshLayerMapper.cpp
MshModel.cpp
MshQualitySelectionDialog.cpp
MshTabWidget.cpp
@@ -98,7 +97,6 @@ set( HEADERS
GeoTreeItem.h
ModelTreeItem.h
MshItem.h
- MshLayerMapper.h
ProcessItem.h
)
diff --git a/Gui/DataView/CondFromRasterDialog.cpp b/Gui/DataView/CondFromRasterDialog.cpp
index 538e1b3c4a097a1f5eb63be2893040324f18f781..5c1da975ae053353c3ad963571f73f1c1d0eae13 100644
--- a/Gui/DataView/CondFromRasterDialog.cpp
+++ b/Gui/DataView/CondFromRasterDialog.cpp
@@ -52,15 +52,15 @@ void CondFromRasterDialog::on_selectButton_pressed()
QString geotiffExtension("");
#endif
QString fileName = QFileDialog::getOpenFileName(this, "Select raster file",
- settings.value("lastOpenedConditionsFileDirectory").toString(), QString(
- "Raster files (*.asc *.grd);;") .arg(geotiffExtension));
+ settings.value("lastOpenedRasterFileDirectory").toString(),
+ QString("Raster files (*.asc *.grd);;").arg(geotiffExtension));
if (!fileName.isEmpty())
{
this->rasterEdit->setText(fileName);
- QDir dir = QDir(fileName);
- settings.setValue("lastOpenedConditionsFileDirectory", dir.absolutePath());
+ QFileInfo fi(fileName);
+ settings.setValue("lastOpenedRasterFileDirectory", fi.absolutePath());
}
}
diff --git a/Gui/DataView/MeshLayerEditDialog.cpp b/Gui/DataView/MeshLayerEditDialog.cpp
index 9db1c4b653222f6239492fc371493e418f5b60e6..5692a07ce325ae833b948273d37df2874fcb4cab 100644
--- a/Gui/DataView/MeshLayerEditDialog.cpp
+++ b/Gui/DataView/MeshLayerEditDialog.cpp
@@ -15,11 +15,15 @@
#include "MeshLayerEditDialog.h"
+// ThirdParty/logog
+#include "logog/include/logog.hpp"
+
#include "OGSError.h"
#include "StringTools.h"
#include "Mesh.h"
#include <QCheckBox>
+#include <QFileInfo>
#include <QFileDialog>
#include <QPushButton>
#include <QSettings>
@@ -188,7 +192,7 @@ void MeshLayerEditDialog::accept()
const std::string imgPath ( this->_edits[0]->text().toStdString() );
const double noDataReplacementValue = this->_noDataReplacementEdit->text().toDouble();
if (!imgPath.empty())
- result = MshLayerMapper::LayerMapping(new_mesh, imgPath, nLayers, 0, noDataReplacementValue);
+ result = MeshLayerMapper::LayerMapping(*new_mesh, imgPath, nLayers, 0, noDataReplacementValue);
}
else
{
@@ -197,12 +201,11 @@ void MeshLayerEditDialog::accept()
{
// "100" is just a default size to have any value for extruding 2D elements.
// The actual mapping based on a raster file will be performed later.
- float thickness = (_use_rasters) ? 100 : (this->_edits[i]->text().toFloat());
- if (thickness > std::numeric_limits<float>::epsilon())
- layer_thickness.push_back(thickness);
+ const float thickness = (_use_rasters) ? 100 : (this->_edits[i]->text().toFloat());
+ layer_thickness.push_back(thickness);
}
- new_mesh = MshLayerMapper::CreateLayers(_msh, layer_thickness);
+ new_mesh = MeshLayerMapper::CreateLayers(*_msh, layer_thickness);
if (_use_rasters)
{
@@ -212,13 +215,13 @@ void MeshLayerEditDialog::accept()
const double noDataReplacement = (i==0) ? 0.0 : -9999.0;
if (!imgPath.empty())
{
- result = MshLayerMapper::LayerMapping(new_mesh, imgPath, nLayers, i, noDataReplacement);
+ result = MeshLayerMapper::LayerMapping(*new_mesh, imgPath, nLayers, i, noDataReplacement);
if (result==0) break;
}
}
if (this->_edits[0]->text().length()>0)
{
- MeshLib::Mesh* final_mesh = MshLayerMapper::blendLayersWithSurface(new_mesh, nLayers, this->_edits[0]->text().toStdString());
+ MeshLib::Mesh* final_mesh = MeshLayerMapper::blendLayersWithSurface(*new_mesh, nLayers, this->_edits[0]->text().toStdString());
delete new_mesh;
new_mesh = final_mesh;
}
@@ -250,9 +253,9 @@ void MeshLayerEditDialog::getFileName()
QPushButton* button = dynamic_cast<QPushButton*>(this->sender());
QSettings settings;
QString filename = QFileDialog::getOpenFileName(this, "Select raster file to open",
- settings.value("lastOpenedFileDirectory").toString(),
+ settings.value("lastOpenedRasterFileDirectory").toString(),
"ASCII raster files (*.asc);;All files (* *.*)");
_fileButtonMap[button]->setText(filename);
- QDir dir = QDir(filename);
- settings.setValue("lastOpenedFileDirectory", dir.absolutePath());
+ QFileInfo fi(filename);
+ settings.setValue("lastOpenedRasterFileDirectory", fi.absolutePath());
}
diff --git a/Gui/DataView/MeshLayerEditDialog.h b/Gui/DataView/MeshLayerEditDialog.h
index 1b02632fb182eced09000282a8ec44f8d92e767e..89fabac92366e00ed6a5e243e4f2f92d362ffa1c 100644
--- a/Gui/DataView/MeshLayerEditDialog.h
+++ b/Gui/DataView/MeshLayerEditDialog.h
@@ -18,7 +18,7 @@
#include "ui_MeshLayerEdit.h"
#include <QDialog>
-#include "MshLayerMapper.h"
+#include "MeshGenerators/MeshLayerMapper.h"
class QPushButton;
class QCheckBox;
diff --git a/MathLib/MathTools.cpp b/MathLib/MathTools.cpp
index 7de70163cbea2329342798b05848b02c8f505e3c..651fa75369b10a4218753378150339db1694d346 100644
--- a/MathLib/MathTools.cpp
+++ b/MathLib/MathTools.cpp
@@ -86,13 +86,4 @@ double calcTetrahedronVolume(const double* x1, const double* x2, const double* x
+ (x1[2] - x4[2]) * ((x2[0] - x4[0]) * (x3[1] - x4[1]) - (x2[1] - x4[1]) * (x3[0] - x4[0]))) / 6.0;
}
-void MPhi2D(double* vf, double r, double s)
-{
- vf[0] = (1.0 + r) * (1.0 + s);
- vf[1] = (1.0 - r) * (1.0 + s);
- vf[2] = (1.0 - r) * (1.0 - s);
- vf[3] = (1.0 + r) * (1.0 - s);
- for (unsigned i = 0; i < 4; i++)
- vf[i] *= 0.25;
-}
} // namespace
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index f205d25b598c700d82d2edbdf60929117d9c3f22..585d141b14c405b6143eeb70d2ac50a59fc559ab 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -177,8 +177,6 @@ T fastpow (T base, std::size_t exp)
return result;
}
-/// 2D linear interpolation function (TODO adopted from geo_mathlib)
-void MPhi2D(double* vf, double r, double s);
} // namespace
#endif /* MATHTOOLS_H_ */
diff --git a/Gui/DataView/MshLayerMapper.cpp b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
similarity index 50%
rename from Gui/DataView/MshLayerMapper.cpp
rename to MeshLib/MeshGenerators/MeshLayerMapper.cpp
index 24b3731ffdabcfe21321966ed758c29973f0a686..ffaa8c575c49c79e4ee86c3cdcfb420921170ebf 100644
--- a/Gui/DataView/MshLayerMapper.cpp
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
@@ -1,8 +1,8 @@
/**
- * \file
+ * \file MeshLayerMapper.cpp
* \author Karsten Rink
* \date 2010-11-01
- * \brief Implementation of the MshLayerMapper class.
+ * \brief Implementation of the MeshLayerMapper class.
*
* \copyright
* Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
@@ -18,7 +18,7 @@
// ThirdParty/logog
#include "logog/include/logog.hpp"
-#include "MshLayerMapper.h"
+#include "MeshLayerMapper.h"
// GeoLib
#include "Raster.h"
@@ -33,199 +33,182 @@
#include "MathTools.h"
-MeshLib::Mesh* MshLayerMapper::CreateLayers(const MeshLib::Mesh* mesh, const std::vector<float> &thickness)
+MeshLib::Mesh* MeshLayerMapper::CreateLayers(const MeshLib::Mesh &mesh, const std::vector<float> &layer_thickness_vector)
{
- std::size_t nLayers(thickness.size());
- bool throw_error(false);
- for (unsigned i=0; i<nLayers; ++i)
- if (thickness[i]<=0)
- throw_error = true;
- if (nLayers < 1 || mesh->getDimension() != 2)
- throw_error = true;
-
- if (throw_error)
+ std::vector<float> thickness;
+ for (std::size_t i=0; i<layer_thickness_vector.size(); ++i)
+ if (layer_thickness_vector[i] > std::numeric_limits<float>::epsilon())
+ thickness.push_back(layer_thickness_vector[i]);
+ else
+ WARN ("Ignoring layer %d with thickness %f.", i, layer_thickness_vector[i]);
+
+ const std::size_t nLayers(thickness.size());
+ if (nLayers < 1 || mesh.getDimension() != 2)
{
ERR("MshLayerMapper::CreateLayers(): A 2D mesh with nLayers > 0 is required as input.");
return nullptr;
}
- const size_t nNodes = mesh->getNNodes();
+ const 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(),
+ const std::size_t nElems (std::count_if(mesh.getElements().begin(), mesh.getElements().end(),
[](MeshLib::Element const* elem) { return (elem->getDimension() == 2);}));
- const std::vector<MeshLib::Node*> nodes = mesh->getNodes();
- const std::vector<MeshLib::Element*> elems = mesh->getElements();
- std::vector<MeshLib::Node*> new_nodes(nNodes + (nLayers * nNodes));
- std::vector<MeshLib::Element*> new_elems(nElems * nLayers);
- double z_offset(0.0);
+ const std::size_t nOrgElems (mesh.getNElements());
+ const std::vector<MeshLib::Node*> &nodes = mesh.getNodes();
+ const std::vector<MeshLib::Element*> &elems = mesh.getElements();
+ std::vector<MeshLib::Node*> new_nodes;
+ new_nodes.reserve(nNodes + (nLayers * nNodes));
+ std::vector<MeshLib::Element*> new_elems;
+ new_elems.reserve(nElems * nLayers);
+ double z_offset (0.0);
for (unsigned layer_id = 0; layer_id <= nLayers; ++layer_id)
{
// add nodes for new layer
unsigned node_offset (nNodes * layer_id);
- unsigned elem_offset (nElems * (layer_id-1));
- if (layer_id>0) z_offset += thickness[layer_id-1];
+ if (layer_id > 0) z_offset += thickness[layer_id-1];
for (unsigned i = 0; i < nNodes; ++i)
{
const double* coords = nodes[i]->getCoords();
- new_nodes[node_offset+i] = new MeshLib::Node(coords[0], coords[1], coords[2]-z_offset, node_offset+i);
+ new_nodes.push_back (new MeshLib::Node(coords[0], coords[1], coords[2]-z_offset, node_offset+i));
}
// starting with 2nd layer create prism or hex elements connecting the last layer with the current one
+ const unsigned elem_offset (nElems * (layer_id-1));
if (layer_id > 0)
{
- if (thickness[layer_id-1] > 0)
+ node_offset -= nNodes;
+ const unsigned mat_id (nLayers - layer_id);
+
+ for (unsigned i = 0; i < nOrgElems; ++i)
{
- node_offset -= nNodes;
- const unsigned mat_id (nLayers - layer_id);
+ const MeshLib::Element* sfc_elem( elems[i] );
+ if (sfc_elem->getDimension() != 2)
+ {
+ WARN("MshLayerMapper::CreateLayers() - Method can only handle 2D mesh elements.");
+ WARN("Skipping Element %d of type \"%s\".", i, MeshElemType2String(sfc_elem->getGeomType()).c_str());
+ continue;
+ }
+
+ const unsigned nElemNodes(sfc_elem->getNNodes());
+ MeshLib::Node** e_nodes = new MeshLib::Node*[2*nElemNodes];
- // counts the 2d elements (within the layer),
- // used as a part of index computation in new_elems
- std::size_t cnt(0);
- for (unsigned i = 0; i < mesh->getNElements(); ++i)
+ for (unsigned j=0; j<nElemNodes; ++j)
{
- const MeshLib::Element* sfc_elem( elems[i] );
- if (sfc_elem->getDimension() == 2)
- {
- 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[elem_offset+cnt] = new MeshLib::Prism(e_nodes, mat_id);
- else if (sfc_elem->getGeomType() == MeshElemType::QUAD) // extrude quads to hexes
- new_elems[elem_offset+cnt] = new MeshLib::Hex(e_nodes, mat_id);
- cnt++;
- }
- else
- {
- WARN("MshLayerMapper::CreateLayers() - Method can only handle 2D mesh elements.");
- WARN("Skipping Element %d of type \"%s\".", i, MeshElemType2String(sfc_elem->getGeomType()).c_str());
- }
+ 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];
}
- }
- else
- {
- ERR("Error in MshLayerMapper::CreateLayers() - Layer thickness for layer %d is %f (needs to be >0).", (layer_id-1), thickness[layer_id-1]);
- return nullptr;
+ 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("SubsurfaceMesh", new_nodes, new_elems);
}
-int MshLayerMapper::LayerMapping(MeshLib::Mesh* new_mesh, const std::string &rasterfile,
+int MeshLayerMapper::LayerMapping(MeshLib::Mesh &new_mesh, const std::string &rasterfile,
const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue = 0.0)
{
- if (new_mesh == nullptr)
+ if (nLayers < layer_id)
{
- ERR("MshLayerMapper::LayerMapping() - Passed Mesh is NULL.");
+ ERR("MshLayerMapper::LayerMapping() - Mesh has only %d Layers, cannot assign layer %d.", nLayers, layer_id);
return 0;
}
- if (nLayers >= layer_id)
+ const GeoLib::Raster *raster(GeoLib::Raster::getRasterFromASCFile(rasterfile));
+ if (! raster) {
+ ERR("MshLayerMapper::LayerMapping - could not read raster file %s", rasterfile.c_str());
+ return 0;
+ }
+ const double x0(raster->getOrigin()[0]);
+ const double y0(raster->getOrigin()[1]);
+ const double delta(raster->getRasterPixelDistance());
+ const double no_data(raster->getNoDataValue());
+ const std::size_t width(raster->getNCols());
+ const std::size_t height(raster->getNRows());
+ double const*const elevation(raster->begin());
+
+ const std::pair<double, double> xDim(x0, x0 + width * delta); // extension in x-dimension
+ const std::pair<double, double> yDim(y0, y0 + height * 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)
{
- GeoLib::Raster *raster(GeoLib::Raster::getRasterFromASCFile(rasterfile));
- if (! raster) {
- ERR("MshLayerMapper::LayerMapping - could not read raster file %s", rasterfile.c_str());
- return 0;
- }
- const double x0(raster->getOrigin()[0]);
- const double y0(raster->getOrigin()[1]);
- const double delta(raster->getRasterPixelDistance());
- const double no_data(raster->getNoDataValue());
- const std::size_t width(raster->getNCols());
- const std::size_t height(raster->getNRows());
- double const*const elevation(raster->begin());
-
- const std::pair<double, double> xDim(x0, x0 + width * delta); // extension in x-dimension
- const std::pair<double, double> yDim(y0, y0 + height * delta); // extension in y-dimension
+ const double* coords (nodes[i]->getCoords());
- 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)
+ if (!isNodeOnRaster(*nodes[i], xDim, yDim))
{
- const double* coords = nodes[i]->getCoords();
+ // use either default value or elevation from layer above
+ const double new_elevation = (layer_id == 0) ? noDataReplacementValue : (*nodes[i-nNodesPerLayer])[2];
+ nodes[i]->updateCoordinates(coords[0], coords[1], noDataReplacementValue);
+ continue;
+ }
- if (!isNodeOnRaster(*nodes[i], xDim, yDim))
+ // position in raster
+ const double xPos ((coords[0] - xDim.first) / delta);
+ const double yPos ((coords[1] - yDim.first) / delta);
+ // raster cell index
+ const size_t xIdx (static_cast<size_t>(floor(xPos)));
+ const size_t yIdx (static_cast<size_t>(floor(yPos)));
+
+ // weights for bilinear interpolation
+ const double xShift = fabs(xPos-(xIdx+half_delta))/delta;
+ const double yShift = fabs(yPos-(yIdx+half_delta))/delta;
+ std::array<double,4> weight = { (1-xShift)*(1-xShift), xShift*(1-yShift), xShift*yShift, (1-xShift)*yShift };
+
+ // neightbors to include in interpolation
+ const int xShiftIdx = (xPos-xIdx-half_delta>=0) ? 1 : -1;
+ const int yShiftIdx = (yPos-yIdx-half_delta>=0) ? 1 : -1;
+ const std::array<int,4> x_nb = { 0, xShiftIdx, xShiftIdx, 0 };
+ const std::array<int,4> y_nb = { 0, 0, yShiftIdx, yShiftIdx };
+
+ // get pixel values
+ std::array<double,4> pix_val;
+ unsigned no_data_count (0);
+ for (unsigned j=0; j<4; ++j)
+ {
+ pix_val[j] = elevation[(yIdx + y_nb[j])*width + (xIdx + x_nb[j])];
+ if (fabs(pix_val[j] - no_data) < std::numeric_limits<double>::epsilon())
{
- if (layer_id == 0) // use default value
- nodes[i]->updateCoordinates(coords[0], coords[1], noDataReplacementValue);
- else // use z-value from layer above
- nodes[i]->updateCoordinates(coords[0], coords[1], (*nodes[i-nNodesPerLayer])[2]);
- continue;
+ weight[j] = 0;
+ no_data_count++;
}
+ }
- // position in raster
- const double xPos ((coords[0] - xDim.first) / delta);
- const double yPos ((coords[1] - yDim.first) / delta);
- // raster cell index
- const size_t xIdx (static_cast<size_t>(floor(xPos)));
- const size_t yIdx (static_cast<size_t>(floor(yPos)));
-
- // deviation of mesh node from centre of raster cell ( in [-1:1) because it is normalised by delta/2 )
- const double xShift = (xPos-xIdx-half_delta)/half_delta;
- const double yShift = (yPos-yIdx-half_delta)/half_delta;
-
- const int xShiftIdx = static_cast<int>((xShift>=0) ? ceil(xShift) : floor(xShift));
- const int yShiftIdx = static_cast<int>((yShift>=0) ? ceil(yShift) : floor(yShift));
-
- // determining the neighbouring pixels that add weight to the interpolation
- const int x_nb[4] = {0, xShiftIdx, xShiftIdx, 0};
- const int y_nb[4] = {0, 0, yShiftIdx, yShiftIdx};
-
- double locZ[4];
- locZ[0] = elevation[yIdx*width + xIdx];
- if (fabs(locZ[0] - no_data) > std::numeric_limits<double>::epsilon())
- {
- for (unsigned j=1; j<4; ++j)
- {
- locZ[j] = elevation[(yIdx+y_nb[j])*width + (xIdx+x_nb[j])];
- if (fabs(locZ[j] - no_data) < std::numeric_limits<double>::epsilon())
- locZ[j]=locZ[0];
- }
-
- double ome[4];
- double xi = 1-fabs(xShift);
- double eta = 1-fabs(xShift);
- MathLib::MPhi2D(ome, xi, eta);
-
- double z(0.0);
- for(unsigned j = 0; j < 4; ++j)
- z += ome[j] * locZ[j];
-
- nodes[i]->updateCoordinates(coords[0], coords[1], z);
- }
- else
+ // adjust weights if necessary
+ if (no_data_count > 0)
+ {
+ if (no_data_count == 4) // if there is absolutely no data just use the default value
{
- if (layer_id == 0) // use default value
- nodes[i]->updateCoordinates(coords[0], coords[1], noDataReplacementValue);
- else // use z-value from layer above
- nodes[i]->updateCoordinates(coords[0], coords[1], (*nodes[i-nNodesPerLayer])[2]);
+ const double new_elevation = (layer_id == 0) ? noDataReplacementValue : (*nodes[i-nNodesPerLayer])[2];
+ nodes[i]->updateCoordinates(coords[0], coords[1], noDataReplacementValue);
+ continue;
}
+ const double norm = 4/(4-no_data_count);
+ std::for_each(weight.begin(), weight.end(), [&norm](double &val){val*=norm;});
}
- delete raster;
- return 1;
+ // new value
+ double z = MathLib::scalarProduct<double,4>(weight.data(), pix_val.data());
+ nodes[i]->updateCoordinates(coords[0], coords[1], z);
}
- else
- ERR("MshLayerMapper::LayerMapping() - Mesh has only %d Layers, cannot assign layer %d.", nLayers, layer_id);
- return 0;
+
+ delete raster;
+ return 1;
}
-bool MshLayerMapper::isNodeOnRaster(const MeshLib::Node &node,
+bool MeshLayerMapper::isNodeOnRaster(const MeshLib::Node &node,
const std::pair<double, double> &xDim,
const std::pair<double, double> &yDim)
{
@@ -235,16 +218,16 @@ bool MshLayerMapper::isNodeOnRaster(const MeshLib::Node &node,
return true;
}
-MeshLib::Mesh* MshLayerMapper::blendLayersWithSurface(MeshLib::Mesh* mesh, const unsigned nLayers, const std::string &dem_raster)
+MeshLib::Mesh* MeshLayerMapper::blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster)
{
// construct surface mesh from DEM
const MathLib::Vector3 dir(0,0,1);
- MeshLib::Mesh* dem = MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh, dir);
- MshLayerMapper::LayerMapping(dem, dem_raster, 0, 0);
- const std::vector<MeshLib::Node*> dem_nodes (dem->getNodes());
+ MeshLib::Mesh* dem = MeshLib::MeshSurfaceExtraction::getMeshSurface(mesh, dir);
+ MeshLayerMapper::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 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);
@@ -268,7 +251,7 @@ MeshLib::Mesh* MshLayerMapper::blendLayersWithSurface(MeshLib::Mesh* mesh, const
// 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 them)
+ // 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;
@@ -292,6 +275,8 @@ MeshLib::Mesh* MshLayerMapper::blendLayersWithSurface(MeshLib::Mesh* mesh, const
}
}
}
+ 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);
@@ -304,8 +289,8 @@ MeshLib::Mesh* MshLayerMapper::blendLayersWithSurface(MeshLib::Mesh* mesh, const
}
// 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();
+ 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)
{
@@ -332,7 +317,7 @@ MeshLib::Mesh* MshLayerMapper::blendLayersWithSurface(MeshLib::Mesh* mesh, const
if (!nodes_below_surface[elem->getNode(i)->getID()])
top_idx = i-3;
- // construct pyrmid element based on missing node
+ // 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()]];
diff --git a/Gui/DataView/MshLayerMapper.h b/MeshLib/MeshGenerators/MeshLayerMapper.h
similarity index 81%
rename from Gui/DataView/MshLayerMapper.h
rename to MeshLib/MeshGenerators/MeshLayerMapper.h
index 6f483272ed3e1d47aed1a48bcb098e41ae1ea056..bcbd9a7da94e38a561fa583aec776ed439c7115d 100644
--- a/Gui/DataView/MshLayerMapper.h
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.h
@@ -1,8 +1,8 @@
/**
- * \file
+ * \file MeshLayerMapper.h
* \author Karsten Rink
* \date 2010-11-01
- * \brief Definition of the MshLayerMapper class.
+ * \brief Definition of the MeshLayerMapper class.
*
* \copyright
* Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
@@ -12,8 +12,8 @@
*
*/
-#ifndef MSHLAYERMAPPER_H
-#define MSHLAYERMAPPER_H
+#ifndef MESHLAYERMAPPER_H
+#define MESHLAYERMAPPER_H
#include <string>
@@ -27,11 +27,11 @@ namespace MeshLib {
/**
* \brief Manipulating and adding layers to an existing mesh
*/
-class MshLayerMapper
+class MeshLayerMapper
{
public:
- MshLayerMapper() {}
- ~MshLayerMapper() {}
+ MeshLayerMapper() {}
+ ~MeshLayerMapper() {}
/**
* Based on a triangle-or quad mesh this method creates a 3D mesh with with a given number of prism- or hex-layers
@@ -40,13 +40,13 @@ public:
* \param thickness The thickness of each of these newly added layers
* \return A mesh with the requested number of layers of prism/hex elements
*/
- static MeshLib::Mesh* CreateLayers(const MeshLib::Mesh* mesh, const std::vector<float> &thickness);
+ static MeshLib::Mesh* CreateLayers(const MeshLib::Mesh &mesh, const std::vector<float> &layer_thickness_vector);
/**
* 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.
*/
- static int LayerMapping(MeshLib::Mesh* msh, const std::string &rasterfile,
+ static int LayerMapping(MeshLib::Mesh &mesh, const std::string &rasterfile,
const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue);
/**
@@ -55,7 +55,7 @@ public:
* 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.
*/
- static MeshLib::Mesh* blendLayersWithSurface(MeshLib::Mesh* mesh, const unsigned nLayers, const std::string &dem_raster);
+ static MeshLib::Mesh* blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster);
private:
/// Checks if the given mesh is within the dimensions given by xDim and yDim.
@@ -64,4 +64,4 @@ private:
const std::pair<double, double> &yDim);
};
-#endif //MSHLAYERMAPPER_H
+#endif //MESHLAYERMAPPER_H
diff --git a/MeshLib/Node.h b/MeshLib/Node.h
index e6000a6cd5d1e06584a308f48b60168c2a6e40c7..8bdb813c93418f5dc983fde85c54f4038c2f28a4 100644
--- a/MeshLib/Node.h
+++ b/MeshLib/Node.h
@@ -24,9 +24,7 @@
#include "Mesh.h"
#include "MeshEditing/removeMeshNodes.h"
#include "MeshSurfaceExtraction.h"
-#ifdef OGS_BUILD_GUI
- #include "../Gui/DataView/MshLayerMapper.h"
-#endif
+#include "MeshGenerators/MeshLayerMapper.h"
namespace MeshLib {
@@ -38,11 +36,10 @@ class Element;
class Node : public GeoLib::PointWithID
{
/* friend functions: */
-#ifdef OGS_BUILD_GUI
- friend int MshLayerMapper::LayerMapping(MeshLib::Mesh* msh, const std::string &rasterfile, const unsigned nLayers,
+ friend int MeshLayerMapper::LayerMapping(MeshLib::Mesh &mesh, const std::string &rasterfile, const unsigned nLayers,
const unsigned layer_id, double noDataReplacementValue);
- friend MeshLib::Mesh* MshLayerMapper::blendLayersWithSurface(MeshLib::Mesh* mesh, const unsigned nLayers, const std::string &dem_raster);
-#endif
+ friend MeshLib::Mesh* MeshLayerMapper::blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster);
+
/* friend classes: */
friend class Mesh;
friend class MeshRevision;