13#include <boost/tokenizer.hpp>
29using Bitset = boost::dynamic_bitset<>;
33 _path(
BaseLib::extractPath(fname).empty()
35 :
BaseLib::extractPath(fname) +
'/'),
37 _double_precision_binary(false),
38 _bin_pnts_in_double_precision(false)
41 std::ifstream in(
_fname.c_str());
49 bool pnts_read(
false);
53 while (std::getline(in, line))
59 if (line.back() ==
'\r')
63 if (line.substr(0, 8) ==
"HEADER {")
67 if (line.substr(0, 32) ==
"GOCAD_ORIGINAL_COORDINATE_SYSTEM")
72 if (line.substr(0, 7) ==
"AXIS_N ")
76 else if (line.substr(0, 12) ==
"POINTS_FILE ")
80 else if (line.substr(0, 9) ==
"PROPERTY ")
85 else if (line.substr(0, 35) ==
"BINARY_POINTS_IN_DOUBLE_PRECISION 1")
89 else if (line.substr(0, 11) ==
"FLAGS_FILE ")
93 else if (line.substr(0, 18) ==
"REGION_FLAGS_FILE ")
97 else if (line.substr(0, 7) ==
"REGION " ||
98 line.substr(0, 13) ==
"MODEL_REGION ")
102 else if (line.substr(0, 12) ==
"MODEL_LAYER ")
106 else if (line.substr(0, 24) ==
"REGION_FLAGS_BIT_LENGTH ")
108 std::istringstream iss(line);
109 std::istream_iterator<std::string> it(iss);
111 std::size_t bit_length = std::atoi(it->c_str());
112 if (
regions.size() != bit_length)
114 ERR(
"{:d} regions read but {:d} expected.\n",
regions.size(),
116 throw std::runtime_error(
117 "Number of read regions differs from expected.\n");
120 else if (line.substr(0, 9) ==
"FACE_SET ")
133 std::stringstream regions_ss;
134 regions_ss <<
regions.size() <<
" regions read:\n";
136 std::ostream_iterator<Gocad::Region>(regions_ss,
"\t"));
137 DBUG(
"{:s}", regions_ss.str());
139 std::stringstream layers_ss;
140 layers_ss <<
layers.size() <<
" layers read:\n";
142 std::ostream_iterator<Gocad::Layer>(layers_ss,
"\n"));
143 DBUG(
"{:s}", layers_ss.str());
145 std::stringstream properties_ss;
147 <<
" properties read:\n";
149 std::ostream_iterator<Gocad::Property>(properties_ss,
"\n"));
150 DBUG(
"{:s}", properties_ss.str());
181 std::vector<MeshLib::Node*> nodes;
182 std::transform(
_nodes.cbegin(),
_nodes.cend(), std::back_inserter(nodes),
184 { return new MeshLib::Node(*node); });
189 DBUG(
"Creating mesh from Gocad SGrid.");
194 DBUG(
"Mesh created.");
202 for (
auto const& name : prop_names)
210 DBUG(
"Adding Gocad property '{:s}' with {:d} values.", name,
211 prop->_property_data.size());
217 ERR(
"Could not create mesh property '{:s}'.", name);
221 pv->resize(prop->_property_data.size());
222 std::copy(prop->_property_data.cbegin(), prop->_property_data.cend(),
230 while (std::getline(in, line))
232 if (line.front() ==
'}')
236 if (line.substr(0, 27) ==
"double_precision_binary: on")
244 "GocadSGridReader::parseHeader(): _double_precision_binary == "
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());
258 std::stringstream ssx(*(it),
259 std::stringstream::in | std::stringstream::out);
262 std::stringstream ssy(*it, std::stringstream::in | std::stringstream::out);
265 std::stringstream ssz(*it, std::stringstream::in | std::stringstream::out);
269 "x_dim = {:d}, y_dim = {:d}, z_dim = {:d} => #nodes = {:d}, #cells = "
276 std::string& result_string)
const
278 boost::char_separator<char> sep(
" \r");
279 boost::tokenizer<boost::char_separator<char>> tok(line, sep);
280 auto it(tok.begin());
282 result_string =
_path + *it;
305 std::istringstream iss(line);
306 std::istream_iterator<std::string> it(iss);
308 if (*it != std::string(
"FACE_SET"))
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 "
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);
322 while (std::getline(in, line) && face_set_id_cnt < n_of_face_set_ids)
324 if (line.back() ==
'\r')
328 boost::char_separator<char> sep(
"\t ");
329 boost::tokenizer<boost::char_separator<char>> tokens(line, sep);
331 for (
auto tok_it = tokens.begin(); tok_it != tokens.end();)
333 auto const id(
static_cast<std::size_t
>(std::atoi(tok_it->c_str())));
335 auto const face_indicator(
336 static_cast<std::size_t
>(std::atoi(tok_it->c_str())));
341 ERR(
"Face set id {:d} is greater than the number of nodes "
349 std::array<std::size_t, 3>
const c(
353 ERR(
"****** i coord {:d} to big for id {:d}.", c[0],
id);
357 ERR(
"****** j coord {:d} to big for id {:d}.", c[1],
id);
361 ERR(
"****** k coord {:d} to big for id {:d}.", c[2],
id);
363 std::size_t
const cell_id(
366 static_cast<double>(face_indicator);
372 if (face_set_id_cnt != n_of_face_set_ids)
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 "
386 split_node->transmitFaceDirections(*
_nodes[
id]);
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);
402 std::vector<block_t> data;
404 std::fill_n(data.data(), blocks, 0);
405 in.read(
reinterpret_cast<char*
>(data.data()), bytes);
407 return Bitset(data.begin(), data.end());
412 std::ifstream in(
_pnts_fname.c_str(), std::ios::binary);
416 throw std::runtime_error(
"Could not open points file.");
425 while (in && k < n * 3)
437 if ((k + 1) % 3 == 0)
439 const std::size_t layer_transition_idx(
445 if (k != n * 3 && !in.eof())
447 ERR(
"Read different number of points. Expected {:d} floats, got "
454 std::vector<Bitset>
const& rf)
457 "GocadSGridReader::mapRegionFlagsToCellProperties region_flags.size: "
482 std::size_t
const cell_id(
487 if (rf[node_id].test(region.bit))
503 std::string
const& fname(property._property_data_fname);
504 if (property._property_data_fname.empty())
506 WARN(
"Empty filename for property {:s}.", property._property_name);
512 "GocadSGridReader::readElementPropertiesBinary(): Read {:d} float "
513 "properties from binary file.",
516 std::transform(float_properties.cbegin(), float_properties.cend(),
517 float_properties.begin(),
519 { return BaseLib::swapEndianness(val); });
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())
526 ERR(
"Reading of element properties file '{:s}' failed.", fname);
533 std::vector<Bitset> result;
538 ERR(
"readRegionFlagsBinary(): Could not open file '{:s}' for input.\n",
552 if (k != n && !in.eof())
554 ERR(
"Read different number of values. Expected {:d}, got {:d}.\n", n,
560 std::vector<MeshLib::Node*>
const& nodes)
const
562 std::vector<MeshLib::Element*> elements;
564 std::array<MeshLib::Node*, 8> element_nodes{};
592 std::ifstream in(
_fname.c_str());
602 std::stringstream ss;
603 while (std::getline(in, line))
605 std::size_t pos(line.find(
"SPLIT "));
606 if (pos != std::string::npos)
608 ss << line.substr(pos + 6, line.size() - (pos + 6));
610 std::array<std::size_t, 3> grid_coords{};
611 ss >> grid_coords[0];
612 ss >> grid_coords[1];
613 ss >> grid_coords[2];
623 std::bitset<8> affected_cells{};
624 for (std::size_t ac = 0; ac < 8; ++ac)
628 affected_cells[ac] = bit !=
'0';
630 const std::size_t layer_transition_index(
631 _nodes[
id]->getLayerTransitionIndex());
634 layer_transition_index));
640 std::vector<MeshLib::Node*>& nodes,
641 std::vector<MeshLib::Element*>
const& elements)
const
645 std::size_t
const new_node_pos(nodes.size());
646 nodes.push_back(
new MeshLib::Node(split_node->data(), new_node_pos));
649 std::array<std::size_t, 3>
const& gc(split_node->getGridCoords());
651 auto const& affected_cells(split_node->getAffectedCells());
659 const std::size_t idx(
663 if (affected_cells[1] && gc[0] > 0 &&
667 const std::size_t idx(
671 if (affected_cells[2] && gc[1] > 0 &&
675 const std::size_t idx(
679 if (affected_cells[3] && gc[0] > 0 && gc[1] > 0 &&
682 const std::size_t idx(
686 if (affected_cells[4] && gc[2] > 0 &&
690 const std::size_t idx(
694 if (affected_cells[5] && gc[0] > 0 && gc[2] > 0 &&
697 const std::size_t idx(
701 if (affected_cells[6] && gc[1] > 0 && gc[2] > 0 &&
704 const std::size_t idx(
708 if (affected_cells[7] && gc[0] > 0 && gc[1] > 0 && gc[2] > 0)
710 const std::size_t idx(
725 std::find(hex_nodes, hex_nodes + 8, node2sub));
727 if (node_pos != hex_nodes + 8)
730 hex_nodes)[std::distance(hex_nodes, node_pos)] = substitute_node;
735 std::string
const& name)
const
740 { return p._property_name == name; }));
750 std::vector<std::string> names;
753 std::back_inserter(names),
759 std::size_t
const face_set_number)
const
761 std::vector<MeshLib::Node*> face_set_nodes;
762 std::vector<MeshLib::Element*> face_set_elements;
764 for (
auto const node :
_nodes)
766 if (node->isMemberOfFaceSet(face_set_number))
773 if (face_set_nodes.empty())
780 if (node->isMemberOfFaceSet(face_set_number))
782 if (node->getAffectedCells()[0])
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,
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
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());
834 ERR(
"Could not create face for node with id {:d}.",
id);
836 std::copy(begin(quad_nodes), end(quad_nodes),
837 back_inserter(face_set_nodes));
Definition of the Hex class.
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition of the Mesh class.
Definition of the Quad class.
bool parse(std::istream &in)
FaceDirection getFaceDirection(std::size_t const face_set_number) const
Gocad::Property const * getProperty(std::string const &name) const
std::vector< Gocad::Property > _property_meta_data_vecs
std::string const & _fname
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::string _region_flags_fname
std::unique_ptr< MeshLib::Mesh > getMesh() const
GocadSGridReader()=delete
void readElementPropertiesBinary()
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 readSplitInformation()
void parseFaceSet(std::string &line, std::istream &in)
void parseFileName(std::string const &line, std::string &result_string) const
bool _bin_pnts_in_double_precision
Gocad::IndexCalculator _index_calculator
std::vector< GocadSplitNode * > _split_nodes
void parseDims(std::string const &line)
std::vector< std::string > getPropertyNames() const
bool _double_precision_binary
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
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)
Bitset readBits(std::ifstream &in, const std::size_t bits)
Layer parseLayer(std::string const &line, std::vector< Region > const ®ions)
Region parseRegion(std::string const &line)
PropertyVector< T > * getOrCreateMeshProperty(Mesh &mesh, std::string const &property_name, MeshItemType const item_type, int const number_of_components)
TemplateElement< MeshLib::HexRule8 > Hex
double _property_no_data_value
std::string _property_data_fname
std::string _property_class_name
std::string _property_unit
std::string _property_name
std::string _property_data_type
std::vector< double > _property_data