diff --git a/FileIO/ImportFileTypes.h b/FileIO/ImportFileTypes.h
index 6b7eb2ee6eba1f6af6c4915cb88831749fd1c4cc..8f51d734d8de2f8aa75cbf2c710783cc13aba8db 100644
--- a/FileIO/ImportFileTypes.h
+++ b/FileIO/ImportFileTypes.h
@@ -84,7 +84,7 @@ public:
else if (t==ImportFileType::SHAPE)
return "ESRI Shape files (*.shp)";
else if (t==ImportFileType::TETGEN)
- return "TetGen node files (*.node *.poly)";
+ return "TetGen node files (*.node *.poly *.smesh)";
else if (t==ImportFileType::VTK)
return "VTK files (*.vtk *.vti *.vtr *.vts *.vtp *.vtu)";
else return "All files (*.*)";
diff --git a/FileIO/TetGenInterface.cpp b/FileIO/TetGenInterface.cpp
index ad009c5e93b724430ddd8995a4ba680cfb4eb6e2..7eb739682956001949c5c481c715d77e5b9b4bf6 100644
--- a/FileIO/TetGenInterface.cpp
+++ b/FileIO/TetGenInterface.cpp
@@ -26,16 +26,20 @@
// FileIO
#include "TetGenInterface.h"
+// MathLib
+#include "TemplatePoint.h"
+
// MeshLib
#include "Mesh.h"
#include "Node.h"
#include "Elements/Element.h"
#include "Elements/Tet.h"
+#include "MeshInformation.h"
namespace FileIO
{
TetGenInterface::TetGenInterface() :
- _zero_based_idx (false)
+ _zero_based_idx (false), _boundary_markers (false)
{
}
@@ -43,14 +47,20 @@ TetGenInterface::~TetGenInterface()
{
}
-bool TetGenInterface::readTetGenPoly (std::string const& poly_fname,
- GeoLib::GEOObjects &geo_objects)
+bool TetGenInterface::readTetGenGeometry (std::string const& geo_fname,
+ GeoLib::GEOObjects &geo_objects)
{
- std::ifstream poly_stream (poly_fname.c_str());
+ std::ifstream poly_stream (geo_fname.c_str());
if (!poly_stream)
{
- ERR ("TetGenInterface::readTetGenPoly failed to open %s", poly_fname.c_str());
+ ERR ("TetGenInterface::readTetGenPoly() failed to open %s", geo_fname.c_str());
+ return false;
+ }
+ std::string ext (BaseLib::getFileExtension(geo_fname));
+ if (ext.compare("poly") && ext.compare("smesh"))
+ {
+ ERR ("TetGenInterface::readTetGenPoly() - unknown file type (only *.poly or *.smesh supported).");
return false;
}
@@ -70,8 +80,17 @@ bool TetGenInterface::readTetGenPoly (std::string const& poly_fname,
points->push_back(new GeoLib::Point(nodes[k]->getCoords()));
delete nodes[k];
}
+ std::string geo_name (BaseLib::extractBaseNameWithoutExtension(geo_fname));
+ geo_objects.addPointVec(points, geo_name);
+ const std::vector<std::size_t> &id_map (geo_objects.getPointVecObj(geo_name)->getIDMap());
+
std::vector<GeoLib::Surface*> *surfaces = new std::vector<GeoLib::Surface*>;
- if (!parseFacets(poly_stream, *surfaces, *points))
+ bool failed (true);
+ if (ext.compare("poly") == 0)
+ failed = !parsePolyFacets(poly_stream, *surfaces, *points, id_map);
+ else
+ failed = !parseSmeshFacets(poly_stream, *surfaces, *points, id_map);
+ if (failed)
{
// remove surfaces read until now but keep the points
for (std::size_t k=0; k<surfaces->size(); k++)
@@ -80,14 +99,12 @@ bool TetGenInterface::readTetGenPoly (std::string const& poly_fname,
surfaces = nullptr;
}
- std::string geo_name (BaseLib::extractBaseNameWithoutExtension(poly_fname));
- geo_objects.addPointVec(points, geo_name);
if (surfaces)
geo_objects.addSurfaceVec(surfaces, geo_name);
return true;
}
-std::size_t TetGenInterface::getNFacets(std::ifstream &input) const
+std::size_t TetGenInterface::getNFacets(std::ifstream &input)
{
std::string line;
while (!input.fail())
@@ -104,15 +121,19 @@ std::size_t TetGenInterface::getNFacets(std::ifstream &input) const
continue;
const std::list<std::string> fields = BaseLib::splitString(line, ' ');
- return BaseLib::str2number<size_t> (*fields.begin());
- // here this line also includes a flag for boundary markers which we ignore for now
+ std::list<std::string>::const_iterator it = fields.begin();
+ const std::size_t nFacets (BaseLib::str2number<size_t> (*it));
+ if (fields.size() > 1)
+ _boundary_markers = (BaseLib::str2number<size_t> (*(++it)) == 0) ? false : true;
+ return nFacets;
}
return false;
}
-bool TetGenInterface::parseFacets(std::ifstream &input,
- std::vector<GeoLib::Surface*> &surfaces,
- std::vector<GeoLib::Point*> &points) const
+bool TetGenInterface::parsePolyFacets(std::ifstream &input,
+ std::vector<GeoLib::Surface*> &surfaces,
+ const std::vector<GeoLib::Point*> &points,
+ const std::vector<std::size_t> &pnt_id_map)
{
const std::size_t nFacets (this->getNFacets(input));
std::size_t nMultPolys (0);
@@ -163,7 +184,7 @@ bool TetGenInterface::parseFacets(std::ifstream &input,
{
GeoLib::Polyline polyline(points);
for (std::size_t j(0); j<nPoints; ++j)
- polyline.addPoint(BaseLib::str2number<std::size_t>(*(++it))-offset);
+ polyline.addPoint(pnt_id_map[BaseLib::str2number<std::size_t>(*(++it))-offset]);
polyline.closePolyline();
surfaces.push_back(GeoLib::Surface::createSurface(polyline));
@@ -188,6 +209,79 @@ bool TetGenInterface::parseFacets(std::ifstream &input,
return false;
}
+bool TetGenInterface::parseSmeshFacets(std::ifstream &input,
+ std::vector<GeoLib::Surface*> &surfaces,
+ const std::vector<GeoLib::Point*> &points,
+ const std::vector<std::size_t> &pnt_id_map)
+{
+ const std::size_t nFacets (this->getNFacets(input));
+ std::string line;
+ surfaces.reserve(nFacets);
+ std::list<std::string>::const_iterator it;
+
+ const unsigned offset = (_zero_based_idx) ? 0 : 1;
+ std::vector<std::size_t> idx_map;
+
+ for (std::size_t k(0); k<nFacets && !input.fail(); k++)
+ {
+ getline (input, line);
+ if (input.fail())
+ {
+ ERR("TetGenInterface::parseFacets(): Error reading facet %d.", k);
+ return false;
+ }
+
+ BaseLib::simplify(line);
+ if (line.empty() || line.compare(0,1,"#") == 0)
+ {
+ k--;
+ continue;
+ }
+
+ // read facets
+ const std::list<std::string> point_fields = BaseLib::splitString(line, ' ');
+ it = point_fields.begin();
+ const std::size_t nPoints = BaseLib::str2number<std::size_t>(*it);
+ if (nPoints != 3)
+ {
+ ERR ("Smesh-files are currently only supported for triangle meshes.");
+ return false;
+ }
+ std::vector<std::size_t> point_ids;
+ const std::size_t point_field_size = (_boundary_markers) ? nPoints+1 : nPoints;
+ if (point_fields.size() > point_field_size)
+ {
+ for (std::size_t j(0); j<nPoints; ++j)
+ point_ids.push_back(pnt_id_map[BaseLib::str2number<std::size_t>(*(++it))-offset]);
+
+ const std::size_t sfc_marker = (_boundary_markers) ? BaseLib::str2number<std::size_t>(*(++it)) : 0;
+ const std::size_t idx = std::find(idx_map.begin(), idx_map.end(), sfc_marker) - idx_map.begin();
+ if (idx >= surfaces.size())
+ {
+ idx_map.push_back(sfc_marker);
+ surfaces.push_back(new GeoLib::Surface(points));
+ }
+ surfaces[idx]->addTriangle(point_ids[0], point_ids[1], point_ids[2]);
+ }
+ else
+ {
+ ERR("TetGenInterface::parseFacets(): Error reading points for facet %d.", k);
+ return false;
+ }
+ }
+ // here the poly-file potentially defines a number of points to mark holes within the volumes defined by the facets, these are ignored for now
+ // here the poly-file potentially defines a number of region attributes, these are ignored for now
+
+ std::size_t nTotalTriangles (0);
+ for (std::size_t i=0; i<surfaces.size(); ++i)
+ nTotalTriangles += surfaces[i]->getNTriangles();
+ if (nTotalTriangles == nFacets)
+ return true;
+
+ ERR ("TetGenInterface::parseFacets(): Number of expected total triangles (%d) does not match number of found triangles (%d).", surfaces.size(), nTotalTriangles);
+ return false;
+}
+
MeshLib::Mesh* TetGenInterface::readTetGenMesh (std::string const& nodes_fname,
std::string const& ele_fname)
{
@@ -234,18 +328,16 @@ bool TetGenInterface::readNodesFromStream (std::ifstream &ins,
{
std::string line;
getline (ins, line);
- size_t pos_beg (line.find_first_not_of(" "));
size_t n_nodes, dim, n_attributes;
bool boundary_markers;
while (!ins.fail())
{
- line = line.substr(pos_beg);
- if (line.compare(0,1,"#") == 0)
+ BaseLib::simplify(line);
+ if (line.empty() || line.compare(0,1,"#") == 0)
{
// this line is a comment - skip
getline (ins, line);
- pos_beg = line.find_first_not_of(" ");
continue;
}
// read header line
@@ -265,35 +357,25 @@ bool TetGenInterface::parseNodesFileHeader(std::string &line,
std::size_t &n_attributes,
bool &boundary_markers) const
{
- std::size_t pos_beg, pos_end;
-
- // number of nodes
- pos_beg = line.find_first_not_of (" ");
- pos_end = line.find_first_of(" ", pos_beg);
- if (pos_beg != std::string::npos && pos_end != std::string::npos)
- n_nodes = BaseLib::str2number<size_t> (line.substr(pos_beg, pos_end - pos_beg));
- else
+ std::list<std::string> pnt_header = BaseLib::splitString(line, ' ');
+ if (pnt_header.empty())
{
ERR("TetGenInterface::parseNodesFileHeader(): could not read number of nodes specified in header.");
return false;
}
- // dimension
- pos_beg = line.find_first_not_of (" ", pos_end);
- pos_end = line.find_first_of(" ", pos_beg);
- dim = BaseLib::str2number<size_t> (line.substr(pos_beg, pos_end - pos_beg));
- // number of attributes
- pos_beg = line.find_first_not_of (" ", pos_end);
- pos_end = line.find_first_of(" ", pos_beg);
- n_attributes = BaseLib::str2number<size_t> (line.substr(pos_beg, pos_end - pos_beg));
- // boundary marker at nodes?
- pos_beg = line.find_first_not_of (" ", pos_end);
- pos_end = line.find_first_of(" ", pos_beg);
- if (pos_end == std::string::npos)
- pos_end = line.size();
- if ((line.substr(pos_beg, pos_end - pos_beg)).compare("1") == 0)
- boundary_markers = true;
- else
+ auto it = pnt_header.begin();
+ n_nodes = BaseLib::str2number<size_t> (*it);
+ dim = (pnt_header.size()==1) ? 3 : BaseLib::str2number<size_t> (*(++it));
+
+ if (pnt_header.size()<4)
+ {
+ n_attributes = 0;
boundary_markers = false;
+ return true;
+ }
+
+ n_attributes = BaseLib::str2number<size_t> (*(++it));
+ boundary_markers = ((++it)->compare("1") == 0) ? true : false;
return true;
}
@@ -365,18 +447,16 @@ bool TetGenInterface::readElementsFromStream(std::ifstream &ins,
{
std::string line;
getline (ins, line);
- std::size_t pos_beg (line.find_first_not_of(" "));
std::size_t n_tets, n_nodes_per_tet;
bool region_attributes;
while (!ins.fail())
{
- line = line.substr(pos_beg);
- if (line.compare(0,1,"#") == 0)
+ BaseLib::simplify(line);
+ if (line.empty() || line.compare(0,1,"#") == 0)
{
// this line is a comment - skip
getline (ins, line);
- pos_beg = line.find_first_not_of(" ");
continue;
}
@@ -500,9 +580,10 @@ bool TetGenInterface::parseElements(std::ifstream& ins,
return true;
}
-bool TetGenInterface::writeTetGenPoly(const std::string &file_name,
- const GeoLib::GEOObjects &geo_objects,
- const std::string &geo_name) const
+bool TetGenInterface::writeTetGenSmesh(const std::string &file_name,
+ const GeoLib::GEOObjects &geo_objects,
+ const std::string &geo_name,
+ const std::vector<GeoLib::PointWithID> &attribute_points) const
{
std::vector<GeoLib::Point*> const*const points = geo_objects.getPointVec(geo_name);
std::vector<GeoLib::Surface*> const*const surfaces = geo_objects.getSurfaceVec(geo_name);
@@ -516,33 +597,100 @@ bool TetGenInterface::writeTetGenPoly(const std::string &file_name,
WARN ("No surfaces found for geometry %s. Writing points only.", geo_name.c_str());
std::ofstream out( file_name.c_str(), std::ios::out );
+ out.precision(std::numeric_limits<double>::digits10);
// the points header
const std::size_t nPoints (points->size());
- out << nPoints << " 3 0 0\n";
+ out << nPoints << " 3\n";
// the point list
for (std::size_t i=0; i<nPoints; ++i)
out << i << " " << (*(*points)[i])[0] << " " << (*(*points)[i])[1] << " " << (*(*points)[i])[2] << "\n";
// the surfaces header
const std::size_t nSurfaces = (surfaces) ? surfaces->size() : 0;
- out << nSurfaces << " 0\n";
- // the facets list
+ std::size_t nTotalTriangles (0);
+ for (std::size_t i=0; i<nSurfaces; ++i)
+ nTotalTriangles += (*surfaces)[i]->getNTriangles();
+ out << nTotalTriangles << " 1\n";
+
for (std::size_t i=0; i<nSurfaces; ++i)
{
- // the number of polys per facet
const std::size_t nTriangles ((*surfaces)[i]->getNTriangles());
- out << nTriangles << "\n";
+ const std::size_t marker (i+1); // must NOT be 0!
// the poly list
for (std::size_t j=0; j<nTriangles; ++j)
{
const GeoLib::Triangle &tri = *(*(*surfaces)[i])[j];
- out << "3 " << tri[0] << " " << tri[1] << " " << tri[2] << "\n";
+ out << "3 " << tri[0] << " " << tri[1] << " " << tri[2] << " " << marker << "\n";
}
}
out << "0\n"; // the polygon holes list
- out << "0\n"; // the region attribues list
+ // the region attributes list
+ if (attribute_points.empty())
+ out << "0\n";
+ else
+ {
+ const std::size_t nAttributePoints (attribute_points.size());
+ out << nAttributePoints << "\n";
+ for (std::size_t i=0; i<nAttributePoints; ++i)
+ out << i+1 << " " << attribute_points[i][0] << " " << attribute_points[i][1] << " " << attribute_points[i][2] << " " << 10*attribute_points[i].getID() << "\n";
+ }
INFO ("TetGenInterface::writeTetGenPoly() - %d points and %d surfaces successfully written.", nPoints, nSurfaces);
out.close();
return true;
}
+bool TetGenInterface::writeTetGenSmesh(const std::string &file_name,
+ const MeshLib::Mesh &mesh,
+ const std::vector<MeshLib::Node> &attribute_points) const
+{
+ if (mesh.getDimension() != 2)
+ return false;
+
+ const std::vector<MeshLib::Node*> &nodes = mesh.getNodes();
+ const std::vector<MeshLib::Element*> &elements = mesh.getElements();
+
+ std::ofstream out( file_name.c_str(), std::ios::out );
+ out.precision(std::numeric_limits<double>::digits10);
+ // the points header
+ const std::size_t nPoints (nodes.size());
+ out << nPoints << " 3\n";
+ // the point list
+ for (std::size_t i=0; i<nPoints; ++i)
+ out << i << " " << (*nodes[i])[0] << " " << (*nodes[i])[1] << " " << (*nodes[i])[2] << "\n";
+
+ // the surfaces header
+ const std::array<unsigned,7> types = MeshInformation::getNumberOfElementTypes(mesh);
+ const unsigned nTotalTriangles (types[1] + (2*types[2]));
+ out << nTotalTriangles << " 1\n";
+
+ const std::size_t nElements (elements.size());
+ unsigned count(0);
+ for (std::size_t i=0; i<nElements; ++i)
+ {
+ count++;
+ if (elements[i]->getGeomType() == MeshElemType::TRIANGLE)
+ out << "3 " << elements[i]->getNodeIndex(0) << " " << elements[i]->getNodeIndex(1) << " " << elements[i]->getNodeIndex(2) << " " << elements[i]->getValue() << " # " << count << "\n";
+ else if (elements[i]->getGeomType() == MeshElemType::QUAD)
+ {
+ out << "3 " << elements[i]->getNodeIndex(0) << " " << elements[i]->getNodeIndex(1) << " " << elements[i]->getNodeIndex(2) << " " << elements[i]->getValue() << " # " << count << "\n";
+ count++;
+ out << "3 " << elements[i]->getNodeIndex(0) << " " << elements[i]->getNodeIndex(2) << " " << elements[i]->getNodeIndex(3) << " " << elements[i]->getValue() << " # " << count << "\n";
+ }
+ }
+ out << "0\n"; // the polygon holes list
+
+ // the region attributes list
+ if (attribute_points.empty())
+ out << "0\n";
+ else
+ {
+ const std::size_t nAttributePoints (attribute_points.size());
+ out << nAttributePoints << "\n";
+ for (std::size_t i=0; i<nAttributePoints; ++i)
+ out << i+1 << " " << attribute_points[i][0] << " " << attribute_points[i][1] << " " << attribute_points[i][2] << " " << 10*attribute_points[i].getID() << "\n";
+ }
+ INFO ("TetGenInterface::writeTetGenPoly() - %d points and %d surfaces successfully written.", nPoints, nElements);
+ out.close();
+ return true;
+}
+
} // end namespace FileIO
diff --git a/FileIO/TetGenInterface.h b/FileIO/TetGenInterface.h
index 72ca2c04618a2dabfe7e2134e9633a2879395eea..2e5a38555f8f9600e4178c0d323e49d222e7255f 100644
--- a/FileIO/TetGenInterface.h
+++ b/FileIO/TetGenInterface.h
@@ -43,13 +43,13 @@ public:
virtual ~TetGenInterface();
/**
- * Method reads the TetGen mesh from node file and element file.
- * @param poly_fname file name of the poly file
+ * Method reads geometry from a TetGen poly or smesh file.
+ * @param geo_fname file name of the poly file
* @param geo_objects where the geometry is written to
* @return on success the method returns true, otherwise it returns false
*/
- bool readTetGenPoly (std::string const& poly_fname,
- GeoLib::GEOObjects &geo_objects);
+ bool readTetGenGeometry (std::string const& geo_fname,
+ GeoLib::GEOObjects &geo_objects);
/**
* Method reads the TetGen mesh from node file and element file.
@@ -61,30 +61,58 @@ public:
std::string const& ele_fname);
/**
- * Writes the geometry of a given name to TetGen poly-file.
- * @param file_name file name of the new poly file
- * @param geo_objects the container for the geometry.
- * @param geo_name the name for the geometry.
+ * Writes the geometry of a given name to TetGen smesh-file.
+ * @param file_name file name of the new smesh file.
+ * @param geo_objects the container for the geometry.
+ * @param geo_name the name for the geometry containing the subsurface boundary representation used for meshing.
+ * @param attribute_points attribute points containing material IDs (if the vector is empty no attributes are written).
* @return returns true on success and false otherwise.
*/
- bool writeTetGenPoly(const std::string &file_name,
- const GeoLib::GEOObjects &geo_objects,
- const std::string &geo_name) const;
+ bool writeTetGenSmesh(const std::string &file_name,
+ const GeoLib::GEOObjects &geo_objects,
+ const std::string &geo_name,
+ const std::vector<GeoLib::PointWithID> &attribute_points) const;
+
+ /**
+ * Writes the geometry of a given name to TetGen smesh-file.
+ * @param file_name file name of the new smesh file.
+ * @param mesh mesh containing the subsurface boundary representation used for meshing.
+ * @param attribute_points attribute points containing material IDs (if the vector is empty no attributes are written).
+ * @return returns true on success and false otherwise.
+ */
+ bool writeTetGenSmesh(const std::string &file_name,
+ const MeshLib::Mesh &mesh,
+ const std::vector<MeshLib::Node> &attribute_points) const;
private:
/// Returns the declared number of facets in the poly file.
- std::size_t getNFacets(std::ifstream &input) const;
+ std::size_t getNFacets(std::ifstream &input);
/**
* Method parses the lines reading the facets from TetGen poly file
- * @param input the input stream (input)
- * @param surfaces the vector of surfaces to be filled (output)
- * @param points the point vector needed for creating surfaces (input)
+ * @param input the input stream (input)
+ * @param surfaces the vector of surfaces to be filled (output)
+ * @param points the point vector needed for creating surfaces (input)
+ * @param pnt_id_map the id map to compensate for possibly changed point ids after adding the point vector to GEOObjects
* @return true, if the facets have been read correctly, false if the method detects an error
*/
- bool parseFacets(std::ifstream &input,
- std::vector<GeoLib::Surface*> &surfaces,
- std::vector<GeoLib::Point*> &points) const;
+ bool parsePolyFacets(std::ifstream &input,
+ std::vector<GeoLib::Surface*> &surfaces,
+ const std::vector<GeoLib::Point*> &points,
+ const std::vector<std::size_t> &pnt_id_map);
+
+ /**
+ * Method parses the lines reading the facets from TetGen smesh file
+ * @param input the input stream (input)
+ * @param surfaces the vector of surfaces to be filled (output)
+ * @param points the point vector needed for creating surfaces (input)
+ * @param pnt_id_map the id map to compensate for possibly changed point ids after adding the point vector to GEOObjects
+ * @return true, if the facets have been read correctly, false if the method detects an error
+ */
+ bool parseSmeshFacets(std::ifstream &input,
+ std::vector<GeoLib::Surface*> &surfaces,
+ const std::vector<GeoLib::Point*> &points,
+ const std::vector<std::size_t> &pnt_id_map);
/**
* Method reads the nodes from stream and stores them in a node vector.
@@ -162,10 +190,12 @@ private:
std::size_t n_nodes_per_tet,
bool region_attribute) const;
- /**
- * the value is true if the indexing is zero based, else false
- */
+
+ /// the value is true if the indexing is zero based, else false
bool _zero_based_idx;
+
+ /// true if boundary markers are set, false otherwise
+ bool _boundary_markers;
};
}
diff --git a/Gui/DataView/MeshLayerEditDialog.cpp b/Gui/DataView/MeshLayerEditDialog.cpp
index 5692a07ce325ae833b948273d37df2874fcb4cab..fc2011ff510686da324a4573323e5f5cb493d651 100644
--- a/Gui/DataView/MeshLayerEditDialog.cpp
+++ b/Gui/DataView/MeshLayerEditDialog.cpp
@@ -12,7 +12,7 @@
*
*/
-
+#include "MeshGenerators/LayeredVolume.h"
#include "MeshLayerEditDialog.h"
// ThirdParty/logog
@@ -22,6 +22,8 @@
#include "StringTools.h"
#include "Mesh.h"
+#include "TetGenInterface.h"
+
#include <QCheckBox>
#include <QFileInfo>
#include <QFileDialog>
@@ -30,21 +32,24 @@
#include <QLineEdit>
#include <QGroupBox>
#include <QRadioButton>
+#include <QHBoxLayout>
#include <QVBoxLayout>
MeshLayerEditDialog::MeshLayerEditDialog(const MeshLib::Mesh* mesh, QDialog* parent)
: QDialog(parent), _msh(mesh), _n_layers(0),
- _nLayerExplanation (new QLabel("(select \"0\" for surface mapping)")),
- _layerEdit (new QLineEdit("0")),
+ _nLayerExplanation (new QLabel(this)),
+ _layerEdit (new QLineEdit("0", this)),
_noDataReplacementEdit(nullptr),
- _nextButton (new QPushButton("Next")),
+ _nextButton (new QPushButton("Next", this)),
_layerBox (nullptr),
_radioButtonBox (nullptr),
- _layerSelectionLayout (new QGridLayout),
+ _ogsMeshButton (nullptr),
+ _layerSelectionLayout (new QGridLayout(_layerBox)),
_use_rasters(true)
{
setupUi(this);
+ _nLayerExplanation->setText("(select \"0\" for surface mapping)");
this->gridLayoutLayerMapping->addWidget(new QLabel("Please specify the number of layers to add:", this), 0, 0);
this->gridLayoutLayerMapping->addWidget(_nLayerExplanation, 1, 0);
this->gridLayoutLayerMapping->addWidget(_layerEdit, 0, 1);
@@ -53,7 +58,6 @@ MeshLayerEditDialog::MeshLayerEditDialog(const MeshLib::Mesh* mesh, QDialog* par
connect(_nextButton, SIGNAL(pressed()), this, SLOT(nextButtonPressed()));
// configure group box + layout
- _layerBox = new QGroupBox;
this->_layerSelectionLayout->setMargin(10);
this->_layerSelectionLayout->setColumnMinimumWidth(2,10);
this->_layerSelectionLayout->setColumnStretch(0, 80);
@@ -63,16 +67,6 @@ MeshLayerEditDialog::MeshLayerEditDialog(const MeshLib::Mesh* mesh, QDialog* par
MeshLayerEditDialog::~MeshLayerEditDialog()
{
- for (int i = 0; i < _edits.size(); ++i)
- delete _edits[i];
-
- delete _nLayerExplanation;
- delete _layerEdit;
- delete _noDataReplacementEdit;
- delete _nextButton;
- delete _radioButtonBox;
- delete _layerSelectionLayout;
- delete _layerBox;
}
void MeshLayerEditDialog::nextButtonPressed()
@@ -82,12 +76,14 @@ void MeshLayerEditDialog::nextButtonPressed()
_nLayerExplanation->setText("");
_n_layers = static_cast<unsigned>(_layerEdit->text().toInt());
- if (_n_layers > 0)
+ if (_n_layers == 0)
+ this->createWithRasters();
+ else
{
QVBoxLayout* _radiobuttonLayout (new QVBoxLayout(_radioButtonBox));
QRadioButton* selectButton1 (new QRadioButton("Add layers based on raster files", _radioButtonBox));
QRadioButton* selectButton2 (new QRadioButton("Add layers with static thickness", _radioButtonBox));
- _radioButtonBox = new QGroupBox;
+ _radioButtonBox = new QGroupBox(this);
_radiobuttonLayout->addWidget(selectButton1);
_radiobuttonLayout->addWidget(selectButton2);
_radioButtonBox->setLayout(_radiobuttonLayout);
@@ -97,10 +93,6 @@ void MeshLayerEditDialog::nextButtonPressed()
connect(selectButton1, SIGNAL(pressed()), this, SLOT(createWithRasters()));
connect(selectButton2, SIGNAL(pressed()), this, SLOT(createStatic()));
}
- else
- this->createWithRasters();
-
-
}
void MeshLayerEditDialog::createWithRasters()
@@ -110,7 +102,8 @@ void MeshLayerEditDialog::createWithRasters()
this->_radioButtonBox->setEnabled(false);
const QString selectText = (_n_layers>0) ?
"Please specify a raster file for mapping each layer:" :
- "Please specify rasterfile for surface mapping:";
+ "Please specify raster file for surface mapping:";
+ this->_layerBox = new QGroupBox(this);
this->_layerBox->setTitle(selectText);
for (unsigned i = 0; i <= _n_layers+1; ++i)
@@ -119,7 +112,7 @@ void MeshLayerEditDialog::createWithRasters()
if (i==0) text = "Surface";
else if (i>_n_layers) text = "Layer" + QString::number(_n_layers) + "-Bottom";
else text="Layer" + QString::number(i) + "-Top";
- QLineEdit* edit (new QLineEdit());
+ QLineEdit* edit (new QLineEdit(this));
QPushButton* button (new QPushButton("...", _layerBox));
this->_edits.push_back(edit);
@@ -144,103 +137,165 @@ void MeshLayerEditDialog::createWithRasters()
}
this->_layerBox->setLayout(this->_layerSelectionLayout);
this->gridLayoutLayerMapping->addWidget(_layerBox, 4, 0, 1, 3);
+ if (this->_n_layers > 0)
+ this->createMeshToolSelection();
}
void MeshLayerEditDialog::createStatic()
{
this->_use_rasters = false;
this->_radioButtonBox->setEnabled(false);
+ this->_layerBox = new QGroupBox(this);
this->_layerBox->setTitle("Please specify a thickness or each layer");
for (unsigned i = 0; i < this->_n_layers; ++i)
{
QString text("Layer" + QString::number(i) + "-Thickness");
- QLineEdit* staticLayerEdit = new QLineEdit("10");
+ QLineEdit* staticLayerEdit = new QLineEdit("10", this);
staticLayerEdit->setValidator(new QDoubleValidator(staticLayerEdit));
_edits.push_back(staticLayerEdit);
this->_layerSelectionLayout->addWidget(new QLabel(text, _layerBox), i, 0);
this->_layerSelectionLayout->addWidget(_edits[i], i, 1);
}
this->_layerBox->setLayout(this->_layerSelectionLayout);
- this->gridLayoutLayerMapping->addWidget(_layerBox, 5, 0, 1, 3);
+ this->gridLayoutLayerMapping->addWidget(_layerBox, 4, 0, 1, 3);
+ this->createMeshToolSelection();
}
-void MeshLayerEditDialog::accept()
+void MeshLayerEditDialog::createMeshToolSelection()
{
- if (this->_edits.size()>0)
+ QGroupBox* meshToolSelectionBox (new QGroupBox(this));
+ meshToolSelectionBox->setTitle("Select output element type");
+ QHBoxLayout* meshToolSelectionLayout (new QHBoxLayout(meshToolSelectionBox));
+ _ogsMeshButton = new QRadioButton("Prisms", meshToolSelectionBox);
+ QRadioButton* tetgenMeshButton = new QRadioButton("Tetrahedra", meshToolSelectionBox);
+ meshToolSelectionLayout->addWidget(_ogsMeshButton);
+ meshToolSelectionLayout->addWidget(tetgenMeshButton);
+ meshToolSelectionBox->setLayout(meshToolSelectionLayout);
+ _ogsMeshButton->setChecked(true);
+ gridLayoutLayerMapping->addWidget(meshToolSelectionBox, 5, 0, 1, 3);
+}
+
+MeshLib::Mesh* MeshLayerEditDialog::createPrismMesh()
+{
+ const unsigned nLayers = _layerEdit->text().toInt();
+ std::vector<float> layer_thickness;
+ for (unsigned i=0; i<nLayers; ++i)
{
- bool all_paths_set (true);
- if ((_n_layers==0) && _use_rasters && (_edits[0]->text().length()==0))
- all_paths_set = false;
- else
+ // "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.
+ const float thickness = (_use_rasters) ? 100 : (this->_edits[i]->text().toFloat());
+ layer_thickness.push_back(thickness);
+ }
+
+ MeshLib::Mesh* new_mesh = MeshLayerMapper::CreateLayers(*_msh, layer_thickness);
+
+ if (_use_rasters)
+ {
+ for (unsigned i=0; i<=nLayers; ++i)
+ {
+ const std::string imgPath ( this->_edits[i+1]->text().toStdString() );
+ const double noDataReplacement = (i==0) ? 0.0 : -9999.0;
+ if (!MeshLayerMapper::LayerMapping(*new_mesh, imgPath, nLayers, i, noDataReplacement))
+ {
+ delete new_mesh;
+ return nullptr;
+ }
+ }
+ if (this->_edits[0]->text().length()>0)
{
- int start_idx = (_use_rasters) ? 1:0;
- for (int i=start_idx; i<_edits.size(); ++i)
- if (_edits[i]->text().length()==0)
- all_paths_set = false;
+ MeshLib::Mesh* final_mesh = MeshLayerMapper::blendLayersWithSurface(*new_mesh, nLayers, this->_edits[0]->text().toStdString());
+ delete new_mesh;
+ new_mesh = final_mesh;
}
+ }
+ return new_mesh;
+}
- if (all_paths_set)
+MeshLib::Mesh* MeshLayerEditDialog::createTetMesh()
+{
+ const unsigned nLayers = _layerEdit->text().toInt();
+ std::vector<std::string> raster_paths(nLayers+1);
+ for (unsigned i=0; i<=nLayers; ++i)
+ raster_paths[i] = this->_edits[i+1]->text().toStdString();
+ LayeredVolume lv;
+ lv.createGeoVolumes(*_msh, raster_paths);
+
+ MeshLib::Mesh* tg_mesh (lv.getMesh());
+
+ QString file_path("");
+ if (tg_mesh)
+ {
+ QSettings settings;
+ QString filename = QFileDialog::getSaveFileName(this, "Write TetGen input file to",
+ settings.value("lastOpenedTetgenFileDirectory").toString(),
+ "TetGen Geometry (*.smesh)");
+ if (!filename.isEmpty())
{
- int result(1);
- const unsigned nLayers = _layerEdit->text().toInt();
- MeshLib::Mesh* new_mesh (NULL);
+ std::vector<MeshLib::Node> tg_attr (lv.getAttributePoints());
+ FileIO::TetGenInterface tetgen_interface;
+ tetgen_interface.writeTetGenSmesh(filename.toStdString(), *tg_mesh, tg_attr);
+ }
+ }
+ return tg_mesh;
+}
- if (nLayers==0)
- {
- new_mesh = new MeshLib::Mesh(*_msh);
- const std::string imgPath ( this->_edits[0]->text().toStdString() );
- const double noDataReplacementValue = this->_noDataReplacementEdit->text().toDouble();
- if (!imgPath.empty())
- result = MeshLayerMapper::LayerMapping(*new_mesh, imgPath, nLayers, 0, noDataReplacementValue);
- }
- else
- {
- std::vector<float> layer_thickness;
- for (unsigned i=0; i<nLayers; ++i)
- {
- // "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.
- const float thickness = (_use_rasters) ? 100 : (this->_edits[i]->text().toFloat());
- layer_thickness.push_back(thickness);
- }
-
- new_mesh = MeshLayerMapper::CreateLayers(*_msh, layer_thickness);
-
- if (_use_rasters)
- {
- for (unsigned i=0; i<=nLayers; ++i)
- {
- const std::string imgPath ( this->_edits[i+1]->text().toStdString() );
- const double noDataReplacement = (i==0) ? 0.0 : -9999.0;
- if (!imgPath.empty())
- {
- result = MeshLayerMapper::LayerMapping(*new_mesh, imgPath, nLayers, i, noDataReplacement);
- if (result==0) break;
- }
- }
- if (this->_edits[0]->text().length()>0)
- {
- MeshLib::Mesh* final_mesh = MeshLayerMapper::blendLayersWithSurface(*new_mesh, nLayers, this->_edits[0]->text().toStdString());
- delete new_mesh;
- new_mesh = final_mesh;
- }
- }
- }
+void MeshLayerEditDialog::accept()
+{
+ if (this->_edits.isEmpty())
+ {
+ OGSError::box("Please specifiy the number and\n type of layers and press \"Next\"");
+ return;
+ }
- if (new_mesh)
- emit mshEditFinished(new_mesh);
+ bool all_paths_set (true);
+ if (_n_layers==0)
+ {
+ if (_edits[0]->text().isEmpty())
+ all_paths_set = false;
+ }
+ else
+ {
+ int start_idx = (_use_rasters) ? 1:0;
+ for (int i=start_idx; i<_edits.size(); ++i)
+ if (_edits[i]->text().isEmpty())
+ all_paths_set = false;
+ }
- if (!new_mesh || result==0)
- OGSError::box("Error creating mesh");
+ if (!all_paths_set)
+ {
+ OGSError::box("Please specifiy raster files for all layers.");
+ return;
+ }
+
+ const unsigned nLayers = _layerEdit->text().toInt();
+ MeshLib::Mesh* new_mesh (NULL);
- this->done(QDialog::Accepted);
+ if (nLayers==0)
+ {
+ new_mesh = new MeshLib::Mesh(*_msh);
+ const std::string imgPath ( this->_edits[0]->text().toStdString() );
+ const double noDataReplacementValue = this->_noDataReplacementEdit->text().toDouble();
+ if (!MeshLayerMapper::LayerMapping(*new_mesh, imgPath, nLayers, 0, noDataReplacementValue))
+ {
+ delete new_mesh;
+ return;
}
+ }
+ else
+ {
+ if (_ogsMeshButton->isChecked())
+ new_mesh = this->createPrismMesh();
else
- OGSError::box("Please specifiy raster files for all layers.");
+ new_mesh = this->createTetMesh();
}
+
+ if (new_mesh)
+ emit mshEditFinished(new_mesh);
else
- OGSError::box("Please specifiy the number and\n type of layers and press \"Next\"");
+ OGSError::box("Error creating mesh");
+
+ this->done(QDialog::Accepted);
}
void MeshLayerEditDialog::reject()
diff --git a/Gui/DataView/MeshLayerEditDialog.h b/Gui/DataView/MeshLayerEditDialog.h
index 89fabac92366e00ed6a5e243e4f2f92d362ffa1c..a4dcde19d5c72264dcef0c791efaa8e339cd9177 100644
--- a/Gui/DataView/MeshLayerEditDialog.h
+++ b/Gui/DataView/MeshLayerEditDialog.h
@@ -44,6 +44,10 @@ public:
~MeshLayerEditDialog(void);
private:
+ void createMeshToolSelection();
+ MeshLib::Mesh* createPrismMesh();
+ MeshLib::Mesh* createTetMesh();
+
const MeshLib::Mesh* _msh;
unsigned _n_layers;
QMap<QPushButton*, QLineEdit*> _fileButtonMap;
@@ -55,6 +59,7 @@ private:
QPushButton* _nextButton;
QGroupBox* _layerBox;
QGroupBox* _radioButtonBox;
+ QRadioButton* _ogsMeshButton;
QGridLayout* _layerSelectionLayout;
bool _use_rasters;
diff --git a/Gui/VtkVis/VtkCompositeColorByHeightFilter.cpp b/Gui/VtkVis/VtkCompositeColorByHeightFilter.cpp
index 1f11aac763050b6f5ab60284bfb494e6e610c2d4..50a9f4b139a807b5f2cf066cb5faa8631c15b4c4 100644
--- a/Gui/VtkVis/VtkCompositeColorByHeightFilter.cpp
+++ b/Gui/VtkVis/VtkCompositeColorByHeightFilter.cpp
@@ -47,14 +47,16 @@ void VtkCompositeColorByHeightFilter::init()
unsigned char a[4] = { 0, 0, 255, 255 }; // blue
unsigned char b[4] = { 0, 255, 0, 255 }; // green
unsigned char c[4] = { 255, 255, 0, 255 }; // yellow
- unsigned char d[4] = { 255, 0, 0, 255 }; // red
+ unsigned char d[4] = { 155, 100, 50, 255 }; // brown
+ unsigned char e[4] = { 255, 0, 0, 255 }; // red
VtkColorLookupTable* ColorLookupTable = heightFilter->GetColorLookupTable();
ColorLookupTable->setInterpolationType(VtkColorLookupTable::LUTType::LINEAR);
- ColorLookupTable->setColor(-35, a);
- ColorLookupTable->setColor(150, b); // green at about 150m
- ColorLookupTable->setColor(450, c); // yellow at about 450m and changing to red from then on
- ColorLookupTable->setColor(800, d);
- ColorLookupTable->SetTableRange(-35, 800);
+ ColorLookupTable->setColor(-50, a);
+ ColorLookupTable->setColor(200, b); // green at about 20m
+ ColorLookupTable->setColor(500, c); // yellow at about 500m and changing to red from then on
+ ColorLookupTable->setColor(1000, d);
+ ColorLookupTable->setColor(2000, e);
+ ColorLookupTable->SetTableRange(-35, 2000);
ColorLookupTable->Build();
// This passes ownership of the ColorLookupTable to VtkVisPointSetItem
diff --git a/Gui/mainwindow.cpp b/Gui/mainwindow.cpp
index 8dade0950b9306867d5d61e43af2aeaba4364ea4..a9c8178bd0310627d58ff8359840c500b7f341f6 100644
--- a/Gui/mainwindow.cpp
+++ b/Gui/mainwindow.cpp
@@ -687,18 +687,17 @@ void MainWindow::loadFile(ImportFileType::type t, const QString &fileName)
}
else if (t == ImportFileType::TETGEN)
{
- if (fi.suffix().toLower().compare("poly") == 0)
+ if (fi.suffix().toLower().compare("poly") == 0 || fi.suffix().toLower().compare("smesh") == 0)
{
FileIO::TetGenInterface tetgen;
- tetgen.readTetGenPoly(fileName.toStdString(), *(_project.getGEOObjects()));
+ tetgen.readTetGenGeometry(fileName.toStdString(), *(_project.getGEOObjects()));
}
else {
settings.setValue("lastOpenedTetgenFileDirectory", QFileInfo(fileName).absolutePath());
- QString element_fname = QFileDialog::getOpenFileName(this, "Select TetGen element file",
- settings.value("lastOpenedTetgenFileDirectory").toString(),
- "TetGen element files (*.ele);;");
+ QString element_fname(fi.path() + "\\" + fi.completeBaseName() + ".ele");
- if (!fileName.isEmpty() && !element_fname.isEmpty()) {
+ if (!fileName.isEmpty())
+ {
FileIO::TetGenInterface tetgen;
MeshLib::Mesh* mesh (tetgen.readTetGenMesh(fileName.toStdString(), element_fname.toStdString()));
if (mesh)
@@ -1029,7 +1028,7 @@ void MainWindow::mapGeometry(const std::string &geo_name)
void MainWindow::convertMeshToGeometry(const MeshLib::Mesh* mesh)
{
- MeshLib::convertMeshToGeo(*mesh, this->_project.getGEOObjects());
+ MeshLib::convertMeshToGeo(*mesh, *this->_project.getGEOObjects());
}
void MainWindow::exportBoreholesToGMS(std::string listName, std::string fileName)
diff --git a/MeshLib/MeshGenerators/LayeredVolume.cpp b/MeshLib/MeshGenerators/LayeredVolume.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f8be11ebddd34a3c0c83d06ff722d7b13d379316
--- /dev/null
+++ b/MeshLib/MeshGenerators/LayeredVolume.cpp
@@ -0,0 +1,232 @@
+/**
+ * \file LayeredVolume.cpp
+ * \author Karsten Rink
+ * \date 2014-04-11
+ * \brief Implementation of the LayeredVolume class.
+ *
+ * \copyright
+ * Copyright (c) 2013, 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 "LayeredVolume.h"
+
+#include <fstream>
+#include <numeric>
+
+#include "Vector3.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(GeoLib::Raster::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)
+{
+ if (mesh.getDimension() != 2)
+ return false;
+
+ _elevation_epsilon = calcEpsilon(*rasters[0], *rasters.back());
+ if (_elevation_epsilon <= 0)
+ return false;
+
+ // remove line elements, only tri + quad remain
+ MeshLib::ElementExtraction ex(mesh);
+ ex.searchByElementType(MeshElemType::LINE);
+ MeshLib::Mesh* mesh_layer (ex.removeMeshElements("MeshLayer"));
+ if (mesh_layer==nullptr)
+ mesh_layer = new MeshLib::Mesh(mesh);
+
+ // map each layer and attach to subsurface mesh
+ const std::size_t nRasters (rasters.size());
+ std::vector<GeoLib::Point> in_region_points(nRasters-1, GeoLib::Point(0,0,0));
+ for (size_t i=0; i<nRasters; ++i)
+ {
+ const double replacement_value = (i==0) ? noDataReplacementValue : _invalid_value;
+ if (!MeshLayerMapper::LayerMapping(*mesh_layer, *rasters[i], 0, 0, replacement_value))
+ {
+ this->cleanUpOnError();
+ return false;
+ }
+ this->addLayerToMesh(*mesh_layer, 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)
+{
+ 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);
+
+ for (std::size_t i=0; i<nNodes; ++i)
+ {
+ if (layer_id > 0 &&
+ ((*layer_nodes[i])[2] == _invalid_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
+ _nodes.push_back(new MeshLib::Node(layer_nodes[i]->getCoords(), _nodes.size()));
+ }
+
+ const std::vector<MeshLib::Element*> &layer_elements (mesh_layer.getElements());
+ for (MeshLib::Element* elem : layer_elements)
+ {
+ if (elem->getGeomType() == MeshElemType::TRIANGLE)
+ {
+ std::array<MeshLib::Node*,3> tri_nodes = { _nodes[node_id_offset+elem->getNodeIndex(0)],
+ _nodes[node_id_offset+elem->getNodeIndex(1)],
+ _nodes[node_id_offset+elem->getNodeIndex(2)] };
+ _elements.push_back(new MeshLib::Tri(tri_nodes, layer_id+1));
+ }
+ else if (elem->getGeomType() == MeshElemType::QUAD)
+ {
+ std::array<MeshLib::Node*,4> quad_nodes = { _nodes[node_id_offset+elem->getNodeIndex(0)],
+ _nodes[node_id_offset+elem->getNodeIndex(1)],
+ _nodes[node_id_offset+elem->getNodeIndex(2)],
+ _nodes[node_id_offset+elem->getNodeIndex(3)] };
+ _elements.push_back(new MeshLib::Quad(quad_nodes, layer_id+1));
+ }
+ }
+}
+
+void LayeredVolume::addLayerBoundaries(const MeshLib::Mesh &layer, std::size_t nLayers)
+{
+ const unsigned nLayerBoundaries (nLayers-1);
+ const std::size_t nNodes (layer.getNNodes());
+ const std::vector<MeshLib::Element*> &layer_elements (layer.getElements());
+ for (MeshLib::Element* elem : layer_elements)
+ {
+ const std::size_t nElemNodes (elem->getNNodes());
+ for (unsigned i=0; i<nElemNodes; ++i)
+ if (elem->getNeighbor(i) == nullptr)
+ for (unsigned j=0; j<nLayerBoundaries; ++j)
+ {
+ const std::size_t offset (j*nNodes);
+ MeshLib::Node* n0 = _nodes[offset + elem->getNodeIndex(i)];
+ MeshLib::Node* n1 = _nodes[offset + elem->getNodeIndex((i+1)%nElemNodes)];
+ MeshLib::Node* n2 = _nodes[offset + nNodes + elem->getNodeIndex((i+1)%nElemNodes)];
+ MeshLib::Node* n3 = _nodes[offset + nNodes + elem->getNodeIndex(i)];
+
+ if (MathLib::Vector3(*n1, *n2).getLength() > std::numeric_limits<double>::epsilon())
+ {
+ const std::array<MeshLib::Node*,3> tri_nodes = { n0, n1, n2 };
+ _elements.push_back(new MeshLib::Tri(tri_nodes, nLayers+j));
+ }
+ if (MathLib::Vector3(*n0, *n3).getLength() > std::numeric_limits<double>::epsilon())
+ {
+ const std::array<MeshLib::Node*,3> tri_nodes = { n0, n2, n3 };
+ _elements.push_back(new MeshLib::Tri(tri_nodes, nLayers+j));
+ }
+ }
+ }
+}
+
+void LayeredVolume::removeCongruentElements(std::size_t nLayers, std::size_t nElementsPerLayer)
+{
+ for (std::size_t i=1; i<nLayers; ++i)
+ {
+ const std::size_t upper_offset ((i-1) * nElementsPerLayer);
+ const std::size_t lower_offset ( i * nElementsPerLayer);
+ for (std::size_t j=0; j<nElementsPerLayer; ++j)
+ {
+ MeshLib::Element const*const high (_elements[upper_offset+j]);
+ MeshLib::Element *low (_elements[lower_offset+j]);
+
+ unsigned count(0);
+ const std::size_t nElemNodes (low->getNNodes());
+ for (std::size_t k=0; k<nElemNodes; ++k)
+ if (high->getNodeIndex(k) == low->getNodeIndex(k))
+ {
+ low->setNode(k, _nodes[high->getNodeIndex(k)]);
+ count++;
+ }
+
+ if (count == nElemNodes)
+ {
+ delete _elements[upper_offset+j];
+ _elements[upper_offset+j] = nullptr;
+ }
+ else
+ {
+ MeshLib::Node attr = high->getCenterOfGravity();
+ _attribute_points.push_back(MeshLib::Node(attr[0], attr[1], (attr[2] + low->getCenterOfGravity()[2])/2.0, low->getValue()));
+ }
+ }
+ }
+ 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
new file mode 100644
index 0000000000000000000000000000000000000000..49deb3ec6ea156fb9a950e3ddde14af08f57892b
--- /dev/null
+++ b/MeshLib/MeshGenerators/LayeredVolume.h
@@ -0,0 +1,97 @@
+/**
+ * \file LayeredVolume.h
+ * \author Karsten Rink
+ * \date 2014-04-11
+ * \brief Definition of the LayeredVolume class
+ *
+ * \copyright
+ * Copyright (c) 2013, 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 LAYEREDVOLUME_H
+#define LAYEREDVOLUME_H
+
+#include <string>
+#include <vector>
+
+#include "Raster.h"
+#include "Node.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
+{
+public:
+ LayeredVolume();
+ ~LayeredVolume() {}
+
+ /**
+ * 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 rasters Containing all the raster-data for the subsurface layers from top to bottom (starting 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<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 top to bottom (starting 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; }
+
+ /// 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);
+
+ /// Creates boundary surfaces between the mapped layers to make the volumes watertight
+ void addLayerBoundaries(const MeshLib::Mesh &layer, std::size_t nLayers);
+
+ /// 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 59c600ef2878f7bcd62adf4b4e1e5815003a382c..fa8788d0bdd46ad71b59c25f41791cfb21b307f8 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.cpp
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
@@ -109,27 +109,35 @@ MeshLib::Mesh* MeshLayerMapper::CreateLayers(const MeshLib::Mesh &mesh, const st
return new MeshLib::Mesh("SubsurfaceMesh", new_nodes, new_elems);
}
-int MeshLayerMapper::LayerMapping(MeshLib::Mesh &new_mesh, const std::string &rasterfile,
- const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue = 0.0)
+bool MeshLayerMapper::LayerMapping(MeshLib::Mesh &new_mesh, const std::string &rasterfile,
+ const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue = 0.0)
+{
+ const GeoLib::Raster *raster(GeoLib::Raster::getRasterFromASCFile(rasterfile));
+ if (! raster) {
+ 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);
+ 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)
{
if (nLayers < layer_id)
{
ERR("MshLayerMapper::LayerMapping() - Mesh has only %d Layers, cannot assign layer %d.", nLayers, layer_id);
- return 0;
+ return false;
}
- 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 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
@@ -201,8 +209,7 @@ int MeshLayerMapper::LayerMapping(MeshLib::Mesh &new_mesh, const std::string &ra
nodes[i]->updateCoordinates(coords[0], coords[1], z);
}
- delete raster;
- return 1;
+ return true;
}
bool MeshLayerMapper::isNodeOnRaster(const MeshLib::Node &node,
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.h b/MeshLib/MeshGenerators/MeshLayerMapper.h
index bcbd9a7da94e38a561fa583aec776ed439c7115d..cb8c2cb6b5c60fc07bf8af1e09efb231dca09dc3 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.h
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.h
@@ -16,6 +16,7 @@
#define MESHLAYERMAPPER_H
#include <string>
+#include "Raster.h"
class QImage;
@@ -46,8 +47,15 @@ public:
* 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 &mesh, const std::string &rasterfile,
- const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue);
+ static bool LayerMapping(MeshLib::Mesh &mesh, const std::string &rasterfile,
+ const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue);
+
+ /**
+ * 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 bool LayerMapping(MeshLib::Mesh &mesh, const GeoLib::Raster &raster,
+ const unsigned nLayers, const unsigned layer_id, double noDataReplacementValue);
/**
* 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.
diff --git a/MeshLib/Node.h b/MeshLib/Node.h
index 8bdb813c93418f5dc983fde85c54f4038c2f28a4..48684bb6d22e1bacbe86980c3f1d3fde082e9c8c 100644
--- a/MeshLib/Node.h
+++ b/MeshLib/Node.h
@@ -36,8 +36,8 @@ class Element;
class Node : public GeoLib::PointWithID
{
/* friend functions: */
- friend int MeshLayerMapper::LayerMapping(MeshLib::Mesh &mesh, const std::string &rasterfile, const unsigned nLayers,
- const unsigned layer_id, double noDataReplacementValue);
+ friend bool MeshLayerMapper::LayerMapping(MeshLib::Mesh &mesh, const GeoLib::Raster &raster, const unsigned nLayers,
+ const unsigned layer_id, double noDataReplacementValue);
friend MeshLib::Mesh* MeshLayerMapper::blendLayersWithSurface(MeshLib::Mesh &mesh, const unsigned nLayers, const std::string &dem_raster);
/* friend classes: */
diff --git a/MeshLib/convertMeshToGeo.cpp b/MeshLib/convertMeshToGeo.cpp
index 4475e919267032392ddeb7e54595c6eca39cdbcc..8642338dccfa69000ef18a006b6dc1d2e035910f 100644
--- a/MeshLib/convertMeshToGeo.cpp
+++ b/MeshLib/convertMeshToGeo.cpp
@@ -24,10 +24,11 @@
#include "Elements/Tri.h"
#include "Elements/Quad.h"
#include "Node.h"
+#include "MeshInformation.h"
namespace MeshLib {
-bool convertMeshToGeo(const MeshLib::Mesh &mesh, GeoLib::GEOObjects* geo_objects)
+bool convertMeshToGeo(const MeshLib::Mesh &mesh, GeoLib::GEOObjects &geo_objects, double eps)
{
if (mesh.getDimension() != 2)
{
@@ -44,31 +45,38 @@ bool convertMeshToGeo(const MeshLib::Mesh &mesh, GeoLib::GEOObjects* geo_objects
(*points)[i] = new GeoLib::Point(static_cast<GeoLib::Point>(*nodes[i]));
std::string mesh_name (mesh.getName());
- geo_objects->addPointVec(points, mesh_name);
- const std::vector<std::size_t> id_map (geo_objects->getPointVecObj(mesh_name)->getIDMap());
+ geo_objects.addPointVec(points, mesh_name, nullptr, eps);
+ const std::vector<std::size_t> id_map (geo_objects.getPointVecObj(mesh_name)->getIDMap());
// elements to surface triangles conversion
+ const std::pair<unsigned, unsigned> bounds (MeshInformation::getValueBounds(mesh));
+ const unsigned nMatGroups(bounds.second-bounds.first+1);
+ std::vector<GeoLib::Surface*> *sfcs = new std::vector<GeoLib::Surface*>;
+ sfcs->reserve(nMatGroups);
+ for (unsigned i=0; i<nMatGroups; ++i)
+ sfcs->push_back(new GeoLib::Surface(*points));
+
const std::vector<MeshLib::Element*> &elements = mesh.getElements();
- GeoLib::Surface* sfc = new GeoLib::Surface(*points);
const std::size_t nElems (mesh.getNElements());
for (unsigned i=0; i<nElems; ++i)
{
MeshLib::Element* e (elements[i]);
if (e->getGeomType() == MeshElemType::TRIANGLE)
- sfc->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(1)], id_map[e->getNodeIndex(2)]);
+ (*sfcs)[e->getValue()-bounds.first]->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(1)], id_map[e->getNodeIndex(2)]);
if (e->getGeomType() == MeshElemType::QUAD)
{
- sfc->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(1)], id_map[e->getNodeIndex(2)]);
- sfc->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(2)], id_map[e->getNodeIndex(3)]);
+ (*sfcs)[e->getValue()-bounds.first]->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(1)], id_map[e->getNodeIndex(2)]);
+ (*sfcs)[e->getValue()-bounds.first]->addTriangle(id_map[e->getNodeIndex(0)], id_map[e->getNodeIndex(2)], id_map[e->getNodeIndex(3)]);
}
// all other element types are ignored (i.e. lines)
}
- std::vector<GeoLib::Surface*> *sfcs = new std::vector<GeoLib::Surface*>(1);
- (*sfcs)[0] = sfc;
+ std::for_each(sfcs->begin(), sfcs->end(), [](GeoLib::Surface* sfc){ if (sfc->getNTriangles()==0) delete sfc; sfc = nullptr;});
+ auto sfcs_end = std::remove(sfcs->begin(), sfcs->end(), nullptr);
+ sfcs->erase(sfcs_end, sfcs->end());
- geo_objects->addSurfaceVec(sfcs, mesh_name);
+ geo_objects.addSurfaceVec(sfcs, mesh_name);
return true;
}
diff --git a/MeshLib/convertMeshToGeo.h b/MeshLib/convertMeshToGeo.h
index 2953460aa6526ff203b35999e2994dcf18f9d50b..54883619b265cd5e7ad90472e0a7292ca0892bd3 100644
--- a/MeshLib/convertMeshToGeo.h
+++ b/MeshLib/convertMeshToGeo.h
@@ -15,6 +15,8 @@
#ifndef CONVERTMESHTOGEO_H_
#define CONVERTMESHTOGEO_H_
+#include <limits>
+
namespace GeoLib
{
class GEOObjects;
@@ -32,7 +34,7 @@ namespace MeshLib
* converted to geometric triangles, quads are split into two triangles, all other elements
* are ignored.
*/
- bool convertMeshToGeo(const MeshLib::Mesh &mesh, GeoLib::GEOObjects* geo_objects);
+ bool convertMeshToGeo(const MeshLib::Mesh &mesh, GeoLib::GEOObjects &geo_objects, double eps = std::numeric_limits<double>::epsilon());
} // namespace