From 1d2c4a69066d27fda2c6c2b782be79c86f1e0c31 Mon Sep 17 00:00:00 2001
From: Norihiro Watanabe <norihiro.watanabe@ufz.de>
Date: Tue, 20 May 2014 18:35:14 +0200
Subject: [PATCH] add FEFLOEInterface class in FileIO
---
FileIO/CMakeLists.txt | 1 +
FileIO/FEFLOWInterface.cpp | 406 +++++++++++++++++++++++++++++++++++++
FileIO/FEFLOWInterface.h | 145 +++++++++++++
3 files changed, 552 insertions(+)
create mode 100644 FileIO/FEFLOWInterface.cpp
create mode 100644 FileIO/FEFLOWInterface.h
diff --git a/FileIO/CMakeLists.txt b/FileIO/CMakeLists.txt
index ee5e300aac9..00b4bac1400 100644
--- a/FileIO/CMakeLists.txt
+++ b/FileIO/CMakeLists.txt
@@ -27,6 +27,7 @@ SET (SOURCES ${SOURCES} ${SOURCES_BASE_XML})
GET_SOURCE_FILES(SOURCES_BOOST_XML XmlIO/Boost)
SET (SOURCES ${SOURCES} ${SOURCES_BOOST_XML})
IF (QT4_FOUND)
+ SET ( SOURCES ${SOURCES} FEFLOWInterface.h FEFLOWInterface.cpp)
GET_SOURCE_FILES(SOURCES_QT_XML XmlIO/Qt)
SET ( SOURCES ${SOURCES} ${SOURCES_QT_XML})
ENDIF (QT4_FOUND)
diff --git a/FileIO/FEFLOWInterface.cpp b/FileIO/FEFLOWInterface.cpp
new file mode 100644
index 00000000000..02de780c226
--- /dev/null
+++ b/FileIO/FEFLOWInterface.cpp
@@ -0,0 +1,406 @@
+/**
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "FEFLOWInterface.h"
+
+#include "logog/include/logog.hpp"
+#include <QtXml>
+
+#include "StringTools.h"
+#include "FileTools.h"
+#include "Polygon.h"
+#include "Mesh.h"
+#include "Node.h"
+#include "Elements/Line.h"
+#include "Elements/Hex.h"
+#include "Elements/Prism.h"
+#include "Elements/Quad.h"
+#include "Elements/Tet.h"
+#include "Elements/Tri.h"
+
+namespace FileIO
+{
+
+MeshLib::Mesh* FEFLOWInterface::readFEFLOWFile(const std::string &filename)
+{
+ std::ifstream in(filename.c_str());
+ if (!in)
+ {
+ ERR("FEFLOWInterface::readFEFLOWFile(): Could not open file %s.", filename.c_str());
+ return nullptr;
+ }
+
+ FEM_CLASS fem_class;
+ FEM_DIM fem_dim;
+ std::vector<GeoLib::Point*>* points = NULL;
+ std::vector<GeoLib::Polyline*>* lines = NULL;
+
+ bool isXZplane = false;
+
+ std::vector<MeshLib::Node*> vec_nodes;
+ std::vector<MeshLib::Element*> vec_elements;
+
+ std::string line_string;
+ std::stringstream line_stream;
+ while (!in.eof())
+ {
+ getline(in, line_string);
+ //....................................................................
+ // CLASS
+ if (line_string.find("CLASS") != std::string::npos) // the version number follows afterward, e.g. CLASS (v.5.313)
+ {
+ getline(in, line_string);
+ line_stream.str(line_string);
+ // problem class, time mode, problem orientation, dimension, nr. layers for 3D, saturation switch, precision of results, precision of coordinates
+ line_stream >> fem_class.problem_class >> fem_class.time_mode >> fem_class.orientation >> fem_class.dimension >> fem_class.n_layers3d;
+ line_stream.clear();
+ }
+ //....................................................................
+ // DIMENS
+ else if (line_string.compare("DIMENS") == 0)
+ {
+ // DIMENS
+ getline(in, line_string);
+ line_stream.str(line_string);
+ line_stream >> fem_dim.n_nodes >> fem_dim.n_elements >> fem_dim.n_nodes_of_element >> std::ws;
+ // create node pointers with dummy coordinates to create element objects.
+ // True coordinates are set later in COOR and ELEV_I.
+ double dummy_coords[3] = {};
+ vec_nodes.reserve(fem_dim.n_nodes);
+ for (size_t i=0; i<fem_dim.n_nodes; i++)
+ vec_nodes.push_back(new MeshLib::Node(dummy_coords, i));
+ line_stream.clear();
+ }
+ //....................................................................
+ // NODE (node index for elements)
+ else if (line_string.compare("NODE") == 0)
+ {
+ assert(!vec_nodes.empty());
+
+ MeshElemType eleType = MeshElemType::INVALID;
+ if (fem_dim.n_nodes_of_element == 2)
+ eleType = MeshElemType::LINE;
+ else if (fem_dim.n_nodes_of_element == 3)
+ eleType = MeshElemType::TRIANGLE;
+ else if (fem_dim.n_nodes_of_element == 4 && fem_class.dimension == 2)
+ eleType = MeshElemType::TRIANGLE;
+ else if (fem_dim.n_nodes_of_element == 4 && fem_class.dimension == 3)
+ eleType = MeshElemType::TETRAHEDRON;
+ else if (fem_dim.n_nodes_of_element == 6 && fem_class.dimension == 3)
+ eleType = MeshElemType::PRISM;
+ else if (fem_dim.n_nodes_of_element == 8 && fem_class.dimension == 3)
+ eleType = MeshElemType::HEXAHEDRON;
+
+ if (eleType == MeshElemType::INVALID) {
+ ERR("FEFLOWInterface::readFEFLOWFile(): Unsupported element type with the number of node = %d and dim = %d", fem_dim.n_nodes_of_element, fem_class.dimension);
+ return nullptr;
+ }
+
+ vec_elements.reserve(fem_dim.n_elements);
+ for (size_t i=0; i<fem_dim.n_elements; i++)
+ {
+ getline(in, line_string);
+ vec_elements.push_back(readElement(fem_dim, eleType, line_string, vec_nodes));
+ }
+ }
+ //....................................................................
+ // COOR
+ else if (line_string.compare("COOR") == 0)
+ {
+ readNodeCoordinates(in, fem_class, fem_dim, vec_nodes);
+ }
+ //....................................................................
+ // ELEV_I
+ else if (line_string.compare("ELEV_I") == 0)
+ {
+ if (fem_class.dimension == 2)
+ continue;
+ readELEV(in, fem_class, fem_dim, vec_nodes);
+ }
+ //....................................................................
+ // EXTENTS
+ else if (line_string.compare("EXTENTS") == 0)
+ {
+ }
+ //....................................................................
+ // GRAVITY
+ else if (line_string.compare("GRAVITY") == 0)
+ {
+ getline(in, line_string);
+ line_stream.str(line_string);
+ double vec[3] = { };
+ line_stream >> vec[0] >> vec[1] >> vec[2];
+ if (vec[0] == 0.0 && vec[1] == -1.0 && vec[2] == 0.0)
+ // x-z plane
+ isXZplane = true;
+ line_stream.clear();
+ }
+ //....................................................................
+ // SUPERMESH
+ else if (line_string.compare("SUPERMESH") == 0)
+ {
+ readSuperMesh(in, fem_class, &points, &lines);
+ }
+ //....................................................................
+ }
+ in.close();
+
+ if (lines && lines->size() > 1)
+ {
+ this->setMaterialID(vec_elements, lines);
+ }
+
+ if (isXZplane)
+ {
+ for (auto* nod : vec_nodes)
+ {
+ (*nod)[2] = (*nod)[1];
+ (*nod)[1] = 0.0;
+ }
+ if (points)
+ {
+ for (auto* pt : *points)
+ {
+ (*pt)[2] = (*pt)[1];
+ (*pt)[1] = .0;
+ }
+ }
+ }
+
+ std::string project_name(BaseLib::extractBaseNameWithoutExtension(filename));
+ if (_geoObjects && points)
+ _geoObjects->addPointVec(points, project_name);
+ if (_geoObjects && lines)
+ _geoObjects->addPolylineVec(lines, project_name);
+
+ return new MeshLib::Mesh(project_name, vec_nodes, vec_elements);
+}
+
+
+void FEFLOWInterface::readNodeCoordinates(std::ifstream &in, const FEM_CLASS &fem_class, const FEM_DIM &fem_dim, std::vector<MeshLib::Node*> &vec_nodes)
+{
+ const size_t no_nodes_per_layer = (fem_class.dimension == 2) ? fem_dim.n_nodes : fem_dim.n_nodes / (fem_class.n_layers3d + 1);
+ const size_t n_lines = no_nodes_per_layer / 12 + 1;
+ const size_t n_layers = (fem_class.dimension == 3) ? fem_class.n_layers3d + 1 : 1;
+ std::string line_string;
+ std::stringstream line_stream;
+ double x;
+ char dummy_char; // for comma(,)
+ // x, y
+ for (unsigned k = 0; k < 2; k++)
+ {
+ // each line
+ for (size_t i = 0; i < n_lines; i++)
+ {
+ getline(in, line_string);
+ line_stream.str(line_string);
+ // maximum 12 columns
+ for (unsigned j = 0; j < 12; j++)
+ {
+ if (i * 12 + j >= no_nodes_per_layer)
+ break;
+ line_stream >> x >> dummy_char;
+ for (size_t l = 0; l < n_layers; l++)
+ {
+ size_t n = i * 12 + l * no_nodes_per_layer + j;
+ MeshLib::Node* m_nod = vec_nodes[n];
+ if (k == 0)
+ (*m_nod)[0] = x;
+ else
+ (*m_nod)[1] = x;
+ }
+ }
+ line_stream.clear();
+ }
+ }
+}
+
+void FEFLOWInterface::readELEV(std::ifstream &in, const FEM_CLASS &fem_class, const FEM_DIM &fem_dim, std::vector<MeshLib::Node*> &vec_nodes)
+{
+ const size_t no_nodes_per_layer = fem_dim.n_nodes / (fem_class.n_layers3d + 1);
+ double z = .0;
+ size_t n0 = 0;
+ std::string line_string;
+ std::stringstream line_stream;
+ for (size_t l = 0; l < fem_class.n_layers3d + 1; l++)
+ {
+ if (l > 0)
+ getline(in, line_string);
+
+ getline(in, line_string);
+ line_stream.str(line_string);
+ line_stream >> z >> n0;
+ line_stream.clear();
+
+ for (size_t i = 0; i < no_nodes_per_layer; i++)
+ {
+ size_t n = n0 - 1 + i + l * no_nodes_per_layer;
+ (*vec_nodes[n])[2] = z;
+ }
+ }
+}
+
+MeshLib::Element* FEFLOWInterface::readElement(const FEM_DIM &fem_dim, const MeshElemType elem_type, const std::string& line, const std::vector<MeshLib::Node*> &nodes)
+{
+ std::stringstream ss(line);
+ std::string elem_type_str("");
+
+ unsigned idx[8];
+ for (size_t i = 0; i < fem_dim.n_nodes_of_element; ++i)
+ ss >> idx[i];
+ MeshLib::Node** ele_nodes = new MeshLib::Node*[fem_dim.n_nodes_of_element];
+ for (unsigned k(0); k < fem_dim.n_nodes_of_element; ++k)
+ ele_nodes[k] = nodes[idx[k]-1];
+
+ switch (elem_type)
+ {
+ case MeshElemType::LINE:
+ return new MeshLib::Line(ele_nodes);
+ case MeshElemType::TRIANGLE:
+ return new MeshLib::Tri(ele_nodes);
+ case MeshElemType::QUAD:
+ return new MeshLib::Quad(ele_nodes);
+ case MeshElemType::TETRAHEDRON:
+ return new MeshLib::Tet(ele_nodes);
+ case MeshElemType::HEXAHEDRON:
+ return new MeshLib::Hex(ele_nodes);
+ case MeshElemType::PRISM:
+ return new MeshLib::Prism(ele_nodes);
+ default:
+ assert(false);
+ return nullptr;
+ }
+}
+
+void FEFLOWInterface::readPoints(QDomElement &nodesEle, const std::string &tag, int dim, std::vector<GeoLib::Point*> &points)
+{
+ QDomElement xmlEle = nodesEle.firstChildElement(QString::fromStdString(tag));
+ if (xmlEle.isNull())
+ return;
+ QString str_pt_list1 = xmlEle.text().simplified();
+ std::istringstream ss(str_pt_list1.toStdString());
+ while (!ss.eof())
+ {
+ int pt_id = 0;
+ double pt_xyz[3] =
+ { };
+ ss >> pt_id;
+ for (int i = 0; i < dim; i++)
+ ss >> pt_xyz[i];
+ GeoLib::Point* pnt = new GeoLib::Point(pt_xyz[0], pt_xyz[1], pt_xyz[2]); //id?
+ points[pt_id - 1] = pnt;
+ }
+}
+
+//
+void FEFLOWInterface::readSuperMesh(std::ifstream &in, const FEM_CLASS &fem_class, std::vector<GeoLib::Point*>** p_points, std::vector<GeoLib::Polyline*>** p_lines)
+{
+ // get XML strings
+ std::ostringstream oss;
+ std::string line_string;
+ while (true)
+ {
+ getline(in, line_string);
+ BaseLib::trim(line_string);
+ oss << line_string << "\n";
+ if (line_string.find("</supermesh>") != std::string::npos)
+ break;
+ }
+// std::cout << oss.str();
+ QString strXML(oss.str().c_str());
+ //qDebug() << strXML;
+
+ // convert string to XML
+ QDomDocument doc;
+ if (!doc.setContent(strXML))
+ {
+ ERR("FEFLOWInterface::readSuperMesh(): Illegal XML format error");
+ return;
+ }
+
+ // get geometry data from XML
+ QDomElement docElem = doc.documentElement(); // #supermesh
+ // #nodes
+ *p_points = new std::vector<GeoLib::Point*>();
+ std::vector<GeoLib::Point*>* points = *p_points;
+ QDomElement nodesEle = docElem.firstChildElement("nodes");
+ if (nodesEle.isNull())
+ return;
+ {
+ QString str = nodesEle.attribute("count");
+ const long n_points = str.toLong();
+ points->resize(n_points);
+ //fixed
+ readPoints(nodesEle, "fixed", fem_class.dimension, *points);
+ readPoints(nodesEle, "linear", fem_class.dimension, *points);
+ readPoints(nodesEle, "parabolic", fem_class.dimension, *points);
+ }
+
+ // #polygons
+ *p_lines = new std::vector<GeoLib::Polyline*>();
+ std::vector<GeoLib::Polyline*>* lines = *p_lines;
+ QDomElement polygonsEle = docElem.firstChildElement("polygons");
+ if (polygonsEle.isNull())
+ return;
+ {
+ QDomNode child = polygonsEle.firstChild();
+ while (!child.isNull())
+ {
+ if (child.nodeName() != "polygon")
+ {
+ child = child.nextSibling();
+ continue;
+ }
+ QDomElement xmlEle = child.firstChildElement("nodes");
+ if (xmlEle.isNull())
+ continue;
+ QString str = xmlEle.attribute("count");
+ const size_t n_points = str.toLong();
+ QString str_ptId_list = xmlEle.text().simplified();
+ {
+ GeoLib::Polyline* line = new GeoLib::Polyline(*points);
+ lines->push_back(line);
+ std::istringstream ss(str_ptId_list.toStdString());
+ for (size_t i = 0; i < n_points; i++)
+ {
+ int pt_id = 0;
+ ss >> pt_id;
+ line->addPoint(pt_id - 1);
+ }
+ line->addPoint(line->getPointID(0));
+ }
+ child = child.nextSibling();
+ }
+ }
+}
+
+void FEFLOWInterface::setMaterialID(std::vector<MeshLib::Element*> &vec_elements, std::vector<GeoLib::Polyline*>* lines)
+{
+ for (size_t i = 0; i < vec_elements.size(); i++)
+ {
+ MeshLib::Element* e = vec_elements[i];
+ MeshLib::Node gpt = e->getCenterOfGravity();
+ size_t matId = 0;
+ for (size_t j = 0; j < lines->size(); j++)
+ {
+ GeoLib::Polyline* poly = (*lines)[j];
+ if (!poly->isClosed())
+ continue;
+
+ GeoLib::Polygon polygon(*poly, true);
+ if (polygon.isPntInPolygon(gpt[0], gpt[1], gpt[2]))
+ {
+ matId = j;
+ break;
+ }
+ }
+ e->setValue(matId);
+ }
+}
+
+} // end namespace FileIO
diff --git a/FileIO/FEFLOWInterface.h b/FileIO/FEFLOWInterface.h
new file mode 100644
index 00000000000..ca89a6593d1
--- /dev/null
+++ b/FileIO/FEFLOWInterface.h
@@ -0,0 +1,145 @@
+/**
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#ifndef FEFLOWINTERFACE_H_
+#define FEFLOWINTERFACE_H_
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+
+#include "GEOObjects.h"
+#include "MeshEnums.h"
+
+class QDomElement;
+class QString;
+
+namespace MeshLib
+{
+class Mesh;
+class Element;
+class Node;
+}
+
+namespace FileIO
+{
+
+/**
+ * Read FEFLOW files into OGS data structure
+ */
+class FEFLOWInterface
+{
+public:
+ /// Constructor
+ explicit FEFLOWInterface(GeoLib::GEOObjects* obj)
+ : _geoObjects(obj)
+ {
+ }
+
+ /**
+ * read a FEFLOW Model file (*.fem) in ASCII format (Version 5.4)
+ *
+ * This function reads mesh data in addition to geometry data given in Supermesh.
+ *
+ * @param filename FEFLOW file name
+ * @return a pointer to a created OGS mesh
+ */
+ MeshLib::Mesh* readFEFLOWFile(const std::string &filename);
+
+private:
+ // CLASS
+ struct FEM_CLASS
+ {
+ unsigned problem_class;
+ unsigned time_mode;
+ unsigned orientation;
+ unsigned dimension;
+ unsigned n_layers3d;
+ unsigned saturation_flag;
+ unsigned save_fsize_rreal;
+ unsigned save_fsize_creal;
+
+ FEM_CLASS()
+ {
+ problem_class = 0;
+ time_mode = 0;
+ orientation = 0;
+ dimension = 0;
+ n_layers3d = 0;
+ saturation_flag = 0;
+ save_fsize_rreal = 0;
+ save_fsize_creal = 0;
+ }
+ };
+
+ // DIMENSION
+ struct FEM_DIM
+ {
+ size_t n_nodes;
+ size_t n_elements;
+ unsigned n_nodes_of_element;
+ unsigned n_steps;
+ unsigned icrank;
+ unsigned upwind;
+ size_t obs;
+ unsigned optim;
+ unsigned aquifer_type;
+ unsigned nwca;
+ size_t np_cor;
+ unsigned adaptive_mesh;
+ unsigned sp_fem_pcs_id;
+ unsigned sorption_type;
+ unsigned reaction_type;
+ unsigned dispersion_type;
+
+ FEM_DIM()
+ {
+ n_nodes = 0;
+ n_elements = 0;
+ n_nodes_of_element = 0;
+ n_steps = 0;
+ icrank = 0;
+ upwind = 0;
+ obs = 0;
+ optim = 0;
+ aquifer_type = 0;
+ nwca = 0;
+ np_cor = 0;
+ adaptive_mesh = 0;
+ sp_fem_pcs_id = 0;
+ sorption_type = 0;
+ reaction_type = 0;
+ dispersion_type = 0;
+ }
+ };
+
+ /// read node index and create a mesh element
+ MeshLib::Element* readElement(const FEM_DIM &fem_dim, const MeshElemType elem_type, const std::string& line, const std::vector<MeshLib::Node*> &nodes);
+
+ /// read node coordinates
+ void readNodeCoordinates(std::ifstream &in, const FEM_CLASS &fem_class, const FEM_DIM &fem_dim, std::vector<MeshLib::Node*> &nodes);
+
+ /// read elevation data
+ void readELEV(std::ifstream &in, const FEM_CLASS &fem_class, const FEM_DIM &fem_dim, std::vector<MeshLib::Node*> &vec_nodes);
+
+ /// read Supermesh data
+ void readSuperMesh(std::ifstream &feflow_file, const FEM_CLASS &fem_class, std::vector<GeoLib::Point*>** points, std::vector<GeoLib::Polyline*>** lines);
+
+ //// read point data in Supermesh
+ void readPoints(QDomElement &nodesEle, const std::string &tag, int dim, std::vector<GeoLib::Point*> &points);
+
+ /// set element material IDs
+ void setMaterialID(std::vector<MeshLib::Element*> &vec_elements, std::vector<GeoLib::Polyline*>* lines);
+
+ //// Geometric objects
+ GeoLib::GEOObjects* _geoObjects;
+};
+} // end namespace FileIO
+
+#endif /* FEFLOWINTERFACE_H_ */
--
GitLab