13 #include <boost/tokenizer.hpp>
28 using Bitset = boost::dynamic_bitset<>;
34 _double_precision_binary(false),
35 _bin_pnts_in_double_precision(false)
38 std::ifstream in(
_fname.c_str());
46 bool pnts_read(
false);
50 while (std::getline(in, line))
56 if (line.back() ==
'\r')
60 if (line.substr(0, 8) ==
"HEADER {")
64 if (line.substr(0, 32) ==
"GOCAD_ORIGINAL_COORDINATE_SYSTEM")
69 if (line.substr(0, 7) ==
"AXIS_N ")
73 else if (line.substr(0, 12) ==
"POINTS_FILE ")
77 else if (line.substr(0, 9) ==
"PROPERTY ")
82 else if (line.substr(0, 35) ==
"BINARY_POINTS_IN_DOUBLE_PRECISION 1")
86 else if (line.substr(0, 11) ==
"FLAGS_FILE ")
90 else if (line.substr(0, 18) ==
"REGION_FLAGS_FILE ")
94 else if (line.substr(0, 7) ==
"REGION " ||
95 line.substr(0, 13) ==
"MODEL_REGION ")
99 else if (line.substr(0, 12) ==
"MODEL_LAYER ")
103 else if (line.substr(0, 24) ==
"REGION_FLAGS_BIT_LENGTH ")
105 std::istringstream iss(line);
106 std::istream_iterator<std::string> it(iss);
108 std::size_t bit_length = std::atoi(it->c_str());
109 if (
regions.size() != bit_length)
111 ERR(
"{:d} regions read but {:d} expected.\n",
regions.size(),
113 throw std::runtime_error(
114 "Number of read regions differs from expected.\n");
117 else if (line.substr(0, 9) ==
"FACE_SET ")
130 std::stringstream regions_ss;
131 regions_ss <<
regions.size() <<
" regions read:\n";
133 std::ostream_iterator<Gocad::Region>(regions_ss,
"\t"));
134 DBUG(
"{:s}", regions_ss.str());
136 std::stringstream layers_ss;
137 layers_ss <<
layers.size() <<
" layers read:\n";
139 std::ostream_iterator<Gocad::Layer>(layers_ss,
"\n"));
140 DBUG(
"{:s}", layers_ss.str());
142 std::stringstream properties_ss;
144 <<
" properties read:\n";
146 std::ostream_iterator<Gocad::Property>(properties_ss,
"\n"));
147 DBUG(
"{:s}", properties_ss.str());
178 std::vector<MeshLib::Node*> nodes;
179 std::transform(
_nodes.cbegin(),
_nodes.cend(), std::back_inserter(nodes),
181 { return new MeshLib::Node(*node); });
186 DBUG(
"Creating mesh from Gocad SGrid.");
190 DBUG(
"Mesh created.");
198 for (
auto const&
name : prop_names)
206 DBUG(
"Adding Gocad property '{:s}' with {:d} values.",
name,
207 prop->_property_data.size());
209 auto pv = MeshLib::getOrCreateMeshProperty<double>(
213 ERR(
"Could not create mesh property '{:s}'.",
name);
217 pv->resize(prop->_property_data.size());
218 std::copy(prop->_property_data.cbegin(), prop->_property_data.cend(),
226 while (std::getline(in, line))
228 if (line.front() ==
'}')
232 if (line.substr(0, 27) ==
"double_precision_binary: on")
240 "GocadSGridReader::parseHeader(): _double_precision_binary == "
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());
254 std::stringstream ssx(*(it),
255 std::stringstream::in | std::stringstream::out);
258 std::stringstream ssy(*it, std::stringstream::in | std::stringstream::out);
261 std::stringstream ssz(*it, std::stringstream::in | std::stringstream::out);
265 "x_dim = {:d}, y_dim = {:d}, z_dim = {:d} => #nodes = {:d}, #cells = "
272 std::string& result_string)
const
274 boost::char_separator<char> sep(
" \r");
275 boost::tokenizer<boost::char_separator<char>> tok(line, sep);
276 auto it(tok.begin());
278 result_string =
_path + *it;
301 std::istringstream iss(line);
302 std::istream_iterator<std::string> it(iss);
304 if (*it != std::string(
"FACE_SET"))
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 "
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);
318 while (std::getline(in, line) && face_set_id_cnt < n_of_face_set_ids)
320 if (line.back() ==
'\r')
324 boost::char_separator<char> sep(
"\t ");
325 boost::tokenizer<boost::char_separator<char>> tokens(line, sep);
327 for (
auto tok_it = tokens.begin(); tok_it != tokens.end();)
329 auto const id(
static_cast<std::size_t
>(std::atoi(tok_it->c_str())));
331 auto const face_indicator(
332 static_cast<std::size_t
>(std::atoi(tok_it->c_str())));
337 ERR(
"Face set id {:d} is greater than the number of nodes "
345 std::array<std::size_t, 3>
const c(
349 ERR(
"****** i coord {:d} to big for id {:d}.",
c[0],
id);
353 ERR(
"****** j coord {:d} to big for id {:d}.",
c[1],
id);
357 ERR(
"****** k coord {:d} to big for id {:d}.",
c[2],
id);
359 std::size_t
const cell_id(
362 static_cast<double>(face_indicator);
368 if (face_set_id_cnt != n_of_face_set_ids)
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 "
382 split_node->transmitFaceDirections(*
_nodes[
id]);
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);
398 std::vector<block_t> data;
400 std::fill_n(data.data(), blocks, 0);
401 in.read(
reinterpret_cast<char*
>(data.data()), bytes);
403 return Bitset(data.begin(), data.end());
408 std::ifstream in(
_pnts_fname.c_str(), std::ios::binary);
412 throw std::runtime_error(
"Could not open points file.");
421 while (in && k < n * 3)
433 if ((k + 1) % 3 == 0)
435 const std::size_t layer_transition_idx(
441 if (k != n * 3 && !in.eof())
443 ERR(
"Read different number of points. Expected {:d} floats, got "
450 std::vector<Bitset>
const& rf)
453 "GocadSGridReader::mapRegionFlagsToCellProperties region_flags.size: "
478 std::size_t
const cell_id(
483 if (rf[node_id].test(region.bit))
499 std::string
const& fname(property._property_data_fname);
500 if (property._property_data_fname.empty())
502 WARN(
"Empty filename for property {:s}.", property._property_name);
505 std::vector<float> float_properties =
508 "GocadSGridReader::readElementPropertiesBinary(): Read {:d} float "
509 "properties from binary file.",
512 std::transform(float_properties.cbegin(), float_properties.cend(),
513 float_properties.begin(),
515 { return BaseLib::swapEndianness(val); });
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())
522 ERR(
"Reading of element properties file '{:s}' failed.", fname);
529 std::vector<Bitset> result;
534 ERR(
"readRegionFlagsBinary(): Could not open file '{:s}' for input.\n",
548 if (k != n && !in.eof())
550 ERR(
"Read different number of values. Expected {:d}, got {:d}.\n", n,
556 std::vector<MeshLib::Node*>
const& nodes)
const
558 std::vector<MeshLib::Element*> elements;
560 std::array<MeshLib::Node*, 8> element_nodes{};
588 std::ifstream in(
_fname.c_str());
598 std::stringstream ss;
599 while (std::getline(in, line))
601 std::size_t pos(line.find(
"SPLIT "));
602 if (pos != std::string::npos)
604 ss << line.substr(pos + 6, line.size() - (pos + 6));
606 std::array<std::size_t, 3> grid_coords{};
607 ss >> grid_coords[0];
608 ss >> grid_coords[1];
609 ss >> grid_coords[2];
619 std::bitset<8> affected_cells{};
620 for (std::size_t ac = 0; ac < 8; ++ac)
624 affected_cells[ac] = bit !=
'0';
626 const std::size_t layer_transition_index(
627 _nodes[
id]->getLayerTransitionIndex());
630 layer_transition_index));
636 std::vector<MeshLib::Node*>& nodes,
637 std::vector<MeshLib::Element*>
const& elements)
const
641 std::size_t
const new_node_pos(nodes.size());
646 std::array<std::size_t, 3>
const& gc(split_node->getGridCoords());
648 auto const& affected_cells(split_node->getAffectedCells());
656 const std::size_t idx(
660 if (affected_cells[1] && gc[0] > 0 &&
664 const std::size_t idx(
668 if (affected_cells[2] && gc[1] > 0 &&
672 const std::size_t idx(
676 if (affected_cells[3] && gc[0] > 0 && gc[1] > 0 &&
679 const std::size_t idx(
683 if (affected_cells[4] && gc[2] > 0 &&
687 const std::size_t idx(
691 if (affected_cells[5] && gc[0] > 0 && gc[2] > 0 &&
694 const std::size_t idx(
698 if (affected_cells[6] && gc[1] > 0 && gc[2] > 0 &&
701 const std::size_t idx(
705 if (affected_cells[7] && gc[0] > 0 && gc[1] > 0 && gc[2] > 0)
707 const std::size_t idx(
722 std::find(hex_nodes, hex_nodes + 8, node2sub));
724 if (node_pos != hex_nodes + 8)
727 hex_nodes)[std::distance(hex_nodes, node_pos)] = substitute_node;
732 std::string
const&
name)
const
737 { return p._property_name == name; }));
747 std::vector<std::string> names;
750 std::back_inserter(names),
756 std::size_t
const face_set_number)
const
758 std::vector<MeshLib::Node*> face_set_nodes;
759 std::vector<MeshLib::Element*> face_set_elements;
761 for (
auto const node :
_nodes)
763 if (node->isMemberOfFaceSet(face_set_number))
770 if (face_set_nodes.empty())
777 if (node->isMemberOfFaceSet(face_set_number))
779 if (node->getAffectedCells()[0])
787 std::string
const mesh_name(
"FaceSet-" + std::to_string(face_set_number));
788 return std::make_unique<MeshLib::Mesh>(mesh_name, face_set_nodes,
793 GocadNode* face_set_node, std::size_t face_set_number,
794 std::vector<MeshLib::Node*>& face_set_nodes,
795 std::vector<MeshLib::Element*>& face_set_elements)
const
797 std::array<MeshLib::Node*, 4> quad_nodes{};
798 quad_nodes[0] =
new GocadNode(*face_set_node);
799 const std::size_t id(face_set_node->
getID());
830 ERR(
"Could not create face for node with id {:d}.",
id);
832 std::copy(begin(quad_nodes), end(quad_nodes),
833 back_inserter(face_set_nodes));
Definition of the Hex class.
void ERR(char const *fmt, Args const &... args)
void DBUG(char const *fmt, Args const &... args)
void WARN(char const *fmt, Args const &... 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
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
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
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.
std::string extractPath(std::string const &pathname)
std::string extractBaseNameWithoutExtension(std::string const &pathname)
double swapEndianness(double const &v)
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)
boost::dynamic_bitset<> Bitset
void copy(PETScVector const &x, PETScVector &y)
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