13#include <boost/tokenizer.hpp>
29using Bitset = boost::dynamic_bitset<>;
33 _path(
BaseLib::extractPath(fname)),
35 _double_precision_binary(false),
36 _bin_pnts_in_double_precision(false)
39 std::ifstream in(
_fname.c_str());
47 bool pnts_read(
false);
51 while (std::getline(in, line))
57 if (line.back() ==
'\r')
61 if (line.substr(0, 8) ==
"HEADER {")
65 if (line.substr(0, 32) ==
"GOCAD_ORIGINAL_COORDINATE_SYSTEM")
70 if (line.substr(0, 7) ==
"AXIS_N ")
74 else if (line.substr(0, 12) ==
"POINTS_FILE ")
78 else if (line.substr(0, 9) ==
"PROPERTY ")
83 else if (line.substr(0, 35) ==
"BINARY_POINTS_IN_DOUBLE_PRECISION 1")
87 else if (line.substr(0, 11) ==
"FLAGS_FILE ")
91 else if (line.substr(0, 18) ==
"REGION_FLAGS_FILE ")
95 else if (line.substr(0, 7) ==
"REGION " ||
96 line.substr(0, 13) ==
"MODEL_REGION ")
100 else if (line.substr(0, 12) ==
"MODEL_LAYER ")
104 else if (line.substr(0, 24) ==
"REGION_FLAGS_BIT_LENGTH ")
106 std::istringstream iss(line);
107 std::istream_iterator<std::string> it(iss);
109 std::size_t bit_length = std::atoi(it->c_str());
110 if (
regions.size() != bit_length)
112 ERR(
"{:d} regions read but {:d} expected.\n",
regions.size(),
114 throw std::runtime_error(
115 "Number of read regions differs from expected.\n");
118 else if (line.substr(0, 9) ==
"FACE_SET ")
131 std::stringstream regions_ss;
132 regions_ss <<
regions.size() <<
" regions read:\n";
134 std::ostream_iterator<Gocad::Region>(regions_ss,
"\t"));
135 DBUG(
"{:s}", regions_ss.str());
137 std::stringstream layers_ss;
138 layers_ss <<
layers.size() <<
" layers read:\n";
140 std::ostream_iterator<Gocad::Layer>(layers_ss,
"\n"));
141 DBUG(
"{:s}", layers_ss.str());
143 std::stringstream properties_ss;
145 <<
" properties read:\n";
147 std::ostream_iterator<Gocad::Property>(properties_ss,
"\n"));
148 DBUG(
"{:s}", properties_ss.str());
179 std::vector<MeshLib::Node*> nodes;
180 std::transform(
_nodes.cbegin(),
_nodes.cend(), std::back_inserter(nodes),
182 { return new MeshLib::Node(*node); });
187 DBUG(
"Creating mesh from Gocad SGrid.");
192 DBUG(
"Mesh created.");
200 for (
auto const& name : prop_names)
208 DBUG(
"Adding Gocad property '{:s}' with {:d} values.", name,
209 prop->_property_data.size());
211 auto pv = MeshLib::getOrCreateMeshProperty<double>(
215 ERR(
"Could not create mesh property '{:s}'.", name);
219 pv->resize(prop->_property_data.size());
220 std::copy(prop->_property_data.cbegin(), prop->_property_data.cend(),
228 while (std::getline(in, line))
230 if (line.front() ==
'}')
234 if (line.substr(0, 27) ==
"double_precision_binary: on")
242 "GocadSGridReader::parseHeader(): _double_precision_binary == "
249 std::size_t x_dim(0);
250 std::size_t y_dim(0);
251 std::size_t z_dim(0);
252 boost::tokenizer<> tok(line);
253 auto it(tok.begin());
256 std::stringstream ssx(*(it),
257 std::stringstream::in | std::stringstream::out);
260 std::stringstream ssy(*it, std::stringstream::in | std::stringstream::out);
263 std::stringstream ssz(*it, std::stringstream::in | std::stringstream::out);
267 "x_dim = {:d}, y_dim = {:d}, z_dim = {:d} => #nodes = {:d}, #cells = "
274 std::string& result_string)
const
276 boost::char_separator<char> sep(
" \r");
277 boost::tokenizer<boost::char_separator<char>> tok(line, sep);
278 auto it(tok.begin());
280 result_string =
_path + *it;
303 std::istringstream iss(line);
304 std::istream_iterator<std::string> it(iss);
306 if (*it != std::string(
"FACE_SET"))
308 ERR(
"Expected FACE_SET keyword but '{:s}' found.", it->c_str());
309 throw std::runtime_error(
310 "In GocadSGridReader::parseFaceSet() expected FACE_SET keyword not "
316 auto const n_of_face_set_ids(
317 static_cast<std::size_t
>(std::atoi(it->c_str())));
318 std::size_t face_set_id_cnt(0);
320 while (std::getline(in, line) && face_set_id_cnt < n_of_face_set_ids)
322 if (line.back() ==
'\r')
326 boost::char_separator<char> sep(
"\t ");
327 boost::tokenizer<boost::char_separator<char>> tokens(line, sep);
329 for (
auto tok_it = tokens.begin(); tok_it != tokens.end();)
331 auto const id(
static_cast<std::size_t
>(std::atoi(tok_it->c_str())));
333 auto const face_indicator(
334 static_cast<std::size_t
>(std::atoi(tok_it->c_str())));
339 ERR(
"Face set id {:d} is greater than the number of nodes "
347 std::array<std::size_t, 3>
const c(
351 ERR(
"****** i coord {:d} to big for id {:d}.", c[0],
id);
355 ERR(
"****** j coord {:d} to big for id {:d}.", c[1],
id);
359 ERR(
"****** k coord {:d} to big for id {:d}.", c[2],
id);
361 std::size_t
const cell_id(
364 static_cast<double>(face_indicator);
370 if (face_set_id_cnt != n_of_face_set_ids)
372 ERR(
"Expected {:d} number of face set ids, read {:d}.",
373 n_of_face_set_ids, face_set_id_cnt);
374 throw std::runtime_error(
375 "Expected number of face set points does not match number of read "
384 split_node->transmitFaceDirections(*
_nodes[
id]);
395 using block_t = Bitset::block_type;
396 auto const bytes =
static_cast<std::size_t
>(std::ceil(bits / 8.));
397 std::size_t
const blocks =
398 bytes + 1 <
sizeof(block_t) ? 1 : (bytes + 1) /
sizeof(block_t);
400 std::vector<block_t> data;
402 std::fill_n(data.data(), blocks, 0);
403 in.read(
reinterpret_cast<char*
>(data.data()), bytes);
405 return Bitset(data.begin(), data.end());
410 std::ifstream in(
_pnts_fname.c_str(), std::ios::binary);
414 throw std::runtime_error(
"Could not open points file.");
423 while (in && k < n * 3)
435 if ((k + 1) % 3 == 0)
437 const std::size_t layer_transition_idx(
443 if (k != n * 3 && !in.eof())
445 ERR(
"Read different number of points. Expected {:d} floats, got "
452 std::vector<Bitset>
const& rf)
455 "GocadSGridReader::mapRegionFlagsToCellProperties region_flags.size: "
480 std::size_t
const cell_id(
485 if (rf[node_id].test(region.bit))
501 std::string
const& fname(property._property_data_fname);
502 if (property._property_data_fname.empty())
504 WARN(
"Empty filename for property {:s}.", property._property_name);
510 "GocadSGridReader::readElementPropertiesBinary(): Read {:d} float "
511 "properties from binary file.",
514 std::transform(float_properties.cbegin(), float_properties.cend(),
515 float_properties.begin(),
517 { return BaseLib::swapEndianness(val); });
519 property._property_data.resize(float_properties.size());
520 std::copy(float_properties.begin(), float_properties.end(),
521 property._property_data.begin());
522 if (property._property_data.empty())
524 ERR(
"Reading of element properties file '{:s}' failed.", fname);
531 std::vector<Bitset> result;
536 ERR(
"readRegionFlagsBinary(): Could not open file '{:s}' for input.\n",
550 if (k != n && !in.eof())
552 ERR(
"Read different number of values. Expected {:d}, got {:d}.\n", n,
558 std::vector<MeshLib::Node*>
const& nodes)
const
560 std::vector<MeshLib::Element*> elements;
562 std::array<MeshLib::Node*, 8> element_nodes{};
590 std::ifstream in(
_fname.c_str());
600 std::stringstream ss;
601 while (std::getline(in, line))
603 std::size_t pos(line.find(
"SPLIT "));
604 if (pos != std::string::npos)
606 ss << line.substr(pos + 6, line.size() - (pos + 6));
608 std::array<std::size_t, 3> grid_coords{};
609 ss >> grid_coords[0];
610 ss >> grid_coords[1];
611 ss >> grid_coords[2];
621 std::bitset<8> affected_cells{};
622 for (std::size_t ac = 0; ac < 8; ++ac)
626 affected_cells[ac] = bit !=
'0';
628 const std::size_t layer_transition_index(
629 _nodes[
id]->getLayerTransitionIndex());
632 layer_transition_index));
638 std::vector<MeshLib::Node*>& nodes,
639 std::vector<MeshLib::Element*>
const& elements)
const
643 std::size_t
const new_node_pos(nodes.size());
644 nodes.push_back(
new MeshLib::Node(split_node->data(), new_node_pos));
647 std::array<std::size_t, 3>
const& gc(split_node->getGridCoords());
649 auto const& affected_cells(split_node->getAffectedCells());
657 const std::size_t idx(
661 if (affected_cells[1] && gc[0] > 0 &&
665 const std::size_t idx(
669 if (affected_cells[2] && gc[1] > 0 &&
673 const std::size_t idx(
677 if (affected_cells[3] && gc[0] > 0 && gc[1] > 0 &&
680 const std::size_t idx(
684 if (affected_cells[4] && gc[2] > 0 &&
688 const std::size_t idx(
692 if (affected_cells[5] && gc[0] > 0 && gc[2] > 0 &&
695 const std::size_t idx(
699 if (affected_cells[6] && gc[1] > 0 && gc[2] > 0 &&
702 const std::size_t idx(
706 if (affected_cells[7] && gc[0] > 0 && gc[1] > 0 && gc[2] > 0)
708 const std::size_t idx(
723 std::find(hex_nodes, hex_nodes + 8, node2sub));
725 if (node_pos != hex_nodes + 8)
728 hex_nodes)[std::distance(hex_nodes, node_pos)] = substitute_node;
733 std::string
const& name)
const
738 { return p._property_name == name; }));
748 std::vector<std::string> names;
751 std::back_inserter(names),
757 std::size_t
const face_set_number)
const
759 std::vector<MeshLib::Node*> face_set_nodes;
760 std::vector<MeshLib::Element*> face_set_elements;
762 for (
auto const node :
_nodes)
764 if (node->isMemberOfFaceSet(face_set_number))
771 if (face_set_nodes.empty())
778 if (node->isMemberOfFaceSet(face_set_number))
780 if (node->getAffectedCells()[0])
788 std::string
const mesh_name(
"FaceSet-" + std::to_string(face_set_number));
789 return std::make_unique<MeshLib::Mesh>(
790 mesh_name, face_set_nodes, face_set_elements,
795 GocadNode* face_set_node, std::size_t face_set_number,
796 std::vector<MeshLib::Node*>& face_set_nodes,
797 std::vector<MeshLib::Element*>& face_set_elements)
const
799 std::array<MeshLib::Node*, 4> quad_nodes{};
800 quad_nodes[0] =
new GocadNode(*face_set_node);
801 const std::size_t id(face_set_node->
getID());
832 ERR(
"Could not create face for node with id {:d}.",
id);
834 std::copy(begin(quad_nodes), end(quad_nodes),
835 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)
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