OGS
GocadSGridReader.cpp
Go to the documentation of this file.
1
10#include "GocadSGridReader.h"
11
12#include <algorithm>
13#include <boost/tokenizer.hpp>
14#include <cstddef>
15#include <fstream>
16#include <iostream>
17#include <sstream>
18
19#include "BaseLib/FileTools.h"
22#include "MeshLib/Mesh.h"
23
24namespace FileIO
25{
26namespace Gocad
27{
28using Bitset = boost::dynamic_bitset<>;
29
30GocadSGridReader::GocadSGridReader(std::string const& fname)
31 : _fname(fname),
32 _path(BaseLib::extractPath(fname)),
33 _n_face_sets(0),
34 _double_precision_binary(false),
35 _bin_pnts_in_double_precision(false)
36{
37 // check if file exists
38 std::ifstream in(_fname.c_str());
39 if (!in)
40 {
41 ERR("Could not open '{:s}'.", _fname);
42 in.close();
43 return;
44 }
45
46 bool pnts_read(false);
47
48 // read information in the stratigraphic grid file
49 std::string line;
50 while (std::getline(in, line))
51 {
52 if (line.empty()) // skip empty lines
53 {
54 continue;
55 }
56 if (line.back() == '\r') // check dos line ending
57 {
58 line.pop_back();
59 }
60 if (line.substr(0, 8) == "HEADER {")
61 {
62 parseHeader(in);
63 }
64 if (line.substr(0, 32) == "GOCAD_ORIGINAL_COORDINATE_SYSTEM")
65 {
67 continue;
68 }
69 if (line.substr(0, 7) == "AXIS_N ")
70 {
71 parseDims(line);
72 }
73 else if (line.substr(0, 12) == "POINTS_FILE ")
74 {
76 }
77 else if (line.substr(0, 9) == "PROPERTY ")
78 {
81 }
82 else if (line.substr(0, 35) == "BINARY_POINTS_IN_DOUBLE_PRECISION 1")
83 {
85 }
86 else if (line.substr(0, 11) == "FLAGS_FILE ")
87 {
89 }
90 else if (line.substr(0, 18) == "REGION_FLAGS_FILE ")
91 {
93 }
94 else if (line.substr(0, 7) == "REGION " ||
95 line.substr(0, 13) == "MODEL_REGION ")
96 {
97 regions.push_back(Gocad::parseRegion(line));
98 }
99 else if (line.substr(0, 12) == "MODEL_LAYER ")
100 {
101 layers.push_back(Gocad::parseLayer(line, regions));
102 }
103 else if (line.substr(0, 24) == "REGION_FLAGS_BIT_LENGTH ")
104 {
105 std::istringstream iss(line);
106 std::istream_iterator<std::string> it(iss);
107 it++;
108 std::size_t bit_length = std::atoi(it->c_str());
109 if (regions.size() != bit_length)
110 {
111 ERR("{:d} regions read but {:d} expected.\n", regions.size(),
112 bit_length);
113 throw std::runtime_error(
114 "Number of read regions differs from expected.\n");
115 }
116 }
117 else if (line.substr(0, 9) == "FACE_SET ")
118 {
119 // first read the points
120 if (!pnts_read)
121 {
123 pnts_read = true;
124 }
125 parseFaceSet(line, in);
126 }
127 }
128
129#ifndef NDEBUG
130 std::stringstream regions_ss;
131 regions_ss << regions.size() << " regions read:\n";
132 std::copy(regions.begin(), regions.end(),
133 std::ostream_iterator<Gocad::Region>(regions_ss, "\t"));
134 DBUG("{:s}", regions_ss.str());
135
136 std::stringstream layers_ss;
137 layers_ss << layers.size() << " layers read:\n";
138 std::copy(layers.begin(), layers.end(),
139 std::ostream_iterator<Gocad::Layer>(layers_ss, "\n"));
140 DBUG("{:s}", layers_ss.str());
141
142 std::stringstream properties_ss;
143 properties_ss << "meta data for " << _property_meta_data_vecs.size()
144 << " properties read:\n";
146 std::ostream_iterator<Gocad::Property>(properties_ss, "\n"));
147 DBUG("{:s}", properties_ss.str());
148#endif
149
150 // if not done already read the points
151 if (!pnts_read)
152 {
154 }
156 std::vector<Bitset> region_flags = readRegionFlagsBinary();
157 mapRegionFlagsToCellProperties(region_flags);
158
160
161 in.close();
162}
163
165{
166 for (auto node : _nodes)
167 {
168 delete node;
169 }
170 for (auto node : _split_nodes)
171 {
172 delete node;
173 }
174}
175
176std::unique_ptr<MeshLib::Mesh> GocadSGridReader::getMesh() const
177{
178 std::vector<MeshLib::Node*> nodes;
179 std::transform(_nodes.cbegin(), _nodes.cend(), std::back_inserter(nodes),
180 [](MeshLib::Node const* const node)
181 { return new MeshLib::Node(*node); });
182
183 std::vector<MeshLib::Element*> elements(createElements(nodes));
184 applySplitInformation(nodes, elements);
185
186 DBUG("Creating mesh from Gocad SGrid.");
187 std::unique_ptr<MeshLib::Mesh> mesh(new MeshLib::Mesh(
190 DBUG("Mesh created.");
191
192 return mesh;
193}
194
196{
197 std::vector<std::string> const& prop_names(getPropertyNames());
198 for (auto const& name : prop_names)
199 {
200 auto prop = getProperty(name);
201 if (!prop)
202 {
203 continue;
204 }
205
206 DBUG("Adding Gocad property '{:s}' with {:d} values.", name,
207 prop->_property_data.size());
208
209 auto pv = MeshLib::getOrCreateMeshProperty<double>(
211 if (pv == nullptr)
212 {
213 ERR("Could not create mesh property '{:s}'.", name);
214 continue;
215 }
216
217 pv->resize(prop->_property_data.size());
218 std::copy(prop->_property_data.cbegin(), prop->_property_data.cend(),
219 pv->begin());
220 }
221}
222
223void GocadSGridReader::parseHeader(std::istream& in)
224{
225 std::string line;
226 while (std::getline(in, line))
227 {
228 if (line.front() == '}')
229 {
230 return;
231 }
232 if (line.substr(0, 27) == "double_precision_binary: on")
233 {
235 }
236 }
238 {
239 DBUG(
240 "GocadSGridReader::parseHeader(): _double_precision_binary == "
241 "true.");
242 }
243}
244
245void GocadSGridReader::parseDims(std::string const& line)
246{
247 std::size_t x_dim(0);
248 std::size_t y_dim(0);
249 std::size_t z_dim(0);
250 boost::tokenizer<> tok(line);
251 auto it(tok.begin());
252 it++; // overread token "AXIS"
253 it++; // overread "N"
254 std::stringstream ssx(*(it),
255 std::stringstream::in | std::stringstream::out);
256 ssx >> x_dim;
257 it++;
258 std::stringstream ssy(*it, std::stringstream::in | std::stringstream::out);
259 ssy >> y_dim;
260 it++;
261 std::stringstream ssz(*it, std::stringstream::in | std::stringstream::out);
262 ssz >> z_dim;
263 _index_calculator = Gocad::IndexCalculator(x_dim, y_dim, z_dim);
264 DBUG(
265 "x_dim = {:d}, y_dim = {:d}, z_dim = {:d} => #nodes = {:d}, #cells = "
266 "{:d}",
267 x_dim, y_dim, z_dim, _index_calculator._n_nodes,
269}
270
271void GocadSGridReader::parseFileName(std::string const& line,
272 std::string& result_string) const
273{
274 boost::char_separator<char> sep(" \r");
275 boost::tokenizer<boost::char_separator<char>> tok(line, sep);
276 auto it(tok.begin());
277 ++it; // overread POINTS_FILE or FLAGS_FILE or REGION_FLAGS_FILE
278 result_string = _path + *it;
279}
280
285void GocadSGridReader::parseFaceSet(std::string& line, std::istream& in)
286{
287 // create and initialize a Gocad::Property object for storing face set data
288 Gocad::Property face_set_property;
289 face_set_property._property_id = _n_face_sets;
290 face_set_property._property_name = "FaceSet";
291 face_set_property._property_class_name = "FaceSetData";
292 face_set_property._property_unit = "unitless";
293 face_set_property._property_data_type = "double";
294 face_set_property._property_data_fname = "";
295 face_set_property._property_no_data_value = -1.0;
296 face_set_property._property_data.resize(_index_calculator._n_cells);
297 std::fill(face_set_property._property_data.begin(),
298 face_set_property._property_data.end(),
299 face_set_property._property_no_data_value);
300
301 std::istringstream iss(line);
302 std::istream_iterator<std::string> it(iss);
303 // Check first word is FACE_SET
304 if (*it != std::string("FACE_SET"))
305 {
306 ERR("Expected FACE_SET keyword but '{:s}' found.", it->c_str());
307 throw std::runtime_error(
308 "In GocadSGridReader::parseFaceSet() expected FACE_SET keyword not "
309 "found.");
310 }
311 ++it;
312 face_set_property._property_name += *it;
313 ++it;
314 auto const n_of_face_set_ids(
315 static_cast<std::size_t>(std::atoi(it->c_str())));
316 std::size_t face_set_id_cnt(0);
317
318 while (std::getline(in, line) && face_set_id_cnt < n_of_face_set_ids)
319 {
320 if (line.back() == '\r')
321 {
322 line.pop_back();
323 }
324 boost::char_separator<char> sep("\t ");
325 boost::tokenizer<boost::char_separator<char>> tokens(line, sep);
326
327 for (auto tok_it = tokens.begin(); tok_it != tokens.end();)
328 {
329 auto const id(static_cast<std::size_t>(std::atoi(tok_it->c_str())));
330 tok_it++;
331 auto const face_indicator(
332 static_cast<std::size_t>(std::atoi(tok_it->c_str())));
333 tok_it++;
334
335 if (id >= _index_calculator._n_nodes)
336 {
337 ERR("Face set id {:d} is greater than the number of nodes "
338 "({:d}).",
340 }
341 else
342 {
343 static_cast<GocadNode*>(_nodes[id])
344 ->setFaceSet(_n_face_sets, face_indicator);
345 std::array<std::size_t, 3> const c(
347 if (c[0] >= _index_calculator._x_dim - 1)
348 {
349 ERR("****** i coord {:d} to big for id {:d}.", c[0], id);
350 }
351 if (c[1] >= _index_calculator._y_dim - 1)
352 {
353 ERR("****** j coord {:d} to big for id {:d}.", c[1], id);
354 }
355 if (c[2] >= _index_calculator._z_dim - 1)
356 {
357 ERR("****** k coord {:d} to big for id {:d}.", c[2], id);
358 }
359 std::size_t const cell_id(
360 _index_calculator.getCellIdx(c[0], c[1], c[2]));
361 face_set_property._property_data[cell_id] =
362 static_cast<double>(face_indicator);
363 }
364 face_set_id_cnt++;
365 }
366 }
367
368 if (face_set_id_cnt != n_of_face_set_ids)
369 {
370 ERR("Expected {:d} number of face set ids, read {:d}.",
371 n_of_face_set_ids, face_set_id_cnt);
372 throw std::runtime_error(
373 "Expected number of face set points does not match number of read "
374 "points.");
375 }
376 _n_face_sets++;
377
378 // pre condition: split nodes are read already
379 for (auto split_node : _split_nodes)
380 {
381 std::size_t const id(_index_calculator(split_node->getGridCoords()));
382 split_node->transmitFaceDirections(*_nodes[id]);
383 }
384
385 _property_meta_data_vecs.push_back(face_set_property);
386}
387
388// Reads given number of bits (rounded up to next byte) into a bitset.
389// Used for reading region information which can be represented by some
390// number of bits.
391Bitset readBits(std::ifstream& in, const std::size_t bits)
392{
393 using block_t = Bitset::block_type;
394 auto const bytes = static_cast<std::size_t>(std::ceil(bits / 8.));
395 std::size_t const blocks =
396 bytes + 1 < sizeof(block_t) ? 1 : (bytes + 1) / sizeof(block_t);
397
398 std::vector<block_t> data;
399 data.resize(blocks);
400 std::fill_n(data.data(), blocks, 0);
401 in.read(reinterpret_cast<char*>(data.data()), bytes);
402
403 return Bitset(data.begin(), data.end());
404}
405
407{
408 std::ifstream in(_pnts_fname.c_str(), std::ios::binary);
409 if (!in)
410 {
411 ERR("Could not open points file '{:s}'.", _pnts_fname);
412 throw std::runtime_error("Could not open points file.");
413 }
414
415 std::size_t const n = _index_calculator._n_nodes;
416 _nodes.resize(n);
417
418 double coords[3];
419
420 std::size_t k = 0;
421 while (in && k < n * 3)
422 {
424 {
425 coords[k % 3] =
426 BaseLib::swapEndianness(BaseLib::readBinaryValue<double>(in));
427 }
428 else
429 {
430 coords[k % 3] =
431 BaseLib::swapEndianness(BaseLib::readBinaryValue<float>(in));
432 }
433 if ((k + 1) % 3 == 0)
434 {
435 const std::size_t layer_transition_idx(
437 _nodes[k / 3] = new GocadNode(coords, k / 3, layer_transition_idx);
438 }
439 k++;
440 }
441 if (k != n * 3 && !in.eof())
442 {
443 ERR("Read different number of points. Expected {:d} floats, got "
444 "{:d}.\n",
445 n * 3, k);
446 }
447}
448
450 std::vector<Bitset> const& rf)
451{
452 DBUG(
453 "GocadSGridReader::mapRegionFlagsToCellProperties region_flags.size: "
454 "{:d}",
455 rf.size());
456
457 Gocad::Property region_flags;
458 region_flags._property_id = 0;
459 region_flags._property_name = "RegionFlags";
460 region_flags._property_class_name = "RegionFlags";
461 region_flags._property_unit = "unitless";
462 region_flags._property_data_type = "int";
463 region_flags._property_data_fname = "";
464 region_flags._property_no_data_value = -1;
465 std::size_t const n = _index_calculator._n_cells;
466 region_flags._property_data.resize(n);
467 std::fill(region_flags._property_data.begin(),
468 region_flags._property_data.end(), -1);
469
470 // region flags are stored in each node ijk and give the region index for
471 // the ijk-th cell.
472 for (std::size_t i(0); i < _index_calculator._x_dim - 1; i++)
473 {
474 for (std::size_t j(0); j < _index_calculator._y_dim - 1; j++)
475 {
476 for (std::size_t k(0); k < _index_calculator._z_dim - 1; k++)
477 {
478 std::size_t const cell_id(
480 std::size_t const node_id(_index_calculator({i, j, k}));
481 for (auto& region : regions)
482 {
483 if (rf[node_id].test(region.bit))
484 {
485 region_flags._property_data[cell_id] += region.bit + 1;
486 }
487 }
488 }
489 }
490 }
491
492 _property_meta_data_vecs.push_back(region_flags);
493}
494
496{
497 for (auto& property : _property_meta_data_vecs)
498 {
499 std::string const& fname(property._property_data_fname);
500 if (property._property_data_fname.empty())
501 {
502 WARN("Empty filename for property {:s}.", property._property_name);
503 continue;
504 }
505 std::vector<float> float_properties =
506 BaseLib::readBinaryArray<float>(fname, _index_calculator._n_cells);
507 DBUG(
508 "GocadSGridReader::readElementPropertiesBinary(): Read {:d} float "
509 "properties from binary file.",
511
512 std::transform(float_properties.cbegin(), float_properties.cend(),
513 float_properties.begin(),
514 [](float const& val)
515 { return BaseLib::swapEndianness(val); });
516
517 property._property_data.resize(float_properties.size());
518 std::copy(float_properties.begin(), float_properties.end(),
519 property._property_data.begin());
520 if (property._property_data.empty())
521 {
522 ERR("Reading of element properties file '{:s}' failed.", fname);
523 }
524 }
525}
526
528{
529 std::vector<Bitset> result;
530
531 std::ifstream in(_region_flags_fname.c_str());
532 if (!in)
533 {
534 ERR("readRegionFlagsBinary(): Could not open file '{:s}' for input.\n",
536 in.close();
537 return result;
538 }
539
540 std::size_t const n = _index_calculator._n_nodes;
541 result.resize(n);
542
543 std::size_t k = 0;
544 while (in && k < n)
545 {
546 result[k++] = readBits(in, regions.size());
547 }
548 if (k != n && !in.eof())
549 {
550 ERR("Read different number of values. Expected {:d}, got {:d}.\n", n,
551 k);
552 }
553 return result;
554}
555std::vector<MeshLib::Element*> GocadSGridReader::createElements(
556 std::vector<MeshLib::Node*> const& nodes) const
557{
558 std::vector<MeshLib::Element*> elements;
559 elements.resize(_index_calculator._n_cells);
560 std::array<MeshLib::Node*, 8> element_nodes{};
561 std::size_t cnt(0);
562 for (std::size_t k(0); k < _index_calculator._z_dim - 1; k++)
563 {
564 for (std::size_t j(0); j < _index_calculator._y_dim - 1; j++)
565 {
566 for (std::size_t i(0); i < _index_calculator._x_dim - 1; i++)
567 {
568 element_nodes[0] = nodes[_index_calculator({i, j, k})];
569 element_nodes[1] = nodes[_index_calculator({i + 1, j, k})];
570 element_nodes[2] = nodes[_index_calculator({i + 1, j + 1, k})];
571 element_nodes[3] = nodes[_index_calculator({i, j + 1, k})];
572 element_nodes[4] = nodes[_index_calculator({i, j, k + 1})];
573 element_nodes[5] = nodes[_index_calculator({i + 1, j, k + 1})];
574 element_nodes[6] =
575 nodes[_index_calculator({i + 1, j + 1, k + 1})];
576 element_nodes[7] = nodes[_index_calculator({i, j + 1, k + 1})];
577 elements[cnt] = new MeshLib::Hex(
578 element_nodes, _index_calculator.getCellIdx(i, j, k));
579 cnt++;
580 }
581 }
582 }
583 return elements;
584}
585
587{
588 std::ifstream in(_fname.c_str());
589 if (!in)
590 {
591 ERR("Could not open '{:s}'.", _fname);
592 in.close();
593 return;
594 }
595
596 // read split information from the stratigraphic grid file
597 std::string line;
598 std::stringstream ss;
599 while (std::getline(in, line))
600 {
601 std::size_t pos(line.find("SPLIT "));
602 if (pos != std::string::npos)
603 {
604 ss << line.substr(pos + 6, line.size() - (pos + 6));
605 // read position in grid
606 std::array<std::size_t, 3> grid_coords{};
607 ss >> grid_coords[0];
608 ss >> grid_coords[1];
609 ss >> grid_coords[2];
610 // read coordinates for the split node
611 double coords[3];
612 ss >> coords[0];
613 ss >> coords[1];
614 ss >> coords[2];
615 // read the id
616 std::size_t id;
617 ss >> id;
618 // read the affected cells
619 std::bitset<8> affected_cells{};
620 for (std::size_t ac = 0; ac < 8; ++ac)
621 {
622 char bit;
623 ss >> bit;
624 affected_cells[ac] = bit != '0';
625 }
626 const std::size_t layer_transition_index(
627 _nodes[id]->getLayerTransitionIndex());
628 _split_nodes.push_back(new GocadSplitNode(coords, id, grid_coords,
629 affected_cells,
630 layer_transition_index));
631 }
632 }
633}
634
636 std::vector<MeshLib::Node*>& nodes,
637 std::vector<MeshLib::Element*> const& elements) const
638{
639 for (auto split_node : _split_nodes)
640 {
641 std::size_t const new_node_pos(nodes.size());
642 nodes.push_back(new MeshLib::Node(split_node->data(), new_node_pos));
643
644 // get grid coordinates
645 std::array<std::size_t, 3> const& gc(split_node->getGridCoords());
646 // get affected cells
647 auto const& affected_cells(split_node->getAffectedCells());
648 // get mesh node to substitute in elements
649 MeshLib::Node const* const node2sub(nodes[_index_calculator(gc)]);
650
651 if (affected_cells[0] && gc[0] < _index_calculator._x_dim - 1 &&
652 gc[1] < _index_calculator._y_dim - 1 &&
653 gc[2] < _index_calculator._z_dim - 1)
654 {
655 const std::size_t idx(
656 _index_calculator.getCellIdx(gc[0], gc[1], gc[2]));
657 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
658 }
659 if (affected_cells[1] && gc[0] > 0 &&
660 gc[1] < _index_calculator._y_dim - 1 &&
661 gc[2] < _index_calculator._z_dim - 1)
662 {
663 const std::size_t idx(
664 _index_calculator.getCellIdx(gc[0] - 1, gc[1], gc[2]));
665 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
666 }
667 if (affected_cells[2] && gc[1] > 0 &&
668 gc[0] < _index_calculator._x_dim - 1 &&
669 gc[2] < _index_calculator._z_dim - 1)
670 {
671 const std::size_t idx(
672 _index_calculator.getCellIdx(gc[0], gc[1] - 1, gc[2]));
673 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
674 }
675 if (affected_cells[3] && gc[0] > 0 && gc[1] > 0 &&
676 gc[2] < _index_calculator._z_dim - 1)
677 {
678 const std::size_t idx(
679 _index_calculator.getCellIdx(gc[0] - 1, gc[1] - 1, gc[2]));
680 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
681 }
682 if (affected_cells[4] && gc[2] > 0 &&
683 gc[0] < _index_calculator._x_dim - 1 &&
684 gc[1] < _index_calculator._y_dim - 1)
685 {
686 const std::size_t idx(
687 _index_calculator.getCellIdx(gc[0], gc[1], gc[2] - 1));
688 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
689 }
690 if (affected_cells[5] && gc[0] > 0 && gc[2] > 0 &&
691 gc[1] < _index_calculator._y_dim - 1)
692 {
693 const std::size_t idx(
694 _index_calculator.getCellIdx(gc[0] - 1, gc[1], gc[2] - 1));
695 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
696 }
697 if (affected_cells[6] && gc[1] > 0 && gc[2] > 0 &&
698 gc[0] < _index_calculator._x_dim - 1)
699 {
700 const std::size_t idx(
701 _index_calculator.getCellIdx(gc[0], gc[1] - 1, gc[2] - 1));
702 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
703 }
704 if (affected_cells[7] && gc[0] > 0 && gc[1] > 0 && gc[2] > 0)
705 {
706 const std::size_t idx(
707 _index_calculator.getCellIdx(gc[0] - 1, gc[1] - 1, gc[2] - 1));
708 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
709 }
710 }
711}
712
714 MeshLib::Node const* node2sub,
715 MeshLib::Node* substitute_node)
716{
717 // get the node pointers of the cell
718 MeshLib::Node* const* hex_nodes(hex->getNodes());
719 // search for the position the split node will be set to
720 MeshLib::Node* const* node_pos(
721 std::find(hex_nodes, hex_nodes + 8, node2sub));
722 // set the split node instead of the node2sub
723 if (node_pos != hex_nodes + 8)
724 {
725 const_cast<MeshLib::Node**>(
726 hex_nodes)[std::distance(hex_nodes, node_pos)] = substitute_node;
727 }
728}
729
731 std::string const& name) const
732{
733 auto const it(std::find_if(_property_meta_data_vecs.begin(),
735 [&name](Gocad::Property const& p)
736 { return p._property_name == name; }));
737 if (it == _property_meta_data_vecs.end())
738 {
739 return nullptr;
740 }
741 return &*it;
742}
743
744std::vector<std::string> GocadSGridReader::getPropertyNames() const
745{
746 std::vector<std::string> names;
747 std::transform(_property_meta_data_vecs.begin(),
749 std::back_inserter(names),
750 [](Gocad::Property const& p) { return p._property_name; });
751 return names;
752}
753
754std::unique_ptr<MeshLib::Mesh> GocadSGridReader::getFaceSetMesh(
755 std::size_t const face_set_number) const
756{
757 std::vector<MeshLib::Node*> face_set_nodes;
758 std::vector<MeshLib::Element*> face_set_elements;
759
760 for (auto const node : _nodes)
761 {
762 if (node->isMemberOfFaceSet(face_set_number))
763 {
764 addFaceSetQuad(node, face_set_number, face_set_nodes,
765 face_set_elements);
766 }
767 }
768
769 if (face_set_nodes.empty())
770 {
771 return nullptr;
772 }
773
774 for (auto const node : _split_nodes)
775 {
776 if (node->isMemberOfFaceSet(face_set_number))
777 {
778 if (node->getAffectedCells()[0])
779 {
780 addFaceSetQuad(node, face_set_number, face_set_nodes,
781 face_set_elements);
782 }
783 }
784 }
785
786 std::string const mesh_name("FaceSet-" + std::to_string(face_set_number));
787 return std::make_unique<MeshLib::Mesh>(mesh_name, face_set_nodes,
788 face_set_elements);
789}
790
792 GocadNode* face_set_node, std::size_t face_set_number,
793 std::vector<MeshLib::Node*>& face_set_nodes,
794 std::vector<MeshLib::Element*>& face_set_elements) const
795{
796 std::array<MeshLib::Node*, 4> quad_nodes{};
797 quad_nodes[0] = new GocadNode(*face_set_node);
798 const std::size_t id(face_set_node->getID());
799 std::array<std::size_t, 3> const c(_index_calculator.getCoordsForID(id));
800
801 const FaceDirection dir(face_set_node->getFaceDirection(face_set_number));
802 switch (dir)
803 {
804 case FaceDirection::U:
805 quad_nodes[1] = new GocadNode(
806 *_nodes[_index_calculator({c[0], c[1] + 1, c[2]})]);
807 quad_nodes[2] = new GocadNode(
808 *_nodes[_index_calculator({c[0], c[1] + 1, c[2] + 1})]);
809 quad_nodes[3] = new GocadNode(
810 *_nodes[_index_calculator({c[0], c[1], c[2] + 1})]);
811 break;
812 case FaceDirection::V:
813 quad_nodes[1] = new GocadNode(
814 *_nodes[_index_calculator({c[0] + 1, c[1], c[2]})]);
815 quad_nodes[2] = new GocadNode(
816 *_nodes[_index_calculator({c[0] + 1, c[1], c[2] + 1})]);
817 quad_nodes[3] = new GocadNode(
818 *_nodes[_index_calculator({c[0], c[1], c[2] + 1})]);
819 break;
820 case FaceDirection::W:
821 quad_nodes[1] = new GocadNode(
822 *_nodes[_index_calculator({c[0] + 1, c[1], c[2]})]);
823 quad_nodes[2] = new GocadNode(
824 *_nodes[_index_calculator({c[0] + 1, c[1] + 1, c[2]})]);
825 quad_nodes[3] = new GocadNode(
826 *_nodes[_index_calculator({c[0], c[1] + 1, c[2]})]);
827 break;
828 default:
829 ERR("Could not create face for node with id {:d}.", id);
830 }
831 std::copy(begin(quad_nodes), end(quad_nodes),
832 back_inserter(face_set_nodes));
833 face_set_elements.push_back(new MeshLib::Quad(quad_nodes));
834}
835
836} // end namespace Gocad
837} // end namespace FileIO
Filename manipulation routines.
Definition of the Hex class.
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:44
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:29
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:39
Definition of the Mesh class.
Definition of the Quad class.
FaceDirection getFaceDirection(std::size_t const face_set_number) const
Definition: GocadNode.h:90
Gocad::Property const * getProperty(std::string const &name) const
std::vector< Gocad::Property > _property_meta_data_vecs
std::vector< MeshLib::Element * > createElements(std::vector< MeshLib::Node * > const &nodes) const
std::unique_ptr< MeshLib::Mesh > getFaceSetMesh(std::size_t const face_set_number) const
std::vector< GocadNode * > _nodes
std::vector< Gocad::Region > regions
std::unique_ptr< MeshLib::Mesh > getMesh() const
std::vector< Bitset > readRegionFlagsBinary() const
void addGocadPropertiesToMesh(MeshLib::Mesh &mesh) const
void applySplitInformation(std::vector< MeshLib::Node * > &nodes, std::vector< MeshLib::Element * > const &elements) const
void mapRegionFlagsToCellProperties(std::vector< Bitset > const &rf)
Gocad::CoordinateSystem _coordinate_system
void parseFaceSet(std::string &line, std::istream &in)
void parseFileName(std::string const &line, std::string &result_string) const
Gocad::IndexCalculator _index_calculator
std::vector< GocadSplitNode * > _split_nodes
static void modifyElement(MeshLib::Element *hex, MeshLib::Node const *node2sub, MeshLib::Node *substitute_node)
void parseDims(std::string const &line)
std::vector< std::string > getPropertyNames() const
void addFaceSetQuad(GocadNode *face_set_node, std::size_t face_set_number, std::vector< MeshLib::Node * > &face_set_nodes, std::vector< MeshLib::Element * > &face_set_elements) const
void parseHeader(std::istream &in)
std::vector< Gocad::Layer > layers
std::size_t getCellIdx(std::size_t u, std::size_t v, std::size_t w) const
std::array< std::size_t, 3 > getCoordsForID(std::size_t id) const
std::size_t getID() const
Definition: Point3dWithID.h:62
virtual Node *const * getNodes() const =0
Get array of element nodes.
std::string extractPath(std::string const &pathname)
Definition: FileTools.cpp:207
std::string extractBaseNameWithoutExtension(std::string const &pathname)
Definition: FileTools.cpp:180
double swapEndianness(double const &v)
Definition: FileTools.cpp:147
Property parseGocadPropertyMetaData(std::string &line, std::istream &in, std::string const &path)
Definition: Property.cpp:30
Bitset readBits(std::ifstream &in, const std::size_t bits)
Layer parseLayer(std::string const &line, std::vector< Region > const &regions)
Definition: Layer.cpp:29
Region parseRegion(std::string const &line)
Definition: Region.cpp:26
boost::dynamic_bitset<> Bitset
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:37
TemplateElement< MeshLib::HexRule8 > Hex
Definition: Hex.h:25
double _property_no_data_value
Definition: Property.h:29
std::string _property_data_fname
Definition: Property.h:28
std::string _property_class_name
Definition: Property.h:25
std::string _property_unit
Definition: Property.h:26
std::string _property_name
Definition: Property.h:24
std::size_t _property_id
Definition: Property.h:23
std::string _property_data_type
Definition: Property.h:27
std::vector< double > _property_data
Definition: Property.h:43