Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_seq_get_vtk.cc
Go to the documentation of this file.
1
21//
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
34namespace rheolef {
35using 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//
45static
46string
47skip_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}
62template <class T>
63idiststream&
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");
80 size_type nnod;
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
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
static variant_type last_variant_by_dimension(size_type dim)
static variant_type first_variant_by_dimension(size_type dim)
variant_type variant() const
#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)")
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
idiststream & geo_get_vtk(idiststream &ips, geo_basic< T, sequential > &omega)
template idiststream & geo_get_vtk< Float >(idiststream &, geo_basic< Float, sequential > &)
size_t vtk_cell_type2variant(size_t vtk_cell_type)
size_t nv2vtk_cell_type(size_t map_dim, size_t nv)
STL namespace.
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