Flow123d  master-f44eb46
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 #include "mesh/bc_mesh.hh"
33 
34 
35 
36 using namespace std;
37 
38 
40 : BaseMeshReader(file_name)
41 {
42  tok_.set_comment_pattern( "#");
43  data_section_name_ = "$ElementData";
44  has_compatible_mesh_ = false;
46 }
47 
48 
49 
50 GmshMeshReader::~GmshMeshReader() // Tokenizer close the file automatically
51 {}
52 
53 
54 
56  using namespace boost;
57  unsigned int n_nodes;
58  MessageOut() << "- Reading nodes...";
59  tok_.set_position( Tokenizer::Position() );
60 
61  if (! tok_.skip_to("$Nodes")) THROW(ExcMissingSection() << EI_Section("$Nodes") << EI_GMSHFile(tok_.f_name()) );
62  try {
63  tok_.next_line(false);
64  n_nodes = lexical_cast<unsigned int> (*tok_);
65  mesh->init_node_vector( n_nodes );
66  if (n_nodes == 0) THROW( ExcZeroNodes() << EI_Position(tok_.position_msg()) );
67  ++tok_; // end of line
68 
69  for (unsigned int i = 0; i < n_nodes; ++i) {
70  tok_.next_line();
71 
72  unsigned int id = lexical_cast<unsigned int> (*tok_); ++tok_; // node id
73  arma::vec3 coords; // node coordinates
74  coords(0) = lexical_cast<double> (*tok_); ++tok_;
75  coords(1) = lexical_cast<double> (*tok_); ++tok_;
76  coords(2) = lexical_cast<double> (*tok_); ++tok_;
77  ++tok_; // skip mesh size parameter
78 
79  mesh->add_node(id, coords);
80  }
81 
82  } catch (bad_lexical_cast &) {
83  THROW(ExcWrongFormat() << EI_Type("number") << EI_TokenizerMsg(tok_.position_msg()) << EI_MeshFile(tok_.f_name()) );
84  }
85  MessageOut().fmt("... {} nodes read. \n", n_nodes);
86 }
87 
88 
90  using namespace boost;
91  MessageOut() << "- Reading elements...";
92 
93  if (! tok_.skip_to("$Elements")) THROW(ExcMissingSection() << EI_Section("$Elements") << EI_GMSHFile(tok_.f_name()) );
94  try {
95  tok_.next_line(false);
96  unsigned int n_elements = lexical_cast<unsigned int> (*tok_);
97  if (n_elements == 0) THROW( ExcZeroElements() << EI_Position(tok_.position_msg()) );
98  ++tok_; // end of line
99 
100  std::vector<unsigned int> node_ids; //node_ids of elements
101  node_ids.resize(4); // maximal count of nodes
102 
103  mesh->init_element_vector(n_elements);
104 
105  for (unsigned int i = 0; i < n_elements; ++i) {
106  tok_.next_line();
107  unsigned int id = lexical_cast<unsigned int>(*tok_); ++tok_;
108 
109  //get element type: supported:
110  // 1 Line (2 nodes)
111  // 2 Triangle (3 nodes)
112  // 4 Tetrahedron (4 nodes)
113  // 15 Point (1 node)
114  unsigned int type = lexical_cast<unsigned int>(*tok_); ++tok_;
115  unsigned int dim;
116  switch (type) {
117  case 1:
118  dim = 1;
119  break;
120  case 2:
121  dim = 2;
122  break;
123  case 4:
124  dim = 3;
125  break;
126  case 15:
127  dim = 0;
128  break;
129  default:
130  dim = 0;
131  THROW(ExcUnsupportedType() << EI_ElementId(id) << EI_ElementType(type) << EI_GMSHFile(tok_.f_name()) );
132  break;
133  }
134 
135  //get number of tags (at least 2)
136  unsigned int n_tags = lexical_cast<unsigned int>(*tok_);
137  if (n_tags < 2) THROW( ExcTooManyElementTags() << EI_ElementId(id) << EI_Position(tok_.position_msg()) );
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  MessageOut().fmt("... {} bulk elements, {} boundary elements. \n", mesh->n_elements(), mesh->bc_mesh()->n_elements());
160 }
161 
162 
163 
165  ASSERT_PTR(mesh).error("Argument mesh is NULL.\n");
166 
167  using namespace boost;
168 
169  if (! tok_.skip_to("$PhysicalNames", "$Nodes") ) return;
170  try {
171  tok_.next_line(false);
172  unsigned int n_physicals = lexical_cast<unsigned int> (*tok_);
173  ++tok_; // end of line
174 
175  for (unsigned int i = 0; i < n_physicals; ++i) {
176  tok_.next_line();
177  // format of one line:
178  // dim physical-id physical-name
179 
180  unsigned int dim = lexical_cast<unsigned int>(*tok_); ++tok_;
181  unsigned int id = lexical_cast<unsigned int>(*tok_); ++tok_;
182  string name = *tok_; ++tok_;
183  mesh->add_physical_name( dim, id, name );
184  }
185 
186  } catch (bad_lexical_cast &) {
187  THROW(ExcWrongFormat() << EI_Type("number") << EI_TokenizerMsg(tok_.position_msg()) << EI_MeshFile(tok_.f_name()) );
188  }
189 
190 }
191 
192 
193 // Is assumed to be called just after tok.skip_to("..")
194 // reads the header from the tokenizer @p tok and return it as the second parameter
196  using namespace boost;
197  try {
198  // string tags
199  tok_.next_line(false);
200  unsigned int n_str = lexical_cast<unsigned int>(*tok_); ++tok_;
201  head.field_name="";
202  head.interpolation_scheme = "";
203  if (n_str > 0) {
204  tok_.next_line(); n_str--;
205  head.field_name= *tok_; ++tok_; // unquoted by tokenizer if needed
206  }
207  if (n_str > 0) {
208  tok_.next_line(); n_str--;
209  head.interpolation_scheme = *tok_; ++tok_;
210  }
211  for(;n_str>0;n_str--) tok_.next_line(false); // skip possible remaining tags
212 
213  //real tags
214  tok_.next_line();
215  unsigned int n_real = lexical_cast<unsigned int>(*tok_); ++tok_;
216  head.time=0.0;
217  if (n_real>0) {
218  tok_.next_line(); n_real--;
219  head.time=lexical_cast<double>(*tok_); ++tok_;
220  }
221  for(;n_real>0;n_real--) tok_.next_line(false);
222 
223  // int tags
224  tok_.next_line();
225  unsigned int n_int = lexical_cast<unsigned int>(*tok_); ++tok_;
226  head.time_index=0;
227  head.n_components=1;
228  head.n_entities=0;
229  head.partition_index=0;
230  if (n_int>0) {
231  tok_.next_line(); n_int--;
232  head.time_index=lexical_cast<unsigned int>(*tok_); ++tok_;
233  }
234  if (n_int>0) {
235  tok_.next_line(); n_int--;
236  head.n_components=lexical_cast<unsigned int>(*tok_); ++tok_;
237  }
238  if (n_int>0) {
239  tok_.next_line(); n_int--;
240  head.n_entities=lexical_cast<unsigned int>(*tok_); ++tok_;
241  }
242  for(;n_int>0;n_int--) tok_.next_line(false);
243  head.position = tok_.get_position();
244  head.discretization = OutputTime::DiscreteSpace::ELEM_DATA;
245  } catch (bad_lexical_cast &) {
246  THROW(ExcWrongFormat() << EI_Type("$ElementData header") << EI_TokenizerMsg(tok_.position_msg()) << EI_MeshFile(tok_.f_name()) );
247  }
248 }
249 
250 
251 
253  static int imax = std::numeric_limits<int>::max();
254  unsigned int id, i_row;
255  unsigned int n_bulk_read = 0, n_bdr_read = 0;
256  std::vector<int> bulk_el_ids = this->get_element_ids(false); // bulk
257  bulk_el_ids.push_back( imax ); // put 'save' item at the end of vector
258  vector<int>::const_iterator bulk_id_iter = bulk_el_ids.begin();
259  std::vector<int> bdr_el_ids = this->get_element_ids(true); // boundary
260  bdr_el_ids.push_back( imax ); // put 'save' item at the end of vector
261  vector<int>::const_iterator bdr_id_iter = bdr_el_ids.begin();
262 
263  // read @p data buffer as we have correct header with already passed time
264  // we assume that @p data buffer is big enough
265  tok_.set_position(header.position);
266 
267  // read data
268  for (i_row = 0; i_row < header.n_entities; ++i_row)
269  try {
270  tok_.next_line();
271  id = boost::lexical_cast<unsigned int>(*tok_); ++tok_;
272 
273  while ( std::min(*bulk_id_iter, *bdr_id_iter) < (int)id) { // skip initialization of some rows in data if ID is missing
274  if (*bulk_id_iter < *bdr_id_iter) ++bulk_id_iter;
275  else ++bdr_id_iter;
276  }
277 
278  if (*bulk_id_iter == (int)id) {
279  // bulk
280  data_cache.read_ascii_data(tok_, header.n_components, (bulk_id_iter - bulk_el_ids.begin()) );
281  ++n_bulk_read; ++bulk_id_iter;
282  } else if (*bdr_id_iter == (int)id) {
283  // boundary
284  unsigned int bdr_shift = data_cache.get_boundary_begin();
285  data_cache.read_ascii_data(tok_, header.n_components, (bdr_id_iter - bdr_el_ids.begin() + bdr_shift) );
286  ++n_bdr_read; ++bdr_id_iter;
287  } else {
288  if ( (*bulk_id_iter != imax) | (*bdr_id_iter != imax) )
289  WarningOut().fmt("In file '{}', '$ElementData' section for field '{}', time: {}.\nData ID {} is not in order. Skipping rest of data.\n",
290  tok_.f_name(), header.field_name, header.time, id);
291  break;
292  }
293 
294  } catch (boost::bad_lexical_cast &) {
295  THROW(ExcWrongFormat() << EI_Type("$ElementData line") << EI_TokenizerMsg(tok_.position_msg())
296  << EI_MeshFile(tok_.f_name()) );
297  }
298  // possibly skip remaining lines after break
299  while (i_row < header.n_entities) tok_.next_line(false), ++i_row;
300 
301  LogOut().fmt("time: {}; {} bulk and {} boundary entities of field {} read.\n",
302  header.time, n_bulk_read, n_bdr_read, header.field_name);
303 }
304 
305 
306 
308 {
309  header_table_.clear();
310  MeshDataHeader header;
311  while ( !tok_.eof() ) {
312  if ( tok_.skip_to("$ElementData") ) {
313  read_data_header(header);
314  HeaderTable::iterator it = header_table_.find(header.field_name);
315 
316  if (it == header_table_.end()) { // field doesn't exists, insert new vector to map
318  vec.push_back(header);
319  header_table_[header.field_name]=vec;
320  } else if ( header.time <= it->second.back().time ) { // time is in wrong order. can't be add
321  WarningOut().fmt("Wrong time order: field '{}', time '{}', file '{}'. Skipping this '$ElementData' section.\n",
322  header.field_name, header.time, tok_.f_name() );
323  } else { // add new time step
324  it->second.push_back(header);
325  }
326  }
327  }
328 
329  tok_.set_position( Tokenizer::Position() );
330 }
331 
332 
333 
335 {
336  // check discretization, only type element_data or undefined is supported
337  if (header_query.discretization != OutputTime::DiscreteSpace::ELEM_DATA) {
338  if (header_query.discretization != OutputTime::DiscreteSpace::UNDEFINED && header_query.discretization != OutputTime::DiscreteSpace::NATIVE_DATA) {
339  WarningOut().fmt(
340  "Unsupported discretization for field '{}', time: {} and GMSH format.\nType 'ELEM_DATA' of discretization will be used.\n",
341  header_query.field_name, header_query.time);
342  }
343  header_query.discretization = OutputTime::DiscreteSpace::ELEM_DATA;
344  }
345 
346  HeaderTable::iterator table_it = header_table_.find(header_query.field_name);
347 
348  if (table_it == header_table_.end()) {
349  // no data found
350  THROW( ExcFieldNameNotFound() << EI_FieldName(header_query.field_name) << EI_MeshFile(tok_.f_name()));
351  }
352 
353  auto comp = [](double t, const MeshDataHeader &a) {
354  return t * (1.0 + 2.0*numeric_limits<double>::epsilon()) < a.time;
355  };
356 
357  std::vector<MeshDataHeader>::iterator headers_it = std::upper_bound(table_it->second.begin(),
358  table_it->second.end(),
359  header_query.time,
360  comp);
361 
362  if (headers_it == table_it->second.begin()) {
363  THROW( ExcFieldNameNotFound() << EI_FieldName(header_query.field_name)
364  << EI_MeshFile(tok_.f_name()) << EI_Time(header_query.time));
365  }
366 
367  --headers_it;
368  return *headers_it;
369 }
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
Definition: asserts.hh:341
Tokenizer tok_
Tokenizer used for reading ASCII file format.
std::string data_section_name_
Store name of field data section specify for type of mesh file.
std::vector< int > const & get_element_ids(bool boundary_domain)
unsigned int get_boundary_begin() const
Returns start position of boundary data in cache.
virtual void read_ascii_data(Tokenizer &tok, unsigned int n_components, unsigned int i_row)=0
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
void read_element_data(ElementDataCacheBase &data_cache, MeshDataHeader header) override
GmshMeshReader(const FilePath &file_name)
void read_data_header(MeshDataHeader &head)
MeshDataHeader & find_header(HeaderQuery &header_query) override
void read_nodes(Mesh *mesh)
void read_elements(Mesh *mesh)
virtual ~GmshMeshReader()
HeaderTable header_table_
Table with data of ElementData headers.
void make_header_table() override
void read_physical_names(Mesh *mesh) override
void init_element_vector(unsigned int size)
Initialize element_vec_, set size and reset counters of boundary and bulk elements.
Definition: mesh.cc:1142
void init_node_vector(unsigned int size)
Initialize node_vec_, set size.
Definition: mesh.cc:1152
unsigned int n_elements() const
Definition: mesh.h:111
Definition: mesh.h:362
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:1081
BCMesh * bc_mesh() const override
Implement MeshBase::bc_mesh(), getter of boundary mesh.
Definition: mesh.h:567
void add_node(unsigned int node_id, arma::vec3 coords)
Add new node of given id and coordinates to mesh.
Definition: mesh.cc:1086
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:1094
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:278
#define MessageOut()
Macro defining 'message' record of log.
Definition: logger.hh:275
#define LogOut()
Macro defining 'log' record of log.
Definition: logger.hh:281
ArmaVec< double, N > vec
Definition: armor.hh:885
std::string field_name
Name of field.
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 and PVD reader)
unsigned int partition_index
Currently d'ont used.
unsigned int n_components
Number of values on one row.
std::string interpolation_scheme
Currently d'ont used.
unsigned int n_entities
Number of rows.
OutputTime::DiscreteSpace discretization
Flag marks input discretization of data of VTK file.
Tokenizer::Position position
Position of data in mesh file.
double time
Time of field data (used only for GMSH reader)
unsigned int time_index
Currently d'ont used.
std::string field_name
Name of field.