From fb03f6ab2a3a675c1ab0c8562bdf6f967a2fe2e2 Mon Sep 17 00:00:00 2001
From: Karsten Rink <karsten.rink@ufz.de>
Date: Wed, 29 Oct 2014 15:53:24 +0100
Subject: [PATCH] subdivided gmsh reading method a bit and added handling of
 physical names tag

---
 FileIO/GMSHInterface.cpp | 216 +++++++++++++++++++++------------------
 FileIO/GMSHInterface.h   |   7 +-
 2 files changed, 122 insertions(+), 101 deletions(-)

diff --git a/FileIO/GMSHInterface.cpp b/FileIO/GMSHInterface.cpp
index 912351f5c81..54e98f4c56c 100644
--- a/FileIO/GMSHInterface.cpp
+++ b/FileIO/GMSHInterface.cpp
@@ -107,126 +107,67 @@ MeshLib::Mesh* GMSHInterface::readGMSHMesh(std::string const& fname)
 	std::vector<MeshLib::Node*> nodes;
 	std::vector<MeshLib::Element*> elements;
 
-	if (line.find("$MeshFormat") != std::string::npos)
+	if (line.find("$MeshFormat") == std::string::npos)
 	{
-		getline(in, line); // version-number file-type data-size
-		if (line.substr(0,3).compare("2.2") != 0) {
-			WARN("Wrong gmsh file format version.");
-			return nullptr;
-		}
-		if (line.substr(4,1).compare("0") != 0) {
-			WARN("Currently reading gmsh binary file type is not supported.");
-			return nullptr;
-		}
-		getline(in, line); //$EndMeshFormat
-		getline(in, line); //$Nodes Keywords
+		in.close();
+		WARN ("No GMSH file format recognized.");
+		return nullptr;
+	}
+
+	getline(in, line); // version-number file-type data-size
+	if (line.substr(0,3).compare("2.2") != 0) {
+		WARN("Wrong gmsh file format version.");
+		return nullptr;
+	}
+
+	if (line.substr(4,1).compare("0") != 0) {
+		WARN("Currently reading gmsh binary file type is not supported.");
+		return nullptr;
+	}
+	getline(in, line); //$EndMeshFormat
 
-		size_t n_nodes(0);
-		size_t n_elements(0);
-		while (line.find("$EndElements") == std::string::npos)
+	std::map<unsigned, unsigned> id_map;
+	while (line.find("$EndElements") == std::string::npos)
+	{
+		// Node data
+		getline(in, line); //$Nodes Keywords
+		if (line.find("$Nodes") != std::string::npos)
 		{
-			// Node data
+			size_t n_nodes(0);
 			long id;
 			double x, y, z;
 			in >> n_nodes >> std::ws;
 			nodes.resize(n_nodes);
-			std::map<unsigned, unsigned> id_map;
 			for (size_t i = 0; i < n_nodes; i++) {
 				in >> id >> x >> y >> z >> std::ws;
 				id_map.insert(std::map<unsigned, unsigned>::value_type(id, i));
 				nodes[i] = new MeshLib::Node(x,y,z,id);
 			}
 			getline(in, line); // End Node keyword $EndNodes
+		}
 
-			// Element data
-			getline(in, line); // Element keyword $Elements
+		// Element data
+		if (line.find("$Elements") != std::string::npos)
+		{
+			size_t n_elements(0);
 			in >> n_elements >> std::ws; // number-of-elements
 			elements.reserve(n_elements);
-			unsigned idx, type, n_tags, dummy, mat_id;
 			for (size_t i = 0; i < n_elements; i++)
 			{
-				MeshLib::Element* elem (NULL);
-				std::vector<unsigned> node_ids;
-				std::vector<MeshLib::Node*> elem_nodes;
-				in >> idx >> type >> n_tags >> dummy >> mat_id;
-
-				// skip tags
-				for (size_t j = 2; j < n_tags; j++)
-					in >> dummy;
-
-				switch (type)
-				{
-				case 1: {
-					readNodeIDs(in, 2, node_ids, id_map);
-					// edge_nodes array will be deleted from Line object
-					MeshLib::Node** edge_nodes = new MeshLib::Node*[2];
-					edge_nodes[0] = nodes[node_ids[0]];
-					edge_nodes[1] = nodes[node_ids[1]];
-					elem = new MeshLib::Line(edge_nodes, 0);
-					break;
-				}
-				case 2: {
-					readNodeIDs(in, 3, node_ids, id_map);
-					MeshLib::Node** tri_nodes = new MeshLib::Node*[3];
-					tri_nodes[0] = nodes[node_ids[2]];
-					tri_nodes[1] = nodes[node_ids[1]];
-					tri_nodes[2] = nodes[node_ids[0]];
-					elem = new MeshLib::Tri(tri_nodes, mat_id);
-					break;
-				}
-				case 3: {
-					readNodeIDs(in, 4, node_ids, id_map);
-					MeshLib::Node** quad_nodes = new MeshLib::Node*[4];
-					for (unsigned k(0); k < 4; k++)
-						quad_nodes[k] = nodes[node_ids[k]];
-					elem = new MeshLib::Quad(quad_nodes, mat_id);
-					break;
-				}
-				case 4: {
-					readNodeIDs(in, 4, node_ids, id_map);
-					MeshLib::Node** tet_nodes = new MeshLib::Node*[5];
-					for (unsigned k(0); k < 4; k++)
-						tet_nodes[k] = nodes[node_ids[k]];
-					elem = new MeshLib::Tet(tet_nodes, mat_id);
-					break;
-				}
-				case 5: {
-					readNodeIDs(in, 8, node_ids, id_map);
-					MeshLib::Node** hex_nodes = new MeshLib::Node*[8];
-					for (unsigned k(0); k < 8; k++)
-						hex_nodes[k] = nodes[node_ids[k]];
-					elem = new MeshLib::Hex(hex_nodes, mat_id);
-					break;
-				}
-				case 6: {
-					readNodeIDs(in, 6, node_ids, id_map);
-					MeshLib::Node** prism_nodes = new MeshLib::Node*[6];
-					for (unsigned k(0); k < 6; k++)
-						prism_nodes[k] = nodes[node_ids[k]];
-					elem = new MeshLib::Prism(prism_nodes, mat_id);
-					break;
-				}
-				case 7: {
-					readNodeIDs(in, 5, node_ids, id_map);
-					MeshLib::Node** pyramid_nodes = new MeshLib::Node*[5];
-					for (unsigned k(0); k < 5; k++)
-						pyramid_nodes[k] = nodes[node_ids[k]];
-					elem = new MeshLib::Pyramid(pyramid_nodes, mat_id);
-					break;
-				}
-				case 15:
-					in >> dummy; // skip rest of line
-					continue;
-					break;
-				default:
-						WARN("GMSHInterface::readGMSHMesh(): Unknown element type %d.", type);
-				}
-				in >> std::ws;
+				MeshLib::Element* elem (readElement(in, nodes, id_map));
 
-				if (type > 0 && type < 8)
+				if (elem)
 					elements.push_back(elem);
 			}
+			getline(in, line); // END keyword
+		}
 
+		if (line.find("PhysicalNames") != std::string::npos)
+		{
+			size_t n_lines(0);
+			in >> n_lines >> std::ws; // number-of-lines
+			for (size_t i = 0; i < n_lines; i++)
+				getline(in, line);
 			getline(in, line); // END keyword
 		}
 	}
@@ -243,16 +184,91 @@ MeshLib::Mesh* GMSHInterface::readGMSHMesh(std::string const& fname)
 void GMSHInterface::readNodeIDs(std::ifstream &in,
                                 unsigned n_nodes,
                                 std::vector<unsigned> &node_ids,
-                                std::map<unsigned, unsigned> &id_map)
+                                std::map<unsigned, unsigned> const& id_map)
 {
 	unsigned idx;
 	for (unsigned i = 0; i < n_nodes; i++)
 	{
 		in >> idx;
-		node_ids.push_back(id_map[idx]);
+		node_ids.push_back(id_map.at(idx));
+	}
+}
+
+MeshLib::Element* GMSHInterface::readElement(std::ifstream &in, std::vector<MeshLib::Node*> const& nodes, std::map<unsigned, unsigned> const& id_map)
+{	
+	unsigned idx, type, n_tags, dummy, mat_id;
+	MeshLib::Element* elem (nullptr);
+	std::vector<unsigned> node_ids;
+	std::vector<MeshLib::Node*> elem_nodes;
+	in >> idx >> type >> n_tags >> dummy >> mat_id;
+
+	// skip tags
+	for (size_t j = 2; j < n_tags; j++)
+		in >> dummy;
+
+	switch (type)
+	{
+	case 1: {
+		readNodeIDs(in, 2, node_ids, id_map);
+		// edge_nodes array will be deleted from Line object
+		MeshLib::Node** edge_nodes = new MeshLib::Node*[2];
+		edge_nodes[0] = nodes[node_ids[0]];
+		edge_nodes[1] = nodes[node_ids[1]];
+		return new MeshLib::Line(edge_nodes, 0);
+	}
+	case 2: {
+		readNodeIDs(in, 3, node_ids, id_map);
+		MeshLib::Node** tri_nodes = new MeshLib::Node*[3];
+		tri_nodes[0] = nodes[node_ids[2]];
+		tri_nodes[1] = nodes[node_ids[1]];
+		tri_nodes[2] = nodes[node_ids[0]];
+		return new MeshLib::Tri(tri_nodes, mat_id);
+	}
+	case 3: {
+		readNodeIDs(in, 4, node_ids, id_map);
+		MeshLib::Node** quad_nodes = new MeshLib::Node*[4];
+		for (unsigned k(0); k < 4; k++)
+			quad_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Quad(quad_nodes, mat_id);
+	}
+	case 4: {
+		readNodeIDs(in, 4, node_ids, id_map);
+		MeshLib::Node** tet_nodes = new MeshLib::Node*[5];
+		for (unsigned k(0); k < 4; k++)
+			tet_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Tet(tet_nodes, mat_id);
+	}
+	case 5: {
+		readNodeIDs(in, 8, node_ids, id_map);
+		MeshLib::Node** hex_nodes = new MeshLib::Node*[8];
+		for (unsigned k(0); k < 8; k++)
+			hex_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Hex(hex_nodes, mat_id);
 	}
+	case 6: {
+		readNodeIDs(in, 6, node_ids, id_map);
+		MeshLib::Node** prism_nodes = new MeshLib::Node*[6];
+		for (unsigned k(0); k < 6; k++)
+			prism_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Prism(prism_nodes, mat_id);
+	}
+	case 7: {
+		readNodeIDs(in, 5, node_ids, id_map);
+		MeshLib::Node** pyramid_nodes = new MeshLib::Node*[5];
+		for (unsigned k(0); k < 5; k++)
+			pyramid_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Pyramid(pyramid_nodes, mat_id);
+	}
+	case 15:
+		in >> dummy; // skip rest of line
+		break;
+	default:
+		WARN("GMSHInterface::readGMSHMesh(): Unknown element type %d.", type);
+	}
+	return nullptr;
 }
 
+
 bool GMSHInterface::write()
 {
 	_out << "// GMSH input file created by OpenGeoSys " << BaseLib::BuildInfo::ogs_version_and_persons;
diff --git a/FileIO/GMSHInterface.h b/FileIO/GMSHInterface.h
index 880ff631787..ddaa0aff0bf 100644
--- a/FileIO/GMSHInterface.h
+++ b/FileIO/GMSHInterface.h
@@ -30,6 +30,8 @@
 
 namespace MeshLib {
 	class Mesh;
+	class Element;
+	class Node;
 }
 
 namespace FileIO
@@ -87,6 +89,9 @@ protected:
 	bool write();
 
 private:
+	/// Reads a mesh element from the input stream
+	static MeshLib::Element* readElement(std::ifstream &in, std::vector<MeshLib::Node*> const& nodes, std::map<unsigned, unsigned> const& id_map);
+
 	/**
 	 * 1. get and merge data from _geo_objs
 	 * 2. compute topological hierarchy
@@ -94,7 +99,7 @@ private:
 	 */
 	void writeGMSHInputFile(std::ostream & out);
 
-	static void readNodeIDs(std::ifstream &in, unsigned n_nodes, std::vector<unsigned> &node_ids, std::map<unsigned, unsigned> &id_map);
+	static void readNodeIDs(std::ifstream &in, unsigned n_nodes, std::vector<unsigned> &node_ids, std::map<unsigned, unsigned> const& id_map);
 
 	void writePoints(std::ostream& out) const;
 
-- 
GitLab