Loading [MathJax]/extensions/tex2jax.js
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 pv->assign(prop->_property_data);
221 }
222}
223
224void GocadSGridReader::parseHeader(std::istream& in)
225{
226 std::string line;
227 while (std::getline(in, line))
228 {
229 if (line.front() == '}')
230 {
231 return;
232 }
233 if (line.substr(0, 27) == "double_precision_binary: on")
234 {
236 }
237 }
239 {
240 DBUG(
241 "GocadSGridReader::parseHeader(): _double_precision_binary == "
242 "true.");
243 }
244}
245
246void GocadSGridReader::parseDims(std::string const& line)
247{
248 std::size_t x_dim(0);
249 std::size_t y_dim(0);
250 std::size_t z_dim(0);
251 boost::tokenizer<> tok(line);
252 auto it(tok.begin());
253 it++; // overread token "AXIS"
254 it++; // overread "N"
255 std::stringstream ssx(*(it),
256 std::stringstream::in | std::stringstream::out);
257 ssx >> x_dim;
258 it++;
259 std::stringstream ssy(*it, std::stringstream::in | std::stringstream::out);
260 ssy >> y_dim;
261 it++;
262 std::stringstream ssz(*it, std::stringstream::in | std::stringstream::out);
263 ssz >> z_dim;
264 _index_calculator = Gocad::IndexCalculator(x_dim, y_dim, z_dim);
265 DBUG(
266 "x_dim = {:d}, y_dim = {:d}, z_dim = {:d} => #nodes = {:d}, #cells = "
267 "{:d}",
268 x_dim, y_dim, z_dim, _index_calculator._n_nodes,
270}
271
272void GocadSGridReader::parseFileName(std::string const& line,
273 std::string& result_string) const
274{
275 boost::char_separator<char> sep(" \r");
276 boost::tokenizer<boost::char_separator<char>> tok(line, sep);
277 auto it(tok.begin());
278 ++it; // overread POINTS_FILE or FLAGS_FILE or REGION_FLAGS_FILE
279 result_string = _path + *it;
280}
281
286void GocadSGridReader::parseFaceSet(std::string& line, std::istream& in)
287{
288 // create and initialize a Gocad::Property object for storing face set data
289 Gocad::Property face_set_property;
290 face_set_property._property_id = _n_face_sets;
291 face_set_property._property_name = "FaceSet";
292 face_set_property._property_class_name = "FaceSetData";
293 face_set_property._property_unit = "unitless";
294 face_set_property._property_data_type = "double";
295 face_set_property._property_data_fname = "";
296 face_set_property._property_no_data_value = -1.0;
297 face_set_property._property_data.resize(_index_calculator._n_cells);
298 std::fill(face_set_property._property_data.begin(),
299 face_set_property._property_data.end(),
300 face_set_property._property_no_data_value);
301
302 std::istringstream iss(line);
303 std::istream_iterator<std::string> it(iss);
304 // Check first word is FACE_SET
305 if (*it != std::string("FACE_SET"))
306 {
307 ERR("Expected FACE_SET keyword but '{:s}' found.", it->c_str());
308 throw std::runtime_error(
309 "In GocadSGridReader::parseFaceSet() expected FACE_SET keyword not "
310 "found.");
311 }
312 ++it;
313 face_set_property._property_name += *it;
314 ++it;
315 auto const n_of_face_set_ids(
316 static_cast<std::size_t>(std::atoi(it->c_str())));
317 std::size_t face_set_id_cnt(0);
318
319 while (std::getline(in, line) && face_set_id_cnt < n_of_face_set_ids)
320 {
321 if (line.back() == '\r')
322 {
323 line.pop_back();
324 }
325 boost::char_separator<char> sep("\t ");
326 boost::tokenizer<boost::char_separator<char>> tokens(line, sep);
327
328 for (auto tok_it = tokens.begin(); tok_it != tokens.end();)
329 {
330 auto const id(static_cast<std::size_t>(std::atoi(tok_it->c_str())));
331 tok_it++;
332 auto const face_indicator(
333 static_cast<std::size_t>(std::atoi(tok_it->c_str())));
334 tok_it++;
335
336 if (id >= _index_calculator._n_nodes)
337 {
338 ERR("Face set id {:d} is greater than the number of nodes "
339 "({:d}).",
341 }
342 else
343 {
344 static_cast<GocadNode*>(_nodes[id])
345 ->setFaceSet(_n_face_sets, face_indicator);
346 std::array<std::size_t, 3> const c(
348 if (c[0] >= _index_calculator._x_dim - 1)
349 {
350 ERR("****** i coord {:d} to big for id {:d}.", c[0], id);
351 }
352 if (c[1] >= _index_calculator._y_dim - 1)
353 {
354 ERR("****** j coord {:d} to big for id {:d}.", c[1], id);
355 }
356 if (c[2] >= _index_calculator._z_dim - 1)
357 {
358 ERR("****** k coord {:d} to big for id {:d}.", c[2], id);
359 }
360 std::size_t const cell_id(
361 _index_calculator.getCellIdx(c[0], c[1], c[2]));
362 face_set_property._property_data[cell_id] =
363 static_cast<double>(face_indicator);
364 }
365 face_set_id_cnt++;
366 }
367 }
368
369 if (face_set_id_cnt != n_of_face_set_ids)
370 {
371 ERR("Expected {:d} number of face set ids, read {:d}.",
372 n_of_face_set_ids, face_set_id_cnt);
373 throw std::runtime_error(
374 "Expected number of face set points does not match number of read "
375 "points.");
376 }
377 _n_face_sets++;
378
379 // pre condition: split nodes are read already
380 for (auto split_node : _split_nodes)
381 {
382 std::size_t const id(_index_calculator(split_node->getGridCoords()));
383 split_node->transmitFaceDirections(*_nodes[id]);
384 }
385
386 _property_meta_data_vecs.push_back(face_set_property);
387}
388
389// Reads given number of bits (rounded up to next byte) into a bitset.
390// Used for reading region information which can be represented by some
391// number of bits.
392Bitset readBits(std::ifstream& in, const std::size_t bits)
393{
394 using block_t = Bitset::block_type;
395 auto const bytes = static_cast<std::size_t>(std::ceil(bits / 8.));
396 std::size_t const blocks =
397 bytes + 1 < sizeof(block_t) ? 1 : (bytes + 1) / sizeof(block_t);
398
399 std::vector<block_t> data;
400 data.resize(blocks);
401 std::fill_n(data.data(), blocks, 0);
402 in.read(reinterpret_cast<char*>(data.data()), bytes);
403
404 return Bitset(data.begin(), data.end());
405}
406
408{
409 std::ifstream in(_pnts_fname.c_str(), std::ios::binary);
410 if (!in)
411 {
412 ERR("Could not open points file '{:s}'.", _pnts_fname);
413 throw std::runtime_error("Could not open points file.");
414 }
415
416 std::size_t const n = _index_calculator._n_nodes;
417 _nodes.resize(n);
418
419 double coords[3];
420
421 std::size_t k = 0;
422 while (in && k < n * 3)
423 {
425 {
426 coords[k % 3] =
428 }
429 else
430 {
431 coords[k % 3] =
433 }
434 if ((k + 1) % 3 == 0)
435 {
436 const std::size_t layer_transition_idx(
438 _nodes[k / 3] = new GocadNode(coords, k / 3, layer_transition_idx);
439 }
440 k++;
441 }
442 if (k != n * 3 && !in.eof())
443 {
444 ERR("Read different number of points. Expected {:d} floats, got "
445 "{:d}.\n",
446 n * 3, k);
447 }
448}
449
451 std::vector<Bitset> const& rf)
452{
453 DBUG(
454 "GocadSGridReader::mapRegionFlagsToCellProperties region_flags.size: "
455 "{:d}",
456 rf.size());
457
458 Gocad::Property region_flags;
459 region_flags._property_id = 0;
460 region_flags._property_name = "RegionFlags";
461 region_flags._property_class_name = "RegionFlags";
462 region_flags._property_unit = "unitless";
463 region_flags._property_data_type = "int";
464 region_flags._property_data_fname = "";
465 region_flags._property_no_data_value = -1;
466 std::size_t const n = _index_calculator._n_cells;
467 region_flags._property_data.resize(n);
468 std::fill(region_flags._property_data.begin(),
469 region_flags._property_data.end(), -1);
470
471 // region flags are stored in each node ijk and give the region index for
472 // the ijk-th cell.
473 for (std::size_t i(0); i < _index_calculator._x_dim - 1; i++)
474 {
475 for (std::size_t j(0); j < _index_calculator._y_dim - 1; j++)
476 {
477 for (std::size_t k(0); k < _index_calculator._z_dim - 1; k++)
478 {
479 std::size_t const cell_id(
481 std::size_t const node_id(_index_calculator({i, j, k}));
482 for (auto& region : regions)
483 {
484 if (rf[node_id].test(region.bit))
485 {
486 region_flags._property_data[cell_id] += region.bit + 1;
487 }
488 }
489 }
490 }
491 }
492
493 _property_meta_data_vecs.push_back(region_flags);
494}
495
497{
498 for (auto& property : _property_meta_data_vecs)
499 {
500 std::string const& fname(property._property_data_fname);
501 if (property._property_data_fname.empty())
502 {
503 WARN("Empty filename for property {:s}.", property._property_name);
504 continue;
505 }
506 std::vector<float> float_properties = BaseLib::readBinaryVector<float>(
507 fname, 0, _index_calculator._n_cells);
508 DBUG(
509 "GocadSGridReader::readElementPropertiesBinary(): Read {:d} float "
510 "properties from binary file.",
512
513 std::transform(float_properties.cbegin(), float_properties.cend(),
514 float_properties.begin(),
515 [](float const& val)
516 { return BaseLib::swapEndianness(val); });
517
518 property._property_data.resize(float_properties.size());
519 std::copy(float_properties.begin(), float_properties.end(),
520 property._property_data.begin());
521 if (property._property_data.empty())
522 {
523 ERR("Reading of element properties file '{:s}' failed.", fname);
524 }
525 }
526}
527
529{
530 std::vector<Bitset> result;
531
532 std::ifstream in(_region_flags_fname.c_str());
533 if (!in)
534 {
535 ERR("readRegionFlagsBinary(): Could not open file '{:s}' for input.\n",
537 in.close();
538 return result;
539 }
540
541 std::size_t const n = _index_calculator._n_nodes;
542 result.resize(n);
543
544 std::size_t k = 0;
545 while (in && k < n)
546 {
547 result[k++] = readBits(in, regions.size());
548 }
549 if (k != n && !in.eof())
550 {
551 ERR("Read different number of values. Expected {:d}, got {:d}.\n", n,
552 k);
553 }
554 return result;
555}
556std::vector<MeshLib::Element*> GocadSGridReader::createElements(
557 std::vector<MeshLib::Node*> const& nodes) const
558{
559 std::vector<MeshLib::Element*> elements;
560 elements.resize(_index_calculator._n_cells);
561 std::array<MeshLib::Node*, 8> element_nodes{};
562 std::size_t cnt(0);
563 for (std::size_t k(0); k < _index_calculator._z_dim - 1; k++)
564 {
565 for (std::size_t j(0); j < _index_calculator._y_dim - 1; j++)
566 {
567 for (std::size_t i(0); i < _index_calculator._x_dim - 1; i++)
568 {
569 element_nodes[0] = nodes[_index_calculator({i, j, k})];
570 element_nodes[1] = nodes[_index_calculator({i + 1, j, k})];
571 element_nodes[2] = nodes[_index_calculator({i + 1, j + 1, k})];
572 element_nodes[3] = nodes[_index_calculator({i, j + 1, k})];
573 element_nodes[4] = nodes[_index_calculator({i, j, k + 1})];
574 element_nodes[5] = nodes[_index_calculator({i + 1, j, k + 1})];
575 element_nodes[6] =
576 nodes[_index_calculator({i + 1, j + 1, k + 1})];
577 element_nodes[7] = nodes[_index_calculator({i, j + 1, k + 1})];
578 elements[cnt] = new MeshLib::Hex(
579 element_nodes, _index_calculator.getCellIdx(i, j, k));
580 cnt++;
581 }
582 }
583 }
584 return elements;
585}
586
588{
589 std::ifstream in(_fname.c_str());
590 if (!in)
591 {
592 ERR("Could not open '{:s}'.", _fname);
593 in.close();
594 return;
595 }
596
597 // read split information from the stratigraphic grid file
598 std::string line;
599 std::stringstream ss;
600 while (std::getline(in, line))
601 {
602 std::size_t pos(line.find("SPLIT "));
603 if (pos != std::string::npos)
604 {
605 ss << line.substr(pos + 6, line.size() - (pos + 6));
606 // read position in grid
607 std::array<std::size_t, 3> grid_coords{};
608 ss >> grid_coords[0];
609 ss >> grid_coords[1];
610 ss >> grid_coords[2];
611 // read coordinates for the split node
612 double coords[3];
613 ss >> coords[0];
614 ss >> coords[1];
615 ss >> coords[2];
616 // read the id
617 std::size_t id;
618 ss >> id;
619 // read the affected cells
620 std::bitset<8> affected_cells{};
621 for (std::size_t ac = 0; ac < 8; ++ac)
622 {
623 char bit;
624 ss >> bit;
625 affected_cells[ac] = bit != '0';
626 }
627 const std::size_t layer_transition_index(
628 _nodes[id]->getLayerTransitionIndex());
629 _split_nodes.push_back(new GocadSplitNode(coords, id, grid_coords,
630 affected_cells,
631 layer_transition_index));
632 }
633 }
634}
635
637 std::vector<MeshLib::Node*>& nodes,
638 std::vector<MeshLib::Element*> const& elements) const
639{
640 for (auto split_node : _split_nodes)
641 {
642 std::size_t const new_node_pos(nodes.size());
643 nodes.push_back(new MeshLib::Node(split_node->data(), new_node_pos));
644
645 // get grid coordinates
646 std::array<std::size_t, 3> const& gc(split_node->getGridCoords());
647 // get affected cells
648 auto const& affected_cells(split_node->getAffectedCells());
649 // get mesh node to substitute in elements
650 MeshLib::Node const* const node2sub(nodes[_index_calculator(gc)]);
651
652 if (affected_cells[0] && gc[0] < _index_calculator._x_dim - 1 &&
653 gc[1] < _index_calculator._y_dim - 1 &&
654 gc[2] < _index_calculator._z_dim - 1)
655 {
656 const std::size_t idx(
657 _index_calculator.getCellIdx(gc[0], gc[1], gc[2]));
658 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
659 }
660 if (affected_cells[1] && gc[0] > 0 &&
661 gc[1] < _index_calculator._y_dim - 1 &&
662 gc[2] < _index_calculator._z_dim - 1)
663 {
664 const std::size_t idx(
665 _index_calculator.getCellIdx(gc[0] - 1, gc[1], gc[2]));
666 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
667 }
668 if (affected_cells[2] && gc[1] > 0 &&
669 gc[0] < _index_calculator._x_dim - 1 &&
670 gc[2] < _index_calculator._z_dim - 1)
671 {
672 const std::size_t idx(
673 _index_calculator.getCellIdx(gc[0], gc[1] - 1, gc[2]));
674 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
675 }
676 if (affected_cells[3] && gc[0] > 0 && gc[1] > 0 &&
677 gc[2] < _index_calculator._z_dim - 1)
678 {
679 const std::size_t idx(
680 _index_calculator.getCellIdx(gc[0] - 1, gc[1] - 1, gc[2]));
681 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
682 }
683 if (affected_cells[4] && gc[2] > 0 &&
684 gc[0] < _index_calculator._x_dim - 1 &&
685 gc[1] < _index_calculator._y_dim - 1)
686 {
687 const std::size_t idx(
688 _index_calculator.getCellIdx(gc[0], gc[1], gc[2] - 1));
689 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
690 }
691 if (affected_cells[5] && gc[0] > 0 && gc[2] > 0 &&
692 gc[1] < _index_calculator._y_dim - 1)
693 {
694 const std::size_t idx(
695 _index_calculator.getCellIdx(gc[0] - 1, gc[1], gc[2] - 1));
696 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
697 }
698 if (affected_cells[6] && gc[1] > 0 && gc[2] > 0 &&
699 gc[0] < _index_calculator._x_dim - 1)
700 {
701 const std::size_t idx(
702 _index_calculator.getCellIdx(gc[0], gc[1] - 1, gc[2] - 1));
703 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
704 }
705 if (affected_cells[7] && gc[0] > 0 && gc[1] > 0 && gc[2] > 0)
706 {
707 const std::size_t idx(
708 _index_calculator.getCellIdx(gc[0] - 1, gc[1] - 1, gc[2] - 1));
709 modifyElement(elements[idx], node2sub, nodes[new_node_pos]);
710 }
711 }
712}
713
715 MeshLib::Node const* const node2sub,
716 MeshLib::Node* const substitute_node)
717{
718 // get the node pointers of the cell
719 MeshLib::Node* const* hex_nodes(hex->getNodes());
720 // search for the position the split node will be set to
721 MeshLib::Node* const* node_pos(
722 std::find(hex_nodes, hex_nodes + 8, node2sub));
723 // set the split node instead of the node2sub
724 if (node_pos != hex_nodes + 8)
725 {
726 const_cast<MeshLib::Node**>(
727 hex_nodes)[std::distance(hex_nodes, node_pos)] = substitute_node;
728 }
729}
730
732 std::string const& name) const
733{
734 auto const it(std::find_if(_property_meta_data_vecs.begin(),
736 [&name](Gocad::Property const& p)
737 { return p._property_name == name; }));
738 if (it == _property_meta_data_vecs.end())
739 {
740 return nullptr;
741 }
742 return &*it;
743}
744
745std::vector<std::string> GocadSGridReader::getPropertyNames() const
746{
747 std::vector<std::string> names;
748 std::transform(_property_meta_data_vecs.begin(),
750 std::back_inserter(names),
751 [](Gocad::Property const& p) { return p._property_name; });
752 return names;
753}
754
755std::unique_ptr<MeshLib::Mesh> GocadSGridReader::getFaceSetMesh(
756 std::size_t const face_set_number) const
757{
758 std::vector<MeshLib::Node*> face_set_nodes;
759 std::vector<MeshLib::Element*> face_set_elements;
760
761 for (auto const node : _nodes)
762 {
763 if (node->isMemberOfFaceSet(face_set_number))
764 {
765 addFaceSetQuad(node, face_set_number, face_set_nodes,
766 face_set_elements);
767 }
768 }
769
770 if (face_set_nodes.empty())
771 {
772 return nullptr;
773 }
774
775 for (auto const node : _split_nodes)
776 {
777 if (node->isMemberOfFaceSet(face_set_number))
778 {
779 if (node->getAffectedCells()[0])
780 {
781 addFaceSetQuad(node, face_set_number, face_set_nodes,
782 face_set_elements);
783 }
784 }
785 }
786
787 std::string const mesh_name("FaceSet-" + std::to_string(face_set_number));
788 return std::make_unique<MeshLib::Mesh>(
789 mesh_name, face_set_nodes, face_set_elements,
790 true /* compute_element_neighbors */);
791}
792
794 GocadNode* face_set_node, std::size_t face_set_number,
795 std::vector<MeshLib::Node*>& face_set_nodes,
796 std::vector<MeshLib::Element*>& face_set_elements) const
797{
798 std::array<MeshLib::Node*, 4> quad_nodes{};
799 quad_nodes[0] = new GocadNode(*face_set_node);
800 const std::size_t id(face_set_node->getID());
801 std::array<std::size_t, 3> const c(_index_calculator.getCoordsForID(id));
802
803 const FaceDirection dir(face_set_node->getFaceDirection(face_set_number));
804 switch (dir)
805 {
806 case FaceDirection::U:
807 quad_nodes[1] = new GocadNode(
808 *_nodes[_index_calculator({c[0], c[1] + 1, c[2]})]);
809 quad_nodes[2] = new GocadNode(
810 *_nodes[_index_calculator({c[0], c[1] + 1, c[2] + 1})]);
811 quad_nodes[3] = new GocadNode(
812 *_nodes[_index_calculator({c[0], c[1], c[2] + 1})]);
813 break;
814 case FaceDirection::V:
815 quad_nodes[1] = new GocadNode(
816 *_nodes[_index_calculator({c[0] + 1, c[1], c[2]})]);
817 quad_nodes[2] = new GocadNode(
818 *_nodes[_index_calculator({c[0] + 1, c[1], c[2] + 1})]);
819 quad_nodes[3] = new GocadNode(
820 *_nodes[_index_calculator({c[0], c[1], c[2] + 1})]);
821 break;
822 case FaceDirection::W:
823 quad_nodes[1] = new GocadNode(
824 *_nodes[_index_calculator({c[0] + 1, c[1], c[2]})]);
825 quad_nodes[2] = new GocadNode(
826 *_nodes[_index_calculator({c[0] + 1, c[1] + 1, c[2]})]);
827 quad_nodes[3] = new GocadNode(
828 *_nodes[_index_calculator({c[0], c[1] + 1, c[2]})]);
829 break;
830 default:
831 ERR("Could not create face for node with id {:d}.", id);
832 }
833 std::copy(begin(quad_nodes), end(quad_nodes),
834 back_inserter(face_set_nodes));
835 face_set_elements.push_back(new MeshLib::Quad(quad_nodes));
836}
837
838} // end namespace Gocad
839} // 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