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