OpenLB 1.7
Loading...
Searching...
No Matches
vtiReader.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2015 Mathias J. Krause, Benjamin Förster
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22*/
23
24
25#ifndef VTI_READER_HH
26#define VTI_READER_HH
27
28#include <stdlib.h>
29#include <iostream>
30#include <fstream>
31#include <sstream>
32#include <iomanip>
33#include <math.h>
34#include <cassert>
35
36#include "geometry/cuboid2D.h"
37#include "vtiReader.h"
39
40namespace olb {
41
42//template< typename T, typename BaseType> class SuperData3D;
43
44/* ------------------ BaseVTIreader -----------------*/
45
46template<typename T>
47BaseVTIreader<T>::BaseVTIreader( const std::string& fName, int dim,
48 std::string dName, const std::string class_name )
49 : clout(class_name), _dim(dim), _size(0), _origin(_dim, 0), _extent(_dim, 0),
50 _delta(0), _xmlReader(fName), _nCuboids(0)
51{
52 // Read WholeExtent (2 * _dim ints) from XML File and calculate _extent
53 std::vector<int> wholeExtend = readExtent(&_xmlReader["ImageData"], "WholeExtent");
54 _extent = getNbNodes(wholeExtend);
55
56 // Read _delta
57 std::stringstream stream_val_1(_xmlReader["ImageData"].getAttribute("Spacing"));
58 stream_val_1 >> _delta;
59
60 // Read _origin
61 std::stringstream stream_val_2(_xmlReader["ImageData"].getAttribute("Origin"));
62 for (auto& origin_i : _origin) {
63 stream_val_2 >> origin_i;
64 }
65
66 // Read _size and count cuboids (_nCuboids)
67 for ( auto& piece : _xmlReader["ImageData"] ) {
68 if (piece->getName() == "Piece") {
69 // Read _size from first dName PointData Tag
70 if ( _nCuboids == 0 ) {
71 for (auto &dataArray : (*piece)["PointData"]) {
72 if (dataArray->getAttribute("Name") == dName && dataArray->getName() == "DataArray") {
73 _size = this->getSize(*dataArray);
74 }
75 }
76 }
77
78 _nCuboids++;
79 }
80 }
81}
82
83template<typename T>
85{
86 clout << "Information on VTIreader Data:" << std::endl;
87 clout << "Origin: ";
88 for (auto& origin_i : _origin) {
89 clout << origin_i << " ";
90 }
91 clout << std::endl;
92
93 clout << "Extend: ";
94 for (auto& extend_i : _extent) {
95 clout << extend_i << " ";
96 }
97 clout << std::endl;
98
99 clout << "Spacing: " << _delta << std::endl;
100}
101
102template<typename T>
103std::vector<int> BaseVTIreader<T>::readExtent(const XMLreader* reader, std::string extAttrName)
104{
105 // An extent is in the form of four (or six) integers "x0 x1 y0 y1 z0 z1"
106 // each representing a node number
107 std::stringstream extstr(reader->getAttribute(extAttrName));
108 std::vector<int> extents;
109 int tmp;
110 for (int i = 0; i < 2 * _dim; ++i) {
111 extstr >> tmp;
112 extents.push_back(tmp);
113 }
114 return extents;
115}
116
117template<typename T>
118std::vector<int> BaseVTIreader<T>::getNbNodes(std::vector<int>& extents)
119{
120 // Convert 4D (or 6D) extents vector into 2D (3D) extent vector
121 std::vector<int> nNodes;
122 for ( int i = 0; i < _dim; i++ ) {
123 nNodes.push_back(extents[ 2*i + 1 ] - extents[ 2* i ] + 1);
124 }
125 return nNodes;
126}
127
128template<typename T>
130{
131 // read the NumberOfComponents-Attribute (VTI standard) and return as integer
132 int size = std::atoi((tag.getAttribute("NumberOfComponents")).c_str());
133 if (size == 0) {
134 clout << "VTI READER: NumberOfComponents zero or not given!" << std::endl;
135 clout << "EXAMPLE: <DataArray Name='physVelocity' NumberOfComponents='3'" << std::endl;
136 exit(1);
137 }
138 return size;
139}
140
141
142/* -------------------- BaseVTIreader3D --------------------- */
143
144template<typename T, typename BaseType>
145BaseVTIreader3D<T,BaseType>::BaseVTIreader3D( const std::string& fName, std::string dName,
146 const std::string class_name)
147 : BaseVTIreader<T>(fName, 3, dName, class_name)
148{
149}
150
151template<typename T, typename BaseType>
153{
154 if (piece->getName() == "Piece") {
155 std::vector<int> extents = this->readExtent(piece, "Extent");
156 std::vector<int> extent = this->getNbNodes(extents);
157 // int extents[i] is node number => multiply with _delta to get coordinate
158 cuboid.init(extents[0] * this->_delta,
159 extents[2] * this->_delta,
160 extents[4] * this->_delta,
161 this->_delta,
162 extent[0],
163 extent[1],
164 extent[2]);
165 }
166}
167
168template<typename T, typename BaseType>
170 const XMLreader& pieceTag, const std::string dName)
171{
172 /*** This is the main data reader method. All other methods use this one for accessing data ***/
173
174 // Calculate number of cells in each direction
175 std::vector<int> extents = this->readExtent(&pieceTag, "Extent");
176 std::vector<int> extent = this->getNbNodes(extents);
177
178 // Iterate through all <DataArray> tags and take the one with the given Name attribute
179 for (auto & dataArray : pieceTag["PointData"]) {
180 if (dataArray->getAttribute("Name") == dName && dataArray->getName() == "DataArray") {
181 std::string data_str;
182 if (dataArray->read(data_str)) {
183 std::stringstream stream_val(data_str);
184
185 // Careful: respect ordering in VTI File
186 for (int iz = 0; iz < blockData.getNz(); iz++) {
187 for (int iy = 0; iy < blockData.getNy(); iy++) {
188 for (int ix = 0; ix < blockData.getNx(); ix++) {
189 for (unsigned iSize=0; iSize < blockData.getSize(); iSize++) {
190 BaseType tmp;
191 stream_val >> tmp;
192 // write tmp into blockData
193 blockData.get({ix, iy, iz}, iSize) = tmp;
194 }
195 }
196 }
197 }
198 }
199
200 // set correct blockData Name (GenericF)
201 // blockData.getName() = dName;
202
203 return true;
204 }
205 }
206 return false;
207}
208
209
210/* ---------------- BlockVTIreader3D -------------------*/
211
212template<typename T,typename BaseType>
213BlockVTIreader3D<T,BaseType>::BlockVTIreader3D(const std::string& fName, const std::string& dName )
214 : BaseVTIreader3D<T,BaseType>(fName, dName, "BlockVTIreader3D"),
215 _cuboid(this->_origin, this->_delta, this->_extent),
216 _blockData(_cuboid, 0, this->_size)
217{
218 // Only read the first <Piece> tag in the XML file
219 this->readBlockData(_blockData, this->_xmlReader["ImageData"]["Piece"], dName);
220}
221
222
223template<typename T,typename BaseType>
228
229template<typename T,typename BaseType>
234
235
236/* --------------- SuperVTIreader3D ---------------------*/
237
238template<typename T,typename BaseType>
240{
241 delete _loadBalancer;
242 delete _cGeometry;
243 delete _superData;
244}
245
246template<typename T,typename BaseType>
247SuperVTIreader3D<T,BaseType>::SuperVTIreader3D(const std::string& fName, const std::string dName )
248 : BaseVTIreader3D<T,BaseType>(fName, dName, "SuperVTIreader3D")
249{
250 this->clout << "Start reading \"" << fName << "\"... "
251 << "(" << this->_nCuboids << " cuboids)" << std::endl;
252
253 // Create CuboidGeometry
254 _cGeometry = new CuboidGeometry3D<T> (this->_origin, this->_delta, this->_extent);
255
256 this->clout << "* Reading Cuboid Geometry..." << std::endl;
257
258 // Fill CuboidGeometry
259 readCuboidGeometry();
260
261 _cGeometry->printExtended();
262
263 // Create LoadBalancer
264 _loadBalancer = new HeuristicLoadBalancer<T> (*_cGeometry);
265
266 // Create SuperData (allocation of the data, this->_size is already known!)
267 _superData = new SuperData<3,T,BaseType> ( *_cGeometry, *_loadBalancer, 2, this->_size);
268
269 this->clout << "* Reading BlockData..." << std::endl;
270
271 // Fill data objects
272 readSuperData(dName);
273
274 this->clout << "VTI Reader finished." << std::endl;
275}
276
277template<typename T,typename BaseType>
279{
280 //_cGeometry->clearCuboids();
281 for ( auto& piece : this->_xmlReader["ImageData"] ) {
282 if (piece->getName() == "Piece") {
283 std::vector<int> extents = this->readExtent(piece, "Extent");
284 std::vector<int> extent = this->getNbNodes(extents);
285 // int extent[i] is node number => multiply with _delta to get coordinate
286 Cuboid3D<T> cuboid(extents[0] * this->_delta,
287 extents[2] * this->_delta,
288 extents[4] * this->_delta,
289 this->_delta,
290 extent[0],
291 extent[1],
292 extent[2]);
293 _cGeometry->add(cuboid);
294 }
295 }
296}
297
298template<typename T,typename BaseType>
299void SuperVTIreader3D<T,BaseType>::readSuperData(const std::string dName)
300{
301 int counter = 0;
302 // Iterate over all <Piece> tags
303 for (auto & piece : this->_xmlReader["ImageData"]) {
304 if (piece->getName() == "Piece") {
305 this->readBlockData(_superData->getBlock(counter), *piece, dName);
306 counter++;
307 }
308 }
309}
310
311template<typename T,typename BaseType>
316
317template<typename T,typename BaseType>
322
323template<typename T,typename BaseType>
328
329
330
331
332
333
334// ************************************ old 2D ********************************************** //
335
336/*
337template<typename T>
338void VTIreader2D<T>::printInfo()
339{
340 clout << "Origin: " << _x0 << " " << _y0 << std::endl;
341 clout << "Extend: " << _x << " " << _y << std::endl;
342 clout << "Spacing: " << _delta << std::endl;
343}
344
345template<typename T>
346VTIreader2D<T>::VTIreader2D(const std::string& fName ) : XMLreader(fName)
347{
348 // get _x, _y, _z
349 int x, y, z;
350 std::stringstream stream_val_0((*this)["ImageData"].getAttribute("WholeExtent"));
351 stream_val_0 >> x >> _x >> y >> _y >> z >> _z;
352 _x = _x-x+1;
353 _y = _y-y+1;
354 _z = _z-z+1;
355
356 // get _delta
357 std::stringstream stream_val_1((*this)["ImageData"].getAttribute("Spacing"));
358 stream_val_1 >> _delta;
359
360 std::stringstream stream_val_2((*this)["ImageData"].getAttribute("Origin"));
361 stream_val_2 >> _x0 >> _y0 >> _z0;
362}
363
364template<typename T>
365VTIreader2D<T>::~VTIreader2D() {}
366
367template<typename T>
368void VTIreader2D<T>::getCuboid(Cuboid2D<T>& cuboid)
369{
370 cuboid.init(_x0, _y0, _delta, _x, _y);
371}
372
373template<typename T>
374void VTIreader2D<T>::getCuboids(std::vector<Cuboid2D<T>* >& cuboids)
375{
376 std::vector<XMLreader*>::const_iterator it;
377 int i = 0;
378 for (it = (*this)["ImageData"].begin(); it != (*this)["ImageData"].end(); it++) {
379 if ((*it)->getName() == "Piece") {
380 ++i;
381 std::stringstream extstr((*it)->getAttribute("Extent"));
382 int extents[6];
383 for (int i = 0; i < 6; ++i) {
384 extstr >> extents[i];
385 }
386 Cuboid2D<T>* cuboid = new Cuboid2D<T>(extents[0], extents[2], _delta,
387 extents[1]-extents[0]+1,
388 extents[3]-extents[2]+1);
389 cuboids.push_back(cuboid);
390 }
391 }
392}
393*/
394/*
395template<typename T>
396void VTIreader2D<T>::getScalarMultiPieceData(std::vector<const ScalarFieldBase2D<T>* >& bases, const std::string dName)
397{
398 std::vector<XMLreader*>::const_iterator it;
399 for (it = (*this)["ImageData"].begin(); it != (*this)["ImageData"].end(); it++) {
400 if ((*it)->getName() == "Piece") {
401 std::stringstream extstr((*it)->getAttribute("Extent"));
402 int extents[6];
403 for (int i = 0; i < 6; i++) {
404 extstr >> extents[i];
405 }
406 std::vector<XMLreader*>::const_iterator it2;
407 for (it2 = (*it)->operator[]("PointData").begin(); it2 != (*it)->operator[]("PointData").end(); it2++) {
408 if ((*it2)->getAttribute("Name") == dName && (*it2)->getName() == "DataArray") {
409 ScalarField2D<T>* tmp = new ScalarField2D<T>(extents[1]-extents[0]+1, extents[3]-extents[2]+1);
410 tmp->construct();
411 std::string data_str;
412 if ((*it2)->read(data_str)) {
413 std::stringstream stream_val(data_str);
414 for (int iz = 0; iz < extents[5]-extents[4]+1; iz++) {
415 for (int iy = 0; iy < extents[3]-extents[2]+1; iy++) {
416 for (int ix = 0; ix < extents[1]-extents[0]+1; ix++) {
417 T tmp2;
418 stream_val >> tmp2;
419 tmp->get(ix, iy) = tmp2;
420 }
421 }
422 }
423 }
424 bases.push_back(tmp);
425 }
426 }
427 }
428 }
429}
430
431template<typename T>
432void VTIreader2D<T>::getVectorMultiPieceData(std::vector<const TensorFieldBase2D<T, 2>* >& bases, const std::string dName)
433{
434 std::vector<XMLreader*>::const_iterator it;
435 for (it = (*this)["ImageData"].begin(); it != (*this)["ImageData"].end(); it++) {
436 if ((*it)->getName() == "Piece") {
437 std::stringstream extstr((*it)->getAttribute("Extent"));
438 int extents[6];
439 for (int i=0; i<6; i++) {
440 extstr >> extents[i];
441 }
442 std::vector<XMLreader*>::const_iterator it2;
443 for (it2 = (*it)->operator[]("PointData").begin(); it2 != (*it)->operator[]("PointData").end(); it2++) {
444 if ((*it2)->getAttribute("Name") == dName && (*it2)->getName() == "DataArray") {
445 TensorField2D<T, 2>* tmp = new TensorField2D<T, 2>(extents[1]-extents[0]+1, extents[3]-extents[2]+1);
446 tmp->construct();
447 std::string data_str;
448 if ((*it2)->read(data_str)) {
449 std::stringstream stream_val(data_str);
450 for (int iz = 0; iz < extents[5]-extents[4]+1; iz++) {
451 for (int iy = 0; iy < extents[3]-extents[2]+1; iy++) {
452 for (int ix = 0; ix < extents[1]-extents[0]+1; ix++) {
453 T tmpval;
454 for (int i=0; i < 2; i++) {
455 stream_val >> tmpval;
456 tmp->get(ix, iy)[i] = tmpval;
457 }
458 stream_val >> tmpval;
459 }
460 }
461 }
462 }
463 bases.push_back(tmp);
464 }
465 }
466 }
467 }
468}
469
470template<typename T>
471bool VTIreader2D<T>::getScalarData(ScalarField2D<T>* base, const std::string dName)
472{
473 std::vector<XMLreader*>::const_iterator it;
474 for (it = (*this)["ImageData"]["Piece"]["PointData"].begin(); it != (*this)["ImageData"]["Piece"]["PointData"].end(); it++) {
475 if ((*it)->getAttribute("Name") == dName && (*it)->getName() == "DataArray") {
476 ScalarField2D<T>* tmp = new ScalarField2D<T>(_x, _y);
477 tmp->construct();
478 std::string data_str;
479 if ((*it)->read(data_str)) {
480 std::stringstream stream_val(data_str);
481 for (int iz = 0; iz < _z; iz++) {
482 for (int iy = 0; iy < _y; iy++) {
483 for (int ix = 0; ix < _x; ix++) {
484 T tmp2;
485 stream_val >> tmp2;
486 tmp->get(ix, iy) = tmp2;
487 }
488 }
489 }
490 }
491 base->swap(*tmp);
492 delete tmp;
493 return true;
494 }
495 }
496 return false;
497}
498
499template<typename T>
500bool VTIreader2D<T>::getVectorData(TensorField2D<T, 2>* base, const std::string dName)
501{
502 std::vector<XMLreader*>::const_iterator it;
503 for (it = (*this)["ImageData"]["Piece"]["PointData"].begin(); it != (*this)["ImageData"]["Piece"]["PointData"].end(); it++) {
504 if ((*it)->getAttribute("Name") == dName && (*it)->getName() == "DataArray") {
505 TensorField2D<T, 2>* tmp = new TensorField2D<T, 2>(_x, _y);
506 tmp->construct();
507 std::string data_str;
508 if ((*it)->read(data_str)) {
509 std::stringstream stream_val(data_str);
510 for (int iz = 0; iz < _z; iz++) {
511 for (int iy = 0; iy < _y; iy++) {
512 for (int ix = 0; ix < _x; ix++) {
513 T tmpval;
514 for (int i = 0; i < 2; i++) {
515 stream_val >> tmpval;
516 tmp->get(ix, iy)[i] = tmpval;
517 }
518 stream_val >> tmpval;
519 }
520 }
521 }
522 }
523 base->swap(*tmp);
524 delete tmp;
525 return true;
526 }
527 }
528 return false;
529}
530*/
531
532
533}
534#endif
BaseVTIreader3D(const std::string &fName, std::string dName, const std::string class_name="BaseVTIreader3D")
Definition vtiReader.hh:145
bool readBlockData(BlockData< 3, T, BaseType > &blockData, const XMLreader &pointDataTag, const std::string dName)
Reads from pointDataTag and fills blockData.
Definition vtiReader.hh:169
void readCuboid(Cuboid3D< T > &cuboid, XMLreader *piece)
Reads cuboid from piece node.
Definition vtiReader.hh:152
std::vector< int > getNbNodes(std::vector< int > &extents)
Converts 4D (or 6D) extents vector into 2D (3D) nb_nodes vector.
Definition vtiReader.hh:118
int _size
Size of Data Field.
Definition vtiReader.h:90
std::vector< T > _origin
Definition vtiReader.h:92
OstreamManager clout
Definition vtiReader.h:86
std::vector< int > readExtent(const XMLreader *reader, std::string extAttrName)
Reads Extent from extAttrName from XML Tag and returns as vector.
Definition vtiReader.hh:103
std::vector< int > _extent
Definition vtiReader.h:94
XMLreader _xmlReader
Definition vtiReader.h:96
BaseVTIreader(const std::string &fName, int dim, std::string dName, const std::string class_name="BaseVTIreader")
Definition vtiReader.hh:47
int getSize(const XMLreader &tag)
Reads size from XML tag (attribute "NumberOfComponents")
Definition vtiReader.hh:129
unsigned getSize() const
Definition blockData.hh:118
U & get(std::size_t iCell, int iD=0)
Definition blockData.hh:94
int getNy() const
Read only access to block height.
int getNx() const
Read only access to block width.
int getNz() const
Read only access to block height.
Cuboid3D< T > & getCuboid()
Definition vtiReader.hh:230
BlockVTIreader3D(const std::string &fName, const std::string &dName)
Definition vtiReader.hh:213
BlockData< 3, T, BaseType > _blockData
Definition vtiReader.h:133
BlockData< 3, T, BaseType > & getBlockData()
Definition vtiReader.hh:224
A regular single 3D cuboid is the basic component of a 3D cuboid structure which defines the grid.
Definition cuboid3D.h:58
void init(T globPosX, T globPosY, T globPosZ, T delta, int nX, int nY, int nZ)
Initializes the cuboid.
Definition cuboid3D.hh:111
A cuboid geometry represents a voxel mesh.
Constructs a load balancer from a given cuboid geometry using a heurist.
Base class for all LoadBalancer.
~SuperVTIreader3D() override
Definition vtiReader.hh:239
LoadBalancer< T > & getLoadBalancer()
Definition vtiReader.hh:324
SuperData< 3, T, BaseType > & getSuperData()
Definition vtiReader.hh:312
SuperVTIreader3D(const std::string &fName, const std::string dName)
Definition vtiReader.hh:247
CuboidGeometry3D< T > & getCuboidGeometry()
Definition vtiReader.hh:318
std::string getAttribute(const std::string &aName) const
Definition xmlReader.h:541
std::string getName() const
return the name of the element
Definition xmlReader.h:402
The description of a single 2D cuboid – header file.
Top level namespace for all of OpenLB.
typename util::BaseTypeHelper< T >::type BaseType
Definition baseType.h:59
The VTI reader is able to read from VTI files and create and fill corresponding data structures.