Rheolef  7.2
an efficient C++ finite element environment
geo_seq_get_vtk.cc
Go to the documentation of this file.
1 //
22 // geo: vtk input
23 //
24 // author: Pierre.Saramito@imag.fr
25 //
26 // 5 march 2012
27 //
28 #include "rheolef/geo.h"
29 #include "rheolef/rheostream.h"
30 #include "rheolef/iorheo.h"
31 
32 #include "vtk_cell_type.h"
33 
34 namespace rheolef {
35 using namespace std;
36 
37 // skip optional line, e.g.
38 // "OFFSETS vtktypeint64"
39 // "CONNECTIVITY vtktypeint64"
40 // that do not start by a digit
41 // if a digit is founded, returns an empty string
42 // otherwise, returns the first keyword, e.g. "OFFSETS" or "CONNECTIVITY"
43 // and skip until the end of line
44 //
45 static
46 string
47 skip_optional_line (istream& is)
48 {
49  is >> ws;
50  char c = is.peek();
51  if (isdigit(c)){
52  return string("");
53  }
54  string key;
55  is >> key;
56  c = is.peek();
57  while (c != '\n' && is.good()) {
58  is.get(c);
59  }
60  return key;
61 }
62 template <class T>
63 idiststream&
65 {
66  using namespace std;
68  typedef typename geo_basic<T,sequential>::node_type node_type;
69 
70  check_macro (ips.good(), "bad input stream for vtk");
71  istream& is = ips.is();
72  geo_header hdr;
73  hdr.order = 1;
74  // ------------------------------------------------------------------------
75  // 1) load the coordinates
76  // POINTS <np> float
77  // {xi yi zi} i=0..np-1
78  // ------------------------------------------------------------------------
79  check_macro (scatch(is,"POINTS"), "unexpected vtk input file");
81  is >> nnod;
82  if (nnod == 0) {
83  warning_macro("empty vtk mesh file");
84  return ips;
85  }
87  scatch (is,"float");
88  node.get_values (ips);
89  check_macro (ips.good(), "bad input stream for vtk");
90  // ------------------------------------------------------------------------
91  // 1.2) compute dimension
92  // ------------------------------------------------------------------------
93  point xmin, xmax;
94  for (size_type k = 0; k < 3; ++k) {
95  xmin[k] = std::numeric_limits<T>::max();
96  xmax[k] = std::numeric_limits<T>::min();
97  }
98  for (size_t i = 0, n = node.size(); i < n; ++i) {
99  for (size_type k = 0; k < 3; ++k) {
100  xmin[k] = std::min (node[i][k], xmin[k]);
101  xmax[k] = std::max (node[i][k], xmax[k]);
102  }
103  }
104  hdr.dimension = 3;
105  for (size_type k = 0; k < 3; ++k) {
106  if (fabs(xmax[k] - xmin[k]) < std::numeric_limits<T>::epsilon())
107  hdr.dimension--;
108  }
109  // ------------------------------------------------------------------------
110  // 2) load the domain connectivity
111  // LINES <ne> <ncs>
112  // 2 {s1 s2} j=0..ne-1
113  // or
114  // POLYGONS <ne> <ncs>
115  // p {s1 .. sp} j=0..ne-1
116  // ------------------------------------------------------------------------
117  // 2.1) get connectivity header
118  string mark;
119  is >> ws >> mark;
120  if (mark == "VERTICES") { // skip unused paragraph in polydata v3.0
121  size_type n, n2, d, idx;
122  is >> n >> n2;
123  for (size_type i = 0; i < n; i++) {
124  is >> d >> idx;
125  }
126  }
127  do {
128  if (mark == "POLYGONS" || mark == "LINES" || mark == "CELLS") break;
129  } while (is >> ws >> mark);
130  bool have_vtk_cell_type = false;
131  string mesh_fmt = mark;
132  if (mesh_fmt == "LINES") hdr.map_dimension = 1;
133  else if (mesh_fmt == "POLYGONS") hdr.map_dimension = 2;
134  else if (mesh_fmt == "CELLS") { hdr.map_dimension = 0; have_vtk_cell_type = true; } // yet unknown
135  else error_macro ("geo: unexpected `" << mesh_fmt << "' in vtk input."
136  << " Expect LINES, POLYGONS or CELLS.");
137  // note: when mesh_fmt == "CELLS", we will also get CELL_TYPES and deduce the map_dimension
138  size_type ne, ncs;
139  is >> ne >> ncs;
140  // skip optional line: "OFFSETS vtktypeint64"
141  mark = skip_optional_line (is);
142  // 2.2) get all connectivity in a flat trash array
143  std::vector<size_type> cell_trash (ncs);
144  std::vector<size_type> vtk_cell_type (ne);
145  typedef geo_element_auto<> element_t;
146  typedef disarray<element_t,sequential> disarray_t;
147  disarray_t elt (ne);
148  if (mark == "") {
149  // old vtk polydata format
150  for (size_type ics = 0; ics < ncs; ics++) {
151  is >> cell_trash [ics];
152  }
153  // skip optional line: "CONNECTIVITY vtktypeint64"
154  if (have_vtk_cell_type) {
155  is >> ws >> mark;
156  check_macro (mark == "CELL_TYPES", "geo: unexpected `" << mark << "' in vtk input."
157  << " Expect CELL_TYPES.");
158  size_type ne2;
159  is >> ne2;
160  check_macro (ne == ne2, "geo: unexpected CELL_TYPES size=`" << ne2 << "' , expect size="<<ne);
161  for (size_type ie = 0; ie < ne; ie++) {
162  is >> vtk_cell_type [ie];
163  }
164  }
165  for (size_type ie = 0, ics = 0; ie < ne; ie++) {
166  size_type nloc = cell_trash [ics];
167  ics++;
168  size_type variant = (!have_vtk_cell_type) ?
170  vtk_cell_type2variant (vtk_cell_type [ie]);
171  elt[ie].reset (variant, hdr.order);
172  for (size_type iloc = 0 ; iloc < nloc; iloc++, ics++) {
173  elt[ie][iloc] = cell_trash [ics];
174  }
175  hdr.dis_size_by_dimension [elt[ie].dimension()]++;
176  hdr.dis_size_by_variant [elt[ie].variant()]++;
177  hdr.map_dimension = std::max (hdr.map_dimension, elt[ie].dimension());
178  }
179  hdr.dis_size_by_dimension [0] = nnod;
180  hdr.dis_size_by_variant [0] = nnod;
181  } else {
182  // new vtk polydata format
183  for (size_t count_key = 0; mark == "OFFSETS" || mark == "CONNECTIVITY"; ++count_key) {
184  if (mark == "OFFSETS") {
185  size_type prec_offset = 0;
186  is >> prec_offset;
187  if (ne > 0) ne--; // first offset is zero
188  for (size_type ie = 0; ie < ne; ie++) {
189  size_type offset;
190  is >> offset;
191  size_type cell_nv = size_type(offset - prec_offset);
192  vtk_cell_type [ie] = nv2vtk_cell_type (hdr.map_dimension, cell_nv);
193  prec_offset = offset;
194  }
195  } else { // CONNECTIVITY
196  for (size_type ics = 0; ics < ncs; ics++) {
197  is >> cell_trash [ics];
198  }
199  }
200  if (count_key == 1) break;
201  // else : read next keyword
202  mark = skip_optional_line (is);
203  }
204  for (size_type ie = 0, ics = 0; ie < ne; ie++) {
205  size_type variant = vtk_cell_type2variant (vtk_cell_type [ie]);
206  elt[ie].reset (variant, hdr.order);
207  size_type nloc = elt [ie].size();
208  for (size_type iloc = 0 ; iloc < nloc; iloc++, ics++) {
209  check_macro (ics < ncs, "unexpected ics="<<ics<<" out of range [0:"<<ncs<<"[");
210  elt[ie][iloc] = cell_trash [ics];
211  }
212  hdr.dis_size_by_dimension [elt[ie].dimension()]++;
213  hdr.dis_size_by_variant [elt[ie].variant()]++;
214  hdr.map_dimension = std::max (hdr.map_dimension, elt[ie].dimension());
215  }
216  hdr.dis_size_by_dimension [0] = nnod;
217  hdr.dis_size_by_variant [0] = nnod;
218  }
219  // ------------------------------------------------------------------------
220  // 2) split elements by variants
221  // ------------------------------------------------------------------------
222  std::array<disarray_t, reference_element::max_variant> tmp_geo_element;
223  if (hdr.map_dimension > 0) {
226  geo_element_auto<> init_val (variant, hdr.order);
227  tmp_geo_element [variant] = disarray_t (hdr.dis_size_by_variant [variant], init_val);
228  }
229  std::array<size_type, 4> counter;
230  counter.fill(0);
231  for (size_type ie = 0; ie < ne; ie++) {
232  size_type variant = elt[ie].variant();
233  size_type ige = counter[variant];
234  tmp_geo_element [variant] [ige] = elt[ie];
235  counter[variant]++;
236  }
237  }
238  // ------------------------------------------------------------------------
239  // 3) build & upgrade
240  // ------------------------------------------------------------------------
241  omega.build_from_data (hdr, node, tmp_geo_element, true);
242  return ips;
243 }
244 // ----------------------------------------------------------------------------
245 // instanciation in library
246 // ----------------------------------------------------------------------------
248 
249 }// namespace rheolef
field::size_type size_type
Definition: branch.cc:430
see the point page for the full documentation
see the disarray page for the full documentation
Definition: disarray.h:497
void build_from_data(const geo_header &hdr, const disarray< node_type, sequential > &node, std::array< disarray< geo_element_auto<>, sequential >, reference_element::max_variant > &tmp_geo_element, bool do_upgrade)
generic mesh with rerefence counting
Definition: geo.h:1089
idiststream: see the diststream page for the full documentation
Definition: diststream.h:336
std::istream & is()
Definition: diststream.h:400
bool good() const
Definition: diststream.cc:124
static variant_type last_variant_by_dimension(size_type dim)
static variant_type first_variant_by_dimension(size_type dim)
variant_type variant() const
size_t size_type
Definition: basis_get.cc:76
#define error_macro(message)
Definition: dis_macros.h:49
#define warning_macro(message)
Definition: dis_macros.h:53
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
This file is part of Rheolef.
bool scatch(std::istream &in, const std::string &ch, bool full_match=true)
scatch: see the rheostream page for the full documentation
Definition: scatch.icc:44
size_t vtk_cell_type2variant(size_t vtk_cell_type)
size_t nv2vtk_cell_type(size_t map_dim, size_t nv)
idiststream & geo_get_vtk(idiststream &ips, geo_basic< T, sequential > &omega)
template idiststream & geo_get_vtk< Float >(idiststream &, geo_basic< Float, sequential > &)
size_type map_dimension
Definition: geo_header.h:40
size_type dis_size_by_dimension[4]
Definition: geo_header.h:44
size_type dimension
Definition: geo_header.h:39
size_type dis_size_by_variant[reference_element::max_variant]
Definition: geo_header.h:43
Float epsilon