Rheolef  7.1
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 template <class T>
38 idiststream&
39 geo_get_vtk (idiststream& ips, geo_basic<T,sequential>& omega)
40 {
41  using namespace std;
43  typedef typename geo_basic<T,sequential>::node_type node_type;
44 
45  check_macro (ips.good(), "bad input stream for vtk");
46  istream& is = ips.is();
47  geo_header hdr;
48  hdr.order = 1;
49  // ------------------------------------------------------------------------
50  // 1) load the coordinates
51  // POINTS <np> float
52  // {xi yi zi} i=0..np-1
53  // ------------------------------------------------------------------------
54  check_macro (scatch(is,"POINTS"), "unexpected vtk input file");
56  is >> nnod;
57  if (nnod == 0) {
58  warning_macro("empty vtk mesh file");
59  return ips;
60  }
62  scatch (is,"float");
63  node.get_values (ips);
64  check_macro (ips.good(), "bad input stream for vtk");
65  // ------------------------------------------------------------------------
66  // 1.2) compute dimension
67  // ------------------------------------------------------------------------
68  point xmin, xmax;
69  for (size_type k = 0; k < 3; ++k) {
70  xmin[k] = std::numeric_limits<T>::max();
71  xmax[k] = std::numeric_limits<T>::min();
72  }
73  for (size_t i = 0, n = node.size(); i < n; ++i) {
74  for (size_type k = 0; k < 3; ++k) {
75  xmin[k] = std::min (node[i][k], xmin[k]);
76  xmax[k] = std::max (node[i][k], xmax[k]);
77  }
78  }
79  hdr.dimension = 3;
80  for (size_type k = 0; k < 3; ++k) {
81  if (fabs(xmax[k] - xmin[k]) < std::numeric_limits<T>::epsilon())
82  hdr.dimension--;
83  }
84  // ------------------------------------------------------------------------
85  // 2) load the domain connectivity
86  // LINES <ne> <ncs>
87  // 2 {s1 s2} j=0..ne-1
88  // or
89  // POLYGONS <ne> <ncs>
90  // p {s1 .. sp} j=0..ne-1
91  // ------------------------------------------------------------------------
92  // 2.1) get connectivity header
93  string mark;
94  is >> ws >> mark;
95  if (mark == "VERTICES") { // skip unused paragraph in polydata v3.0
96  size_type n, n2, d, idx;
97  is >> n >> n2;
98  for (size_type i = 0; i < n; i++) {
99  is >> d >> idx;
100  }
101  }
102  do {
103  if (mark == "POLYGONS" || mark == "LINES" || mark == "CELLS") break;
104  } while (is >> ws >> mark);
105  bool have_vtk_cell_type = false;
106  string mesh_fmt = mark;
107  if (mesh_fmt == "LINES") hdr.map_dimension = 1;
108  else if (mesh_fmt == "POLYGONS") hdr.map_dimension = 2;
109  else if (mesh_fmt == "CELLS") { hdr.map_dimension = 0; have_vtk_cell_type = true; } // yet unknown
110  else error_macro ("geo: unexpected `" << mesh_fmt << "' in vtk input."
111  << " Expect LINES, POLYGONS or CELLS.");
112  // note: when mesh_fmt == "CELLS", we will also get CELL_TYPES and deduce the map_dimension
113  size_type ne, ncs;
114  is >> ne >> ncs;
115  //
116  // 2.2) get all connectivity in a flat trash array
117  std::vector<size_type> cell_trash (ncs);
118  for (size_type ics = 0; ics < ncs; ics++) {
119  is >> cell_trash [ics];
120  }
121  std::vector<size_type> vtk_cell_type (ne);
122  if (have_vtk_cell_type) {
123  is >> ws >> mark;
124  check_macro (mark == "CELL_TYPES", "geo: unexpected `" << mark << "' in vtk input."
125  << " Expect CELL_TYPES.");
126  size_type ne2;
127  is >> ne2;
128  check_macro (ne == ne2, "geo: unexpected CELL_TYPES size=`" << ne2 << "' , expect size="<<ne);
129  for (size_type ie = 0; ie < ne; ie++) {
130  is >> vtk_cell_type [ie];
131  }
132  }
133  typedef geo_element_auto<> element_t;
134  typedef disarray<element_t,sequential> disarray_t;
135  disarray_t elt (ne);
136  for (size_type ie = 0, ics = 0; ie < ne; ie++) {
137  size_type nloc = cell_trash [ics];
138  ics++;
139  size_type variant = (!have_vtk_cell_type) ?
141  vtk_cell_type2variant (vtk_cell_type [ie]);
142  elt[ie].reset (variant, hdr.order);
143  for (size_type iloc = 0 ; iloc < nloc; iloc++, ics++) {
144  elt[ie][iloc] = cell_trash [ics];
145  }
146  hdr.dis_size_by_dimension [elt[ie].dimension()]++;
147  hdr.dis_size_by_variant [elt[ie].variant()]++;
148  hdr.map_dimension = std::max (hdr.map_dimension, elt[ie].dimension());
149  }
150  hdr.dis_size_by_dimension [0] = nnod;
151  hdr.dis_size_by_variant [0] = nnod;
152  // ------------------------------------------------------------------------
153  // 2) split elements by variants
154  // ------------------------------------------------------------------------
155  std::array<disarray_t, reference_element::max_variant> tmp_geo_element;
156  if (hdr.map_dimension > 0) {
159  geo_element_auto<> init_val (variant, hdr.order);
160  tmp_geo_element [variant] = disarray_t (hdr.dis_size_by_variant [variant], init_val);
161  }
162  std::array<size_type, 4> counter;
163  counter.fill(0);
164  for (size_type ie = 0; ie < ne; ie++) {
165  size_type variant = elt[ie].variant();
166  size_type ige = counter[variant];
167  tmp_geo_element [variant] [ige] = elt[ie];
168  counter[variant]++;
169  }
170  }
171  // ------------------------------------------------------------------------
172  // 3) build & upgrade
173  // ------------------------------------------------------------------------
174  omega.build_from_data (hdr, node, tmp_geo_element, true);
175  return ips;
176 }
177 // ----------------------------------------------------------------------------
178 // instanciation in library
179 // ----------------------------------------------------------------------------
180 template idiststream& geo_get_vtk<Float> (idiststream&, geo_basic<Float,sequential>&);
181 
182 }// namespace rheolef
field::size_type size_type
Definition: branch.cc:425
see the point page for the full documentation
see the disarray page for the full documentation
Definition: disarray.h:459
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
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:52
size_t vtk_cell_type2variant(size_t vtk_cell_type)
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