VtkReader.h 9.38 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/


/** \file VtkReader.h */

#ifndef AMDIS_VTKREADER_DETAIL_H
#define AMDIS_VTKREADER_DETAIL_H

// need some extension-methods/libs (pugixml, nanoflann)
#ifdef HAVE_EXTENSIONS

#include <cstring>
31
32
33
34
35
36
37
#include <boost/archive/iterators/binary_from_base64.hpp>
#include <boost/archive/iterators/transform_width.hpp>
#include <boost/iostreams/filtering_streambuf.hpp>
#include <boost/iostreams/copy.hpp>
#ifdef HAVE_COMPRESSION
#include <boost/iostreams/filter/zlib.hpp>
#endif
38
39
40
41
#include "DOFVector.h"
#include "Mesh.h"
#include "boost/filesystem.hpp"
#include "pugixml.hpp"
42
// extensions {
43
44
#include "kdtree_nanoflann.h"
#include "VectorOperations.h"
45
// }
46

47
48
#define BLOCKSIZE 32768

49
50
51
52
53
54
55
56
namespace AMDiS 
{
  namespace io 
  {
    namespace VtkReader 
    {
      namespace detail
      {
57
58
59
60
61
62
	typedef boost::archive::iterators::transform_width<
	          boost::archive::iterators::binary_from_base64<
	              std::string::const_iterator
	          >, 8, 6 
	        > base64_binary;
	
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
	inline void valueVector2type(std::vector<double> p, 
				    double& value)
	{
	  TEST_EXIT(p.size() != 0)("Not enough data for assignment!\n");
	  value = p[0];
	}

	
	inline void valueVector2type(std::vector<double> p, 
				    WorldVector<double>& value)
	{
	  TEST_EXIT(static_cast<int>(p.size()) == Global::getGeo(WORLD))("Not enough data for assignment!\n");
	  for (int i = 0; i < Global::getGeo(WORLD); i++)
	    value[i] = p[i];
	}

	
	inline void valueVector2type(std::vector<double> p, 
				    std::vector<double>& value)
	{
	  TEST_EXIT(p.size() != 0)("Not enough data for assignment!\n");
	  for (size_t i = 0; i < p.size(); i++)
	    value.push_back(p[i]);
	}
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
	
	inline std::string base64ToStr(std::string base64)
	{
	  using namespace std;
	  unsigned int paddChars = count(base64.begin(), base64.end(), '=');
	  std::replace(base64.begin(),base64.end(),'=','A'); 
	  string result(base64_binary(base64.begin()), base64_binary(base64.end()));
	  result.erase(result.end()-paddChars,result.end());
	  return result;
	}

#ifdef HAVE_COMPRESSION
	inline std::string decompress(std::string text)
	{
	  std::stringstream tmp1, tmp2;
	  tmp1.str(text);
	  boost::iostreams::filtering_streambuf<boost::iostreams::input> in;
	  in.push(boost::iostreams::zlib_decompressor());
	  in.push(tmp1);
	  boost::iostreams::copy(in, tmp2);
	  return tmp2.str();
	}
#endif

111
	inline std::string getInnerDataArray(std::string& input, bool zlib, bool base64)
112
113
114
115
116
	{
	  FUNCNAME("VtkReader::detail::getInnerDataArray()");
	  
	  using namespace std;
	  
117
118
	  char* tmp = NULL;
	  int* ptr = NULL;
119
120
121
122
123
	  int nBytes = 0;
	  string header = "", body = "", data = "";
	  
	  if(zlib) {
#ifdef HAVE_COMPRESSION
Siqi Ling's avatar
Siqi Ling committed
124
	    string s = input.substr(0, 8);
125
126
127
	    if (base64)
	      s = detail::base64ToStr(s);
	    tmp =  const_cast<char*>(input.c_str());
128
129
130
	    ptr = reinterpret_cast<int*>(tmp);
	    
	    int nBlocks = *ptr;
131
132
	    int headerSize = (base64) ? (((4 * nBlocks + 12) % 3) > 0) ? 4 * ((4 * nBlocks + 12) / 3 + 1) : 4 * ((4 * nBlocks + 12) / 3)
			    : 4 * nBlocks + 12;
133
134
	    header = input.substr(0, headerSize);
	    body = input.substr(headerSize);
135
136
137
138
	    if (base64) {
	      header = detail::base64ToStr(header);
	      body = detail::base64ToStr(body);
	    }
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
	    
	    int blockSize, finalSize, offset = 0;
	    tmp =  const_cast<char*>(header.c_str());
	    ptr = reinterpret_cast<int*>(tmp);
	    ptr++;
	    blockSize = *ptr++;
	    finalSize = *ptr++;
	    nBytes = (nBlocks - 1) * blockSize + finalSize;
	    
	    for(int i = 0; i < nBlocks; i++, ptr++) {
	      string subBody = body.substr(offset, *ptr);
	      data += decompress(subBody);
	      offset += *ptr;
	    }
#else
            ERROR_EXIT("HAVE_COMPRESSION OFF. VtkReader cannot read APPENDED_COMPRESSED vtu files.\n");
#endif
	  } else {
157
 	    header = (base64) ? detail::base64ToStr(input) : input;
158
159
160
161
162
163
164
165
166
167
168
169
	    tmp =  const_cast<char*>(header.c_str());
	    ptr = reinterpret_cast<int*>(tmp);
	    nBytes = *ptr;
	    data = header.substr(4);
	  }
	  TEST_EXIT((unsigned)nBytes == data.length())("Compressed DataArray is wrong. nBytes: %i and data length: %i.\n", nBytes, data.length());
	  return data;
	}
	
	inline void binary2pointList(std::string& input,
				     std::string type,
			             bool zlib,
170
				     bool base64,
171
172
173
174
175
				     std::vector<WorldVector<double> >& pointList)
	{
	  int dow = Global::getGeo(WORLD);
	  int dowExplicit = 3;
	  
176
	  std::string inner = getInnerDataArray(input, zlib, base64);
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
	  int nBytes = inner.length();
	  char* tmp = const_cast<char*>(inner.c_str());
	  
	  pointList.clear();

	  if(type == "Float32") {
	    int size = nBytes / 4;
	    float* ptr = reinterpret_cast<float*>(tmp); 
	    
	    for(int i = 0; i < size;) {
	      WorldVector<double> p;
	      for (int j = 0; j < dowExplicit && i < size; i++, j++, ptr++)
		if (j < dow)
		  p[j] = boost::lexical_cast<double>(*ptr);
	      pointList.push_back(p);
	    }
	  }
	  else if(type == "Float64") {
	    int size = nBytes / 8;
	    double* ptr = reinterpret_cast<double*>(tmp); 
	    
	    for(int i = 0; i < size;) {
	      WorldVector<double> p;
	      for (int j = 0; j < dowExplicit && i < size; i++, j++, ptr++)
		if (j < dow)
		  p[j] = *ptr;
	      pointList.push_back(p);
	    }
	  } else {
	    ERROR_EXIT("Currently only supports Float32 and Float64 data type.\n");
	  }
	}
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235

	
	inline void string2pointList(std::string& input, 
				    std::vector<WorldVector<double> >& pointList)
	{
	  int dow = Global::getGeo(WORLD);
	  int dowExplicit = 3;

	  std::vector<std::string> valueStringList;
	  std::stringstream ss(input);
	  std::string buf;
	  while (ss >> buf)
	    valueStringList.push_back(buf);

	  pointList.clear();

	  std::vector<std::string>::iterator it;
	  size_t j = 0;
	  for (it = valueStringList.begin(); it != valueStringList.end();j++) {
	    WorldVector<double> p;
	    for (int i = 0; i < dowExplicit && it != valueStringList.end(); i++, it++)
	      if (i < dow)
		p[i] = boost::lexical_cast<double>(*it);
	    pointList.push_back(p);
	  }
	}

236
237
238
239
	template<typename T>
	void binary2valueList(std::string& input, 
			            std::string type,
				    bool zlib,
240
				    bool base64,
241
242
243
244
245
246
247
				    std::vector<T>& valueList, 
				    int numComponent = 1, 
				    int numComponentMax = -1)
	{
	  if (numComponentMax < 0)
	    numComponentMax = numComponent;
	  
248
	  std::string inner = getInnerDataArray(input, zlib, base64);
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
	  int nBytes = inner.length();
	  char* tmp = const_cast<char*>(inner.c_str());
	  
	  valueList.clear();
	  
	  if(type == "Float32") {
	    int size = nBytes / 4;
	    float* ptr = reinterpret_cast<float*>(tmp); 
	    
	    for(int i = 0; i < size;) {
	      std::vector<double> p;
	      for (int j = 0; j < numComponentMax && i < size; i++, j++, ptr++) 
		if (j < numComponent)
   		  p.push_back(boost::lexical_cast<double>(*ptr));
		
	      T value; valueVector2type(p, value);
	      valueList.push_back(value);
	    }
	  } else if(type == "Float64") {
	    int size = nBytes / 8;
	    double* ptr = reinterpret_cast<double*>(tmp); 
	    
	    for(int i = 0; i < size;) {
	      std::vector<double> p;
	      for (int j = 0; j < numComponentMax && i < size; i++, j++, ptr++) 
		if (j < numComponent) 
		  p.push_back(*ptr);
	      T value; valueVector2type(p, value);
	      valueList.push_back(value);
	    }
	  }
	  else {
	    ERROR_EXIT("Currently only supports Float32 and Float64 data type.\n");
	  }
	}
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
	
	template<typename T>
	void string2valueList(std::string& input, 
				    std::vector<T>& valueList, 
				    int numComponent = 1, 
				    int numComponentMax = -1)
	{
	  if (numComponentMax < 0)
	    numComponentMax = numComponent;

	  std::vector<std::string> valueStringList;
	  std::stringstream ss(input);
	  std::string buf;
	  while (ss >> buf)
	    valueStringList.push_back(buf);

	  valueList.clear();

	  std::vector<std::string>::iterator it;
	  for (it = valueStringList.begin(); it != valueStringList.end();) {
	    std::vector<double> p;
	    for (int i = 0; i < numComponentMax && it != valueStringList.end(); i++, it++)
	      if (i < numComponent)
		p.push_back(boost::lexical_cast<double>(*it));
308
	      
309
310
311
312
313
314
315
	    T value; valueVector2type(p, value);
	    valueList.push_back(value);
	  }
	}


	// find point in mesh using KD-tree structure
316
	inline size_t getNearestIndex(extensions::KD_Tree& tree, 
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
				      WorldVector<double>& x)
	{
	  // do a knn search
	  const size_t num_nnp = 1;
	  std::vector<size_t> ret_indexes(num_nnp);
	  std::vector<double> out_dists_sqr(num_nnp);

	  nanoflann::KNNResultSet<double> resultSet(num_nnp);
	  resultSet.init(&ret_indexes[0], &out_dists_sqr[0] );

	  tree.index->findNeighbors(resultSet, x.begin(), nanoflann::SearchParams(10));
	  return ret_indexes[0];
	}

      } // end namespace detail    
    } // end namespace VtkReader
  } // end namespace io
} // end namespace AMDiS

#endif // HAVE_EXTENSIONS

#endif // AMDIS_VTKREADER_DETAIL_H