swh:1:snp:6088ab52ef49920e01e3f334cdf4d5d6c8a822b9
Tip revision: 1566f1eaf7b7c7ad4fd639f451b6da4eb54298f2 authored by Tobias Meisel on 06 September 2021, 11:34:58 UTC
[MeL/IO] HDF/XDMF: New multimesh output tests
[MeL/IO] HDF/XDMF: New multimesh output tests
Tip revision: 1566f1e
GMSInterface.cpp
/**
* \file
* \author Karsten Rink
* \date 2010-06-08
* \brief Implementation of the GMSInterface class.
*
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
* @file GMSInterface.cpp
* @date 2010-06-08
* @author Karsten Rink
*/
#include "GMSInterface.h"
#include <fstream>
#include <list>
#include "BaseLib/FileTools.h"
#include "BaseLib/Logging.h"
#include "BaseLib/StringTools.h"
#include "GeoLib/StationBorehole.h"
#include "MeshLib/Elements/Prism.h"
#include "MeshLib/Elements/Pyramid.h"
#include "MeshLib/Elements/Tet.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/MeshEnums.h"
#include "MeshLib/Node.h"
namespace
{
template <typename It>
std::array<double, 3> parsePointCoordinates(It& it)
{
return {std::strtod((++it)->c_str(), nullptr),
std::strtod((++it)->c_str(), nullptr),
std::strtod((++it)->c_str(), nullptr)};
}
} // namespace
namespace FileIO
{
int GMSInterface::readBoreholesFromGMS(std::vector<GeoLib::Point*>* boreholes,
const std::string& filename)
{
std::ifstream in(filename.c_str());
if (!in.is_open())
{
ERR("GMSInterface::readBoreholeFromGMS(): Could not open file {:s}.",
filename);
return 0;
}
double depth(-9999.0);
std::string line;
std::string cName;
std::string sName;
std::list<std::string>::const_iterator it;
GeoLib::StationBorehole* newBorehole = nullptr;
/* skipping first line because it contains field names */
std::getline(in, line);
/* read all stations */
while (std::getline(in, line))
{
std::list<std::string> fields = BaseLib::splitString(line, '\t');
if (fields.size() >= 5)
{
if (*fields.begin() == cName) // add new layer
{
it = fields.begin();
auto const pnt = parsePointCoordinates(it);
// 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)
{
if (newBorehole == nullptr)
OGS_FATAL("Trying to access a nullptr.");
newBorehole->addSoilLayer(pnt[0], pnt[1], pnt[2], sName);
sName = (*(++it));
depth = pnt[2];
}
else
{
WARN(
"GMSInterface::readBoreholeFromGMS(): Skipped layer "
"'{:s}' in borehole '{:s}' because of thickness 0.0.",
sName, cName);
}
}
else // add new borehole
{
if (newBorehole != nullptr)
{
newBorehole->setDepth((*newBorehole)[2] - depth);
boreholes->push_back(newBorehole);
}
cName = *fields.begin();
it = fields.begin();
auto const pnt = parsePointCoordinates(it);
sName = (*(++it));
newBorehole = GeoLib::StationBorehole::createStation(
cName, pnt[0], pnt[1], pnt[2], 0);
depth = pnt[2];
}
}
else
{
ERR("GMSInterface::readBoreholeFromGMS(): Error reading format.");
}
}
// write the last borehole from the file
if (newBorehole != nullptr)
{
newBorehole->setDepth((*newBorehole)[2] - depth);
boreholes->push_back(newBorehole);
}
in.close();
if (boreholes->empty())
{
return 0;
}
return 1;
}
void GMSInterface::writeBoreholesToGMS(
const std::vector<GeoLib::Point*>* stations, const std::string& filename)
{
std::ofstream out(filename.c_str(), std::ios::out);
// write header
out << "name"
<< "\t" << std::fixed << "X"
<< "\t"
<< "Y"
<< "\t"
<< "Z"
<< "\t"
<< "soilID"
<< "\n";
for (auto station_as_point : *stations)
{
auto* station = static_cast<GeoLib::StationBorehole*>(station_as_point);
std::vector<GeoLib::Point*> profile = station->getProfile();
std::vector<std::string> soilNames = station->getSoilNames();
// std::size_t idx = 0;
std::string current_soil_name;
std::size_t nLayers = profile.size();
for (std::size_t i = 1; i < nLayers; i++)
{
if ((i > 1) && (soilNames[i] == soilNames[i - 1]))
{
continue;
}
current_soil_name = soilNames[i];
out << station->getName() << "\t" << std::fixed
<< (*(profile[i - 1]))[0] << "\t" << (*(profile[i - 1]))[1]
<< "\t" << (*(profile[i - 1]))[2] << "\t"
<< current_soil_name /*idx*/ << "\n";
}
out << station->getName() << "\t" << std::fixed
<< (*(profile[nLayers - 1]))[0] << "\t"
<< (*(profile[nLayers - 1]))[1] << "\t"
<< (*(profile[nLayers - 1]))[2] << "\t" << current_soil_name
<< "\n"; // this line marks the end of the borehole
}
out.close();
}
MeshLib::Mesh* GMSInterface::readGMS3DMMesh(const std::string& filename)
{
std::string line;
std::ifstream in(filename.c_str());
if (!in.is_open())
{
ERR("GMSInterface::readGMS3DMMesh(): Could not open file {:s}.",
filename);
return nullptr;
}
// Read data from file
std::getline(in, line); // "MESH3D"
if (line != "MESH3D")
{
ERR("GMSInterface::readGMS3DMMesh(): Could not read expected file "
"header.");
return nullptr;
}
INFO("GMSInterface::readGMS3DMMesh(): Read GMS-3DM mesh.");
std::vector<MeshLib::Node*> nodes;
std::vector<MeshLib::Element*> elements;
std::vector<int> mat_ids;
std::map<unsigned, unsigned> id_map;
// elements are listed before nodes in 3dm-format, therefore
// traverse file twice and read first nodes and then elements
std::string dummy;
unsigned id(0);
unsigned count(0);
double x[3];
// read nodes
while (std::getline(in, line))
{
if (line[0] == 'N') // "ND" for Node
{
std::stringstream str(line);
str >> dummy >> id >> x[0] >> x[1] >> x[2];
auto* node = new MeshLib::Node(x, id);
id_map.insert(std::pair<unsigned, unsigned>(id, count++));
nodes.push_back(node);
}
}
in.close();
// NOTE: Element types E8H (Hex), E4Q (Quad), E3T (Tri) are not implemented
// yet read elements
in.open(filename.c_str());
std::getline(in, line); // "MESH3D"
unsigned node_idx[6];
int mat_id(0);
while (std::getline(in, line))
{
std::string element_id(line.substr(0, 3));
std::stringstream str(line);
if (element_id == "E6W") // Prism
{
str >> dummy >> id >> node_idx[0] >> node_idx[1] >> node_idx[2] >>
node_idx[3] >> node_idx[4] >> node_idx[5] >> mat_id;
auto** prism_nodes = new MeshLib::Node*[6];
for (unsigned k(0); k < 6; k++)
{
prism_nodes[k] = nodes[id_map.find(node_idx[k])->second];
}
elements.push_back(new MeshLib::Prism(prism_nodes));
mat_ids.push_back(mat_id);
}
else if (element_id == "E4T") // Tet
{
str >> dummy >> id >> node_idx[0] >> node_idx[1] >> node_idx[2] >>
node_idx[3] >> mat_id;
auto** tet_nodes = new MeshLib::Node*[4];
for (unsigned k(0); k < 4; k++)
{
tet_nodes[k] = nodes[id_map.find(node_idx[k])->second];
}
elements.push_back(new MeshLib::Tet(tet_nodes));
mat_ids.push_back(mat_id);
}
else if ((element_id == "E4P") ||
(element_id ==
"E5P")) // Pyramid (both do exist for some reason)
{
str >> dummy >> id >> node_idx[0] >> node_idx[1] >> node_idx[2] >>
node_idx[3] >> node_idx[4] >> mat_id;
auto** pyramid_nodes = new MeshLib::Node*[5];
for (unsigned k(0); k < 5; k++)
{
pyramid_nodes[k] = nodes[id_map.find(node_idx[k])->second];
}
elements.push_back(new MeshLib::Pyramid(pyramid_nodes));
mat_ids.push_back(mat_id);
}
else if (element_id == "ND ")
{ // Node
continue; // skip because nodes have already been read
}
else // default
{
WARN(
"GMSInterface::readGMS3DMMesh() - Element type '{:s}' not "
"recognised.",
element_id);
return nullptr;
}
}
in.close();
INFO("GMSInterface::readGMS3DMMesh(): finished.");
const std::string mesh_name =
BaseLib::extractBaseNameWithoutExtension(filename);
MeshLib::Properties properties;
if (mat_ids.size() == elements.size())
{
auto* const opt_pv = properties.createNewPropertyVector<int>(
"MaterialIDs", MeshLib::MeshItemType::Cell);
if (!opt_pv)
{
ERR("Could not create PropertyVector for material ids.");
BaseLib::cleanupVectorElements(nodes, elements);
return nullptr;
}
opt_pv->reserve(mat_ids.size());
std::copy(mat_ids.cbegin(), mat_ids.cend(),
std::back_inserter(*opt_pv));
}
else
{
ERR("Ignoring Material IDs information (does not match number of "
"elements).");
}
return new MeshLib::Mesh(mesh_name, nodes, elements, properties);
}
} // namespace FileIO