OGS
GMSInterface.cpp
Go to the documentation of this file.
1
18#include "GMSInterface.h"
19
20#include <fstream>
21#include <list>
22
23#include "BaseLib/FileTools.h"
24#include "BaseLib/Logging.h"
25#include "BaseLib/StringTools.h"
30#include "MeshLib/Mesh.h"
31#include "MeshLib/MeshEnums.h"
32#include "MeshLib/Node.h"
33
34namespace
35{
36template <typename It>
37std::array<double, 3> parsePointCoordinates(It& it)
38{
39 return {std::strtod((++it)->c_str(), nullptr),
40 std::strtod((++it)->c_str(), nullptr),
41 std::strtod((++it)->c_str(), nullptr)};
42}
43} // namespace
44
45namespace FileIO
46{
47int GMSInterface::readBoreholesFromGMS(std::vector<GeoLib::Point*>& boreholes,
48 const std::string& filename)
49{
50 std::ifstream in(filename.c_str());
51 if (!in.is_open())
52 {
53 ERR("GMSInterface::readBoreholeFromGMS(): Could not open file {:s}.",
54 filename);
55 return 0;
56 }
57
58 double depth(-9999.0);
59 std::string line;
60 std::string cName;
61 std::string sName;
62 std::list<std::string>::const_iterator it;
63 GeoLib::StationBorehole* newBorehole = nullptr;
64
65 /* skipping first line because it contains field names */
66 std::getline(in, line);
67
68 /* read all stations */
69 while (std::getline(in, line))
70 {
71 std::list<std::string> fields = BaseLib::splitString(line, '\t');
72
73 if (fields.size() >= 5)
74 {
75 if (*fields.begin() == cName) // add new layer
76 {
77 it = fields.begin();
78 auto const pnt = parsePointCoordinates(it);
79
80 // check if current layer has a thickness of 0.0.
81 // if so skip it since it will mess with the vtk-visualisation
82 // later on!
83 if (pnt[2] != depth)
84 {
85 if (newBorehole == nullptr)
86 OGS_FATAL("Trying to access a nullptr.");
87 newBorehole->addSoilLayer(pnt[0], pnt[1], pnt[2], sName);
88 sName = (*(++it));
89 depth = pnt[2];
90 }
91 else
92 {
93 WARN(
94 "GMSInterface::readBoreholeFromGMS(): Skipped layer "
95 "'{:s}' in borehole '{:s}' because of thickness 0.0.",
96 sName, cName);
97 }
98 }
99 else // add new borehole
100 {
101 if (newBorehole != nullptr)
102 {
103 newBorehole->setDepth((*newBorehole)[2] - depth);
104 boreholes.push_back(newBorehole);
105 }
106 cName = *fields.begin();
107 it = fields.begin();
108 auto const pnt = parsePointCoordinates(it);
109 sName = (*(++it));
111 cName, pnt[0], pnt[1], pnt[2], 0);
112 depth = pnt[2];
113 }
114 }
115 else
116 {
117 ERR("GMSInterface::readBoreholeFromGMS(): Error reading format.");
118 }
119 }
120 // write the last borehole from the file
121 if (newBorehole != nullptr)
122 {
123 newBorehole->setDepth((*newBorehole)[2] - depth);
124 boreholes.push_back(newBorehole);
125 }
126 in.close();
127
128 if (boreholes.empty())
129 {
130 return 0;
131 }
132 return 1;
133}
134
136 const std::vector<GeoLib::Point*>* stations, const std::string& filename)
137{
138 std::ofstream out(filename.c_str(), std::ios::out);
139
140 // write header
141 out << "name"
142 << "\t" << std::fixed << "X"
143 << "\t"
144 << "Y"
145 << "\t"
146 << "Z"
147 << "\t"
148 << "soilID"
149 << "\n";
150
151 for (auto station_as_point : *stations)
152 {
153 auto* station = static_cast<GeoLib::StationBorehole*>(station_as_point);
154 std::vector<GeoLib::Point*> profile = station->getProfile();
155 std::vector<std::string> soilNames = station->getSoilNames();
156 // std::size_t idx = 0;
157 std::string current_soil_name;
158
159 std::size_t nLayers = profile.size();
160 for (std::size_t i = 1; i < nLayers; i++)
161 {
162 if ((i > 1) && (soilNames[i] == soilNames[i - 1]))
163 {
164 continue;
165 }
166 current_soil_name = soilNames[i];
167
168 out << station->getName() << "\t" << std::fixed
169 << (*(profile[i - 1]))[0] << "\t" << (*(profile[i - 1]))[1]
170 << "\t" << (*(profile[i - 1]))[2] << "\t"
171 << current_soil_name /*idx*/ << "\n";
172 }
173 out << station->getName() << "\t" << std::fixed
174 << (*(profile[nLayers - 1]))[0] << "\t"
175 << (*(profile[nLayers - 1]))[1] << "\t"
176 << (*(profile[nLayers - 1]))[2] << "\t" << current_soil_name
177 << "\n"; // this line marks the end of the borehole
178 }
179
180 out.close();
181}
182
183MeshLib::Mesh* GMSInterface::readGMS3DMMesh(const std::string& filename)
184{
185 std::string line;
186
187 std::ifstream in(filename.c_str());
188 if (!in.is_open())
189 {
190 ERR("GMSInterface::readGMS3DMMesh(): Could not open file {:s}.",
191 filename);
192 return nullptr;
193 }
194
195 // Read data from file
196 std::getline(in, line); // "MESH3D"
197 if (line != "MESH3D")
198 {
199 ERR("GMSInterface::readGMS3DMMesh(): Could not read expected file "
200 "header.");
201 return nullptr;
202 }
203
204 INFO("GMSInterface::readGMS3DMMesh(): Read GMS-3DM mesh.");
205 std::vector<MeshLib::Node*> nodes;
206 std::vector<MeshLib::Element*> elements;
207 std::vector<int> mat_ids;
208 std::map<unsigned, unsigned> id_map;
209
210 // elements are listed before nodes in 3dm-format, therefore
211 // traverse file twice and read first nodes and then elements
212 std::string dummy;
213 unsigned id(0);
214 unsigned count(0);
215 double x[3];
216 // read nodes
217 while (std::getline(in, line))
218 {
219 if (line[0] == 'N') // "ND" for Node
220 {
221 std::stringstream str(line);
222 str >> dummy >> id >> x[0] >> x[1] >> x[2];
223 auto* node = new MeshLib::Node(x, id);
224 id_map.insert(std::pair<unsigned, unsigned>(id, count++));
225 nodes.push_back(node);
226 }
227 }
228 in.close();
229
230 // NOTE: Element types E8H (Hex), E4Q (Quad), E3T (Tri) are not implemented
231 // yet read elements
232 in.open(filename.c_str());
233 std::getline(in, line); // "MESH3D"
234 unsigned node_idx[6];
235 int mat_id(0);
236 while (std::getline(in, line))
237 {
238 std::string element_id(line.substr(0, 3));
239 std::stringstream str(line);
240
241 if (element_id == "E6W") // Prism
242 {
243 str >> dummy >> id >> node_idx[0] >> node_idx[1] >> node_idx[2] >>
244 node_idx[3] >> node_idx[4] >> node_idx[5] >> mat_id;
245 auto** prism_nodes = new MeshLib::Node*[6];
246 for (unsigned k(0); k < 6; k++)
247 {
248 prism_nodes[k] = nodes[id_map.find(node_idx[k])->second];
249 }
250 elements.push_back(new MeshLib::Prism(prism_nodes));
251 mat_ids.push_back(mat_id);
252 }
253 else if (element_id == "E4T") // Tet
254 {
255 str >> dummy >> id >> node_idx[0] >> node_idx[1] >> node_idx[2] >>
256 node_idx[3] >> mat_id;
257 auto** tet_nodes = new MeshLib::Node*[4];
258 for (unsigned k(0); k < 4; k++)
259 {
260 tet_nodes[k] = nodes[id_map.find(node_idx[k])->second];
261 }
262 elements.push_back(new MeshLib::Tet(tet_nodes));
263 mat_ids.push_back(mat_id);
264 }
265 else if ((element_id == "E4P") ||
266 (element_id ==
267 "E5P")) // Pyramid (both do exist for some reason)
268 {
269 str >> dummy >> id >> node_idx[0] >> node_idx[1] >> node_idx[2] >>
270 node_idx[3] >> node_idx[4] >> mat_id;
271 auto** pyramid_nodes = new MeshLib::Node*[5];
272 for (unsigned k(0); k < 5; k++)
273 {
274 pyramid_nodes[k] = nodes[id_map.find(node_idx[k])->second];
275 }
276 elements.push_back(new MeshLib::Pyramid(pyramid_nodes));
277 mat_ids.push_back(mat_id);
278 }
279 else if (element_id == "ND ")
280 { // Node
281 continue; // skip because nodes have already been read
282 }
283 else // default
284 {
285 WARN(
286 "GMSInterface::readGMS3DMMesh() - Element type '{:s}' not "
287 "recognised.",
288 element_id);
289 return nullptr;
290 }
291 }
292
293 in.close();
294 INFO("GMSInterface::readGMS3DMMesh(): finished.");
295
296 const std::string mesh_name =
298 MeshLib::Properties properties;
299 if (mat_ids.size() == elements.size())
300 {
301 auto* const opt_pv = properties.createNewPropertyVector<int>(
302 "MaterialIDs", MeshLib::MeshItemType::Cell);
303 if (!opt_pv)
304 {
305 ERR("Could not create PropertyVector for material ids.");
306 BaseLib::cleanupVectorElements(nodes, elements);
307 return nullptr;
308 }
309 opt_pv->reserve(mat_ids.size());
310 std::copy(mat_ids.cbegin(), mat_ids.cend(),
311 std::back_inserter(*opt_pv));
312 }
313 else
314 {
315 ERR("Ignoring Material IDs information (does not match number of "
316 "elements).");
317 }
318 return new MeshLib::Mesh(mesh_name, nodes, elements, properties);
319}
320
321} // namespace FileIO
#define OGS_FATAL(...)
Definition: Error.h:26
Filename manipulation routines.
Definition of the GMSInterface class.
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:34
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:44
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:39
Definition of mesh-related Enumerations.
Definition of the Mesh class.
Definition of the Node class.
Definition of the Prism class.
Definition of the Pyramid class.
Definition of the StationBorehole class.
Definition of string helper functions.
Definition of the Tet class.
static MeshLib::Mesh * readGMS3DMMesh(const std::string &filename)
Reads a GMS *.3dm file and converts it to an CFEMesh.
static int readBoreholesFromGMS(std::vector< GeoLib::Point * > &boreholes, const std::string &filename)
Imports borehole data from a file in GMS-format.
static void writeBoreholesToGMS(const std::vector< GeoLib::Point * > *stations, const std::string &filename)
A borehole as a geometric object.
static StationBorehole * createStation(const std::string &name, double x, double y, double z, double depth, const std::string &date="")
Creates a new borehole object based on the given parameters.
void addSoilLayer(double thickness, const std::string &soil_name)
Add a soil layer to the boreholes stratigraphy.
const std::vector< Point * > & getProfile() const
void setDepth(double depth)
Sets the depth of the borehole.
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
Definition: Properties.h:36
PropertyVector< T > * createNewPropertyVector(std::string const &name, MeshItemType mesh_item_type, std::size_t n_components=1)
void cleanupVectorElements(std::vector< T * > &items)
Definition: Algorithm.h:300
std::string extractBaseNameWithoutExtension(std::string const &pathname)
Definition: FileTools.cpp:180
std::vector< std::string > splitString(std::string const &str)
Definition: StringTools.cpp:28
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:37
std::array< double, 3 > parsePointCoordinates(It &it)