current state of integration for LB (fixing netcdf include error)

FemLib/BoundaryCondition.cpp

+/**
+ * \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());
+	}
+}

FemLib/BoundaryCondition.h

+/**
+ * \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

FemLib/InitialCondition.cpp

+/**
+ * \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());
+}

FemLib/InitialCondition.h

+/**
+ * \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

FemLib/SourceTerm.cpp

+/**
+ * \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();
+	}
+}

FemLib/SourceTerm.h

+/**
+ * \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

FileIO/GMSInterface.cpp

+/**
+ * \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);
+}

FileIO/GMSInterface.h

+/**
+ * \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_ */

FileIO/Gmsh2GeoIO.cpp

+/*
+ * 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

FileIO/Gmsh2GeoIO.h

+/*
+ * 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_ */

FileIO/MeshIO/GMSHAdaptiveMeshDensity.cpp

+/*
+ * 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

FileIO/MeshIO/GMSHAdaptiveMeshDensity.h

+/*
+ * 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_ */

FileIO/MeshIO/GMSHFixedMeshDensity.cpp

+/*
+ * 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

FileIO/MeshIO/GMSHFixedMeshDensity.h

+/*
+ * 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_ */

FileIO/MeshIO/GMSHInterface.cpp

+/*
+ * 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

FileIO/MeshIO/GMSHInterface.h

+/*
+ * 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_ */

FileIO/MeshIO/GMSHLine.cpp

+/*
+ * 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

FileIO/MeshIO/GMSHLine.h

+/*
+ * 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_ */

FileIO/MeshIO/GMSHLineLoop.cpp

+/*
+ * 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

FileIO/MeshIO/GMSHLineLoop.h

+/*
+ * 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_ */