Flow123d  JS_before_hm-1002-gafa1d04
msh_gmshreader.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file msh_gmshreader.cc
15  * @ingroup mesh
16  * @brief
17  * @author dalibor
18  */
19 
20 #include <istream>
21 #include <string>
22 #include <limits>
23 
24 #include "msh_gmshreader.h"
26 
27 #include "system/system.hh"
28 #include "system/tokenizer.hh"
29 #include "boost/lexical_cast.hpp"
30 
31 #include "mesh/mesh.h"
32 
33 
34 
35 using namespace std;
36 
37 
39 : BaseMeshReader(file_name)
40 {
41  tok_.set_comment_pattern( "#");
42  data_section_name_ = "$ElementData";
43  has_compatible_mesh_ = false;
45 }
46 
47 
48 
49 GmshMeshReader::~GmshMeshReader() // Tokenizer close the file automatically
50 {}
51 
52 
53 
55  using namespace boost;
56  unsigned int n_nodes;
57  MessageOut() << "- Reading nodes...";
58  tok_.set_position( Tokenizer::Position() );
59 
60  if (! tok_.skip_to("$Nodes")) THROW(ExcMissingSection() << EI_Section("$Nodes") << EI_GMSHFile(tok_.f_name()) );
61  try {
62  tok_.next_line(false);
63  n_nodes = lexical_cast<unsigned int> (*tok_);
64  mesh->init_node_vector( n_nodes );
65  INPUT_CHECK( n_nodes > 0, "Zero number of nodes, %s.\n", tok_.position_msg().c_str() );
66  ++tok_; // end of line
67 
68  for (unsigned int i = 0; i < n_nodes; ++i) {
69  tok_.next_line();
70 
71  unsigned int id = lexical_cast<unsigned int> (*tok_); ++tok_; // node id
72  arma::vec3 coords; // node coordinates
73  coords(0) = lexical_cast<double> (*tok_); ++tok_;
74  coords(1) = lexical_cast<double> (*tok_); ++tok_;
75  coords(2) = lexical_cast<double> (*tok_); ++tok_;
76  ++tok_; // skip mesh size parameter
77 
78  mesh->add_node(id, coords);
79  }
80 
81  } catch (bad_lexical_cast &) {
82  THROW(ExcWrongFormat() << EI_Type("number") << EI_TokenizerMsg(tok_.position_msg()) << EI_MeshFile(tok_.f_name()) );
83  }
84  MessageOut().fmt("... {} nodes read. \n", n_nodes);
85 }
86 
87 
89  using namespace boost;
90  MessageOut() << "- Reading elements...";
91 
92  if (! tok_.skip_to("$Elements")) THROW(ExcMissingSection() << EI_Section("$Elements") << EI_GMSHFile(tok_.f_name()) );
93  try {
94  tok_.next_line(false);
95  unsigned int n_elements = lexical_cast<unsigned int> (*tok_);
96  INPUT_CHECK( n_elements > 0, "Zero number of elements, %s.\n", tok_.position_msg().c_str());
97  ++tok_; // end of line
98 
99  std::vector<unsigned int> node_ids; //node_ids of elements
100  node_ids.resize(4); // maximal count of nodes
101 
102  mesh->init_element_vector(n_elements);
103 
104  for (unsigned int i = 0; i < n_elements; ++i) {
105  tok_.next_line();
106  unsigned int id = lexical_cast<unsigned int>(*tok_); ++tok_;
107 
108  //get element type: supported:
109  // 1 Line (2 nodes)
110  // 2 Triangle (3 nodes)
111  // 4 Tetrahedron (4 nodes)
112  // 15 Point (1 node)
113  unsigned int type = lexical_cast<unsigned int>(*tok_); ++tok_;
114  unsigned int dim;
115  switch (type) {
116  case 1:
117  dim = 1;
118  break;
119  case 2:
120  dim = 2;
121  break;
122  case 4:
123  dim = 3;
124  break;
125  case 15:
126  dim = 0;
127  break;
128  default:
129  dim = 0;
130  THROW(ExcUnsupportedType() << EI_ElementId(id) << EI_ElementType(type) << EI_GMSHFile(tok_.f_name()) );
131  break;
132  }
133 
134  //get number of tags (at least 2)
135  unsigned int n_tags = lexical_cast<unsigned int>(*tok_);
136  INPUT_CHECK(n_tags >= 2, "At least two element tags have to be defined for element with id=%d, %s.\n",
137  id, tok_.position_msg().c_str());
138  ++tok_;
139 
140  //get tags 1 and 2
141  unsigned int region_id = lexical_cast<unsigned int>(*tok_); ++tok_; // region_id
142  lexical_cast<unsigned int>(*tok_); ++tok_; // GMSH region number, we do not store this
143  //get remaining tags
144  unsigned int partition_id = 0;
145  if (n_tags > 2) { partition_id = lexical_cast<unsigned int>(*tok_); ++tok_; } // save partition number from the new GMSH format
146  for (unsigned int ti = 3; ti < n_tags; ti++) ++tok_; //skip remaining tags
147 
148  for (unsigned int ni=0; ni<dim+1; ++ni) { // read node ids
149  node_ids[ni] = lexical_cast<unsigned int>(*tok_);
150  ++tok_;
151  }
152  mesh->add_element(id, dim, region_id, partition_id, node_ids);
153  }
154 
155  } catch (bad_lexical_cast &) {
156  THROW(ExcWrongFormat() << EI_Type("number") << EI_TokenizerMsg(tok_.position_msg()) << EI_MeshFile(tok_.f_name()) );
157  }
158 
159  mesh->create_boundary_elements();
160  MessageOut().fmt("... {} bulk elements, {} boundary elements. \n", mesh->n_elements(), mesh->n_elements(true));
161 }
162 
163 
164 
166  ASSERT(mesh).error("Argument mesh is NULL.\n");
167 
168  using namespace boost;
169 
170  if (! tok_.skip_to("$PhysicalNames", "$Nodes") ) return;
171  try {
172  tok_.next_line(false);
173  unsigned int n_physicals = lexical_cast<unsigned int> (*tok_);
174  ++tok_; // end of line
175 
176  for (unsigned int i = 0; i < n_physicals; ++i) {
177  tok_.next_line();
178  // format of one line:
179  // dim physical-id physical-name
180 
181  unsigned int dim = lexical_cast<unsigned int>(*tok_); ++tok_;
182  unsigned int id = lexical_cast<unsigned int>(*tok_); ++tok_;
183  string name = *tok_; ++tok_;
184  mesh->add_physical_name( dim, id, name );
185  }
186 
187  } catch (bad_lexical_cast &) {
188  THROW(ExcWrongFormat() << EI_Type("number") << EI_TokenizerMsg(tok_.position_msg()) << EI_MeshFile(tok_.f_name()) );
189  }
190 
191 }
192 
193 
194 // Is assumed to be called just after tok.skip_to("..")
195 // reads the header from the tokenizer @p tok and return it as the second parameter
197  using namespace boost;
198  try {
199  // string tags
200  tok_.next_line(false);
201  unsigned int n_str = lexical_cast<unsigned int>(*tok_); ++tok_;
202  head.field_name="";
203  head.interpolation_scheme = "";
204  if (n_str > 0) {
205  tok_.next_line(); n_str--;
206  head.field_name= *tok_; ++tok_; // unquoted by tokenizer if needed
207  }
208  if (n_str > 0) {
209  tok_.next_line(); n_str--;
210  head.interpolation_scheme = *tok_; ++tok_;
211  }
212  for(;n_str>0;n_str--) tok_.next_line(false); // skip possible remaining tags
213 
214  //real tags
215  tok_.next_line();
216  unsigned int n_real = lexical_cast<unsigned int>(*tok_); ++tok_;
217  head.time=0.0;
218  if (n_real>0) {
219  tok_.next_line(); n_real--;
220  head.time=lexical_cast<double>(*tok_); ++tok_;
221  }
222  for(;n_real>0;n_real--) tok_.next_line(false);
223 
224  // int tags
225  tok_.next_line();
226  unsigned int n_int = lexical_cast<unsigned int>(*tok_); ++tok_;
227  head.time_index=0;
228  head.n_components=1;
229  head.n_entities=0;
230  head.partition_index=0;
231  if (n_int>0) {
232  tok_.next_line(); n_int--;
233  head.time_index=lexical_cast<unsigned int>(*tok_); ++tok_;
234  }
235  if (n_int>0) {
236  tok_.next_line(); n_int--;
237  head.n_components=lexical_cast<unsigned int>(*tok_); ++tok_;
238  }
239  if (n_int>0) {
240  tok_.next_line(); n_int--;
241  head.n_entities=lexical_cast<unsigned int>(*tok_); ++tok_;
242  }
243  for(;n_int>0;n_int--) tok_.next_line(false);
244  head.position = tok_.get_position();
245  head.discretization = OutputTime::DiscreteSpace::ELEM_DATA;
246  } catch (bad_lexical_cast &) {
247  THROW(ExcWrongFormat() << EI_Type("$ElementData header") << EI_TokenizerMsg(tok_.position_msg()) << EI_MeshFile(tok_.f_name()) );
248  }
249 }
250 
251 
252 
253 void GmshMeshReader::read_element_data(ElementDataCacheBase &data_cache, MeshDataHeader actual_header, unsigned int n_components,
254  bool boundary_domain) {
255  unsigned int id, i_row;
256  unsigned int n_read = 0;
257  std::vector<int> const & el_ids = this->get_element_vector(boundary_domain);
258  vector<int>::const_iterator id_iter = el_ids.begin();
259 
260  // read @p data buffer as we have correct header with already passed time
261  // we assume that @p data buffer is big enough
262  tok_.set_position(actual_header.position);
263 
264  // read data
265  for (i_row = 0; i_row < actual_header.n_entities; ++i_row)
266  try {
267  tok_.next_line();
268  id = boost::lexical_cast<unsigned int>(*tok_); ++tok_;
269  //skip_element = false;
270  while (id_iter != el_ids.end() && *id_iter < (int)id) {
271  ++id_iter; // skip initialization of some rows in data if ID is missing
272  }
273  if (id_iter == el_ids.end()) {
274  WarningOut().fmt("In file '{}', '$ElementData' section for field '{}', time: {}.\nData ID {} not found or is not in order. Skipping rest of data.\n",
275  tok_.f_name(), actual_header.field_name, actual_header.time, id);
276  break;
277  }
278  // save data from the line if ID was found
279  if (*id_iter == (int)id) {
280  data_cache.read_ascii_data(tok_, n_components, (id_iter - el_ids.begin()) );
281  n_read++;
282  }
283  // skip the line if ID on the line < actual ID in the map el_ids
284  } catch (boost::bad_lexical_cast &) {
285  THROW(ExcWrongFormat() << EI_Type("$ElementData line") << EI_TokenizerMsg(tok_.position_msg())
286  << EI_MeshFile(tok_.f_name()) );
287  }
288  // possibly skip remaining lines after break
289  while (i_row < actual_header.n_entities) tok_.next_line(false), ++i_row;
290 
291  LogOut().fmt("time: {}; {} entities of field {} read.\n",
292  actual_header.time, n_read, actual_header.field_name);
293 }
294 
295 
296 
298 {
299  header_table_.clear();
300  MeshDataHeader header;
301  while ( !tok_.eof() ) {
302  if ( tok_.skip_to("$ElementData") ) {
303  read_data_header(header);
304  HeaderTable::iterator it = header_table_.find(header.field_name);
305 
306  if (it == header_table_.end()) { // field doesn't exists, insert new vector to map
308  vec.push_back(header);
309  header_table_[header.field_name]=vec;
310  } else if ( header.time <= it->second.back().time ) { // time is in wrong order. can't be add
311  WarningOut().fmt("Wrong time order: field '{}', time '{}', file '{}'. Skipping this '$ElementData' section.\n",
312  header.field_name, header.time, tok_.f_name() );
313  } else { // add new time step
314  it->second.push_back(header);
315  }
316  }
317  }
318 
319  tok_.set_position( Tokenizer::Position() );
320 }
321 
322 
323 
325 {
326  // check discretization, only type element_data or undefined is supported
327  if (header_query.discretization != OutputTime::DiscreteSpace::ELEM_DATA) {
328  if (header_query.discretization != OutputTime::DiscreteSpace::UNDEFINED && header_query.discretization != OutputTime::DiscreteSpace::NATIVE_DATA) {
329  WarningOut().fmt(
330  "Unsupported discretization for field '{}', time: {} and GMSH format.\nType 'ELEM_DATA' of discretization will be used.\n",
331  header_query.field_name, header_query.time);
332  }
333  header_query.discretization = OutputTime::DiscreteSpace::ELEM_DATA;
334  }
335 
336  HeaderTable::iterator table_it = header_table_.find(header_query.field_name);
337 
338  if (table_it == header_table_.end()) {
339  // no data found
340  THROW( ExcFieldNameNotFound() << EI_FieldName(header_query.field_name) << EI_MeshFile(tok_.f_name()));
341  }
342 
343  auto comp = [](double t, const MeshDataHeader &a) {
344  return t * (1.0 + 2.0*numeric_limits<double>::epsilon()) < a.time;
345  };
346 
347  std::vector<MeshDataHeader>::iterator headers_it = std::upper_bound(table_it->second.begin(),
348  table_it->second.end(),
349  header_query.time,
350  comp);
351 
352  if (headers_it == table_it->second.begin()) {
353  THROW( ExcFieldNameNotFound() << EI_FieldName(header_query.field_name)
354  << EI_MeshFile(tok_.f_name()) << EI_Time(header_query.time));
355  }
356 
357  --headers_it;
358  actual_header_ = *headers_it;
359  return actual_header_;
360 }
ArmaVec< double, N > vec
Definition: armor.hh:861
void make_header_table() override
void read_nodes(Mesh *mesh)
void read_data_header(MeshDataHeader &head)
std::string interpolation_scheme
Currently d&#39;ont used.
void init_element_vector(unsigned int size)
Initialize element_vec_, set size and reset counters of boundary and bulk elements.
Definition: mesh.cc:1127
void add_node(unsigned int node_id, arma::vec3 coords)
Add new node of given id and coordinates to mesh.
Definition: mesh.cc:1068
#define MessageOut()
Macro defining &#39;message&#39; record of log.
Definition: logger.hh:255
std::string field_name
Name of field.
virtual ~GmshMeshReader()
Definition: mesh.h:78
void read_physical_names(Mesh *mesh) override
std::string data_section_name_
Store name of field data section specify for type of mesh file.
GmshMeshReader(const FilePath &file_name)
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library...
Definition: asserts.hh:347
void add_physical_name(unsigned int dim, unsigned int id, std::string name)
Add new node of given id and coordinates to mesh.
Definition: mesh.cc:1063
#define LogOut()
Macro defining &#39;log&#39; record of log.
Definition: logger.hh:261
unsigned int partition_index
Currently d&#39;ont used.
unsigned int time_index
Currently d&#39;ont used.
unsigned int n_entities
Number of rows.
virtual void read_ascii_data(Tokenizer &tok, unsigned int n_components, unsigned int i_row)=0
OutputTime::DiscreteSpace discretization
Flag marks input discretization of data of VTK file.
void read_element_data(ElementDataCacheBase &data_cache, MeshDataHeader actual_header, unsigned int n_components, bool boundary_domain) override
std::string field_name
Name of field.
#define INPUT_CHECK(i,...)
Debugging macros.
Definition: global_defs.h:51
void create_boundary_elements()
Create boundary elements from data of temporary structure, this method MUST be call after read mesh f...
Definition: mesh.cc:1261
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
void read_elements(Mesh *mesh)
HeaderTable header_table_
Table with data of ElementData headers.
Tokenizer::Position position
Position of data in mesh file.
virtual unsigned int n_elements(bool boundary=false) const
Returns count of boundary or bulk elements.
Definition: mesh.h:344
OutputTime::DiscreteSpace discretization
Flag determinate type of discrete of Field (typically is used for native data of VTK) ...
double time
Time of field data (used only for GMSH reader)
unsigned int n_components
Number of values on one row.
double time
Time of field data (used only for GMSH and PVD reader)
MeshDataHeader actual_header_
Header of actual loaded data.
#define WarningOut()
Macro defining &#39;warning&#39; record of log.
Definition: logger.hh:258
void add_element(unsigned int elm_id, unsigned int dim, unsigned int region_id, unsigned int partition_id, std::vector< unsigned int > node_ids)
Add new element of given id to mesh.
Definition: mesh.cc:1075
MeshDataHeader & find_header(HeaderQuery &header_query) override
Tokenizer tok_
Tokenizer used for reading ASCII file format.
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
std::vector< int > const & get_element_vector(bool boundary_domain)
void init_node_vector(unsigned int size)
Initialize node_vec_, set size.
Definition: mesh.cc:1139