From d3b24c03ecdfb8006310031d004a2288a54cee57 Mon Sep 17 00:00:00 2001
From: Karsten Rink <karsten.rink@ufz.de>
Date: Tue, 21 Aug 2012 10:51:22 +0200
Subject: [PATCH] current state of integration for LB (fixing netcdf include
error)
---
FemLib/BoundaryCondition.cpp | 26 ++
FemLib/BoundaryCondition.h | 33 ++
FemLib/InitialCondition.cpp | 18 +
FemLib/InitialCondition.h | 27 ++
FemLib/SourceTerm.cpp | 62 +++
FemLib/SourceTerm.h | 37 ++
FileIO/GMSInterface.cpp | 314 +++++++++++++
FileIO/GMSInterface.h | 64 +++
FileIO/Gmsh2GeoIO.cpp | 123 +++++
FileIO/Gmsh2GeoIO.h | 33 ++
FileIO/MeshIO/GMSHAdaptiveMeshDensity.cpp | 151 +++++++
FileIO/MeshIO/GMSHAdaptiveMeshDensity.h | 49 ++
FileIO/MeshIO/GMSHFixedMeshDensity.cpp | 30 ++
FileIO/MeshIO/GMSHFixedMeshDensity.h | 28 ++
FileIO/MeshIO/GMSHInterface.cpp | 272 +++++++++++
FileIO/MeshIO/GMSHInterface.h | 110 +++++
FileIO/MeshIO/GMSHLine.cpp | 30 ++
FileIO/MeshIO/GMSHLine.h | 29 ++
FileIO/MeshIO/GMSHLineLoop.cpp | 47 ++
FileIO/MeshIO/GMSHLineLoop.h | 34 ++
FileIO/MeshIO/GMSHMeshDensityStrategy.h | 33 ++
FileIO/MeshIO/GMSHNoMeshDensity.h | 34 ++
FileIO/MeshIO/GMSHPoint.cpp | 38 ++
FileIO/MeshIO/GMSHPoint.h | 30 ++
FileIO/MeshIO/GMSHPolygonTree.cpp | 301 ++++++++++++
FileIO/MeshIO/GMSHPolygonTree.h | 91 ++++
FileIO/MeshIO/TetGenInterface.cpp | 528 ++++++++++++++++++++++
FileIO/MeshIO/TetGenInterface.h | 119 +++++
FileIO/NetCDFInterface.cpp | 204 +++++++++
FileIO/NetCDFInterface.h | 43 ++
FileIO/PetrelInterface.cpp | 234 ++++++++++
FileIO/PetrelInterface.h | 40 ++
Gui/DataView/CondFromRasterDialog.cpp | 21 +-
Gui/DataView/CondFromRasterDialog.h | 4 +-
Gui/DataView/DirectConditionGenerator.h | 3 +-
Gui/DataView/FEMConditionSetupDialog.cpp | 8 +-
Gui/DataView/FEMConditionSetupDialog.h | 6 +-
Gui/mainwindow.cpp | 14 +-
scripts/cmake/Find.cmake | 2 +-
39 files changed, 3241 insertions(+), 29 deletions(-)
create mode 100644 FemLib/BoundaryCondition.cpp
create mode 100644 FemLib/BoundaryCondition.h
create mode 100644 FemLib/InitialCondition.cpp
create mode 100644 FemLib/InitialCondition.h
create mode 100644 FemLib/SourceTerm.cpp
create mode 100644 FemLib/SourceTerm.h
create mode 100644 FileIO/GMSInterface.cpp
create mode 100644 FileIO/GMSInterface.h
create mode 100644 FileIO/Gmsh2GeoIO.cpp
create mode 100644 FileIO/Gmsh2GeoIO.h
create mode 100644 FileIO/MeshIO/GMSHAdaptiveMeshDensity.cpp
create mode 100644 FileIO/MeshIO/GMSHAdaptiveMeshDensity.h
create mode 100644 FileIO/MeshIO/GMSHFixedMeshDensity.cpp
create mode 100644 FileIO/MeshIO/GMSHFixedMeshDensity.h
create mode 100644 FileIO/MeshIO/GMSHInterface.cpp
create mode 100644 FileIO/MeshIO/GMSHInterface.h
create mode 100644 FileIO/MeshIO/GMSHLine.cpp
create mode 100644 FileIO/MeshIO/GMSHLine.h
create mode 100644 FileIO/MeshIO/GMSHLineLoop.cpp
create mode 100644 FileIO/MeshIO/GMSHLineLoop.h
create mode 100644 FileIO/MeshIO/GMSHMeshDensityStrategy.h
create mode 100644 FileIO/MeshIO/GMSHNoMeshDensity.h
create mode 100644 FileIO/MeshIO/GMSHPoint.cpp
create mode 100644 FileIO/MeshIO/GMSHPoint.h
create mode 100644 FileIO/MeshIO/GMSHPolygonTree.cpp
create mode 100644 FileIO/MeshIO/GMSHPolygonTree.h
create mode 100644 FileIO/MeshIO/TetGenInterface.cpp
create mode 100644 FileIO/MeshIO/TetGenInterface.h
create mode 100644 FileIO/NetCDFInterface.cpp
create mode 100644 FileIO/NetCDFInterface.h
create mode 100644 FileIO/PetrelInterface.cpp
create mode 100644 FileIO/PetrelInterface.h
diff --git a/FemLib/BoundaryCondition.cpp b/FemLib/BoundaryCondition.cpp
new file mode 100644
index 00000000000..00060e94552
--- /dev/null
+++ b/FemLib/BoundaryCondition.cpp
@@ -0,0 +1,26 @@
+/**
+ * \file BoundaryCondition.cpp
+ * 2011/08/30 KR inital implementation
+ *
+ */
+
+#include "BoundaryCondition.h"
+#include "rf_bc_new.h"
+
+BoundaryCondition::BoundaryCondition(const CBoundaryCondition &bc, const std::string &geometry_name)
+ : FEMCondition(geometry_name, bc.getProcessType(), bc.getProcessPrimaryVariable(),
+ bc.getGeoType(), bc.getGeoName(),
+ bc.getProcessDistributionType(), FEMCondition::BOUNDARY_CONDITION)
+{
+ if (this->getProcessDistributionType() == FiniteElement::CONSTANT ||
+ this->getProcessDistributionType() == FiniteElement::CONSTANT_NEUMANN)
+ this->setConstantDisValue(bc.getGeoNodeValue());
+ else if (this->getProcessDistributionType() == FiniteElement::LINEAR ||
+ this->getProcessDistributionType() == FiniteElement::LINEAR_NEUMANN)
+ {
+ const std::vector<int> bc_nodes(bc.getPointsWithDistribedBC());
+ std::vector<size_t> dis_nodes(bc_nodes.size());
+ for (size_t i=0; i<dis_nodes.size(); i++) dis_nodes[i] = static_cast<size_t>(bc_nodes[i]);
+ this->setDisValues(dis_nodes, bc.getDistribedBC());
+ }
+}
diff --git a/FemLib/BoundaryCondition.h b/FemLib/BoundaryCondition.h
new file mode 100644
index 00000000000..1fcd0115700
--- /dev/null
+++ b/FemLib/BoundaryCondition.h
@@ -0,0 +1,33 @@
+/**
+ * \file BoundaryCondition.h
+ * 2011/08/30 KR inital implementation
+ *
+ */
+
+#ifndef BOUNDARYCONDITION_H
+#define BOUNDARYCONDITION_H
+
+#include "FEMCondition.h"
+
+/**
+ * \brief Adapter class for handling boundary conditions in the user Interface
+ * \sa FEMCondition
+ */
+class BoundaryCondition : public FEMCondition
+{
+public:
+ BoundaryCondition(const std::string &geometry_name)
+ : FEMCondition(geometry_name, FEMCondition::BOUNDARY_CONDITION), _tim_type(0) {};
+ BoundaryCondition(const CBoundaryCondition &bc, const std::string &geometry_name);
+ BoundaryCondition(const FEMCondition &cond)
+ : FEMCondition(cond, FEMCondition::BOUNDARY_CONDITION) {};
+ ~BoundaryCondition() {}
+
+ size_t getTimType() const {return _tim_type; }
+ void setTimType(size_t value) { _tim_type = value; }
+
+private:
+ size_t _tim_type;
+};
+
+#endif //BOUNDARYCONDITION_H
diff --git a/FemLib/InitialCondition.cpp b/FemLib/InitialCondition.cpp
new file mode 100644
index 00000000000..06e3a315d60
--- /dev/null
+++ b/FemLib/InitialCondition.cpp
@@ -0,0 +1,18 @@
+/**
+ * \file InitialCondition.cpp
+ * 2011/08/30 KR inital implementation
+ *
+ */
+
+#include "InitialCondition.h"
+#include "rf_ic_new.h"
+
+InitialCondition::InitialCondition(const CInitialCondition &ic, const std::string &geometry_name)
+ : FEMCondition(geometry_name, ic.getProcessType(), ic.getProcessPrimaryVariable(),
+ ic.getGeoType(),
+ (ic.getGeoType() == GEOLIB::GEODOMAIN) ? "Domain" : ic.getGeoName(),
+ ic.getProcessDistributionType(), FEMCondition::INITIAL_CONDITION)
+{
+ if (this->getProcessDistributionType() == FiniteElement::CONSTANT)
+ this->setConstantDisValue(ic.getGeoNodeValue());
+}
diff --git a/FemLib/InitialCondition.h b/FemLib/InitialCondition.h
new file mode 100644
index 00000000000..6825edc59c4
--- /dev/null
+++ b/FemLib/InitialCondition.h
@@ -0,0 +1,27 @@
+/**
+ * \file InitialCondition.h
+ * 2011/08/30 KR inital implementation
+ *
+ */
+
+#ifndef INITIALCONDITION_H
+#define INITIALCONDITION_H
+
+#include "FEMCondition.h"
+
+/**
+ * \brief Adapter class for handling boundary conditions in the user Interface
+ * \sa FEMCondition
+ */
+class InitialCondition : public FEMCondition
+{
+public:
+ InitialCondition(const std::string &geometry_name)
+ : FEMCondition(geometry_name, FEMCondition::INITIAL_CONDITION) {};
+ InitialCondition(const CInitialCondition &ic, const std::string &geometry_name);
+ InitialCondition(const FEMCondition &cond)
+ : FEMCondition(cond, FEMCondition::INITIAL_CONDITION) {};
+ ~InitialCondition() {}
+};
+
+#endif //INITIALCONDITION_H
diff --git a/FemLib/SourceTerm.cpp b/FemLib/SourceTerm.cpp
new file mode 100644
index 00000000000..a5e194190a8
--- /dev/null
+++ b/FemLib/SourceTerm.cpp
@@ -0,0 +1,62 @@
+/**
+ * \file SourceTerm.cpp
+ * 2011/08/30 KR inital implementation
+ *
+ */
+
+#include "SourceTerm.h"
+#include "rf_st_new.h"
+
+SourceTerm::SourceTerm(const CSourceTerm &st, const std::string &geometry_name)
+ : FEMCondition(geometry_name, st.getProcessType(), st.getProcessPrimaryVariable(),
+ st.getGeoType(), st.getGeoName(),
+ st.getProcessDistributionType(), FEMCondition::SOURCE_TERM)
+{
+ if (this->getProcessDistributionType() == FiniteElement::CONSTANT ||
+ this->getProcessDistributionType() == FiniteElement::CONSTANT_NEUMANN)
+ this->setConstantDisValue(st.getGeoNodeValue());
+ else if (this->getProcessDistributionType() == FiniteElement::LINEAR ||
+ this->getProcessDistributionType() == FiniteElement::LINEAR_NEUMANN)
+ {
+ const std::vector<int> st_nodes(st.getPointsWithDistribedST());
+ std::vector<size_t> dis_nodes(st_nodes.size());
+ for (size_t i=0; i<dis_nodes.size(); i++) dis_nodes[i] = static_cast<size_t>(st_nodes[i]);
+ this->setDisValues(dis_nodes, st.getDistribedST());
+ }
+ else if (this->getProcessDistributionType() == FiniteElement::DIRECT)
+ {
+ //this->_direct_file_name = st.fname;
+ }
+ else
+ std::cout << "Error in SourceTerm() - Unknown Process Distribution Type \"" <<
+ FiniteElement::convertDisTypeToString(st.getProcessDistributionType()) <<
+ "\"..." <<
+ std::endl;
+}
+
+// Legacy function (only required for ascii st-files): reads values for 'direct' source terms
+void SourceTerm::getDirectNodeValues(const std::string &filename,
+ std::vector< std::pair<size_t, double> > &node_values)
+{
+ std::ifstream in(filename.c_str());
+ if (!in.is_open())
+ {
+ std::cout << "Error in getNodeValues() - Could not find file for direct node values..." << std::endl;
+ return;
+ }
+
+ std::stringstream str_in;
+ std::string line("");
+ size_t idx(0);
+ double val(0);
+
+ while ( getline(in, line) )
+ {
+ if (line.find("#STOP") != std::string::npos)
+ return;
+ str_in << line;
+ str_in >> idx >> val;
+ node_values.push_back(std::pair<size_t, double>(idx, val));
+ str_in.clear();
+ }
+}
diff --git a/FemLib/SourceTerm.h b/FemLib/SourceTerm.h
new file mode 100644
index 00000000000..e8fd382f1be
--- /dev/null
+++ b/FemLib/SourceTerm.h
@@ -0,0 +1,37 @@
+/**
+ * \file SourceTerm.h
+ * 2011/08/30 KR inital implementation
+ *
+ */
+
+#ifndef SOURCETERM_H
+#define SOURCETERM_H
+
+#include "FEMCondition.h"
+
+/**
+ * \brief Adapter class for handling boundary conditions in the user Interface
+ * \sa FEMCondition
+ */
+class SourceTerm : public FEMCondition
+{
+public:
+ SourceTerm(const std::string &geometry_name)
+ : FEMCondition(geometry_name, FEMCondition::SOURCE_TERM), _tim_type(0) {}
+ SourceTerm(const CSourceTerm &st, const std::string &geometry_name);
+ SourceTerm(const FEMCondition &cond)
+ : FEMCondition(cond, FEMCondition::SOURCE_TERM) {};
+ ~SourceTerm() {}
+
+ size_t getTimType() const {return _tim_type; }
+ void setTimType(size_t value) { _tim_type = value; }
+
+ // Legacy function (only required for ascii st-files): reads values for 'direct' source terms
+ static void getDirectNodeValues(const std::string &filename,
+ std::vector< std::pair<size_t, double> > &nodes_values);
+
+private:
+ size_t _tim_type;
+};
+
+#endif //SOURCETERM_H
diff --git a/FileIO/GMSInterface.cpp b/FileIO/GMSInterface.cpp
new file mode 100644
index 00000000000..405358fbadd
--- /dev/null
+++ b/FileIO/GMSInterface.cpp
@@ -0,0 +1,314 @@
+/**
+ * \file GMSInterface.cpp
+ * 08/06/2010 KR Initial implementation
+ *
+ */
+
+#include "GMSInterface.h"
+#include "Station.h"
+#include "Mesh.h"
+#include "Node.h"
+#include "Elements/Tet.h"
+#include "Elements/Pyramid.h"
+#include "Elements/Prism.h"
+#include "MSHEnums.h"
+#include "StringTools.h"
+#include <fstream>
+
+int GMSInterface::readBoreholesFromGMS(std::vector<GeoLib::Point*>* boreholes,
+ const std::string &filename)
+{
+ double depth(-9999.0);
+ std::string line(""), cName(""), sName("");
+ std::list<std::string>::const_iterator it;
+ GeoLib::Point* pnt = new GeoLib::Point();
+ GeoLib::StationBorehole* newBorehole = NULL;
+ std::ifstream in( filename.c_str() );
+
+ if (!in.is_open())
+ {
+ std::cout << "GMSInterface::readBoreholeFromGMS() - Could not open file...\n";
+ return 0;
+ }
+
+ /* skipping first line because it contains field names */
+ getline(in, line);
+
+ /* read all stations */
+ while ( getline(in, line) )
+ {
+ std::list<std::string> fields = splitString(line, '\t');
+
+ if (fields.size() >= 5)
+ {
+ if (fields.begin()->compare(cName) == 0) // add new layer
+ {
+ it = fields.begin();
+ (*pnt)[0] = strtod((++it)->c_str(), 0);
+ (*pnt)[1] = strtod((++it)->c_str(), 0);
+ (*pnt)[2] = strtod((++it)->c_str(), 0);
+
+ // check if current layer has a thickness of 0.0.
+ // if so skip it since it will mess with the vtk-visualisation later on!
+ if ((*pnt)[2] != depth)
+ {
+ newBorehole->addSoilLayer((*pnt)[0], (*pnt)[1], (*pnt)[2], sName);
+ sName = (*(++it));
+ depth = (*pnt)[2];
+ }
+ else
+ std::cout << "Warning: Skipped layer \"" << sName << "\" in borehole \""
+ << cName << "\" because of thickness 0.0." << std::endl;
+ }
+ else // add new borehole
+ {
+ if (newBorehole != NULL)
+ {
+ newBorehole->setDepth((*newBorehole)[2] - depth);
+ boreholes->push_back(newBorehole);
+ }
+ cName = *fields.begin();
+ it = fields.begin();
+ (*pnt)[0] = strtod((++it)->c_str(), 0);
+ (*pnt)[1] = strtod((++it)->c_str(), 0);
+ (*pnt)[2] = strtod((++it)->c_str(), 0);
+ sName = (*(++it));
+ newBorehole = GeoLib::StationBorehole::createStation(cName, (*pnt)[0], (*pnt)[1], (*pnt)[2], 0);
+ depth = (*pnt)[2];
+ }
+ }
+ else
+ std::cout <<
+ "GMSInterface::readBoreholeFromGMS() - Error reading format..." <<
+ std::endl;
+ }
+ // write the last borehole from the file
+ if (newBorehole != NULL)
+ {
+ newBorehole->setDepth((*newBorehole)[2] - depth);
+ boreholes->push_back(newBorehole);
+ }
+ in.close();
+
+ if (boreholes->empty())
+ return 0;
+ return 1;
+}
+
+/*
+ // all boreholes to GMS which each borehole in a single file
+ void StationIO::writeBoreholesToGMS(const std::vector<GEOLIB::Point*> *stations)
+ {
+ //std::vector<std::string> soilID(1);
+ std::vector<std::string> soilID = readSoilIDfromFile("d:/BodeTimeline.txt");
+ for (size_t i=0; i<stations->size(); i++)
+ StationIO::writeBoreholeToGMS(static_cast<GEOLIB::StationBorehole*>((*stations)[i]), std::string("Borehole-" + static_cast<GEOLIB::StationBorehole*>((*stations)[i])->getName() + ".txt"), soilID);
+ StationIO::writeSoilIDTable(soilID, "SoilIDReference.txt");
+ }
+ */
+void GMSInterface::writeBoreholesToGMS(const std::vector<GeoLib::Point*>* stations,
+ const std::string &filename)
+{
+ std::ofstream out( filename.c_str(), std::ios::out );
+ size_t idx = 0;
+ std::vector<std::string> soilID = readSoilIDfromFile("d:/BodeTimeline.txt");
+
+ // write header
+ out << "name" << "\t" << std::fixed << "X" << "\t" << "Y" << "\t" << "Z" << "\t" <<
+ "soilID" << std::endl;
+
+ for (size_t j = 0; j < stations->size(); j++)
+ {
+ GeoLib::StationBorehole* station =
+ static_cast<GeoLib::StationBorehole*>((*stations)[j]);
+ std::vector<GeoLib::Point*> profile = station->getProfile();
+ std::vector<std::string> soilNames = station->getSoilNames();
+
+ size_t nLayers = profile.size();
+ for (size_t i = 1; i < nLayers; i++)
+ {
+ if ( (i > 1) && (soilNames[i].compare(soilNames[i - 1]) == 0) )
+ continue;
+ idx = getSoilID(soilID, soilNames[i]);
+
+ out << station->getName() << "\t" << std::fixed <<
+ (*(profile[i - 1]))[0] << "\t"
+ << (*(profile[i - 1]))[1] << "\t" << (*(profile[i - 1]))[2] << "\t"
+ << idx << std::endl;
+ }
+ out << station->getName() << "\t" << std::fixed <<
+ (*(profile[nLayers - 1]))[0] << "\t"
+ << (*(profile[nLayers -
+ 1]))[1] << "\t" << (*(profile[nLayers - 1]))[2] << "\t"
+ << idx << std::endl; // this line marks the end of the borehole
+ }
+
+ out.close();
+ GMSInterface::writeSoilIDTable(soilID, "d:/SoilIDReference.txt");
+}
+
+int GMSInterface::writeBoreholeToGMS(const GeoLib::StationBorehole* station,
+ const std::string &filename,
+ std::vector<std::string> &soilID)
+{
+ std::ofstream out( filename.c_str(), std::ios::out );
+ size_t idx = 0;
+
+ // write header
+ out << "name" << "\t" << std::fixed << "X" << "\t" << "Y" << "\t" << "Z" << "\t" <<
+ "soilID" << std::endl;
+
+ std::vector<GeoLib::Point*> profile = station->getProfile();
+ std::vector<std::string> soilNames = station->getSoilNames();
+
+ // write table
+ size_t nLayers = profile.size();
+ for (size_t i = 1; i < nLayers; i++)
+ {
+ if ( (i > 1) && (soilNames[i].compare(soilNames[i - 1]) == 0) )
+ continue;
+ idx = getSoilID(soilID, soilNames[i]);
+
+ out << station->getName() << "\t" << std::fixed << (*(profile[i - 1]))[0] << "\t"
+ << (*(profile[i - 1]))[1] << "\t" << (*(profile[i - 1]))[2] << "\t"
+ << idx << std::endl;
+ }
+ out << station->getName() << "\t" << std::fixed << (*(profile[nLayers - 1]))[0] << "\t"
+ << (*(profile[nLayers - 1]))[1] << "\t" << (*(profile[nLayers - 1]))[2] << "\t"
+ << idx << std::endl; // this line marks the end of the borehole
+ out.close();
+
+ return 1;
+}
+
+size_t GMSInterface::getSoilID(std::vector<std::string> &soilID, std::string &soilName)
+{
+ for (size_t j = 0; j < soilID.size(); j++)
+ if (soilID[j].compare(soilName) == 0)
+ return j;
+ soilID.push_back(soilName);
+ return soilID.size() - 1;
+}
+
+int GMSInterface::writeSoilIDTable(const std::vector<std::string> &soilID,
+ const std::string &filename)
+{
+ std::ofstream out( filename.c_str(), std::ios::out );
+
+ // write header
+ out << "ID" << "\t" << std::fixed << "Soil name" << std::endl;
+
+ // write table
+ size_t nIDs = soilID.size();
+ for (size_t i = 0; i < nIDs; i++)
+ out << i << "\t" << std::fixed << soilID[i] << "\t" << std::endl;
+ out.close();
+
+ return 1;
+}
+
+std::vector<std::string> GMSInterface::readSoilIDfromFile(const std::string &filename)
+{
+ std::vector<std::string> soilID;
+ std::string line;
+
+ std::ifstream in( filename.c_str() );
+
+ if (in.is_open())
+ while ( getline(in, line) )
+ {
+ trim(line);
+ soilID.push_back(line);
+ }
+ in.close();
+
+ return soilID;
+}
+
+MeshLib::Mesh* GMSInterface::readGMS3DMMesh(std::string filename)
+{
+ std::string line("");
+
+ std::ifstream in(filename.c_str());
+ if (!in.is_open())
+ {
+ std::cout << "GMSInterface::readGMS3DMMesh() - Could not open file..." << std::endl;
+ return NULL;
+ }
+
+ // Read data from file
+ getline(in, line); // "MESH3D"
+ if (line.compare("MESH3D") != 0)
+ {
+ std::cout << "GMSInterface::readGMS3DMMesh() - Could not read expected file header..." << std::endl;
+ return NULL;
+ }
+
+ std::cout << "Read GMS-3DM data...";
+ std::vector<MeshLib::Node*> nodes;
+ std::vector<MeshLib::Element*> elements;
+
+ // elements are listed before nodes in 3dm-format, therefore
+ // traverse file twice and read first nodes and then elements
+ std::string d1,d2;
+ double x[3];
+ // read nodes
+ while ( getline(in, line) )
+ {
+ if (line[0] == 'N') // "ND" for Node
+ {
+ std::stringstream str(line);
+ str >> d1 >> d2 >> x[0] >> x[1] >> x[2];
+ MeshLib::Node* node = new MeshLib::Node(x);
+ nodes.push_back(node);
+ }
+ }
+
+ // NOTE: Element types E8H (Hex), E4Q (Quad), E3T (Tri) are not implemented yet
+ // read elements
+ in.seekg(0, std::ios::beg);
+ unsigned node_idx[6], mat_id;
+ while ( getline(in, line) )
+ {
+ MeshLib::Element* elem (NULL);
+ std::string element_id(line.substr(0,3));
+ std::stringstream str(line);
+
+ if (element_id.compare("E6W") == 0) // Prism
+ {
+ str >> d1 >> d2 >> node_idx[0] >> node_idx[1] >> node_idx[2] >> node_idx[3]
+ >> node_idx[4] >> node_idx[5] >> mat_id;
+ elem = new MeshLib::Prism(nodes[node_idx[0]-1], nodes[node_idx[1]-1], nodes[node_idx[2]-1],
+ nodes[node_idx[3]-1], nodes[node_idx[4]-1], nodes[node_idx[5]-1], mat_id);
+ elements.push_back(elem);
+ }
+ else if (element_id.compare("E4T") == 0) // Tet
+ {
+ str >> d1 >> d2 >> node_idx[0] >> node_idx[1] >> node_idx[2] >> node_idx[3] >> mat_id;
+ elem = new MeshLib::Tet(nodes[node_idx[0]-1], nodes[node_idx[1]-1], nodes[node_idx[2]-1], nodes[node_idx[3]-1], mat_id);
+ elements.push_back(elem);
+ }
+ else if ((element_id.compare("E4P") == 0) || (element_id.compare("E5P") == 0)) // Pyramid (both do exist for some reason)
+ {
+ str >> d1 >> d2 >> node_idx[0] >> node_idx[1] >> node_idx[2] >> node_idx[3] >> node_idx[4] >> mat_id;
+ elem = new MeshLib::Pyramid(nodes[node_idx[0]-1], nodes[node_idx[1]-1], nodes[node_idx[2]-1],
+ nodes[node_idx[3]-1], nodes[node_idx[4]-1], mat_id);
+ elements.push_back(elem);
+ }
+ else if (element_id.compare("ND ") == 0) // Node
+ {
+ continue; // skip because nodes have already been read
+ }
+ else //default
+ {
+ std::cout << "GMSInterface::readGMS3DMMesh() - Element type \"" << element_id << "\"not recognised ..." << std::endl;
+ return NULL;
+ }
+ }
+
+ in.close();
+ std::cout << "finished" << std::endl;
+
+ std::string mesh_name (BaseLib::getFileNameFromPath(filename));
+ return new MeshLib::Mesh(mesh_name, nodes, elements);
+}
diff --git a/FileIO/GMSInterface.h b/FileIO/GMSInterface.h
new file mode 100644
index 00000000000..c40fa921f0e
--- /dev/null
+++ b/FileIO/GMSInterface.h
@@ -0,0 +1,64 @@
+/**
+ * \file GMSInterface.h
+ * 08/06/2010 KR Initial implementation
+ *
+ */
+
+#ifndef GMSINTERFACE_H_
+#define GMSINTERFACE_H_
+
+#include <list>
+#include <vector>
+
+#include "Point.h"
+
+namespace GeoLib {
+ class Station;
+ class StationBorehole;
+}
+
+namespace MeshLib {
+ class Mesh;
+}
+
+/**
+ * \brief Manages the import and export of Aquaveo GMS files into and out of GEOLIB.
+ */
+class GMSInterface
+{
+public:
+ /// Exports borehole data from all boreholes in a list to a file in GMS-format. (Note: there are some hardcoded tmp-files in the method that you might need to change!)
+ static void writeBoreholesToGMS(const std::vector<GeoLib::Point*>* stations,
+ const std::string &filename);
+
+ /// Imports borehole data from a file in GMS-format.
+ static int readBoreholesFromGMS(std::vector<GeoLib::Point*>* boreholes,
+ const std::string &filename);
+
+ /// Exports borehole data from one borehole to a file in GMS-format.
+ static int writeBoreholeToGMS(const GeoLib::StationBorehole* station,
+ const std::string &filename,
+ std::vector<std::string> &soilID);
+
+ /// Writes a file that assigns each soilID-index in the GMS export file a name.
+ static int writeSoilIDTable(const std::vector<std::string> &soilID,
+ const std::string &filename);
+
+ /// Reads a GMS *.3dm file and converts it to an CFEMesh.
+ static MeshLib::Mesh* readGMS3DMMesh(std::string file_name);
+
+private:
+ /**
+ * \brief Reads SoilIDs for Borehole export from an external file
+ *
+ * The method expects a file with the name of one stratigraphic layer at each line. These layers are assigned
+ * ascending IDs, i.e. the first name gets index 0, the second line gets index 1, etc.
+ * \return An array with the names of the stratigraphic layers in which the index for each string equals its ID.
+ */
+ static std::vector<std::string> readSoilIDfromFile(const std::string &filename);
+
+ /// Finds the ID assigned to soilName or creates a new one ( this method is called from writeBoreholeToGMS() )
+ static size_t getSoilID(std::vector<std::string> &soilID, std::string &soilName);
+};
+
+#endif /* GMSINTERFACE_H_ */
diff --git a/FileIO/Gmsh2GeoIO.cpp b/FileIO/Gmsh2GeoIO.cpp
new file mode 100644
index 00000000000..26156d3cab3
--- /dev/null
+++ b/FileIO/Gmsh2GeoIO.cpp
@@ -0,0 +1,123 @@
+/*
+ * Gmsh2GeoIO.cpp
+ *
+ * Created on: Aug 18, 2011
+ * Author: TF
+ */
+
+#include <fstream>
+#include <vector>
+
+// Base
+#include "StringTools.h"
+
+#include "GEOObjects.h"
+#include "Gmsh2GeoIO.h"
+
+namespace FileIO
+{
+void Gmsh2GeoIO::loadMeshAsGeometry (std::string & fname, GeoLib::GEOObjects* geo)
+{
+ // open file
+ std::ifstream ins (fname.c_str());
+ if (!ins)
+ {
+ std::cout << "could not open file " << fname << std::endl;
+ return;
+ }
+
+ std::string line;
+ // read gmsh header
+ getline (ins, line); // $MeshFormat
+ getline (ins, line);
+ getline (ins, line); // $EndMeshFormat
+
+ // read nodes tag
+ getline (ins, line);
+ // read number of nodes
+ getline (ins, line);
+ const size_t n_pnts (str2number<size_t>(line));
+ std::vector<GeoLib::Point*>* pnts (new std::vector<GeoLib::Point*>);
+ for (size_t k(0); k < n_pnts; k++)
+ {
+ getline (ins, line);
+ // parse id
+ size_t pos_beg(0);
+ size_t pos_end (line.find(" "));
+ // the sub string line.substr(pos_beg, pos_end-pos_beg) represents the id
+ // parse x coordinate
+ pos_beg = pos_end + 1;
+ pos_end = line.find(" ", pos_beg);
+ double x (str2number<double>(line.substr(pos_beg, pos_end - pos_beg)));
+ // parse y coordinate
+ pos_beg = pos_end + 1;
+ pos_end = line.find(" ", pos_beg);
+ double y (str2number<double>(line.substr(pos_beg, pos_end - pos_beg)));
+ // parse z coordinate
+ pos_beg = pos_end + 1;
+ pos_end = line.find("\n", pos_beg);
+ double z (str2number<double>(line.substr(pos_beg, pos_end - pos_beg)));
+
+ pnts->push_back (new GeoLib::Point (x,y,z));
+ }
+ // read end nodes tag
+ getline (ins, line);
+
+ geo->addPointVec (pnts, fname);
+
+ std::vector<size_t> const& pnt_id_map (geo->getPointVecObj(fname)->getIDMap());
+ // read element tag
+ getline (ins, line);
+ // read number of elements
+ getline (ins, line);
+ const size_t n_elements (str2number<size_t>(line));
+ GeoLib::Surface* sfc (new GeoLib::Surface (*pnts));
+ for (size_t k(0); k < n_elements; k++)
+ {
+ getline (ins, line);
+ // parse id
+ size_t pos_beg(0);
+ size_t pos_end (line.find(" "));
+ // the sub string line.substr(pos_beg, pos_end-pos_beg) represents the id
+ // parse element type
+ pos_beg = pos_end + 1;
+ pos_end = line.find(" ", pos_beg);
+ size_t ele_type (str2number<size_t>(line.substr(pos_beg, pos_end - pos_beg)));
+ if (ele_type == 2) // read 3 node triangle
+ { // parse number of tags
+ pos_beg = pos_end + 1;
+ pos_end = line.find(" ", pos_beg);
+ const size_t n_tags (str2number<size_t>(line.substr(pos_beg,
+ pos_end - pos_beg)));
+ // (over) read tags
+ for (size_t j(0); j < n_tags; j++)
+ {
+ pos_beg = pos_end + 1;
+ pos_end = line.find(" ", pos_beg);
+ }
+ // parse first id of triangle
+ pos_beg = pos_end + 1;
+ pos_end = line.find(" ", pos_beg);
+ const size_t id0 (str2number<size_t>(line.substr(pos_beg,
+ pos_end - pos_beg)) - 1); // shift -1!
+ // parse second id of triangle
+ pos_beg = pos_end + 1;
+ pos_end = line.find(" ", pos_beg);
+ const size_t id1 (str2number<size_t>(line.substr(pos_beg,
+ pos_end - pos_beg)) - 1); // shift -1!
+ // parse third id of triangle
+ pos_beg = pos_end + 1;
+ pos_end = line.find(" ", pos_beg);
+ const size_t id2 (str2number<size_t>(line.substr(pos_beg,
+ pos_end - pos_beg)) - 1); // shift -1!
+ sfc->addTriangle (pnt_id_map[id0], pnt_id_map[id1], pnt_id_map[id2]);
+ }
+ }
+ // read end element tag
+ getline (ins, line);
+
+ std::vector<GeoLib::Surface*>* sfcs (new std::vector<GeoLib::Surface*>);
+ sfcs->push_back(sfc);
+ geo->addSurfaceVec (sfcs, fname);
+}
+} // end namespace FileIO
diff --git a/FileIO/Gmsh2GeoIO.h b/FileIO/Gmsh2GeoIO.h
new file mode 100644
index 00000000000..cf5163cc9c4
--- /dev/null
+++ b/FileIO/Gmsh2GeoIO.h
@@ -0,0 +1,33 @@
+/*
+ * Gmsh2GeoIO.h
+ *
+ * Created on: Aug 18, 2011
+ * Author: TF
+ */
+
+#ifndef GMSH2GEOIO_H_
+#define GMSH2GEOIO_H_
+
+#include <string>
+
+namespace GeoLib
+{
+class GEOObjects;
+}
+
+namespace FileIO
+{
+class Gmsh2GeoIO
+{
+public:
+ /**
+ * load a surface mesh (gmsh format) as a geometric surface
+ * @param fname file name of the surface mesh
+ * @param geo the object that manages all geometries,
+ * new surface will be put into this container
+ */
+ static void loadMeshAsGeometry (std::string & fname, GeoLib::GEOObjects* geo);
+};
+} // end namespace FileIO
+
+#endif /* GMSH2GEOIO_H_ */
diff --git a/FileIO/MeshIO/GMSHAdaptiveMeshDensity.cpp b/FileIO/MeshIO/GMSHAdaptiveMeshDensity.cpp
new file mode 100644
index 00000000000..af807164764
--- /dev/null
+++ b/FileIO/MeshIO/GMSHAdaptiveMeshDensity.cpp
@@ -0,0 +1,151 @@
+/*
+ * GMSHAdaptiveMeshDensity.cpp
+ *
+ * Created on: Mar 5, 2012
+ * Author: TF
+ */
+
+#include <list>
+
+// FileIO
+#include "MeshIO/GMSHAdaptiveMeshDensity.h"
+
+// GEOLIB
+#include "Polygon.h"
+
+namespace FileIO {
+
+GMSHAdaptiveMeshDensity::GMSHAdaptiveMeshDensity(double pnt_density, double station_density,
+ size_t max_pnts_per_leaf) :
+ _pnt_density(pnt_density), _station_density(station_density),
+ _max_pnts_per_leaf(max_pnts_per_leaf), _quad_tree(NULL)
+{
+}
+
+GMSHAdaptiveMeshDensity::~GMSHAdaptiveMeshDensity()
+{
+ delete _quad_tree;
+}
+
+void GMSHAdaptiveMeshDensity::init(std::vector<GEOLIB::Point const*> const& pnts)
+{
+ // *** QuadTree - determining bounding box
+#ifndef NDEBUG
+ std::cout << "[GMSHAdaptiveMeshDensity::init]" << std::endl;
+ std::cout << "\tcomputing axis aligned bounding box (2D) for quadtree ... " << std::flush;
+#endif
+ GEOLIB::Point min(pnts[0]->getData()), max(pnts[0]->getData());
+ size_t n_pnts(pnts.size());
+ for (size_t k(1); k<n_pnts; k++) {
+ for (size_t j(0); j<2; j++)
+ if ((*(pnts[k]))[j] < min[j]) min[j] = (*(pnts[k]))[j];
+ for (size_t j(0); j<2; j++)
+ if ((*(pnts[k]))[j] > max[j]) max[j] = (*(pnts[k]))[j];
+ }
+ min[2] = 0.0;
+ max[2] = 0.0;
+#ifndef NDEBUG
+ std::cout << "ok" << std::endl;
+#endif
+
+ // *** QuadTree - create object
+#ifndef NDEBUG
+ std::cout << "\tcreating quadtree ... " << std::flush;
+#endif
+ _quad_tree = new GEOLIB::QuadTree<GEOLIB::Point> (min, max, _max_pnts_per_leaf);
+#ifndef NDEBUG
+ std::cout << "ok" << std::endl;
+#endif
+
+ // *** QuadTree - insert points
+ addPoints(pnts);
+}
+
+void GMSHAdaptiveMeshDensity::addPoints(std::vector<GEOLIB::Point const*> const& pnts)
+{
+ // *** QuadTree - insert points
+ const size_t n_pnts(pnts.size());
+#ifndef NDEBUG
+ std::cout << "\tinserting " << n_pnts << " points into quadtree ... " <<
+ std::flush;
+#endif
+ for (size_t k(0); k < n_pnts; k++)
+ _quad_tree->addPoint(pnts[k]);
+#ifndef NDEBUG
+ std::cout << "ok" << std::endl;
+#endif
+ _quad_tree->balance();
+}
+
+double GMSHAdaptiveMeshDensity::getMeshDensityAtPoint(GEOLIB::Point const*const pnt) const
+{
+ GEOLIB::Point ll, ur;
+ _quad_tree->getLeaf(*pnt, ll, ur);
+ return _pnt_density * (ur[0] - ll[0]);
+}
+
+double GMSHAdaptiveMeshDensity::getMeshDensityAtStation(GEOLIB::Point const*const pnt) const
+{
+ GEOLIB::Point ll, ur;
+ _quad_tree->getLeaf(*pnt, ll, ur);
+ return (_station_density * (ur[0] - ll[0]));
+}
+
+void GMSHAdaptiveMeshDensity::getSteinerPoints (std::vector<GEOLIB::Point*> & pnts, size_t additional_levels) const
+{
+ // get Steiner points
+ size_t max_depth(0);
+ _quad_tree->getMaxDepth(max_depth);
+
+ std::list<GEOLIB::QuadTree<GEOLIB::Point>*> leaf_list;
+ _quad_tree->getLeafs(leaf_list);
+
+ for (std::list<GEOLIB::QuadTree<GEOLIB::Point>*>::const_iterator it(leaf_list.begin()); it
+ != leaf_list.end(); it++) {
+ if ((*it)->getPoints().empty()) {
+ // compute point from square
+ GEOLIB::Point ll, ur;
+ (*it)->getSquarePoints(ll, ur);
+ if ((*it)->getDepth() + additional_levels > max_depth) {
+ additional_levels = max_depth - (*it)->getDepth();
+ }
+ const size_t n_pnts_per_quad_dim (MathLib::fastpow(2, additional_levels));
+ const double delta ((ur[0] - ll[0]) / (2 * n_pnts_per_quad_dim));
+ for (size_t i(0); i<n_pnts_per_quad_dim; i++) {
+ for (size_t j(0); j<n_pnts_per_quad_dim; j++) {
+ pnts.push_back(new GEOLIB::Point (ll[0] + (2*i+1) * delta, ll[1] + (2*j+1) * delta, 0.0));
+ }
+ }
+
+ }
+ }
+}
+
+#ifndef NDEBUG
+void GMSHAdaptiveMeshDensity::getQuadTreeGeometry(std::vector<GEOLIB::Point*> &pnts,
+ std::vector<GEOLIB::Polyline*> &plys) const
+{
+ std::list<GEOLIB::QuadTree<GEOLIB::Point>*> leaf_list;
+ _quad_tree->getLeafs(leaf_list);
+
+ for (std::list<GEOLIB::QuadTree<GEOLIB::Point>*>::const_iterator it(leaf_list.begin()); it
+ != leaf_list.end(); it++) {
+ // fetch corner points from leaf
+ GEOLIB::Point *ll(new GEOLIB::Point), *ur(new GEOLIB::Point);
+ (*it)->getSquarePoints(*ll, *ur);
+ size_t pnt_offset (pnts.size());
+ pnts.push_back(ll);
+ pnts.push_back(new GEOLIB::Point((*ur)[0], (*ll)[1], 0.0));
+ pnts.push_back(ur);
+ pnts.push_back(new GEOLIB::Point((*ll)[0], (*ur)[1], 0.0));
+ plys.push_back(new GEOLIB::Polyline(pnts));
+ plys[plys.size()-1]->addPoint(pnt_offset);
+ plys[plys.size()-1]->addPoint(pnt_offset+1);
+ plys[plys.size()-1]->addPoint(pnt_offset+2);
+ plys[plys.size()-1]->addPoint(pnt_offset+3);
+ plys[plys.size()-1]->addPoint(pnt_offset);
+ }
+}
+#endif
+
+} // end namespace FileIO
diff --git a/FileIO/MeshIO/GMSHAdaptiveMeshDensity.h b/FileIO/MeshIO/GMSHAdaptiveMeshDensity.h
new file mode 100644
index 00000000000..0510cbef528
--- /dev/null
+++ b/FileIO/MeshIO/GMSHAdaptiveMeshDensity.h
@@ -0,0 +1,49 @@
+/*
+ * GMSHAdaptiveMeshDensity.h
+ *
+ * Created on: Mar 5, 2012
+ * Author: TF
+ */
+
+#ifndef GMSHADAPTIVEMESHDENSITY_H_
+#define GMSHADAPTIVEMESHDENSITY_H_
+
+// FileIO
+#include "GMSHMeshDensityStrategy.h"
+
+// GEOLIB
+#include "Point.h"
+#include "QuadTree.h"
+#ifndef NDEBUG
+#include "Polyline.h"
+#endif
+
+namespace GEOLIB {
+class Polygon;
+}
+
+namespace FileIO {
+
+class GMSHAdaptiveMeshDensity: public GMSHMeshDensityStrategy {
+public:
+ GMSHAdaptiveMeshDensity(double pnt_density, double station_density, size_t max_pnts_per_leaf);
+ virtual ~GMSHAdaptiveMeshDensity();
+ void init(std::vector<GEOLIB::Point const*> const& pnts);
+ double getMeshDensityAtPoint(GEOLIB::Point const*const pnt) const;
+ void addPoints(std::vector<GEOLIB::Point const*> const& pnts);
+ double getMeshDensityAtStation(GEOLIB::Point const*const) const;
+ void getSteinerPoints (std::vector<GEOLIB::Point*> & pnts, size_t additional_levels = 0) const;
+#ifndef NDEBUG
+ void getQuadTreeGeometry(std::vector<GEOLIB::Point*> &pnts, std::vector<GEOLIB::Polyline*> &plys) const;
+#endif
+
+private:
+ double _pnt_density;
+ double _station_density;
+ size_t _max_pnts_per_leaf;
+ GEOLIB::QuadTree<GEOLIB::Point> *_quad_tree;
+};
+
+} // end namespace FileIO
+
+#endif /* GMSHADAPTIVEMESHDENSITY_H_ */
diff --git a/FileIO/MeshIO/GMSHFixedMeshDensity.cpp b/FileIO/MeshIO/GMSHFixedMeshDensity.cpp
new file mode 100644
index 00000000000..4a7e106da53
--- /dev/null
+++ b/FileIO/MeshIO/GMSHFixedMeshDensity.cpp
@@ -0,0 +1,30 @@
+/*
+ * GMSHFixedMeshDensity.cpp
+ *
+ * Created on: Mar 5, 2012
+ * Author: TF
+ */
+
+#include "MeshIO/GMSHFixedMeshDensity.h"
+
+namespace FileIO {
+
+GMSHFixedMeshDensity::GMSHFixedMeshDensity(double mesh_density) :
+ _mesh_density(mesh_density)
+{
+}
+
+void GMSHFixedMeshDensity::init(std::vector<GEOLIB::Point const*> const& vec)
+{
+ // to avoid a warning here:
+ (void)(vec);
+}
+
+double GMSHFixedMeshDensity::getMeshDensityAtPoint(GEOLIB::Point const*const pnt) const
+{
+ // to avoid a warning here:
+ (void)(const_cast<GEOLIB::Point const*>(pnt));
+ return _mesh_density;
+}
+
+} // end namespace FileIO
diff --git a/FileIO/MeshIO/GMSHFixedMeshDensity.h b/FileIO/MeshIO/GMSHFixedMeshDensity.h
new file mode 100644
index 00000000000..60240b2f7ab
--- /dev/null
+++ b/FileIO/MeshIO/GMSHFixedMeshDensity.h
@@ -0,0 +1,28 @@
+/*
+ * GMSHFixedMeshDensity.h
+ *
+ * Created on: Mar 5, 2012
+ * Author: TF
+ */
+
+#ifndef GMSHFIXEDMESHDENSITY_H_
+#define GMSHFIXEDMESHDENSITY_H_
+
+#include "GMSHMeshDensityStrategy.h"
+
+namespace FileIO {
+
+class GMSHFixedMeshDensity : public GMSHMeshDensityStrategy
+{
+public:
+ GMSHFixedMeshDensity(double mesh_density);
+ void init(std::vector<GEOLIB::Point const*> const& vec);
+ double getMeshDensityAtPoint(GEOLIB::Point const*const) const;
+
+private:
+ double _mesh_density;
+};
+
+}
+
+#endif /* GMSHFIXEDMESHDENSITY_H_ */
diff --git a/FileIO/MeshIO/GMSHInterface.cpp b/FileIO/MeshIO/GMSHInterface.cpp
new file mode 100644
index 00000000000..deb1a8991d7
--- /dev/null
+++ b/FileIO/MeshIO/GMSHInterface.cpp
@@ -0,0 +1,272 @@
+/*
+ * GMSHInterface.cpp
+ *
+ * Created on: Apr 29, 2010
+ * Author: TF
+ */
+
+#include <fstream>
+#include <vector>
+
+// Base
+#include "swap.h"
+#include "Configure.h"
+#include "BuildInfo.h"
+
+// FileIO
+#include "GMSHInterface.h"
+#include "GMSHNoMeshDensity.h"
+#include "GMSHFixedMeshDensity.h"
+#include "GMSHAdaptiveMeshDensity.h"
+
+// GEOLIB
+#include "Point.h"
+#include "Polygon.h"
+#include "Polyline.h"
+#include "QuadTree.h"
+#include "PolylineWithSegmentMarker.h"
+
+// MSH
+#include "msh_elem.h"
+#include "msh_mesh.h"
+
+namespace FileIO {
+GMSHInterface::GMSHInterface(GEOLIB::GEOObjects & geo_objs, bool include_stations_as_constraints,
+ GMSH::MeshDensityAlgorithm mesh_density_algorithm, double param1, double param2,
+ size_t param3, std::vector<std::string>& selected_geometries) :
+ _n_lines(0), _n_plane_sfc(0), _geo_objs(geo_objs), _selected_geometries(selected_geometries),
+ _include_stations_as_constraints(include_stations_as_constraints)
+{
+ switch (mesh_density_algorithm) {
+ case GMSH::NoMeshDensity:
+ _mesh_density_strategy = new GMSHNoMeshDensity;
+ break;
+ case GMSH::FixedMeshDensity:
+ _mesh_density_strategy = new GMSHFixedMeshDensity(param1);
+ break;
+ case GMSH::AdaptiveMeshDensity:
+ _mesh_density_strategy = new GMSHAdaptiveMeshDensity(param1, param2, param3);
+ break;
+ }
+}
+
+bool GMSHInterface::isGMSHMeshFile(const std::string& fname)
+{
+ std::ifstream input(fname.c_str());
+
+ if (!input) {
+ std::cerr << "GMSHInterface::isGMSHMeshFile could not open file " << fname << std::endl;
+ return false;
+ }
+
+ std::string header_first_line;
+ input >> header_first_line;
+ if (header_first_line.find("$MeshFormat") != std::string::npos) {
+ // read version
+ std::string version;
+ getline(input, version);
+ getline(input, version);
+ std::cerr << "found GMSH mesh file version: " << version << std::endl;
+ input.close();
+ return true;
+ }
+
+ return false;
+}
+
+void GMSHInterface::readGMSHMesh(std::string const& fname, MeshLib::CFEMesh* mesh)
+{
+ std::string line;
+ std::ifstream in(fname.c_str(), std::ios::in);
+ getline(in, line); // Node keyword
+
+ if (line.find("$MeshFormat") != std::string::npos) {
+ getline(in, line); // version-number file-type data-size
+ getline(in, line); //$EndMeshFormat
+ getline(in, line); //$Nodes Keywords
+
+ size_t n_nodes(0);
+ size_t n_elements(0);
+ while (line.find("$EndElements") == std::string::npos) {
+ // Node data
+ long id;
+ double x, y, z;
+ in >> n_nodes >> std::ws;
+ for (size_t i = 0; i < n_nodes; i++) {
+ in >> id >> x >> y >> z >> std::ws;
+ mesh->nod_vector.push_back(new MeshLib::CNode(id, x, y, z));
+ }
+ getline(in, line); // End Node keyword $EndNodes
+
+ // Element data
+ getline(in, line); // Element keyword $Elements
+ in >> n_elements >> std::ws; // number-of-elements
+ for (size_t i = 0; i < n_elements; i++) {
+ MeshLib::CElem* elem(new MeshLib::CElem(i));
+ elem->Read(in, 7);
+ if (elem->GetElementType() != MshElemType::INVALID) mesh->ele_vector.push_back(elem);
+ }
+ getline(in, line); // END keyword
+
+ // correct indices TF
+ const size_t n_elements(mesh->ele_vector.size());
+ for (size_t k(0); k < n_elements; k++)
+ mesh->ele_vector[k]->SetIndex(k);
+
+ // ordering nodes and closing gaps TK
+ std::vector<size_t> gmsh_id;
+ size_t counter(0);
+ for (size_t i = 0; i < mesh->nod_vector.size(); i++) {
+ const size_t diff = mesh->nod_vector[i]->GetIndex() - counter;
+ if (diff == 0) {
+ gmsh_id.push_back(i);
+ counter++;
+ } else {
+ for (size_t j = 0; j < diff; j++) {
+ gmsh_id.push_back(i);
+ counter++;
+ }
+ i--;
+ }
+ }
+
+ for (size_t i = 0; i < mesh->ele_vector.size(); i++)
+ for (long j = 0; j < mesh->ele_vector[i]->GetVertexNumber(); j++)
+ mesh->ele_vector[i]->getNodeIndices()[j]
+ = gmsh_id[mesh->ele_vector[i]->GetNodeIndex(j) + 1];
+
+ for (size_t i = 0; i < mesh->nod_vector.size(); i++)
+ mesh->nod_vector[i]->SetIndex(i);
+ // END OF: ordering nodes and closing gaps TK
+ } /*End while*/
+ }
+ in.close();
+}
+
+int GMSHInterface::write(std::ostream& out)
+{
+ out << "// GMSH input file created by OpenGeoSys " << OGS_VERSION << " built on ";
+#ifdef BUILD_TIMESTAMP
+ out << BUILD_TIMESTAMP;
+#endif
+ out << std::endl << std::endl;
+
+ writeGMSHInputFile(out);
+ return 1;
+}
+
+void GMSHInterface::writeGMSHInputFile(std::ostream& out)
+{
+#ifndef NDEBUG
+ std::cerr << "[GMSHInterface::writeGMSHInputFile] get data from GEOObjects ... " << std::flush;
+#endif
+ // *** get and merge data from _geo_objs
+ _gmsh_geo_name = "GMSHGeometry";
+ _geo_objs.mergeGeometries(_selected_geometries, _gmsh_geo_name);
+ std::vector<GEOLIB::Point*> * merged_pnts(const_cast<std::vector<GEOLIB::Point*> *>(_geo_objs.getPointVec(_gmsh_geo_name)));
+ if (! merged_pnts) {
+ std::cerr << "[GMSHInterface::writeGMSHInputFile] did not found any points" << std::endl;
+ return;
+ } else {
+ const size_t n_pnts(merged_pnts->size());
+ for (size_t k(0); k<n_pnts; k++) {
+ (*((*merged_pnts)[k]))[2] = 0.0;
+ }
+ }
+ std::vector<GEOLIB::Polyline*> const* merged_plys(_geo_objs.getPolylineVec(_gmsh_geo_name));
+#ifndef NDEBUG
+ std::cerr << "ok" << std::endl;
+#endif
+
+ // *** compute topological hierarchy of polygons
+ if (merged_plys) {
+ for (std::vector<GEOLIB::Polyline*>::const_iterator it(merged_plys->begin());
+ it!=merged_plys->end(); it++) {
+ if ((*it)->isClosed()) {
+ _polygon_tree_list.push_back(new GMSHPolygonTree(new GEOLIB::Polygon(*(*it), true), NULL, _geo_objs, _gmsh_geo_name, _mesh_density_strategy));
+ }
+ }
+ std::cout << "[GMSHInterface::writeGMSHInputFile] compute topological hierarchy - detected "
+ << _polygon_tree_list.size() << " polygons" << std::endl;
+ GEOLIB::createPolygonTrees<FileIO::GMSHPolygonTree>(_polygon_tree_list);
+ std::cout << "[GMSHInterface::writeGMSHInputFile] compute topological hierarchy - calculated "
+ << _polygon_tree_list.size() << " polygon trees" << std::endl;
+ } else {
+ return;
+ }
+
+ // *** insert stations and polylines (except polygons) in the appropriate object of
+ // class GMSHPolygonTree
+ // *** insert stations
+ const size_t n_geo_names(_selected_geometries.size());
+ for (size_t j(0); j < n_geo_names; j++) {
+ const std::vector<GEOLIB::Point*>* stations (_geo_objs.getStationVec(_selected_geometries[j]));
+ if (stations) {
+ const size_t n_stations(stations->size());
+ for (size_t k(0); k < n_stations; k++) {
+ bool found(false);
+ for (std::list<GMSHPolygonTree*>::iterator it(_polygon_tree_list.begin());
+ it != _polygon_tree_list.end() && !found; it++) {
+ if ((*it)->insertStation((*stations)[k])) {
+ found = true;
+ }
+ }
+ }
+ }
+ }
+ // *** insert polylines
+ const size_t n_plys(merged_plys->size());
+ for (size_t k(0); k<n_plys; k++) {
+ if (! (*merged_plys)[k]->isClosed()) {
+ for (std::list<GMSHPolygonTree*>::iterator it(_polygon_tree_list.begin());
+ it != _polygon_tree_list.end(); it++) {
+ (*it)->insertPolyline(new GEOLIB::PolylineWithSegmentMarker(*(*merged_plys)[k]));
+ }
+ }
+ }
+
+ // *** init mesh density strategies
+ for (std::list<GMSHPolygonTree*>::iterator it(_polygon_tree_list.begin());
+ it != _polygon_tree_list.end(); it++) {
+ (*it)->initMeshDensityStrategy();
+ }
+
+ // *** create GMSH data structures
+ const size_t n_merged_pnts(merged_pnts->size());
+ _gmsh_pnts.resize(n_merged_pnts);
+ for (size_t k(0); k<n_merged_pnts; k++) {
+ _gmsh_pnts[k] = NULL;
+ }
+ for (std::list<GMSHPolygonTree*>::iterator it(_polygon_tree_list.begin());
+ it != _polygon_tree_list.end(); it++) {
+ (*it)->createGMSHPoints(_gmsh_pnts);
+ }
+
+ // *** finally write data :-)
+ writePoints(out);
+ size_t pnt_id_offset(_gmsh_pnts.size());
+ for (std::list<GMSHPolygonTree*>::iterator it(_polygon_tree_list.begin());
+ it != _polygon_tree_list.end(); it++) {
+ (*it)->writeLineLoop(_n_lines, _n_plane_sfc, out);
+ (*it)->writeSubPolygonsAsLineConstraints(_n_lines, _n_plane_sfc-1, out);
+ (*it)->writeLineConstraints(_n_lines, _n_plane_sfc-1, out);
+ (*it)->writeStations(pnt_id_offset, _n_plane_sfc-1, out);
+ (*it)->writeAdditionalPointData(pnt_id_offset, _n_plane_sfc-1, out);
+ }
+
+ _geo_objs.removeSurfaceVec(_gmsh_geo_name);
+ _geo_objs.removePolylineVec(_gmsh_geo_name);
+ _geo_objs.removePointVec(_gmsh_geo_name);
+}
+
+void GMSHInterface::writePoints(std::ostream& out) const
+{
+ const size_t n_gmsh_pnts(_gmsh_pnts.size());
+ for (size_t k(0); k<n_gmsh_pnts; k++) {
+ if (_gmsh_pnts[k]) {
+ out << *(_gmsh_pnts[k]) << std::endl;
+ }
+ }
+}
+
+} // end namespace FileIO
diff --git a/FileIO/MeshIO/GMSHInterface.h b/FileIO/MeshIO/GMSHInterface.h
new file mode 100644
index 00000000000..97d68f34ea7
--- /dev/null
+++ b/FileIO/MeshIO/GMSHInterface.h
@@ -0,0 +1,110 @@
+/*
+ * GMSHInterface.h
+ *
+ * Created on: Apr 29, 2010
+ * Author: TF
+ */
+
+#ifndef GMSHINTERFACE_H_
+#define GMSHINTERFACE_H_
+
+#include <string>
+#include <list>
+
+// FileIO
+#include "Writer.h"
+#include "GMSHPoint.h"
+#include "GMSHPolygonTree.h"
+#include "GMSHMeshDensityStrategy.h"
+
+// GEOLIB
+#include "GEOObjects.h"
+#include "Polygon.h"
+
+namespace MeshLib
+{
+class CFEMesh;
+}
+
+namespace FileIO
+{
+
+namespace GMSH {
+
+enum MeshDensityAlgorithm {
+ NoMeshDensity = 0, //!< do not set the parameter
+ FixedMeshDensity, //!< set the parameter with a fixed value
+ AdaptiveMeshDensity //!< computing the mesh density employing a QuadTree
+};
+
+}
+
+/**
+ * \brief Reads and writes GMSH-files to and from OGS data structures.
+ */
+class GMSHInterface : public Writer
+{
+public:
+
+ /**
+ *
+ * @param geo_objs reference tp instance of class GEOObject that maintains the geometries.
+ * The instance is used for preparation geometries for writing them to the gmsh file format.
+ * @param include_stations_as_constraints switch to enable writing stations as constraints
+ * @param mesh_density_algorithm one of the mesh density algorithms (\@see enum MeshDensityAlgorithm)
+ * @param param1 parameter that can be used for the mesh density algorithm
+ * @param param2 parameter that can be used for the mesh density algorithm
+ * @param param3 parameter that can be used for the mesh density algorithm
+ * @param selected_geometries vector of names of geometries, that should be employed for mesh generation.
+ * @return
+ */
+ GMSHInterface (GEOLIB::GEOObjects & geo_objs,
+ bool include_stations_as_constraints,
+ GMSH::MeshDensityAlgorithm mesh_density_algorithm,
+ double param1, double param2, size_t param3,
+ std::vector<std::string> & selected_geometries);
+
+ /**
+ * checks if there is a GMSH mesh file header
+ * @param fname the file name of the mesh (including the path)
+ * @return true, if the file seems to be a valid GMSH file, else false
+ */
+ static bool isGMSHMeshFile (const std::string& fname);
+ /**
+ * reads a mesh created by GMSH - this implementation is based on the former function GMSH2MSH
+ * @param fname the file name of the mesh (including the path)
+ * @param mesh the new mesh
+ * @return
+ */
+ static void readGMSHMesh (std::string const& fname, MeshLib::CFEMesh* mesh);
+
+protected:
+ int write(std::ostream& stream);
+
+private:
+ /**
+ * 1. get and merge data from _geo_objs
+ * 2. compute topological hierarchy
+ * @param out
+ */
+ void writeGMSHInputFile(std::ostream & out);
+
+ void writePoints(std::ostream& out) const;
+
+ size_t _n_lines;
+ size_t _n_plane_sfc;
+
+ GEOLIB::GEOObjects & _geo_objs;
+ std::vector<std::string>& _selected_geometries;
+ std::string _gmsh_geo_name;
+ std::list<GMSHPolygonTree*> _polygon_tree_list;
+
+ bool _include_stations_as_constraints;
+
+ std::vector<FileIO::GMSHPoint*> _gmsh_pnts;
+
+ GMSHMeshDensityStrategy *_mesh_density_strategy;
+};
+}
+
+#endif /* GMSHINTERFACE_H_ */
diff --git a/FileIO/MeshIO/GMSHLine.cpp b/FileIO/MeshIO/GMSHLine.cpp
new file mode 100644
index 00000000000..fa6cfdcbe32
--- /dev/null
+++ b/FileIO/MeshIO/GMSHLine.cpp
@@ -0,0 +1,30 @@
+/*
+ * GMSHLine.cpp
+ *
+ * Created on: Mar 22, 2012
+ * Author: fischeth
+ */
+
+#include <MeshIO/GMSHLine.h>
+
+namespace FileIO {
+
+GMSHLine::GMSHLine(size_t start_point_id, size_t end_point_id) :
+ _start_pnt_id(start_point_id), _end_pnt_id(end_point_id)
+{}
+
+GMSHLine::~GMSHLine()
+{}
+
+void GMSHLine::write(std::ostream &os, size_t id) const
+{
+ os << "Line(" << id << ") = {" << _start_pnt_id << "," << _end_pnt_id << "};" << std::endl;
+}
+
+void GMSHLine::resetLineData(size_t start_point_id, size_t end_point_id)
+{
+ _start_pnt_id = start_point_id;
+ _end_pnt_id = end_point_id;
+}
+
+} // end namespace FileIO
diff --git a/FileIO/MeshIO/GMSHLine.h b/FileIO/MeshIO/GMSHLine.h
new file mode 100644
index 00000000000..99096aa4d75
--- /dev/null
+++ b/FileIO/MeshIO/GMSHLine.h
@@ -0,0 +1,29 @@
+/*
+ * GMSHLine.h
+ *
+ * Created on: Mar 22, 2012
+ * Author: fischeth
+ */
+
+#ifndef GMSHLINE_H_
+#define GMSHLINE_H_
+
+#include <ostream>
+
+namespace FileIO {
+
+class GMSHLine {
+public:
+ GMSHLine(size_t start_point_id, size_t end_point_id);
+ virtual ~GMSHLine();
+ void write(std::ostream &os, size_t id) const;
+ void resetLineData(size_t start_point_id, size_t end_point_id);
+
+private:
+ size_t _start_pnt_id;
+ size_t _end_pnt_id;
+};
+
+}
+
+#endif /* GMSHLINE_H_ */
diff --git a/FileIO/MeshIO/GMSHLineLoop.cpp b/FileIO/MeshIO/GMSHLineLoop.cpp
new file mode 100644
index 00000000000..eb01789f42a
--- /dev/null
+++ b/FileIO/MeshIO/GMSHLineLoop.cpp
@@ -0,0 +1,47 @@
+/*
+ * GMSHLineLoop.cpp
+ *
+ * Created on: Mar 22, 2012
+ * Author: fischeth
+ */
+
+#include "MeshIO/GMSHLineLoop.h"
+
+namespace FileIO {
+
+GMSHLineLoop::GMSHLineLoop(bool is_sfc) :
+ _is_sfc(is_sfc)
+{}
+
+GMSHLineLoop::~GMSHLineLoop()
+{
+ const size_t n_lines (_lines.size());
+ for (size_t k(0); k<n_lines; k++) {
+ delete _lines[k];
+ }
+}
+
+void GMSHLineLoop::addLine(GMSHLine* line)
+{
+ _lines.push_back(line);
+}
+
+void GMSHLineLoop::write(std::ostream &os, size_t line_offset, size_t sfc_offset) const
+{
+ const size_t n_lines (_lines.size());
+ for (size_t k(0); k<n_lines; k++) {
+ (_lines[k])->write(os, line_offset+k);
+ }
+ os << "Line Loop(" << line_offset+n_lines << ") = {";
+ for (size_t k(0); k < n_lines - 1; k++)
+ os << line_offset + k << ",";
+ os << line_offset + n_lines - 1 << "};" << std::endl;
+
+ if (_is_sfc) {
+ // write plane surface
+ os << "Plane Surface (" << sfc_offset << ") = {" << line_offset+n_lines << "};" << std::endl;
+ }
+
+}
+
+} // end namespace FileIO
diff --git a/FileIO/MeshIO/GMSHLineLoop.h b/FileIO/MeshIO/GMSHLineLoop.h
new file mode 100644
index 00000000000..51dbc4a73b0
--- /dev/null
+++ b/FileIO/MeshIO/GMSHLineLoop.h
@@ -0,0 +1,34 @@
+/*
+ * GMSHLineLoop.h
+ *
+ * Created on: Mar 22, 2012
+ * Author: fischeth
+ */
+
+#ifndef GMSHLINELOOP_H_
+#define GMSHLINELOOP_H_
+
+#include <vector>
+#include <ostream>
+
+#include "GMSHLine.h"
+
+namespace FileIO {
+
+class GMSHLineLoop {
+public:
+ GMSHLineLoop(bool is_sfc=false);
+ virtual ~GMSHLineLoop();
+ void addLine(GMSHLine* line);
+ bool isSurface() const { return _is_sfc; }
+ void setSurface(bool is_sfc) { _is_sfc = is_sfc; }
+ void write(std::ostream &os, size_t offset, size_t sfc_offset = 0) const;
+
+private:
+ std::vector<GMSHLine*> _lines;
+ bool _is_sfc;
+};
+
+}
+
+#endif /* GMSHLINELOOP_H_ */
diff --git a/FileIO/MeshIO/GMSHMeshDensityStrategy.h b/FileIO/MeshIO/GMSHMeshDensityStrategy.h
new file mode 100644
index 00000000000..7b70f6eb075
--- /dev/null
+++ b/FileIO/MeshIO/GMSHMeshDensityStrategy.h
@@ -0,0 +1,33 @@
+/*
+ * GMSHMeshDensityStrategy.h
+ *
+ * Created on: Mar 5, 2012
+ * Author: TF
+ */
+
+#ifndef GMSHMESHDENSITYSTRATEGY_H_
+#define GMSHMESHDENSITYSTRATEGY_H_
+
+#include <ostream>
+#include <vector>
+
+// GEOLIB
+#include "Point.h"
+
+namespace FileIO
+{
+/**
+ * virtual base class GMSHMeshDensityStrategy for classes
+ * GMSHAdaptiveMeshDensity, GMSHFixedMeshDensity and GMSHNoMeshDensity
+ */
+class GMSHMeshDensityStrategy
+{
+public:
+ virtual void init(std::vector<GEOLIB::Point const*> const&) = 0;
+ virtual double getMeshDensityAtPoint(GEOLIB::Point const*const) const = 0;
+};
+
+} // end namespace
+
+
+#endif /* GMSHMESHDENSITYSTRATEGY_H_ */
diff --git a/FileIO/MeshIO/GMSHNoMeshDensity.h b/FileIO/MeshIO/GMSHNoMeshDensity.h
new file mode 100644
index 00000000000..c2ef03c0718
--- /dev/null
+++ b/FileIO/MeshIO/GMSHNoMeshDensity.h
@@ -0,0 +1,34 @@
+/*
+ * GMSHNoMeshDensity.h
+ *
+ * Created on: Mar 5, 2012
+ * Author: fischeth
+ */
+
+#ifndef GMSHNOMESHDENSITY_H_
+#define GMSHNOMESHDENSITY_H_
+
+#include "GMSHMeshDensityStrategy.h"
+
+namespace FileIO {
+
+class GMSHNoMeshDensity: public FileIO::GMSHMeshDensityStrategy {
+public:
+ GMSHNoMeshDensity() {};
+ void init(std::vector<GEOLIB::Point const*> const& vec)
+ {
+ // to avoid a warning here:
+ (void)(vec);
+ }
+
+ double getMeshDensityAtPoint(GEOLIB::Point const*const pnt) const
+ {
+ // to avoid a warning here:
+ (void)(pnt);
+ return 0.0;
+ }
+};
+
+}
+
+#endif /* GMSHNOMESHDENSITY_H_ */
diff --git a/FileIO/MeshIO/GMSHPoint.cpp b/FileIO/MeshIO/GMSHPoint.cpp
new file mode 100644
index 00000000000..e4bc3d4d2de
--- /dev/null
+++ b/FileIO/MeshIO/GMSHPoint.cpp
@@ -0,0 +1,38 @@
+/*
+ * GMSHPoint.cpp
+ *
+ * Created on: Mar 21, 2012
+ * Author: TF
+ */
+
+#include <cmath>
+#include <limits>
+
+#include "MeshIO/GMSHPoint.h"
+
+namespace FileIO {
+
+GMSHPoint::GMSHPoint(GEOLIB::Point const& pnt, size_t id, double mesh_density) :
+ GEOLIB::PointWithID(pnt, id), _mesh_density(mesh_density)
+{}
+
+void GMSHPoint::write(std::ostream &os) const
+{
+ os << "Point(" << _id << ") = {" << _x[0] << "," << _x[1] << ", 0.0";
+ if (fabs(_mesh_density) > std::numeric_limits<double>::epsilon()) {
+ os << ", " << _mesh_density << "};";
+ } else {
+ os << "};";
+ }
+}
+
+GMSHPoint::~GMSHPoint()
+{}
+
+std::ostream& operator<< (std::ostream &os, GMSHPoint const& p)
+{
+ p.write (os);
+ return os;
+}
+
+} // end namespace FileIO
diff --git a/FileIO/MeshIO/GMSHPoint.h b/FileIO/MeshIO/GMSHPoint.h
new file mode 100644
index 00000000000..ee487226320
--- /dev/null
+++ b/FileIO/MeshIO/GMSHPoint.h
@@ -0,0 +1,30 @@
+/*
+ * GMSHPoint.h
+ *
+ * Created on: Mar 21, 2012
+ * Author: TF
+ */
+
+#ifndef GMSHPOINT_H_
+#define GMSHPOINT_H_
+
+// GEOLIB
+#include "PointWithID.h"
+
+namespace FileIO {
+
+class GMSHPoint : public GEOLIB::PointWithID {
+public:
+ GMSHPoint(GEOLIB::Point const& pnt, size_t id, double mesh_density);
+ virtual ~GMSHPoint();
+ void write(std::ostream &os) const;
+private:
+ double _mesh_density;
+};
+
+/** overload the output operator for class GMSHPoint */
+std::ostream& operator<< (std::ostream &os, GMSHPoint const& p);
+
+}
+
+#endif /* GMSHPOINT_H_ */
diff --git a/FileIO/MeshIO/GMSHPolygonTree.cpp b/FileIO/MeshIO/GMSHPolygonTree.cpp
new file mode 100644
index 00000000000..e8f0fc57046
--- /dev/null
+++ b/FileIO/MeshIO/GMSHPolygonTree.cpp
@@ -0,0 +1,301 @@
+/*
+ * GMSHPolygonTree.cpp
+ *
+ * Created on: Mar 27, 2012
+ * Author: fischeth
+ */
+
+#include "MeshIO/GMSHPolygonTree.h"
+
+#include "GMSHNoMeshDensity.h"
+#include "GMSHFixedMeshDensity.h"
+#include "GMSHAdaptiveMeshDensity.h"
+
+namespace FileIO {
+
+GMSHPolygonTree::GMSHPolygonTree(GEOLIB::Polygon* polygon, GMSHPolygonTree* parent,
+ GEOLIB::GEOObjects &geo_objs, std::string const& geo_name,
+ GMSHMeshDensityStrategy * mesh_density_strategy) :
+ GEOLIB::SimplePolygonTree(polygon, parent), _geo_objs(geo_objs), _geo_name(geo_name),
+ _mesh_density_strategy(mesh_density_strategy)
+{}
+
+GMSHPolygonTree::~GMSHPolygonTree()
+{}
+
+bool GMSHPolygonTree::insertStation(GEOLIB::Point const* station)
+{
+
+ if (_node_polygon->isPntInPolygon(*station)) {
+ // try to insert station into the child nodes
+ for (std::list<SimplePolygonTree*>::const_iterator it (_childs.begin());
+ it != _childs.end(); it++) {
+ if (((*it)->getPolygon())->isPntInPolygon (*station)) {
+ return dynamic_cast<GMSHPolygonTree*>((*it))->insertStation (station);
+ }
+ }
+ // station did not fit into child nodes -> insert the station into this node
+ _stations.push_back (station);
+ return true;
+ } else {
+ return false;
+ }
+}
+
+void GMSHPolygonTree::insertPolyline (GEOLIB::PolylineWithSegmentMarker * ply)
+{
+ if (_node_polygon->isPartOfPolylineInPolygon(*ply)) {
+ // check childs
+ for (std::list<SimplePolygonTree*>::const_iterator it (_childs.begin());
+ it != _childs.end(); it++) {
+ dynamic_cast<GMSHPolygonTree*>((*it))->insertPolyline (ply);
+ }
+ _plys.push_back(ply);
+
+ // calculate possible intersection points
+ // pay attention: loop bound is not fix!
+ size_t n_segments (ply->getNumberOfPoints()-1);
+ GEOLIB::Point tmp_pnt;
+ for (size_t k(0); k<n_segments; k++) {
+ if (! ply->isSegmentMarked(k)) {
+ size_t seg_num(0);
+ GEOLIB::Point *intersection_pnt(new GEOLIB::Point);
+ while (_node_polygon->getNextIntersectionPointPolygonLine(*(ply->getPoint(k)), *(ply->getPoint(k+1)), intersection_pnt, seg_num)) {
+ // insert the intersection point to point vector of GEOObjects instance
+ size_t pnt_id;
+ // appendPoints deletes an already existing point
+ if (_geo_objs.appendPoint(intersection_pnt, _geo_name, pnt_id)) {
+ // case: new point
+ // modify the polygon
+ _node_polygon->insertPoint(seg_num+1, pnt_id);
+ // modify the polyline
+ ply->insertPoint(k+1, pnt_id);
+ n_segments++;
+ } else {
+ // case: existing point
+ // check if point id is within the polygon
+ if (! _node_polygon->isPointIDInPolyline(pnt_id)) {
+ _node_polygon->insertPoint(seg_num+1, pnt_id);
+ }
+ // check if point id is in polyline
+ if (! ply->isPointIDInPolyline(pnt_id)) {
+ ply->insertPoint(k+1, pnt_id);
+ n_segments++;
+ }
+ }
+
+ size_t tmp_seg_num(seg_num+1);
+ if (! _node_polygon->getNextIntersectionPointPolygonLine(*(ply->getPoint(k)), *(ply->getPoint(k+1)), &tmp_pnt, tmp_seg_num)) {
+ // check a point of the segment except the end points
+ for (size_t i(0); i<3; i++) {
+ tmp_pnt[i] = ((*(ply->getPoint(k)))[i] + (*(ply->getPoint(k+1)))[i]) / 2;
+ }
+ if (_node_polygon->isPntInPolygon(tmp_pnt)) {
+ ply->markSegment(k, true);
+ // insert line segment as constraint
+ _gmsh_lines_for_constraints.push_back(new GMSHLine(ply->getPointID(k), ply->getPointID(k+1)));
+ }
+ }
+ intersection_pnt = new GEOLIB::Point;
+ seg_num++;
+ }
+
+ // check a point of the segment except the end points
+ for (size_t i(0); i<3; i++) {
+ tmp_pnt[i] = ((*(ply->getPoint(k)))[i] + (*(ply->getPoint(k+1)))[i]) / 2;
+ }
+ if (_node_polygon->isPntInPolygon(tmp_pnt)) {
+ ply->markSegment(k, true);
+ // insert line segment as constraint
+ _gmsh_lines_for_constraints.push_back(new GMSHLine(ply->getPointID(k), ply->getPointID(k+1)));
+ }
+ }
+ }
+ }
+}
+
+void GMSHPolygonTree::initMeshDensityStrategy()
+{
+ if (dynamic_cast<GMSHAdaptiveMeshDensity*> (_mesh_density_strategy)) {
+ // collect points
+ std::vector<GEOLIB::Point const*> pnts;
+ const size_t n_pnts_polygon (_node_polygon->getNumberOfPoints());
+ for (size_t k(0); k<n_pnts_polygon; k++) {
+ pnts.push_back(_node_polygon->getPoint(k));
+ }
+ getPointsFromSubPolygons(pnts);
+
+ const size_t n_plys (_plys.size());
+ for (size_t k(0); k<n_plys; k++) {
+ const size_t n_pnts_in_kth_ply(_plys[k]->getNumberOfPoints());
+ for (size_t j(0); j<n_pnts_in_kth_ply; j++) {
+ pnts.push_back(_plys[k]->getPoint(j));
+ }
+ }
+
+ // give collected points to the mesh density strategy
+ _mesh_density_strategy->init(pnts);
+ // insert constraints
+ dynamic_cast<GMSHAdaptiveMeshDensity*>(_mesh_density_strategy)->addPoints(_stations);
+ std::vector<GEOLIB::Point const*> stations;
+ getStationsInsideSubPolygons(stations);
+ dynamic_cast<GMSHAdaptiveMeshDensity*>(_mesh_density_strategy)->addPoints(stations);
+ }
+}
+
+void GMSHPolygonTree::createGMSHPoints(std::vector<FileIO::GMSHPoint*> & gmsh_pnts) const
+{
+ const size_t n_pnts_polygon (_node_polygon->getNumberOfPoints());
+ for (size_t k(0); k<n_pnts_polygon; k++) {
+ const size_t id (_node_polygon->getPointID(k));
+ GEOLIB::Point const*const pnt(_node_polygon->getPoint(k));
+ gmsh_pnts[id] = new GMSHPoint(*pnt, id, _mesh_density_strategy->getMeshDensityAtPoint(pnt));
+ }
+
+ const size_t n_plys(_plys.size());
+ for (size_t k(0); k<n_plys; k++) {
+ const size_t n_pnts_in_ply(_plys[k]->getNumberOfPoints());
+ for (size_t j(0); j<n_pnts_in_ply; j++) {
+ if (_node_polygon->isPntInPolygon(*(_plys[k]->getPoint(j)))) {
+ const size_t id (_plys[k]->getPointID(j));
+ GEOLIB::Point const*const pnt(_plys[k]->getPoint(j));
+ gmsh_pnts[id] = new GMSHPoint(*pnt, id, _mesh_density_strategy->getMeshDensityAtPoint(pnt));
+ }
+ }
+ }
+
+ // walk through childs
+ for (std::list<SimplePolygonTree*>::const_iterator it (_childs.begin()); it != _childs.end(); it++) {
+ dynamic_cast<GMSHPolygonTree*>((*it))->createGMSHPoints(gmsh_pnts);
+ }
+}
+
+void GMSHPolygonTree::writeLineLoop(size_t &line_offset, size_t &sfc_offset, std::ostream& out) const
+{
+ const size_t n_pnts (_node_polygon->getNumberOfPoints());
+ size_t first_pnt_id(_node_polygon->getPointID(0)), second_pnt_id;
+ for (size_t k(1); k<n_pnts; k++) {
+ second_pnt_id = _node_polygon->getPointID(k);
+ out << "Line(" << line_offset + k-1 << ") = {" << first_pnt_id << "," << second_pnt_id << "};" << std::endl;
+ first_pnt_id = second_pnt_id;
+ }
+ out << "Line Loop(" << line_offset + n_pnts-1 << ") = {";
+ for (size_t k(0); k<n_pnts - 2; k++) {
+ out << line_offset+k << ",";
+ }
+ out << line_offset+n_pnts-2 << "};" << std::endl;
+ out << "Plane Surface(" << sfc_offset << ") = {" << line_offset+n_pnts-1 << "};" << std::endl;
+ line_offset += n_pnts;
+ sfc_offset++;
+}
+
+void GMSHPolygonTree::writeLineConstraints(size_t &line_offset, size_t sfc_number, std::ostream& out) const
+{
+ const size_t n_plys (_plys.size());
+ for (size_t j(0); j<n_plys; j++) {
+ const size_t n_pnts(_plys[j]->getNumberOfPoints());
+ size_t first_pnt_id(_plys[j]->getPointID(0)), second_pnt_id;
+ for (size_t k(1); k<n_pnts; k++) {
+ second_pnt_id = _plys[j]->getPointID(k);
+ if (_plys[j]->isSegmentMarked(k-1) && _node_polygon->isPntInPolygon(*(_plys[j]->getPoint(k)))) {
+ out << "Line(" << line_offset + k-1 << ") = {" << first_pnt_id << "," << second_pnt_id << "};" << std::endl;
+ out << "Line { " << line_offset+k-1 << " } In Surface { " << sfc_number << " };" << std::endl;
+ }
+ first_pnt_id = second_pnt_id;
+ }
+ line_offset += n_pnts;
+ }
+}
+
+void GMSHPolygonTree::writeSubPolygonsAsLineConstraints(size_t &line_offset, size_t sfc_number, std::ostream& out) const
+{
+ for (std::list<SimplePolygonTree*>::const_iterator it (_childs.begin()); it != _childs.end(); it++) {
+ dynamic_cast<GMSHPolygonTree*>((*it))->writeSubPolygonsAsLineConstraints(line_offset, sfc_number, out);
+ }
+
+ if (_parent != NULL) {
+ const size_t n_pnts(_node_polygon->getNumberOfPoints());
+ size_t first_pnt_id(_node_polygon->getPointID(0)), second_pnt_id;
+ for (size_t k(1); k<n_pnts; k++) {
+ second_pnt_id = _node_polygon->getPointID(k);
+ out << "Line(" << line_offset + k-1 << ") = {" << first_pnt_id << "," << second_pnt_id << "};" << std::endl;
+ first_pnt_id = second_pnt_id;
+ out << "Line { " << line_offset+k-1 << " } In Surface { " << sfc_number << " };" << std::endl;
+ }
+ line_offset += n_pnts;
+ }
+
+}
+
+void GMSHPolygonTree::writeStations(size_t & pnt_id_offset, size_t sfc_number, std::ostream& out) const
+{
+ const size_t n_stations(_stations.size());
+ for (size_t k(0); k<n_stations; k++) {
+ out << "Point(" << pnt_id_offset + k << ") = {" << (*(_stations[k]))[0] << "," << (*(_stations[k]))[1] << ", 0.0, ";
+ out << _mesh_density_strategy->getMeshDensityAtPoint(_stations[k]) << "};" << std::endl;
+ out << "Point { " << pnt_id_offset + k << " } In Surface { " << sfc_number << " }; " << std::endl;
+ }
+ pnt_id_offset += n_stations;
+}
+
+void GMSHPolygonTree::writeAdditionalPointData(size_t & pnt_id_offset, size_t sfc_number, std::ostream& out) const
+{
+ if (dynamic_cast<GMSHAdaptiveMeshDensity*>(_mesh_density_strategy)) {
+ std::vector<GEOLIB::Point*> steiner_pnts;
+ dynamic_cast<GMSHAdaptiveMeshDensity*>(_mesh_density_strategy)->getSteinerPoints(steiner_pnts, 0);
+ const size_t n(steiner_pnts.size());
+ for (size_t k(0); k<n; k++) {
+ if (_node_polygon->isPntInPolygon(*(steiner_pnts[k]))) {
+ out << "Point(" << pnt_id_offset + k << ") = {" << (*(steiner_pnts[k]))[0] << "," << (*(steiner_pnts[k]))[1] << ", 0.0, ";
+ out << _mesh_density_strategy->getMeshDensityAtPoint(steiner_pnts[k]) << "};" << std::endl;
+ out << "Point { " << pnt_id_offset + k << " } In Surface { " << sfc_number << " }; " << std::endl;
+ }
+ delete steiner_pnts[k];
+ }
+ pnt_id_offset += n;
+ }
+
+#ifndef NDEBUG
+ if (dynamic_cast<GMSHAdaptiveMeshDensity*>(_mesh_density_strategy)) {
+ std::vector<GEOLIB::Point*> * pnts(new std::vector<GEOLIB::Point*>);
+ std::vector<GEOLIB::Polyline*> *plys(new std::vector<GEOLIB::Polyline*>);
+ dynamic_cast<GMSHAdaptiveMeshDensity*>(_mesh_density_strategy)->getQuadTreeGeometry(*pnts, *plys);
+ std::string quad_tree_geo("QuadTree");
+ _geo_objs.addPointVec(pnts, quad_tree_geo);
+ std::vector<size_t> const& id_map ((_geo_objs.getPointVecObj(quad_tree_geo))->getIDMap());
+ for (size_t k(0); k<plys->size(); k++) {
+ for (size_t j(0); j<(*plys)[k]->getNumberOfPoints(); j++) {
+ ((*plys)[k])->setPointID(j, id_map[((*plys)[k])->getPointID(j)]);
+ }
+ }
+ _geo_objs.addPolylineVec(plys, quad_tree_geo);
+ }
+#endif
+
+}
+
+void GMSHPolygonTree::getPointsFromSubPolygons(std::vector<GEOLIB::Point const*>& pnts)
+{
+ const size_t n_pnts_polygon (_node_polygon->getNumberOfPoints());
+ for (size_t k(0); k<n_pnts_polygon; k++) {
+ pnts.push_back(_node_polygon->getPoint(k));
+ }
+
+ for (std::list<SimplePolygonTree*>::const_iterator it (_childs.begin()); it != _childs.end(); it++) {
+ dynamic_cast<GMSHPolygonTree*>((*it))->getPointsFromSubPolygons(pnts);
+ }
+}
+
+void GMSHPolygonTree::getStationsInsideSubPolygons(std::vector<GEOLIB::Point const*>& stations)
+{
+ const size_t n_stations(_stations.size());
+ for (size_t k(0); k<n_stations; k++) {
+ stations.push_back(_stations[k]);
+ }
+
+ for (std::list<SimplePolygonTree*>::const_iterator it (_childs.begin()); it != _childs.end(); it++) {
+ dynamic_cast<GMSHPolygonTree*>((*it))->getStationsInsideSubPolygons(stations);
+ }
+}
+
+} // end namespace FileIO
diff --git a/FileIO/MeshIO/GMSHPolygonTree.h b/FileIO/MeshIO/GMSHPolygonTree.h
new file mode 100644
index 00000000000..440be4e0993
--- /dev/null
+++ b/FileIO/MeshIO/GMSHPolygonTree.h
@@ -0,0 +1,91 @@
+/*
+ * GMSHPolygonTree.h
+ *
+ * Created on: Mar 27, 2012
+ * Author: fischeth
+ */
+
+#ifndef GMSHPOLYGONTREE_H_
+#define GMSHPOLYGONTREE_H_
+
+#include <vector>
+#include <string>
+
+// FileIO
+#include "GMSHMeshDensityStrategy.h"
+#include "GMSHPoint.h"
+#include "GMSHLine.h"
+
+// GEOLIB
+#include "SimplePolygonTree.h"
+#include "GEOObjects.h"
+#include "PolylineWithSegmentMarker.h"
+
+namespace FileIO {
+
+class GMSHPolygonTree: public GEOLIB::SimplePolygonTree {
+public:
+ GMSHPolygonTree(GEOLIB::Polygon* polygon, GMSHPolygonTree * parent,
+ GEOLIB::GEOObjects &geo_objs, std::string const& geo_name,
+ GMSHMeshDensityStrategy * mesh_density_strategy);
+ virtual ~GMSHPolygonTree();
+
+ /**
+ * If the station point is inside the polygon, the method inserts the station into
+ * the internal vector of stations. This method works recursive!
+ * @param pnt the station point
+ * @return true if the station is inside the polygon
+ */
+ bool insertStation(GEOLIB::Point const* pnt);
+ /**
+ * If at least one (end) point (of a line segment) of the polyline is inside the polygon
+ * the polyline is inserted to the internal vector of polylines.
+ *
+ * Intersection points are inserted into the points vector the polygon and the polyline
+ * are based on. The id of the intersection point is inserted in both the polygon and the
+ * polyline, i.e. the two intersecting line segments are splitt into four line segment.
+ *
+ * Line segments of the polyline that are completely within the polygon are inserted into
+ * the internal vector _gmsh_lines_for_constraints. The childs of this GMSHPolygonTree node
+ * are checked recursively.
+ * @param ply the polyline that should be inserted
+ */
+ void insertPolyline(GEOLIB::PolylineWithSegmentMarker * ply);
+
+ /**
+ * Initialize the mesh density strategy with data. In case of GMSHAdaptiveMeshDensity
+ * an instance of class QuadTree will be set up with the points within the top
+ * level bounding polygon.
+ */
+ void initMeshDensityStrategy();
+
+ /**
+ * Method creates the gmsh point data structures - including the mesh density.
+ * @param gmsh_pnts a vector of pointers to instances of class GMSHPoint
+ */
+ void createGMSHPoints(std::vector<FileIO::GMSHPoint*> & gmsh_pnts) const;
+
+ virtual void writeLineLoop(size_t &line_offset, size_t &sfc_offset, std::ostream& out) const;
+ void writeSubPolygonsAsLineConstraints(size_t &line_offset, size_t sfc_number, std::ostream& out) const;
+ virtual void writeLineConstraints(size_t &line_offset, size_t sfc_number, std::ostream& out) const;
+ void writeStations(size_t & pnt_id_offset, size_t sfc_number, std::ostream& out) const;
+ void writeAdditionalPointData(size_t & pnt_id_offset, size_t sfc_number, std::ostream& out) const;
+
+private:
+ void getPointsFromSubPolygons(std::vector<GEOLIB::Point const*>& pnts);
+ void getStationsInsideSubPolygons(std::vector<GEOLIB::Point const*>& stations);
+ const std::list<SimplePolygonTree*>& getChilds() const;
+ const std::list<GEOLIB::GEOObjects*>& getGeoObjects () const;
+
+ GEOLIB::GEOObjects & _geo_objs;
+ std::string const& _geo_name;
+ std::vector<GEOLIB::Point const*> _stations;
+ std::vector<GEOLIB::PolylineWithSegmentMarker*> _plys;
+ std::vector<FileIO::GMSHLine*> _gmsh_lines_for_constraints;
+
+ GMSHMeshDensityStrategy * _mesh_density_strategy;
+};
+
+}
+
+#endif /* GMSHPOLYGONTREE_H_ */
diff --git a/FileIO/MeshIO/TetGenInterface.cpp b/FileIO/MeshIO/TetGenInterface.cpp
new file mode 100644
index 00000000000..dbd78c8c309
--- /dev/null
+++ b/FileIO/MeshIO/TetGenInterface.cpp
@@ -0,0 +1,528 @@
+/*
+ * TetGenInterface.cpp
+ *
+ * Created on: Sep 12, 2011
+ * Author: TF
+ */
+
+#include <cstddef>
+#include <string>
+
+// FileIO
+#include "MeshIO/TetGenInterface.h"
+
+// Base
+#include "StringTools.h"
+
+// MSH
+#include "msh_elem.h"
+#include "msh_mesh.h"
+#include "msh_node.h"
+
+namespace FileIO
+{
+TetGenInterface::TetGenInterface() :
+ _mesh (NULL), _zero_based_idx (false)
+{
+}
+
+TetGenInterface::~TetGenInterface()
+{
+}
+
+void TetGenInterface::writeTetGenMesh(std::string const& nodes_fname, std::string const& ele_fname,
+ MeshLib::CFEMesh const* const mesh) const
+{
+ writeTetGenNodes(nodes_fname, mesh);
+ writeTetGenElements(ele_fname, mesh);
+}
+
+void TetGenInterface::writeTetGenNodes(std::string const& nodes_fname, MeshLib::CFEMesh const*const mesh) const
+{
+ std::ofstream out(nodes_fname.c_str());
+ if (out) {
+ std::vector<MeshLib::CElem*> const& elements(mesh->getElementVector());
+ const size_t n_elements(elements.size());
+ size_t n_prisms(0);
+ for (size_t k(0); k<n_elements; k++) {
+ if (elements[k]->GetElementType() == MshElemType::PRISM) {
+ n_prisms++;
+ }
+ }
+ const size_t n_nodes(mesh->GetNodesNumber(false));
+ out << n_nodes+4*n_prisms << " 3 0 0" << std::endl;
+ std::vector<MeshLib::CNode*> const& nodes(mesh->getNodeVector());
+ for (size_t k(0); k<n_nodes; k++) {
+ double const*const node(nodes[k]->getData());
+ out << k << " " << node[0] << " " << node[1] << " " << node[2] << std::endl;
+ }
+ // write additional nodes for prisms
+ std::vector<size_t> idxs;
+ for (size_t k(0), idx(n_nodes); k<n_elements; k++) {
+ if (elements[k]->GetElementType() == MshElemType::PRISM) {
+ elements[k]->getNodeIndices(idxs);
+ double const*const n0((nodes[idxs[0]])->getData());
+ double const*const n1((nodes[idxs[1]])->getData());
+ double const*const n2((nodes[idxs[2]])->getData());
+ double const*const n3((nodes[idxs[3]])->getData());
+ double const*const n4((nodes[idxs[4]])->getData());
+ double const*const n5((nodes[idxs[5]])->getData());
+
+ idxs.clear();
+
+ // compute centroid of the prism
+ const double centroid[3] = {(n0[0]+n1[0]+n2[0]+n3[0]+n4[0]+n5[0])/6,
+ (n0[1]+n1[1]+n2[1]+n3[1]+n4[1]+n5[1])/6,
+ (n0[2]+n1[2]+n2[2]+n3[2]+n4[2]+n5[2])/6};
+
+ // center of first rectangle surface
+ const double center0[3] = {(n0[0]+n1[0]+n3[0]+n4[0])/4, (n0[1]+n1[1]+n3[1]+n4[1])/4, (n0[2]+n1[2]+n3[2]+n4[2])/4};
+ // center of second rectangle surface
+ const double center1[3] = {(n1[0]+n2[0]+n4[0]+n5[0])/4, (n1[1]+n2[1]+n4[1]+n5[1])/4, (n1[2]+n2[2]+n4[2]+n5[2])/4};
+ // center of third rectangle surface
+ const double center2[3] = {(n0[0]+n2[0]+n3[0]+n5[0])/4, (n0[1]+n2[1]+n3[1]+n5[1])/4, (n0[2]+n2[2]+n3[2]+n5[2])/4};
+
+ out << idx++ << " " << centroid[0] << " " << centroid[1] << " " << centroid[2] << std::endl;
+ out << idx++ << " " << center0[0] << " " << center0[1] << " " << center0[2] << std::endl;
+ out << idx++ << " " << center1[0] << " " << center1[1] << " " << center1[2] << std::endl;
+ out << idx++ << " " << center2[0] << " " << center2[1] << " " << center2[2] << std::endl;
+ }
+ }
+ out.close();
+ } else {
+ std::cout << "cold not open file for writing nodes" << std::endl;
+ }
+}
+
+void TetGenInterface::writeTetGenElements(std::string const& ele_fname, MeshLib::CFEMesh const*const mesh) const
+{
+ std::ofstream out(ele_fname.c_str());
+ if (out) {
+ std::vector<MeshLib::CElem*> const& elements(mesh->getElementVector());
+ const size_t n_elements(elements.size());
+ // count number of prisms, tetrahedras, hexahedras
+ size_t n_prisms(0), n_tets(0), n_hexs(0);
+ for (size_t k(0); k<n_elements; k++) {
+ switch (elements[k]->GetElementType()) {
+ case MshElemType::PRISM:
+ n_prisms++;
+ break;
+ case MshElemType::TETRAHEDRON:
+ n_tets++;
+ break;
+ case MshElemType::HEXAHEDRON:
+ n_hexs++;
+ break;
+ default:
+ std::cout << "count elements - element type not yet supported" << std::endl;
+ } // end case
+ } // end for
+
+ const size_t n_tetrahedras(n_tets+14*n_prisms);
+ const size_t nodes_offset(mesh->GetNodesNumber(false));
+ size_t cnt_prisms(0);
+ std::vector<size_t> idxs; // node indices
+ out << n_tetrahedras << " 4 0 0" << std::endl;
+ for (size_t k(0), cnt(0); k<n_elements; k++) {
+ elements[k]->getNodeIndices(idxs);
+ switch (elements[k]->GetElementType()) {
+ case MshElemType::PRISM:
+ {
+ out << cnt++ << " " << idxs[0] << " " << idxs[1] << " " << idxs[2] << " " << nodes_offset+cnt_prisms*4 << std::endl;
+ out << cnt++ << " " << idxs[3] << " " << idxs[5] << " " << idxs[4] << " " << nodes_offset+cnt_prisms*4 << std::endl;
+ //
+ out << cnt++ << " " << idxs[0] << " " << nodes_offset+cnt_prisms*4+1 << " " << idxs[1] << " " << nodes_offset+cnt_prisms*4 << std::endl;
+ out << cnt++ << " " << idxs[1] << " " << nodes_offset+cnt_prisms*4+1 << " " << idxs[4] << " " << nodes_offset+cnt_prisms*4 << std::endl;
+ out << cnt++ << " " << idxs[4] << " " << nodes_offset+cnt_prisms*4+1 << " " << idxs[3] << " " << nodes_offset+cnt_prisms*4 << std::endl;
+ out << cnt++ << " " << idxs[3] << " " << nodes_offset+cnt_prisms*4+1 << " " << idxs[0] << " " << nodes_offset+cnt_prisms*4 << std::endl;
+ //
+ out << cnt++ << " " << idxs[1] << " " << nodes_offset+cnt_prisms*4+2 << " " << idxs[2] << " " << nodes_offset+cnt_prisms*4 << std::endl;
+ out << cnt++ << " " << idxs[4] << " " << nodes_offset+cnt_prisms*4+2 << " " << idxs[1] << " " << nodes_offset+cnt_prisms*4 << std::endl;
+ out << cnt++ << " " << idxs[5] << " " << nodes_offset+cnt_prisms*4+2 << " " << idxs[4] << " " << nodes_offset+cnt_prisms*4 << std::endl;
+ out << cnt++ << " " << idxs[2] << " " << nodes_offset+cnt_prisms*4+2 << " " << idxs[5] << " " << nodes_offset+cnt_prisms*4 << std::endl;
+ //
+ out << cnt++ << " " << idxs[2] << " " << nodes_offset+cnt_prisms*4+3 << " " << idxs[0] << " " << nodes_offset+cnt_prisms*4 << std::endl;
+ out << cnt++ << " " << idxs[5] << " " << nodes_offset+cnt_prisms*4+3 << " " << idxs[2] << " " << nodes_offset+cnt_prisms*4 << std::endl;
+ out << cnt++ << " " << idxs[3] << " " << nodes_offset+cnt_prisms*4+3 << " " << idxs[5] << " " << nodes_offset+cnt_prisms*4 << std::endl;
+ out << cnt++ << " " << idxs[0] << " " << nodes_offset+cnt_prisms*4+3 << " " << idxs[3] << " " << nodes_offset+cnt_prisms*4 << std::endl;
+
+ cnt_prisms++;
+
+ break;
+ }
+ case MshElemType::TETRAHEDRON:
+ out << cnt++ << " " << idxs[0] << " " << idxs[1] << " " << idxs[2] << " " << idxs[3] << std::endl;
+ break;
+ case MshElemType::HEXAHEDRON:
+ std::cout << "element type HEXAHEDRON not yet supported" << std::endl;
+ break;
+ case MshElemType::PYRAMID:
+ std::cout << "element type PYRAMID not yet supported" << std::endl;
+ break;
+ case MshElemType::LINE:
+ std::cout << "element type LINE not yet supported" << std::endl;
+ break;
+ case MshElemType::QUAD:
+ std::cout << "element type QUAD not yet supported" << std::endl;
+ break;
+ default:
+ std::cout << "element type not yet supported" << std::endl;
+ } // end case
+ idxs.clear();
+ } // end for
+ out.close();
+ } else {
+ std::cout << "cold not open file for writing elements" << std::endl;
+ }
+}
+
+MeshLib::CFEMesh* TetGenInterface::readTetGenMesh (std::string const& nodes_fname,
+ std::string const& ele_fname)
+{
+ std::ifstream ins_nodes (nodes_fname.c_str());
+ std::ifstream ins_ele (ele_fname.c_str());
+
+ if (!ins_nodes || !ins_ele)
+ {
+ if (!ins_nodes)
+ std::cout << "TetGenInterface::readTetGenMesh failed to open " <<
+ nodes_fname << std::endl;
+ if (!ins_ele)
+ std::cout << "TetGenInterface::readTetGenMesh failed to open " <<
+ ele_fname << std::endl;
+ return NULL;
+ }
+
+ _mesh = new MeshLib::CFEMesh();
+
+ if (!readNodesFromStream (ins_nodes))
+ {
+ delete _mesh;
+ return NULL;
+ }
+
+ _mesh->InitialNodesNumber();
+
+ if (!readElementsFromStream (ins_ele))
+ {
+ delete _mesh;
+ return NULL;
+ }
+
+ return _mesh;
+}
+
+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;
+ bool not_read_header (true);
+
+ while (!ins.fail() && not_read_header)
+ {
+ line = line.substr(pos_beg);
+ if (line.compare(0,1,"#") == 0)
+ {
+ // this line is a comment - skip
+ getline (ins, line);
+ pos_beg = line.find_first_not_of(" ");
+ }
+ else
+ // read header line
+ not_read_header = !parseNodesFileHeader(line,
+ n_nodes,
+ dim,
+ n_attributes,
+ boundary_markers);
+ }
+ if (not_read_header)
+ return false;
+ if (!parseNodes(ins, n_nodes, dim))
+ return false;
+
+ return true;
+}
+
+bool TetGenInterface::parseNodesFileHeader(std::string &line,
+ size_t& n_nodes,
+ size_t& dim,
+ size_t& n_attributes,
+ bool& boundary_markers) const
+{
+ 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 = str2number<size_t> (line.substr(pos_beg, pos_end - pos_beg));
+ else
+ {
+ std::cout <<
+ "TetGenInterface::parseNodesFileHeader could not correct read TetGen mesh header - number of nodes"
+ << std::endl;
+ return false;
+ }
+ // dimension
+ pos_beg = line.find_first_not_of (" ", pos_end);
+ pos_end = line.find_first_of(" ", pos_beg);
+ dim = 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 = 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
+ boundary_markers = false;
+
+ return true;
+}
+
+bool TetGenInterface::parseNodes(std::ifstream& ins, size_t n_nodes, size_t dim)
+{
+ size_t pos_beg, pos_end;
+ std::string line;
+ double* coordinates (static_cast<double*> (alloca (sizeof(double) * dim)));
+
+ for (size_t k(0); k < n_nodes && !ins.fail(); k++)
+ {
+ getline (ins, line);
+ if (!ins.fail())
+ {
+ if (!line.empty())
+ {
+ pos_end = 0;
+ // read id
+ size_t id;
+ pos_beg = line.find_first_not_of(" ", pos_end);
+ pos_end = line.find_first_of(" \n", pos_beg);
+ if (pos_beg != std::string::npos && pos_end != std::string::npos)
+ {
+ id =
+ str2number<size_t>(line.substr(pos_beg, pos_end -
+ pos_beg));
+ if (k == 0 && id == 0)
+ _zero_based_idx = true;
+ }
+ else
+ {
+ std::cout << "error reading id of node " << k <<
+ " in TetGenInterface::parseNodes" << std::endl;
+ return false;
+ }
+ // read coordinates
+ for (size_t i(0); i < dim; i++)
+ {
+ pos_beg = line.find_first_not_of(" ", pos_end);
+ pos_end = line.find_first_of(" \n", pos_beg);
+ if (pos_end == std::string::npos)
+ pos_end = line.size();
+ if (pos_beg != std::string::npos)
+ coordinates[i] =
+ str2number<double>(line.substr(pos_beg,
+ pos_end -
+ pos_beg));
+ else
+ {
+ std::cout << "error reading coordinate " << i <<
+ " of node " << k <<
+ " in TetGenInterface::parseNodes" <<
+ std::endl;
+ return false;
+ }
+ }
+ if (!_zero_based_idx)
+ id--;
+ // since CFEMesh is our friend we can access private data of mesh
+ _mesh->nod_vector.push_back(new MeshLib::CNode(id, coordinates[0],
+ coordinates[1],
+ coordinates[2]));
+ // read attributes and boundary markers ... - at the moment we do not use this information
+ }
+ }
+ else
+ {
+ std::cout << "error reading node " << k <<
+ " in TetGenInterface::parseNodes" << std::endl;
+ return false;
+ }
+ }
+
+ return true;
+}
+
+bool TetGenInterface::readElementsFromStream(std::ifstream &ins)
+{
+ std::string line;
+ getline (ins, line);
+ size_t pos_beg (line.find_first_not_of(" "));
+ size_t n_tets, n_nodes_per_tet;
+ bool region_attributes;
+ bool not_read_header (true);
+
+ while (!ins.fail() && not_read_header)
+ {
+ line = line.substr(pos_beg);
+ if (line.compare(0,1,"#") == 0)
+ {
+ // this line is a comment - skip
+ getline (ins, line);
+ pos_beg = line.find_first_not_of(" ");
+ }
+ else
+ // read header line
+ not_read_header = !parseElementsFileHeader(line,
+ n_tets,
+ n_nodes_per_tet,
+ region_attributes);
+ }
+ if (not_read_header)
+ return false;
+ if (!parseElements(ins, n_tets, n_nodes_per_tet, region_attributes))
+ return false;
+
+ return true;
+}
+
+bool TetGenInterface::parseElementsFileHeader(std::string &line,
+ size_t& n_tets,
+ size_t& n_nodes_per_tet,
+ bool& region_attribute) const
+{
+ size_t pos_beg, pos_end;
+
+ // number of tetrahedras
+ 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_tets = str2number<size_t> (line.substr(pos_beg, pos_end - pos_beg));
+ else
+ {
+ std::cout <<
+ "TetGenInterface::parseElementsFileHeader could not correct read TetGen mesh header - number of tetrahedras"
+ << std::endl;
+ return false;
+ }
+ // nodes per tet - either 4 or 10
+ pos_beg = line.find_first_not_of (" \t", pos_end);
+ pos_end = line.find_first_of(" \t", pos_beg);
+ n_nodes_per_tet = str2number<size_t> (line.substr(pos_beg, pos_end - pos_beg));
+ // region attribute at tetrahedra?
+ pos_beg = line.find_first_not_of (" \t", pos_end);
+ pos_end = line.find_first_of(" \t\n", pos_beg);
+ if (pos_end == std::string::npos)
+ pos_end = line.size();
+ if ((line.substr(pos_beg, pos_end - pos_beg)).compare("1") == 0)
+ region_attribute = true;
+ else
+ region_attribute = false;
+
+ return true;
+}
+
+bool TetGenInterface::parseElements(std::ifstream& ins, size_t n_tets, size_t n_nodes_per_tet,
+ bool region_attribute)
+{
+ size_t pos_beg, pos_end;
+ std::string line;
+ size_t* ids (static_cast<size_t*>(alloca (sizeof (size_t) * n_nodes_per_tet)));
+
+ for (size_t k(0); k < n_tets && !ins.fail(); k++)
+ {
+ getline (ins, line);
+ if (!ins.fail())
+ {
+ if (!line.empty())
+ {
+ pos_end = 0;
+ // read id
+ size_t id;
+ pos_beg = line.find_first_not_of(" ", pos_end);
+ pos_end = line.find_first_of(" \n", pos_beg);
+ if (pos_beg != std::string::npos && pos_end != std::string::npos)
+ id =
+ str2number<size_t>(line.substr(pos_beg, pos_end -
+ pos_beg));
+ else
+ {
+ std::cout << "error reading id of tetrahedra " << k <<
+ " in TetGenInterface::parseElements" << std::endl;
+ return false;
+ }
+ // read node ids
+ for (size_t i(0); i < n_nodes_per_tet; i++)
+ {
+ 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 (pos_beg != std::string::npos && pos_end !=
+ std::string::npos)
+ ids[i] =
+ str2number<size_t>(line.substr(pos_beg,
+ pos_end -
+ pos_beg));
+ else
+ {
+ std::cout << "error reading node " << i <<
+ " of tetrahedra " << k <<
+ " in TetGenInterface::parseElements" <<
+ std::endl;
+ return false;
+ }
+ }
+ if (!_zero_based_idx)
+ {
+ id--;
+ for (size_t i(0); i < n_nodes_per_tet; i++)
+ ids[i]--;
+ }
+ // since CFEMesh is our friend we can access private data of mesh
+ MeshLib::CElem* elem (new MeshLib::CElem(id));
+ elem->setElementProperties (MshElemType::TETRAHEDRON, false);
+ std::vector<MeshLib::CNode*> ele_nodes(n_nodes_per_tet);
+ for (size_t i(0); i < n_nodes_per_tet; i++) {
+ ele_nodes[i] = _mesh->nod_vector[ids[i]];
+ }
+ elem->setNodes (ele_nodes);
+ _mesh->ele_vector.push_back(elem);
+ // read region attribute - this is something like material group
+ if (region_attribute)
+ {
+ 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 (pos_beg != std::string::npos && pos_end !=
+ std::string::npos)
+ elem->setPatchIndex (str2number<int>(line.substr(
+ pos_beg,
+ pos_end
+ - pos_beg)));
+ else
+ {
+ std::cout << "error reading region attribute of tetrahedra " << k
+ << " in TetGenInterface::parseElements" << std::endl;
+ return false;
+ }
+ }
+ }
+ }
+ else
+ {
+ std::cout << "error reading node " << k << " in TetGenInterface::parseElements" << std::endl;
+ return false;
+ }
+ }
+ return true;
+}
+}
diff --git a/FileIO/MeshIO/TetGenInterface.h b/FileIO/MeshIO/TetGenInterface.h
new file mode 100644
index 00000000000..91ca5384001
--- /dev/null
+++ b/FileIO/MeshIO/TetGenInterface.h
@@ -0,0 +1,119 @@
+/*
+ * TetGenInterface.h
+ *
+ * Created on: Sep 12, 2011
+ * Author: TF
+ */
+
+#ifndef TETGENINTERFACE_H_
+#define TETGENINTERFACE_H_
+
+// forward declaration of mesh class
+namespace MeshLib
+{
+class CFEMesh;
+}
+
+namespace FileIO
+{
+/**
+ * class TetGenInterface is used to read meshes created by <a href="http://tetgen.berlios.de/">TetGen</a>
+ */
+class TetGenInterface
+{
+public:
+ TetGenInterface();
+ virtual ~TetGenInterface();
+
+ /**
+ * write a mesh into TetGen mesh file format
+ * @param nodes_fname
+ * @param ele_fname
+ * @param mesh
+ */
+ void writeTetGenMesh (std::string const& nodes_fname, std::string const& ele_fname, MeshLib::CFEMesh const*const mesh) const;
+
+ /**
+ * Method reads the TetGen mesh from node file and element file.
+ * @param nodes_fname file name of the nodes file
+ * @param ele_fname file name of the elements file
+ * @return on success the method returns a (pointer to a) CFEMesh, else the method returns NULL
+ */
+ MeshLib::CFEMesh* readTetGenMesh (std::string const& nodes_fname,
+ std::string const& ele_fname);
+ /** in order to have a direct access to the
+ * data structures for nodes and elements we make
+ * class TetGenInterface a friend of the mesh class
+ */
+ friend class MeshLib::CFEMesh;
+
+private:
+ void writeTetGenNodes(std::string const& nodes_fname, MeshLib::CFEMesh const*const mesh) const;
+ void writeTetGenElements(std::string const& ele_fname, MeshLib::CFEMesh const*const mesh) const;
+ /**
+ * Method reads the nodes from stream and stores them in the node vector of the mesh class.
+ * For this purpose it uses methods parseNodesFileHeader() and parseNodes().
+ * @param input the input stream
+ * @return true, if all information is read, false if the method detects an error
+ */
+ bool readNodesFromStream(std::ifstream &input);
+ /**
+ * Method parses the header of the nodes file created by TetGen
+ * @param line the header is in this string (input)
+ * @param n_nodes number of nodes in the file (output)
+ * @param dim the spatial dimension of the node (output)
+ * @param n_attributes the number of attributes for each node (output)
+ * @param boundary_markers have the nodes boundary information (output)
+ * @return true, if the file header is read, false if the method detects an error
+ */
+ bool parseNodesFileHeader(std::string &line, size_t& n_nodes, size_t& dim,
+ size_t& n_attributes, bool& boundary_markers) const;
+ /**
+ * method parses the lines reading the nodes from TetGen nodes file
+ * @param ins the input stream (input)
+ * @param n_nodes the number of nodes to read (input)
+ * @param dim the spatial dimension of the node (input)
+ * @return true, if the nodes are read, false if the method detects an error
+ */
+ bool parseNodes(std::ifstream& ins, size_t n_nodes, size_t dim);
+
+ /**
+ * Method reads the elements from stream and stores them in the element vector of the mesh class.
+ * For this purpose it uses methods parseElementsFileHeader() and parseElements().
+ * @param input the input stream
+ * @return true, if all information is read, false if the method detects an error
+ */
+ bool readElementsFromStream(std::ifstream &input);
+ /**
+ * Method parses the header of the elements file created by TetGen
+ * @param line
+ * @param n_tets
+ * @param n_nodes_per_tet
+ * @param region_attribute is on output true, if there
+ * @return
+ */
+ bool parseElementsFileHeader(std::string &line, size_t& n_tets, size_t& n_nodes_per_tet,
+ bool& region_attribute) const;
+ /**
+ * Method parses the tetrahedras and put them in the element vector of the mesh class.
+ * @param ins the input stream
+ * @param n_tets the number of tetrahedras that should be read
+ * @param n_nodes_per_tet the number of nodes per tetrahedron
+ * @param region_attribute if region attribute is true, region information is read
+ * @return true, if the tetrahedras are read, false if the method detects an error
+ */
+ bool parseElements(std::ifstream& ins, size_t n_tets, size_t n_nodes_per_tet,
+ bool region_attribute);
+
+ /**
+ * the mesh that is returned if all data is read
+ */
+ MeshLib::CFEMesh* _mesh;
+ /**
+ * the value is true if the indexing is zero based, else false
+ */
+ bool _zero_based_idx;
+};
+}
+
+#endif /* TETGENINTERFACE_H_ */
diff --git a/FileIO/NetCDFInterface.cpp b/FileIO/NetCDFInterface.cpp
new file mode 100644
index 00000000000..8a77acfe049
--- /dev/null
+++ b/FileIO/NetCDFInterface.cpp
@@ -0,0 +1,204 @@
+/**
+ * \file NetCDFInterface.cpp
+ * 29/07/2010 YW Initial implementations
+ */
+#include "NetCDFInterface.h"
+#include <stdio.h>
+#include <netcdf.h>
+
+using namespace GEOLIB;
+using namespace MeshLib;
+
+/* Names of variables. */
+#define DIM_RLAT "rlat"
+#define DIM_RLON "rlon"
+#define LAT_NAME "lat"
+#define LON_NAME "lon"
+/* Length of the original output NetCDF file. */
+#define LEN_ORIGINAL 18
+
+/* Handle errors by printing an error message and exiting with a
+ * non-zero status. */
+#define HANDLE_ERROR(e) {printf("Error: %s\n", nc_strerror(e)); exit(-1); }
+
+namespace FileIO
+{
+void NetCDFInterface::readNetCDFData(std::string &ifname,
+ std::vector<GEOLIB::Point*>* points_vec,
+ GEOLIB::GEOObjects* obj,
+ size_t &NRLAT,
+ size_t &NRLON)
+{
+ int ncid;
+ int rlat_dimid, rlon_dimid;
+ int lat_varid, lon_varid;
+ int choose_varid;
+
+ /* The start and count arrays will tell the netCDF library where to read our data. */
+ static size_t start2[2], count2[2];
+ static size_t start3[3], count3[3];
+ static size_t start4[4], count4[4];
+
+ /* Error handling. */
+ int status;
+
+ /* Open the file. */
+ std::cout << "Open NetCDF file " << ifname << "\n" << std::flush;
+ status = nc_open(ifname.c_str(), NC_NOWRITE, &ncid);
+ if (status != NC_NOERR)
+ HANDLE_ERROR(status);
+
+ /* Get the lengths of dimensions, such as "rlat" and "rlat". */
+ status = nc_inq_dimid(ncid, DIM_RLAT, &rlat_dimid);
+ if (status != NC_NOERR)
+ HANDLE_ERROR(status);
+ status = nc_inq_dimlen(ncid, rlat_dimid, &NRLAT);
+ if (status != NC_NOERR)
+ HANDLE_ERROR(status);
+
+ status = nc_inq_dimid(ncid, DIM_RLON, &rlon_dimid);
+ if (status != NC_NOERR)
+ HANDLE_ERROR(status);
+ status = nc_inq_dimlen(ncid, rlon_dimid, &NRLON);
+ if (status != NC_NOERR)
+ HANDLE_ERROR(status);
+
+ /* Get the choosen variable name from the opened NetCDF file. */
+ std::string nc_fname;
+ size_t indexChar;
+ indexChar = ifname.find_last_of('/');
+ if(indexChar == std::string::npos)
+ nc_fname = ifname;
+ else
+ nc_fname = ifname.substr(indexChar + 1);
+ size_t len_varname = nc_fname.size() - LEN_ORIGINAL;
+ std::string var_name;
+ var_name = nc_fname.substr(0,len_varname);
+
+ /* Get the varids of the netCDF variables, such as longitude, latitude,
+ * temperature or precipitation, and so on. */
+ status = nc_inq_varid(ncid, LON_NAME, &lon_varid);
+ if (status != NC_NOERR)
+ HANDLE_ERROR(status);
+ status = nc_inq_varid(ncid, LAT_NAME, &lat_varid);
+ if (status != NC_NOERR)
+ HANDLE_ERROR(status);
+ status = nc_inq_varid(ncid, var_name.c_str(), &choose_varid);
+ if (status != NC_NOERR)
+ HANDLE_ERROR(status);
+
+ /* Program variables to hold the data we will read. We will only
+ need enough space to hold one timestep of data; one record.*/
+ size_t len_vals = NRLAT * NRLON;
+ float* lat_in = new float[len_vals];
+ float* lon_in = new float[len_vals];
+ float* var_in = new float[len_vals];
+
+ /* Read the data. Since we know the contents of the file we know
+ * that the data arrays in this program are the correct size to
+ * hold one timestep. */
+ count2[0] = NRLAT;
+ count2[1] = NRLON;
+ count3[0] = 1;
+ count3[1] = NRLAT;
+ count3[2] = NRLON;
+ count4[0] = 1;
+ count4[1] = 1;
+ count4[2] = NRLAT;
+ count4[3] = NRLON;
+
+ for (size_t i = 0; i < 2; i++)
+ start2[i] = 0;
+ for (size_t i = 0; i < 3; i++)
+ start3[i] = 0;
+ for (size_t i = 0; i < 4; i++)
+ start4[i] = 0;
+
+ /* Read the above netCDF variables with their corresponding varids. */
+ status = nc_get_vara_float(ncid, lon_varid, start2, count2, lon_in);
+ if (status != NC_NOERR)
+ HANDLE_ERROR(status);
+
+ status = nc_get_vara_float(ncid, lat_varid, start2, count2, lat_in);
+ if (status != NC_NOERR)
+ HANDLE_ERROR(status);
+
+ if (var_name.compare("T_2M") == 0)
+ {
+ status = nc_get_vara_float(ncid, choose_varid, start4, count4, var_in);
+ if (status != NC_NOERR)
+ HANDLE_ERROR(status);
+ }
+ else if (var_name.compare("TOT_PREC") == 0)
+ {
+ status = nc_get_vara_float(ncid, choose_varid, start3, count3, var_in);
+ if (status != NC_NOERR)
+ HANDLE_ERROR(status);
+ }
+ else if (var_name.compare("HSURF") == 0)
+ {
+ status = nc_get_vara_float(ncid, choose_varid, start3, count3, var_in);
+ if (status != NC_NOERR)
+ HANDLE_ERROR(status);
+ }
+
+ /* Close the file. */
+ status = nc_close(ncid);
+ if (status != NC_NOERR)
+ HANDLE_ERROR(status);
+
+ printf("Reading netCDF file successfully!\n");
+
+ for (size_t i = 0; i < len_vals; i++)
+ points_vec->push_back(new Point( lon_in[i], lat_in[i], var_in[i] ));
+
+ delete [] lat_in;
+ delete [] lon_in;
+ delete [] var_in;
+
+ obj->addPointVec(points_vec, ifname);
+}
+
+CFEMesh* NetCDFInterface::createMeshFromPoints(std::vector<GEOLIB::Point*>* points_vec,
+ size_t &NRLAT,
+ size_t &NRLON)
+{
+ std::cout << "Converting points data to mesh with quad elements\n";
+ CFEMesh* mesh = new CFEMesh();
+
+ // setmesh nodes
+ size_t nNodes = points_vec->size();
+ for (size_t i = 0; i < nNodes; i++)
+ {
+ MeshLib::CNode* node = new MeshLib::CNode(i);
+ double coords[3] =
+ {(*(*points_vec)[i])[0], (*(*points_vec)[i])[1], (*(*points_vec)[i])[2]};
+ node->SetCoordinates(coords);
+ mesh->nod_vector.push_back(node);
+ }
+
+ // set mesh elements
+ for (size_t i = 0; i < NRLAT - 1; i++)
+ for (size_t j = 0; j < NRLON - 1; j++)
+ {
+ size_t n_Elems = i * NRLON + j;
+ // TF 11/2011 - using instead the constructor
+// MeshLib::CElem* elem = new MeshLib::CElem();
+// elem->setElementProperties(MshElemType::QUAD);
+// // Assignment for the element nodes
+// elem->nodes_index[0] = (long) (n_Elems);
+// elem->nodes_index[1] = (long) (n_Elems + 1);
+// elem->nodes_index[2] = (long) (n_Elems + NRLON + 1);
+// elem->nodes_index[3] = (long) (n_Elems + NRLON);
+// // Initialize topological properties
+// elem->InitializeMembers();
+ MeshLib::CElem* elem = new MeshLib::CElem(MshElemType::QUAD, n_Elems, n_Elems + 1, n_Elems + NRLON + 1, n_Elems + NRLON, 0);
+
+ mesh->ele_vector.push_back(elem);
+ }
+ // Establish topology of a grid
+ mesh->ConstructGrid();
+ std::cout << "Mesh is generated successfully!\n" << std::endl;
+ return mesh;
+}
+} // end namespace FileIO
diff --git a/FileIO/NetCDFInterface.h b/FileIO/NetCDFInterface.h
new file mode 100644
index 00000000000..dfc3793bc5f
--- /dev/null
+++ b/FileIO/NetCDFInterface.h
@@ -0,0 +1,43 @@
+/**
+ * \file NetCDFInterface.h
+ * 29/07/2010 YW Initial implementation
+ */
+
+#ifndef NetCDFInterface_H
+#define NetCDFInterface_H
+
+#include "GEOObjects.h"
+#include <string>
+#include <vector>
+
+namespace GEOLIB
+{
+class GEOObjects;
+}
+
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace FileIO
+{
+class NetCDFInterface
+{
+public:
+ /// Import climate data from a NetCDF file.
+ static void readNetCDFData(std::string &fname,
+ std::vector<GeoLib::Point*>* points_vec,
+ GEOLIB::GEOObjects* geo_obj,
+ size_t &NRLAT,
+ size_t &NRLON);
+ /// Convert imported point data to an CFEMesh.
+ static MeshLib::Mesh* createMeshFromPoints(std::vector<GeoLib::Point*>* points_vec,
+ size_t &NRLAT,
+ size_t &NRLON);
+
+private:
+};
+} // end namespace
+
+#endif /* NetCDFInterface_H */
diff --git a/FileIO/PetrelInterface.cpp b/FileIO/PetrelInterface.cpp
new file mode 100644
index 00000000000..63771df62ae
--- /dev/null
+++ b/FileIO/PetrelInterface.cpp
@@ -0,0 +1,234 @@
+/*
+ * PetrelInterface.cpp
+ *
+ * Created on: Feb 16, 2010
+ * Author: fischeth
+ */
+
+#include "PetrelInterface.h"
+#include "StringTools.h"
+#include <fstream>
+
+namespace FileIO
+{
+PetrelInterface::PetrelInterface(std::list<std::string> &sfc_fnames,
+ std::list<std::string> &well_path_fnames,
+ std::string &unique_model_name, GEOLIB::GEOObjects* geo_obj) :
+ _unique_name(unique_model_name), pnt_vec(new std::vector<GEOLIB::Point*>),
+ well_vec(new std::vector<GEOLIB::Point*>), ply_vec(new std::vector<
+ GEOLIB::Polyline*>)
+{
+ for (std::list<std::string>::const_iterator it(sfc_fnames.begin()); it
+ != sfc_fnames.end(); ++it)
+ {
+ std::cout << "PetrelInterface::PetrelInterface open surface stream from file " <<
+ *it
+ << " ... " << std::flush;
+ std::ifstream in((*it).c_str());
+ if (in)
+ {
+ std::cout << "done" << std::endl;
+ readPetrelSurface(in);
+ in.close();
+ }
+ else
+ std::cerr << "error opening stream " << std::endl;
+ }
+
+ for (std::list<std::string>::const_iterator it(well_path_fnames.begin()); it
+ != well_path_fnames.end(); ++it)
+ {
+ std::cout << "PetrelInterface::PetrelInterface open well path stream from file "
+ << *it << " ... " << std::flush;
+ std::ifstream in((*it).c_str());
+ if (in)
+ {
+ std::cout << "done" << std::endl;
+ readPetrelWellTrace(in);
+ in.close();
+ }
+ else
+ std::cerr << "error opening stream " << std::endl;
+ }
+
+ // store data in GEOObject
+ geo_obj->addPointVec(pnt_vec, _unique_name);
+ if (well_vec->size() > 0)
+ geo_obj->addStationVec(
+ well_vec,
+ _unique_name);
+ if (ply_vec->size () > 0)
+ geo_obj->addPolylineVec(ply_vec, _unique_name);
+}
+
+PetrelInterface::~PetrelInterface()
+{}
+
+void PetrelInterface::readPetrelSurface(std::istream &in)
+{
+ char buffer[MAX_COLS_PER_ROW];
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ std::string line(buffer);
+
+ if (line.find("# Petrel Points with attributes") != std::string::npos)
+ {
+ // read header
+ // read Version string
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ line = buffer;
+ // read string BEGIN HEADER
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ line = buffer;
+
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ line = buffer;
+ while (line.find("END HEADER") == std::string::npos)
+ {
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ line = buffer;
+ }
+
+ // read points
+ size_t idx(pnt_vec->size());
+ while (in)
+ {
+ pnt_vec->push_back(new GEOLIB::Point);
+ in >> *((*pnt_vec)[idx]);
+ if (!in)
+ {
+ delete (*pnt_vec)[idx];
+ pnt_vec->pop_back();
+ }
+ else
+ idx++;
+ }
+ }
+ else
+ std::cerr << "error reading petrel points: " << line << std::endl;
+}
+
+void PetrelInterface::readPetrelWellTrace(std::istream &in)
+{
+ char buffer[MAX_COLS_PER_ROW];
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ std::string line(buffer);
+
+ if (line.find("# WELL TRACE FROM PETREL") != std::string::npos)
+ {
+ // read header
+ // read well name
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ line = buffer;
+ std::list<std::string> str_list(splitString(line, ' '));
+ std::list<std::string>::const_iterator it(str_list.begin());
+ while (it != str_list.end())
+ std::cout << *it++ << " " << std::flush;
+ std::cout << std::endl;
+
+ // read well head x coordinate
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ line = buffer;
+ str_list = splitString(line, ' ');
+ it = str_list.begin();
+ while (it != str_list.end())
+ std::cout << *it++ << " " << std::flush;
+ std::cout << std::endl;
+ it = (str_list.end())--;
+ it--;
+ char* buf;
+ double well_head_x(strtod((*it).c_str(), &buf));
+
+ // read well head y coordinate
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ line = buffer;
+ str_list = splitString(line, ' ');
+ it = str_list.begin();
+ while (it != str_list.end())
+ std::cout << *it++ << " " << std::flush;
+ std::cout << std::endl;
+ it = (str_list.end())--;
+ it--;
+ double well_head_y(strtod((*it).c_str(), &buf));
+
+ // read well KB
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ line = buffer;
+ str_list = splitString(line, ' ');
+ it = str_list.begin();
+ while (it != str_list.end())
+ std::cout << *it++ << " " << std::flush;
+ std::cout << std::endl;
+ it = (str_list.end())--;
+ it--;
+ double well_kb(strtod((*it).c_str(), &buf));
+
+ std::cout << "PetrelInterface::readPetrelWellTrace: " << well_head_x << "," <<
+ well_head_y << "," << well_kb << std::endl;
+ well_vec->push_back(
+ static_cast<GEOLIB::StationBorehole*> (new GEOLIB::StationBorehole(
+ well_head_x, well_head_y,
+ well_kb)));
+
+ // read well type
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ line = buffer;
+ str_list = splitString(line, ' ');
+ it = str_list.begin();
+ while (it != str_list.end())
+ std::cout << *it++ << " " << std::flush;
+ std::cout << std::endl;
+ std::string type(*((str_list.end())--));
+
+ readPetrelWellTraceData(in);
+ }
+}
+
+void PetrelInterface::readPetrelWellTraceData(std::istream &in)
+{
+ char buffer[MAX_COLS_PER_ROW];
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ std::string line(buffer);
+ std::list<std::string> str_list;
+ std::list<std::string>::const_iterator it;
+
+ // read yet another header lines
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ line = buffer;
+ while (line.substr(0, 1).compare("#") == 0)
+ {
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ line = buffer;
+ }
+
+ // read column information
+ str_list = splitString(line, ' ');
+ it = str_list.begin();
+ while (it != str_list.end())
+ std::cout << *it++ << " " << std::flush;
+ std::cout << std::endl;
+
+ // read points
+ double md, x, y, z, tvd, dx, dy, azim, incl, dls;
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ line = buffer;
+ while (in)
+ {
+ if (line.size() > 1 && line.substr(0, 1).compare("#") != 0)
+ {
+ std::stringstream stream(line);
+ stream >> md;
+ stream >> x >> y >> z;
+ // pnt_vec->push_back (new GEOLIB::Point (x,y,z));
+ static_cast<GEOLIB::StationBorehole*> ((*well_vec)[well_vec->size()
+ - 1])->addSoilLayer(
+ x,
+ y,
+ z,
+ "unknown");
+ stream >> tvd >> dx >> dy >> azim >> incl >> dls;
+ }
+ in.getline(buffer, MAX_COLS_PER_ROW);
+ line = buffer;
+ }
+}
+} // end namespace FileIO
diff --git a/FileIO/PetrelInterface.h b/FileIO/PetrelInterface.h
new file mode 100644
index 00000000000..fa440c64ccb
--- /dev/null
+++ b/FileIO/PetrelInterface.h
@@ -0,0 +1,40 @@
+/*
+ * PetrelInterface.h
+ *
+ * Created on: Feb 16, 2010
+ * Author: fischeth
+ */
+
+#ifndef PETRELIO_H_
+#define PETRELIO_H_
+
+#include "GEOObjects.h"
+#include <iostream>
+#include <list>
+#include <string>
+#include <vector>
+
+namespace FileIO
+{
+class PetrelInterface
+{
+public:
+ PetrelInterface(std::list<std::string> &sfc_fnames,
+ std::list<std::string> &well_path_fnames,
+ std::string &unique_model_name,
+ GeoLib::GEOObjects* obj);
+ virtual ~PetrelInterface();
+
+private:
+ void readPetrelSurface (std::istream &in);
+ void readPetrelWellTrace (std::istream &in);
+ void readPetrelWellTraceData (std::istream &in);
+ std::string _unique_name;
+ std::vector<GeoLib::Point*>* pnt_vec;
+ std::vector<GeoLib::Point*>* well_vec;
+ std::vector<GeoLib::Polyline*>* ply_vec;
+ static const size_t MAX_COLS_PER_ROW = 256;
+};
+} // end namespace FileIO
+
+#endif /* PETRELIO_H_ */
diff --git a/Gui/DataView/CondFromRasterDialog.cpp b/Gui/DataView/CondFromRasterDialog.cpp
index 94247f2a81a..b4462352cca 100644
--- a/Gui/DataView/CondFromRasterDialog.cpp
+++ b/Gui/DataView/CondFromRasterDialog.cpp
@@ -13,8 +13,8 @@
#include "OGSError.h"
#include "StrictDoubleValidator.h"
-CondFromRasterDialog::CondFromRasterDialog(const std::map<std::string, MeshLib::Mesh*> &msh_map, QDialog* parent)
- : QDialog(parent), _msh_map(msh_map)
+CondFromRasterDialog::CondFromRasterDialog(const std::vector<MeshLib::Mesh*> &msh_vec, QDialog* parent)
+ : QDialog(parent), _msh_vec(msh_vec)
{
setupUi(this);
@@ -23,9 +23,8 @@ CondFromRasterDialog::CondFromRasterDialog(const std::map<std::string, MeshLib::
this->scalingEdit->setText("1.0");
this->scalingEdit->setValidator (_scale_validator);
- for (std::map<std::string, MeshLib::Mesh*>::const_iterator it = _msh_map.begin();
- it != _msh_map.end(); ++it)
- this->meshBox->addItem(QString::fromStdString(it->first));
+ for (std::vector<MeshLib::Mesh*>::const_iterator it = _msh_vec.begin(); it != _msh_vec.end(); ++it)
+ this->meshBox->addItem(QString::fromStdString((*it)->getName()));
this->directButton->setChecked(true);
}
@@ -75,9 +74,15 @@ void CondFromRasterDialog::accept()
return;
}
- const MeshLib::Mesh* mesh = (_msh_map.find(mesh_name))->second;
- //std::string direct_node_name(raster_name + ".txt");
+ MeshLib::Mesh* mesh (NULL);
+ for (size_t i=0; i<_msh_vec.size(); i++)
+ if (_msh_vec[i]->getName().compare(mesh_name) == 0)
+ {
+ mesh = _msh_vec[i];
+ break;
+ }
+
if (this->directButton->isChecked())
{
DirectConditionGenerator dcg;
@@ -93,7 +98,7 @@ void CondFromRasterDialog::accept()
}
MeshLib::Mesh* new_mesh = const_cast<MeshLib::Mesh*>(mesh);
DirectConditionGenerator dcg;
- //TODO6: direct_values = dcg.directWithSurfaceIntegration(*new_mesh, raster_name, scaling_factor);
+ direct_values = dcg.directWithSurfaceIntegration(*new_mesh, raster_name, scaling_factor);
//dcg.writeToFile(direct_node_name);
}
diff --git a/Gui/DataView/CondFromRasterDialog.h b/Gui/DataView/CondFromRasterDialog.h
index 571670ed0cf..15e2e9e9dcd 100644
--- a/Gui/DataView/CondFromRasterDialog.h
+++ b/Gui/DataView/CondFromRasterDialog.h
@@ -25,11 +25,11 @@ class CondFromRasterDialog : public QDialog, private Ui_CondFromRaster
Q_OBJECT
public:
- CondFromRasterDialog(const std::map<std::string, MeshLib::Mesh*> &msh_map, QDialog* parent = 0);
+ CondFromRasterDialog(const std::vector<MeshLib::Mesh*> &msh_vec, QDialog* parent = 0);
~CondFromRasterDialog(void);
private:
- const std::map<std::string, MeshLib::Mesh*> _msh_map;
+ const std::vector<MeshLib::Mesh*> _msh_vec;
StrictDoubleValidator* _scale_validator;
private slots:
diff --git a/Gui/DataView/DirectConditionGenerator.h b/Gui/DataView/DirectConditionGenerator.h
index 249268949b8..c347af18440 100644
--- a/Gui/DataView/DirectConditionGenerator.h
+++ b/Gui/DataView/DirectConditionGenerator.h
@@ -21,8 +21,7 @@ public:
const std::vector< std::pair<size_t,double> >& directToSurfaceNodes(const MeshLib::Mesh &mesh, const std::string &filename);
- //TODO6
- //const std::vector< std::pair<size_t,double> >& directWithSurfaceIntegration(MeshLib::Mesh &mesh, const std::string &filename, double scaling);
+ const std::vector< std::pair<size_t,double> >& directWithSurfaceIntegration(MeshLib::Mesh &mesh, const std::string &filename, double scaling);
int writeToFile(const std::string &name) const;
diff --git a/Gui/DataView/FEMConditionSetupDialog.cpp b/Gui/DataView/FEMConditionSetupDialog.cpp
index a9d96b526b6..38f58d9650b 100644
--- a/Gui/DataView/FEMConditionSetupDialog.cpp
+++ b/Gui/DataView/FEMConditionSetupDialog.cpp
@@ -42,7 +42,7 @@ FEMConditionSetupDialog::FEMConditionSetupDialog(const FEMCondition &cond, QDial
setValuesFromCond();
}
-FEMConditionSetupDialog::FEMConditionSetupDialog(const std::string &name, const MeshLib::CFEMesh* mesh, QDialog* parent)
+FEMConditionSetupDialog::FEMConditionSetupDialog(const std::string &name, const MeshLib::Mesh* mesh, QDialog* parent)
: QDialog(parent), _cond(name, FEMCondition::UNSPECIFIED), _set_on_points(false), _combobox(NULL), directButton(NULL),
_mesh(mesh), _first_value_validator(NULL)
{
@@ -255,9 +255,9 @@ void FEMConditionSetupDialog::directButton_pressed()
}
else
{
- std::map<std::string, MeshLib::CFEMesh*> msh_map;
- msh_map.insert( std::pair<std::string, MeshLib::CFEMesh*>(this->_cond.getGeoName(), const_cast<MeshLib::CFEMesh*>(this->_mesh)) );
- CondFromRasterDialog dlg(msh_map);
+ std::vector<MeshLib::Mesh*> msh_vec;
+ msh_vec.push_back( const_cast<MeshLib::Mesh*>(this->_mesh) );
+ CondFromRasterDialog dlg(msh_vec);
//connect(&dlg, SIGNAL(directNodesWritten(std::string)), this, SLOT(direct_path_changed(std::string)));
connect(&dlg, SIGNAL(transmitDisValues(std::vector< std::pair<size_t,double> >)),
this, SLOT(addDisValues(std::vector< std::pair<size_t,double> >)));
diff --git a/Gui/DataView/FEMConditionSetupDialog.h b/Gui/DataView/FEMConditionSetupDialog.h
index 89e6e2bf277..b7d68d0e0cb 100644
--- a/Gui/DataView/FEMConditionSetupDialog.h
+++ b/Gui/DataView/FEMConditionSetupDialog.h
@@ -20,7 +20,7 @@ namespace GeoLib {
}
namespace MeshLib {
- class CFEMesh;
+ class Mesh;
}
/**
@@ -44,7 +44,7 @@ public:
FEMConditionSetupDialog(const FEMCondition &cond, QDialog* parent = 0);
/// Constructor for creating DIRECT FEM conditions on MeshNodes.
- FEMConditionSetupDialog(const std::string &name, const MeshLib::CFEMesh* mesh, QDialog* parent = 0);
+ FEMConditionSetupDialog(const std::string &name, const MeshLib::Mesh* mesh, QDialog* parent = 0);
~FEMConditionSetupDialog(void);
@@ -58,7 +58,7 @@ private:
bool _set_on_points;
QComboBox* _combobox; //needed for on_points & linear conds
QPushButton* directButton; // needed for direct conditions
- const MeshLib::CFEMesh* _mesh; // needed for direct conditions
+ const MeshLib::Mesh* _mesh; // needed for direct conditions
StrictDoubleValidator* _first_value_validator;
private slots:
diff --git a/Gui/mainwindow.cpp b/Gui/mainwindow.cpp
index 3345aff50af..0835ee32d85 100644
--- a/Gui/mainwindow.cpp
+++ b/Gui/mainwindow.cpp
@@ -47,22 +47,18 @@
#include "BoundaryCondition.h"
#include "InitialCondition.h"
#include "SourceTerm.h"
-#include "rf_bc_new.h"
-#include "rf_ic_new.h"
-#include "rf_st_new.h"
-#include "FEMIO/BoundaryConditionIO.h"
// FileIO includes
-#include "FEFLOWInterface.h"
+// TODO6 #include "FEFLOWInterface.h"
#include "GMSInterface.h"
-#include "GeoIO/Gmsh2GeoIO.h"
-#include "GocadInterface.h"
+#include "Gmsh2GeoIO.h"
+// TODO6 #include "GocadInterface.h"
#include "MeshIO/GMSHInterface.h"
#include "MeshIO/TetGenInterface.h"
#include "NetCDFInterface.h"
-#include "OGSIOVer4.h"
+// TODO6 #include "OGSIOVer4.h"
#include "PetrelInterface.h"
-#include "StationIO.h"
+// TODO6 #include "StationIO.h"
#include "XmlIO/XmlCndInterface.h"
#include "XmlIO/XmlGmlInterface.h"
#include "XmlIO/XmlGspInterface.h"
diff --git a/scripts/cmake/Find.cmake b/scripts/cmake/Find.cmake
index 090cd9ebcbb..cb85ac80982 100644
--- a/scripts/cmake/Find.cmake
+++ b/scripts/cmake/Find.cmake
@@ -89,7 +89,7 @@ IF(VTK_NETCDF_FOUND)
ELSE()
SET(NETCDF_CXX TRUE)
FIND_PACKAGE(NetCDF)
- IF(NOT NETCDF_FOUND AND OGS_USE_QT)
+ IF(NOT NETCDF_FOUND AND OGS_BUILD_GUI)
MESSAGE(FATAL_ERROR "NetCDF was not found but is required for OGS_USE_QT!")
ENDIF()
ENDIF()
--
GitLab