00001 #ifndef VIENNAGRID_IO_VTK_READER_HPP
00002 #define VIENNAGRID_IO_VTK_READER_HPP
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00026 #include <fstream>
00027 #include <iostream>
00028 #include <sstream>
00029 #include <vector>
00030 #include <deque>
00031 #include <string>
00032 #include <algorithm>
00033
00034 #include "viennagrid/forwards.h"
00035 #include "viennagrid/point.hpp"
00036 #include "viennagrid/io/vtk_common.hpp"
00037 #include "viennagrid/io/data_accessor.hpp"
00038 #include "viennagrid/io/helper.hpp"
00039 #include "viennagrid/io/xml_tag.hpp"
00040
00041 namespace viennagrid
00042 {
00043 namespace io
00044 {
00045
00050 template <typename DomainType>
00051 class vtk_reader
00052 {
00053 protected:
00054 typedef typename DomainType::config_type DomainConfiguration;
00055
00056 typedef typename DomainConfiguration::numeric_type CoordType;
00057 typedef typename DomainConfiguration::coordinate_system_tag CoordinateSystemTag;
00058 typedef typename DomainConfiguration::cell_tag CellTag;
00059
00060 typedef typename result_of::point<DomainConfiguration>::type PointType;
00061 typedef typename result_of::ncell<DomainConfiguration, 0>::type VertexType;
00062 typedef typename result_of::ncell<DomainConfiguration, CellTag::dim>::type CellType;
00063
00064
00065 typedef typename viennagrid::result_of::ncell_range<DomainType, 0>::type VertexRange;
00066 typedef typename viennagrid::result_of::iterator<VertexRange>::type VertexIterator;
00067
00068 typedef typename viennagrid::result_of::ncell_range<DomainType, 1>::type EdgeRange;
00069 typedef typename viennagrid::result_of::iterator<EdgeRange>::type EdgeIterator;
00070
00071 typedef typename viennagrid::result_of::ncell_range<DomainType, CellTag::dim-1>::type FacetRange;
00072 typedef typename viennagrid::result_of::iterator<FacetRange>::type FacetIterator;
00073
00074 typedef typename viennagrid::result_of::ncell_range<DomainType, CellTag::dim>::type CellRange;
00075 typedef typename viennagrid::result_of::iterator<CellRange>::type CellIterator;
00076
00077 std::ifstream reader;
00078 std::map<PointType, std::size_t, point_less> global_points;
00079 std::map<std::size_t, PointType> global_points_2;
00080 std::deque<std::deque<std::size_t> > local_to_global_map;
00081 std::deque<std::deque<std::size_t> > local_cell_vertices;
00082 std::deque<std::deque<std::size_t> > local_cell_offsets;
00083 std::deque<std::size_t> local_cell_num;
00084
00085
00086 std::deque<std::deque<std::pair<std::string, std::deque<double> > > > local_scalar_vertex_data;
00087 std::deque<std::deque<std::pair<std::string, std::deque<double> > > > local_vector_vertex_data;
00088 std::deque<std::deque<std::pair<std::string, std::deque<double> > > > local_scalar_cell_data;
00089 std::deque<std::deque<std::pair<std::string, std::deque<double> > > > local_vector_cell_data;
00090
00092 void openFile(std::string const & filename)
00093 {
00094 reader.open(filename.c_str());
00095 if(!reader)
00096 {
00097 throw cannot_open_file_exception(filename);
00098 }
00099 }
00100
00102 void closeFile()
00103 {
00104 reader.close();
00105 }
00106
00108 bool lowercase_compare(std::string const & s1, std::string const & s2)
00109 {
00110 std::string s1_lower = s1;
00111 std::transform(s1.begin(), s1.end(), s1_lower.begin(), char_to_lower<>());
00112
00113 std::string s2_lower = s2;
00114 std::transform(s2.begin(), s2.end(), s2_lower.begin(), char_to_lower<>());
00115
00116 return s1_lower == s2_lower;
00117 }
00118
00120 void checkNextToken(std::string const & expectedToken)
00121 {
00122 std::string token;
00123 reader >> token;
00124
00125 if ( !lowercase_compare(token, expectedToken) )
00126 {
00127 std::cerr << "* vtk_reader::operator(): Expected \"" << expectedToken << "\", but got \"" << token << "\"" << std::endl;
00128 std::stringstream ss;
00129 ss << "Expected \"" << expectedToken << "\", but got \"" << token << "\"";
00130 throw bad_file_format_exception("", ss.str());
00131 }
00132 }
00133
00135 void readNodeCoordinates(long nodeNum, long numberOfComponents, std::size_t seg_id)
00136 {
00137 double nodeCoord;
00138 local_to_global_map[seg_id].resize(nodeNum);
00139
00140 for(long i = 0; i < nodeNum; i++)
00141 {
00142 PointType p;
00143
00144 for(long j = 0; j < numberOfComponents; j++)
00145 {
00146 reader >> nodeCoord;
00147 if (j < CoordinateSystemTag::dim)
00148 p[j] = nodeCoord;
00149
00150 }
00151
00152
00153 if (global_points.find(p) == global_points.end())
00154 {
00155 std::size_t new_global_id = global_points.size();
00156
00157 global_points.insert( std::make_pair(p, new_global_id) );
00158 global_points_2.insert( std::make_pair(new_global_id, p) );
00159 local_to_global_map[seg_id][i] = new_global_id;
00160 }
00161 else
00162 {
00163 local_to_global_map[seg_id][i] = global_points[p];
00164
00165 }
00166 }
00167 }
00168
00170 void readCellIndices(std::size_t seg_id)
00171 {
00172 long cellNode = 0;
00173 std::string token;
00174 reader >> token;
00175
00176 while (token != "</DataArray>")
00177 {
00178 assert( strChecker::myIsNumber(token) && "Cell vertex index is not a number!" );
00179
00180 cellNode = atoi(token.c_str());
00181 local_cell_vertices[seg_id].push_back(cellNode);
00182
00183 reader >> token;
00184 }
00185 }
00186
00188 void readOffsets(std::size_t seg_id)
00189 {
00190
00191
00192
00193
00194
00195 std::string token;
00196 long offset = 0;
00197 reader >> token;
00198
00199 while(token != "</DataArray>")
00200 {
00201 assert( strChecker::myIsNumber(token) && "Cell offset is not a number!" );
00202
00203 offset = atoi(token.c_str());
00204 local_cell_offsets[seg_id].push_back(offset);
00205
00206
00207 reader >> token;
00208 }
00209 }
00210
00212 void readTypes()
00213 {
00214
00215
00216
00217
00218
00219 std::string token;
00220 long type = 0;
00221 reader >> token;
00222
00223 while(token != "</DataArray>")
00224 {
00225 assert( strChecker::myIsNumber(token) && "Cell type is not a number!" );
00226
00227 type = atoi(token.c_str());
00228 assert(type == ELEMENT_TAG_TO_VTK_TYPE<CellTag>::value && "Error in VTK reader: Type mismatch!");
00229
00230
00231 reader >> token;
00232 }
00233 }
00234
00236 void readData(std::deque<double> & container)
00237 {
00238 std::string token;
00239 reader >> token;
00240
00241 while (string_to_lower(token) != "</dataarray>")
00242 {
00243 container.push_back( atof(token.c_str()) );
00244 reader >> token;
00245 }
00246 }
00247
00249 template <typename ContainerType, typename NameContainerType>
00250 void readPointCellData(size_t seg_id,
00251 ContainerType & scalar_data,
00252 ContainerType & vector_data,
00253 NameContainerType & data_names_scalar,
00254 NameContainerType & data_names_vector)
00255 {
00256 std::string name;
00257 std::size_t components = 1;
00258
00259 xml_tag<> tag;
00260
00261 tag.parse(reader);
00262
00263 while (tag.name() == "dataarray")
00264 {
00265 tag.check_attribute("name", "");
00266 name = tag.get_value("name");
00267
00268 if (tag.has_attribute("numberofcomponents"))
00269 components = atoi(tag.get_value("numberofcomponents").c_str());
00270
00271
00272 if (components == 1)
00273 {
00274 data_names_scalar.push_back(std::make_pair(seg_id, name));
00275 scalar_data[seg_id].push_back( std::make_pair(name, std::deque<double>()) );
00276 readData(scalar_data[seg_id].back().second);
00277 }
00278 else if (components == 3)
00279 {
00280 data_names_vector.push_back(std::make_pair(seg_id, name));
00281 vector_data[seg_id].push_back( std::make_pair(name, std::deque<double>()) );
00282 readData(vector_data[seg_id].back().second);
00283 }
00284 else
00285 throw bad_file_format_exception("", "Number of components for data invalid!");
00286
00287 tag.parse(reader);
00288 }
00289
00290
00291 if (tag.name() != "/pointdata" && tag.name() != "/celldata")
00292 throw bad_file_format_exception("", "XML Parse error: Expected </PointData> or </CellData>!");
00293
00294 }
00295
00297
00299 void setupVertices(DomainType & domain)
00300 {
00301 for (std::size_t i=0; i<global_points_2.size(); ++i)
00302 {
00303 VertexType v;
00304 v.point() = global_points_2[i];
00305 domain.push_back(v);
00306 }
00307 }
00308
00310 void setupCells(DomainType & domain, std::size_t seg_id)
00311 {
00312
00313
00314
00315
00316
00317
00318
00319
00320 CellType cell;
00321 long numVertices = 0;
00322 long offsetIdx = 0;
00323
00324 std::deque<std::size_t> const & offsets = local_cell_offsets[seg_id];
00325
00326 for (std::size_t i = 0; i < local_cell_num[seg_id]; i++)
00327 {
00328
00329
00330
00331
00332
00333 if(i == 0)
00334 {
00335 offsetIdx = 0;
00336 numVertices = offsets[i];
00337 }
00338 else
00339 {
00340 offsetIdx = offsets[i-1];
00341 numVertices = offsets[i]-offsets[i-1];
00342 }
00343
00344 std::vector<VertexType *> vertices(numVertices);
00345
00346
00347
00348
00349
00350
00351
00352 vtk_to_viennagrid_orientations<CellTag> reorderer;
00353 for (long j = 0; j < numVertices; j++)
00354 {
00355 long reordered_j = reorderer(j);
00356 std::size_t local_index = local_cell_vertices[seg_id][reordered_j + offsetIdx];
00357 std::size_t global_vertex_index = local_to_global_map[seg_id][local_index];
00358 vertices[j] = &(viennagrid::ncells<0>(domain)[global_vertex_index]);
00359
00360 }
00361
00362 cell.vertices(&(vertices[0]));
00363 cell.id(i);
00364 domain.segments()[seg_id].push_back(cell);
00365 }
00366 }
00367
00369 template <typename ContainerType>
00370 void setupDataVertex(DomainType & domain, std::size_t seg_id, ContainerType const & container, std::size_t num_components)
00371 {
00372
00373
00374
00375 bool name_registered = false;
00376 std::size_t data_accessor_index = 0;
00377 if (num_components == 1)
00378 {
00379 for (std::size_t i=0; i<vertex_data_scalar_names.size(); ++i)
00380 {
00381 if (container.first == vertex_data_scalar_names[i])
00382 {
00383 name_registered = true;
00384 data_accessor_index = i;
00385 #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00386 std::cout << "* vtk_reader::operator(): Found registered scalar name " << container.first << std::endl;
00387 #endif
00388 break;
00389 }
00390 }
00391 }
00392 else
00393 {
00394 for (std::size_t i=0; i<vertex_data_vector_names.size(); ++i)
00395 {
00396 if (container.first == vertex_data_vector_names[i])
00397 {
00398 name_registered = true;
00399 data_accessor_index = i;
00400 #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00401 std::cout << "* vtk_reader::operator(): Found registered vertex name " << container.first << std::endl;
00402 #endif
00403 break;
00404 }
00405 }
00406 }
00407
00408
00409
00410
00411 if (name_registered)
00412 {
00413 if (num_components == 1)
00414 {
00415 #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00416 std::cout << "* vtk_reader::operator(): Reading scalar quantity "
00417 << container.first << " using data_accessor to vertices." << std::endl;
00418 #endif
00419 for (std::size_t i=0; i<container.second.size(); ++i)
00420 {
00421 std::size_t global_vertex_id = local_to_global_map[seg_id][i];
00422 VertexType const & vertex = viennagrid::ncells<0>(domain)[global_vertex_id];
00423 vertex_data_scalar[data_accessor_index](vertex, seg_id, 0, (container.second)[i]);
00424 }
00425 }
00426 else
00427 {
00428 #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00429 std::cout << "* vtk_reader::operator(): Reading vector quantity "
00430 << container.first << " using data_accessor to vertices." << std::endl;
00431 #endif
00432 assert( 3 * viennagrid::ncells<0>(domain.segments()[seg_id]).size() == container.second.size() );
00433 for (std::size_t i=0; i<container.second.size() / 3; ++i)
00434 {
00435 std::size_t global_vertex_id = local_to_global_map[seg_id][i];
00436 VertexType const & vertex = viennagrid::ncells<0>(domain)[global_vertex_id];
00437 vertex_data_vector[data_accessor_index](vertex, seg_id, 0, (container.second)[3*i]);
00438 vertex_data_vector[data_accessor_index](vertex, seg_id, 1, (container.second)[3*i+1]);
00439 vertex_data_vector[data_accessor_index](vertex, seg_id, 2, (container.second)[3*i+2]);
00440 }
00441 }
00442 }
00443 else
00444 {
00445 if (num_components == 1)
00446 {
00447 #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00448 std::cout << "* vtk_reader::operator(): Reading scalar quantity "
00449 << container.first << " using viennadata to vertices." << std::endl;
00450 #endif
00451 for (std::size_t i=0; i<container.second.size(); ++i)
00452 {
00453 std::size_t global_vertex_id = local_to_global_map[seg_id][i];
00454 VertexType const & vertex = viennagrid::ncells<0>(domain)[global_vertex_id];
00455 viennadata::access<std::string, double>(container.first)(vertex) = (container.second)[i];
00456 }
00457 }
00458 else
00459 {
00460 #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00461 std::cout << "* vtk_reader::operator(): Reading vector quantity "
00462 << container.first << " using viennadata to vertices." << std::endl;
00463 #endif
00464 assert( 3 * viennagrid::ncells<0>(domain.segments()[seg_id]).size() == container.second.size());
00465 for (std::size_t i=0; i<container.second.size() / 3; ++i)
00466 {
00467 std::size_t global_vertex_id = local_to_global_map[seg_id][i];
00468 VertexType const & vertex = viennagrid::ncells<0>(domain)[global_vertex_id];
00469 viennadata::access<std::string, std::vector<double> >(container.first)(vertex).resize(3);
00470 viennadata::access<std::string,
00471 std::vector<double> >(container.first)(vertex)[0] = (container.second)[3*i];
00472 viennadata::access<std::string,
00473 std::vector<double> >(container.first)(vertex)[1] = (container.second)[3*i+1];
00474 viennadata::access<std::string,
00475 std::vector<double> >(container.first)(vertex)[2] = (container.second)[3*i+2];
00476 }
00477 }
00478 }
00479
00480 }
00481
00483 template <typename ContainerType>
00484 void setupDataCell(DomainType & domain, std::size_t seg_id, ContainerType const & container, std::size_t num_components)
00485 {
00486
00487
00488
00489 bool name_registered = false;
00490 std::size_t data_accessor_index = 0;
00491 if (num_components == 1)
00492 {
00493 for (std::size_t i=0; i<cell_data_scalar_names.size(); ++i)
00494 {
00495 if (container.first == cell_data_scalar_names[i])
00496 {
00497 name_registered = true;
00498 data_accessor_index = i;
00499 #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00500 std::cout << "* vtk_reader::operator(): Found registered cell scalar name " << container.first << std::endl;
00501 #endif
00502 break;
00503 }
00504 }
00505 }
00506 else
00507 {
00508 for (std::size_t i=0; i<cell_data_vector_names.size(); ++i)
00509 {
00510 if (container.first == cell_data_vector_names[i])
00511 {
00512 name_registered = true;
00513 data_accessor_index = i;
00514 #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00515 std::cout << "* vtk_reader::operator(): Found registered cell vector name " << container.first << std::endl;
00516 #endif
00517 break;
00518 }
00519 }
00520 }
00521
00522
00523
00524
00525 if (name_registered)
00526 {
00527 if (num_components == 1)
00528 {
00529 #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00530 std::cout << "* vtk_reader::operator(): Reading scalar quantity "
00531 << container.first << " using data_accessor to cells." << std::endl;
00532 #endif
00533 for (std::size_t i=0; i<container.second.size(); ++i)
00534 {
00535 CellType const & cell = viennagrid::ncells<CoordinateSystemTag::dim>(domain.segments()[seg_id])[i];
00536 cell_data_scalar[data_accessor_index](cell, seg_id, 0, (container.second)[i]);
00537 }
00538 }
00539 else
00540 {
00541 #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00542 std::cout << "* vtk_reader::operator(): Reading vector quantity "
00543 << container.first << " using data_accessor to cells." << std::endl;
00544 #endif
00545 assert( 3 * viennagrid::ncells<CoordinateSystemTag::dim>(domain.segments()[seg_id]).size() == container.second.size());
00546 for (std::size_t i=0; i<container.second.size() / 3; ++i)
00547 {
00548 CellType const & cell = viennagrid::ncells<CoordinateSystemTag::dim>(domain.segments()[seg_id])[i];
00549 cell_data_vector[data_accessor_index](cell, seg_id, 0, (container.second)[3*i]);
00550 cell_data_vector[data_accessor_index](cell, seg_id, 1, (container.second)[3*i+1]);
00551 cell_data_vector[data_accessor_index](cell, seg_id, 2, (container.second)[3*i+2]);
00552 }
00553 }
00554 }
00555 else
00556 {
00557 if (num_components == 1)
00558 {
00559 #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00560 std::cout << "* vtk_reader::operator(): Reading scalar quantity "
00561 << container.first << " using viennadata to cells." << std::endl;
00562 #endif
00563 for (std::size_t i=0; i<container.second.size(); ++i)
00564 {
00565 CellType const & cell = viennagrid::ncells<CoordinateSystemTag::dim>(domain.segments()[seg_id])[i];
00566 viennadata::access<std::string,
00567 double>(container.first)(cell) = (container.second)[i];
00568 }
00569 }
00570 else
00571 {
00572 #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00573 std::cout << "* vtk_reader::operator(): Reading vector quantity "
00574 << container.first << " using viennadata to cells." << std::endl;
00575 #endif
00576 assert( 3 * viennagrid::ncells<CoordinateSystemTag::dim>(domain.segments()[seg_id]).size() == container.second.size());
00577 for (std::size_t i=0; i<container.second.size() / 3; ++i)
00578 {
00579 CellType const & cell = viennagrid::ncells<CoordinateSystemTag::dim>(domain.segments()[seg_id])[i];
00580 viennadata::access<std::string, std::vector<double> >(container.first)(cell).resize(3);
00581 viennadata::access<std::string,
00582 std::vector<double> >(container.first)(cell)[0] = (container.second)[3*i];
00583 viennadata::access<std::string,
00584 std::vector<double> >(container.first)(cell)[1] = (container.second)[3*i+1];
00585 viennadata::access<std::string,
00586 std::vector<double> >(container.first)(cell)[2] = (container.second)[3*i+2];
00587 }
00588 }
00589 }
00590
00591 }
00592
00594 void setupData(DomainType & domain, std::size_t seg_id)
00595 {
00596 for (size_t i=0; i<local_scalar_vertex_data[seg_id].size(); ++i)
00597 {
00598 setupDataVertex(domain, seg_id, local_scalar_vertex_data[seg_id][i], 1);
00599 }
00600
00601 for (size_t i=0; i<local_vector_vertex_data[seg_id].size(); ++i)
00602 {
00603 setupDataVertex(domain, seg_id, local_vector_vertex_data[seg_id][i], 3);
00604 }
00605
00606
00607 for (size_t i=0; i<local_scalar_cell_data[seg_id].size(); ++i)
00608 {
00609 setupDataCell(domain, seg_id, local_scalar_cell_data[seg_id][i], 1);
00610 }
00611
00612 for (size_t i=0; i<local_vector_cell_data[seg_id].size(); ++i)
00613 {
00614 setupDataCell(domain, seg_id, local_vector_cell_data[seg_id][i], 3);
00615 }
00616
00617 }
00618
00620 void parse_vtu_segment(std::string filename, std::size_t seg_id)
00621 {
00622
00623 try
00624 {
00625 openFile(filename);
00626
00627 long nodeNum = 0;
00628 long numberOfComponents = 0;
00629
00630 xml_tag<> tag;
00631
00632 tag.parse(reader);
00633 if (tag.name() != "?xml" && tag.name() != "?xml?")
00634 throw bad_file_format_exception(filename, "Parse error: No opening ?xml tag!");
00635
00636 tag.parse_and_check_name(reader, "vtkfile", filename);
00637 tag.parse_and_check_name(reader, "unstructuredgrid", filename);
00638
00639 tag.parse_and_check_name(reader, "piece", filename);
00640
00641 tag.check_attribute("numberofpoints", filename);
00642
00643 nodeNum = atoi(tag.get_value("numberofpoints").c_str());
00644 #ifdef VIENNAGRID_DEBUG_IO
00645 std::cout << "#Nodes: " << nodeNum << std::endl;
00646 #endif
00647
00648 tag.check_attribute("numberofcells", filename);
00649
00650 local_cell_num[seg_id] = atoi(tag.get_value("numberofcells").c_str());
00651 #ifdef VIENNAGRID_DEBUG_IO
00652 std::cout << "#Cells: " << local_cell_num[seg_id] << std::endl;
00653 #endif
00654
00655 tag.parse_and_check_name(reader, "points", filename);
00656
00657 tag.parse_and_check_name(reader, "dataarray", filename);
00658 tag.check_attribute("numberofcomponents", filename);
00659
00660 numberOfComponents = atoi(tag.get_value("numberofcomponents").c_str());
00661 readNodeCoordinates(nodeNum, numberOfComponents, seg_id);
00662
00663 tag.parse_and_check_name(reader, "/dataarray", filename);
00664 tag.parse_and_check_name(reader, "/points", filename);
00665
00666 tag.parse(reader);
00667 if (tag.name() == "pointdata")
00668 {
00669 readPointCellData(seg_id, local_scalar_vertex_data, local_vector_vertex_data,
00670 vertex_data_scalar_read, vertex_data_vector_read);
00671 tag.parse(reader);
00672 }
00673
00674 if (tag.name() != "cells")
00675 throw bad_file_format_exception(filename, "Parse error: Expected Cells tag!");
00676
00677 for (std::size_t i=0; i<3; ++i)
00678 {
00679 tag.parse_and_check_name(reader, "dataarray", filename);
00680 tag.check_attribute("name", filename);
00681
00682 if (tag.get_value("name") == "connectivity")
00683 readCellIndices(seg_id);
00684 else if (tag.get_value("name") == "offsets")
00685 readOffsets(seg_id);
00686 else if (tag.get_value("name") == "types")
00687 readTypes();
00688 else
00689 throw bad_file_format_exception(filename, "Parse error: <DataArray> is not named 'connectivity', 'offsets' or 'types'!");
00690 }
00691
00692 tag.parse_and_check_name(reader, "/cells", filename);
00693
00694 tag.parse(reader);
00695 if (tag.name() == "celldata")
00696 {
00697 readPointCellData(seg_id, local_scalar_cell_data, local_vector_cell_data,
00698 cell_data_scalar_read, cell_data_vector_read);
00699 tag.parse(reader);
00700 }
00701
00702
00703 if (tag.name() != "/piece")
00704 throw bad_file_format_exception(filename, "Parse error: Expected </Piece> tag!");
00705
00706 tag.parse_and_check_name(reader, "/unstructuredgrid", filename);
00707 tag.parse_and_check_name(reader, "/vtkfile", filename);
00708
00709 closeFile();
00710 }
00711 catch (std::exception const & ex) {
00712 std::cerr << "Problems while reading file " << filename << std::endl;
00713 std::cerr << "what(): " << ex.what() << std::endl;
00714 }
00715
00716 }
00717
00719 void process_vtu(std::string const & filename)
00720 {
00721 local_cell_num.resize(1);
00722 local_cell_offsets.resize(1);
00723 local_cell_vertices.resize(1);
00724 local_to_global_map.resize(1);
00725
00726 local_scalar_vertex_data.resize(1);
00727 local_scalar_cell_data.resize(1);
00728 local_vector_vertex_data.resize(1);
00729 local_vector_cell_data.resize(1);
00730
00731 parse_vtu_segment(filename, 0);
00732 }
00733
00735 void process_pvd(std::string const & filename)
00736 {
00737 std::deque<std::string> filenames;
00738
00739
00740 std::string::size_type pos = filename.rfind("/");
00741 std::string path_to_pvd;
00742 if (pos == std::string::npos)
00743 pos = filename.rfind("\\");
00744
00745 if (pos != std::string::npos)
00746 path_to_pvd = filename.substr(0, pos + 1);
00747
00748 openFile(filename);
00749
00750
00751
00752
00753 xml_tag<> tag;
00754
00755 tag.parse(reader);
00756 if (tag.name() != "?xml" && tag.name() != "?xml?")
00757 throw bad_file_format_exception(filename, "Parse error: No opening <?xml?> tag!");
00758
00759 tag.parse(reader);
00760 if (tag.name() != "vtkfile")
00761 throw bad_file_format_exception(filename, "Parse error: VTKFile tag expected!");
00762
00763 if (!tag.has_attribute("type"))
00764 throw bad_file_format_exception(filename, "Parse error: VTKFile tag has no attribute 'type'!");
00765
00766 if (string_to_lower(tag.get_value("type")) != "collection")
00767 throw bad_file_format_exception(filename, "Parse error: Type-attribute of VTKFile tag is not 'Collection'!");
00768
00769
00770 tag.parse(reader);
00771 if (tag.name() != "collection")
00772 throw bad_file_format_exception(filename, "Parse error: Collection tag expected!");
00773
00774 long seg_index = 0;
00775 while (reader.good())
00776 {
00777 tag.parse(reader);
00778
00779 if (tag.name() == "/collection")
00780 break;
00781
00782 if (tag.name() != "dataset")
00783 throw bad_file_format_exception(filename, "Parse error: DataSet tag expected!");
00784
00785 if (tag.has_attribute("part"))
00786 {
00787 if (atoi(tag.get_value("part").c_str()) != seg_index)
00788 throw bad_file_format_exception(filename, "Parser limitation: Segments must be provided with increasing index!");
00789 }
00790
00791 if (!tag.has_attribute("file"))
00792 throw bad_file_format_exception(filename, "Parse error: DataSet tag has no file attribute!");
00793
00794 filenames.push_back(tag.get_value("file"));
00795
00796 ++seg_index;
00797 }
00798
00799 tag.parse(reader);
00800 if (tag.name() != "/vtkfile")
00801 throw bad_file_format_exception(filename, "Parse error: Closing VTKFile tag expected!");
00802
00803 closeFile();
00804
00805 assert(filenames.size() > 0 && "No segments in pvd-file specified!");
00806
00807
00808
00809
00810 local_cell_num.resize(filenames.size());
00811 local_cell_offsets.resize(filenames.size());
00812 local_cell_vertices.resize(filenames.size());
00813 local_to_global_map.resize(filenames.size());
00814
00815 local_scalar_vertex_data.resize(filenames.size());
00816 local_scalar_cell_data.resize(filenames.size());
00817 local_vector_vertex_data.resize(filenames.size());
00818 local_vector_cell_data.resize(filenames.size());
00819
00820
00821
00822
00823 for (std::size_t i=0; i<filenames.size(); ++i)
00824 {
00825 #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00826 std::cout << "Parsing file " << path_to_pvd + filenames[i] << std::endl;
00827 #endif
00828 parse_vtu_segment(path_to_pvd + filenames[i], i);
00829 }
00830
00831 }
00832
00833
00834 public:
00835
00841 int operator()(DomainType & domain, std::string const & filename)
00842 {
00843 std::string::size_type pos = filename.rfind(".")+1;
00844 std::string extension = filename.substr(pos, filename.size());
00845
00846 if(extension == "vtu")
00847 {
00848 process_vtu(filename);
00849 }
00850 else if(extension == "pvd")
00851 {
00852 process_pvd(filename);
00853 }
00854 else
00855 {
00856 std::cerr << "Error: vtk-reader does not support file extension: " << extension << std::endl;
00857 std::cerr << " the vtk-reader supports: " << std::endl;
00858 std::cerr << " *.vtu -- single-segment domains" << std::endl;
00859 std::cerr << " *.pvd -- multi-segment domains" << std::endl;
00860 return EXIT_FAILURE;
00861 }
00862
00863
00864
00865
00866 setupVertices(domain);
00867 domain.segments().resize(local_cell_num.size());
00868 for (size_t seg_id = 0; seg_id < local_cell_num.size(); ++seg_id)
00869 {
00870 setupCells(domain, seg_id);
00871 setupData(domain, seg_id);
00872 }
00873
00874 return EXIT_SUCCESS;
00875 }
00876
00878 std::vector<std::string> scalar_vertex_data_names(std::size_t segment_id) const
00879 {
00880 std::vector<std::string> ret;
00881 for (std::size_t i=0; i<local_scalar_vertex_data[segment_id].size(); ++i)
00882 ret.push_back(local_scalar_vertex_data[segment_id][i].first);
00883
00884 return ret;
00885 }
00886
00888 std::vector<std::string> vector_vertex_data_names(std::size_t segment_id) const
00889 {
00890 std::vector<std::string> ret;
00891 for (std::size_t i=0; i<local_vector_vertex_data[segment_id].size(); ++i)
00892 ret.push_back(local_vector_vertex_data[segment_id][i].first);
00893
00894 return ret;
00895 }
00896
00898 std::vector<std::string> scalar_cell_data_names(std::size_t segment_id) const
00899 {
00900 std::vector<std::string> ret;
00901 for (std::size_t i=0; i<local_scalar_cell_data[segment_id].size(); ++i)
00902 ret.push_back(local_scalar_cell_data[segment_id][i].first);
00903
00904 return ret;
00905 }
00906
00908 std::vector<std::string> vector_cell_data_names(std::size_t segment_id) const
00909 {
00910 std::vector<std::string> ret;
00911 for (std::size_t i=0; i<local_vector_cell_data[segment_id].size(); ++i)
00912 ret.push_back(local_vector_cell_data[segment_id][i].first);
00913
00914 return ret;
00915 }
00916
00918
00925 template <typename T>
00926 void add_scalar_data_on_vertices(T const & accessor, std::string name)
00927 {
00928 data_accessor_wrapper<VertexType> wrapper(accessor.clone());
00929 vertex_data_scalar.push_back(wrapper);
00930 vertex_data_scalar_names.push_back(name);
00931 }
00932
00939 template <typename T>
00940 void add_vector_data_on_vertices(T const & accessor, std::string name)
00941 {
00942 data_accessor_wrapper<VertexType> wrapper(accessor.clone());
00943 vertex_data_vector.push_back(wrapper);
00944 vertex_data_vector_names.push_back(name);
00945 }
00946
00953 template <typename T>
00954 void add_normal_data_on_vertices(T const & accessor, std::string name)
00955 {
00956 data_accessor_wrapper<VertexType> wrapper(accessor.clone());
00957 vertex_data_normal.push_back(wrapper);
00958 vertex_data_normal_names.push_back(name);
00959 }
00960
00961
00968 template <typename T>
00969 void add_scalar_data_on_cells(T const & accessor, std::string name)
00970 {
00971 data_accessor_wrapper<CellType> wrapper(accessor.clone());
00972 cell_data_scalar.push_back(wrapper);
00973 cell_data_scalar_names.push_back(name);
00974 }
00975
00982 template <typename T>
00983 void add_vector_data_on_cells(T const & accessor, std::string name)
00984 {
00985 data_accessor_wrapper<CellType> wrapper(accessor.clone());
00986 cell_data_vector.push_back(wrapper);
00987 cell_data_vector_names.push_back(name);
00988 }
00989
00990
00997 template <typename T>
00998 void add_normal_data_on_cells(T const & accessor, std::string name)
00999 {
01000 data_accessor_wrapper<CellType> wrapper(accessor.clone());
01001 cell_data_normal.push_back(wrapper);
01002 cell_data_normal_names.push_back(name);
01003 }
01004
01005
01006
01007
01009 std::vector<std::pair<std::size_t, std::string> > const & get_scalar_data_on_vertices() const
01010 {
01011 return vertex_data_scalar_read;
01012 }
01013
01015 std::vector<std::pair<std::size_t, std::string> > const & get_vector_data_on_vertices() const
01016 {
01017 return vertex_data_vector_read;
01018 }
01019
01021 std::vector<std::pair<std::size_t, std::string> > const & get_scalar_data_on_cells() const
01022 {
01023 return cell_data_scalar_read;
01024 }
01025
01027 std::vector<std::pair<std::size_t, std::string> > const & get_vector_data_on_cells() const
01028 {
01029 return cell_data_vector_read;
01030 }
01031
01032 private:
01033 std::vector< data_accessor_wrapper<VertexType> > vertex_data_scalar;
01034 std::vector< std::string > vertex_data_scalar_names;
01035
01036 std::vector< data_accessor_wrapper<VertexType> > vertex_data_vector;
01037 std::vector< std::string > vertex_data_vector_names;
01038
01039 std::vector< data_accessor_wrapper<VertexType> > vertex_data_normal;
01040 std::vector< std::string > vertex_data_normal_names;
01041
01042 std::vector< data_accessor_wrapper<CellType> > cell_data_scalar;
01043 std::vector< std::string > cell_data_scalar_names;
01044
01045 std::vector< data_accessor_wrapper<CellType> > cell_data_vector;
01046 std::vector< std::string > cell_data_vector_names;
01047
01048 std::vector< data_accessor_wrapper<CellType> > cell_data_normal;
01049 std::vector< std::string > cell_data_normal_names;
01050
01051
01052 std::vector<std::pair<std::size_t, std::string> > vertex_data_scalar_read;
01053 std::vector<std::pair<std::size_t, std::string> > vertex_data_vector_read;
01054
01055 std::vector<std::pair<std::size_t, std::string> > cell_data_scalar_read;
01056 std::vector<std::pair<std::size_t, std::string> > cell_data_vector_read;
01057 };
01058
01059
01061 template < typename DomainType >
01062 int import_vtk(DomainType & domain, std::string const & filename)
01063 {
01064 vtk_reader<DomainType> vtk_reader;
01065 return vtk_reader(domain, filename);
01066 }
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01086 template <typename KeyType, typename DataType, typename DomainType>
01087 vtk_reader<DomainType> & add_scalar_data_on_vertices(vtk_reader<DomainType> & reader,
01088 KeyType const & key,
01089 std::string quantity_name)
01090 {
01091 typedef typename DomainType::config_type DomainConfiguration;
01092 typedef typename result_of::ncell<DomainConfiguration, 0>::type VertexType;
01093
01094 data_accessor_wrapper<VertexType> wrapper(new global_scalar_data_accessor<VertexType, KeyType, DataType>(key));
01095 reader.add_scalar_data_on_vertices(wrapper, quantity_name);
01096 return reader;
01097 }
01098
01108 template <typename KeyType, typename DataType, typename DomainType>
01109 vtk_reader<DomainType> & add_scalar_data_on_vertices_per_segment(vtk_reader<DomainType> & reader,
01110 KeyType const & key,
01111 std::string quantity_name)
01112 {
01113 typedef typename DomainType::config_type DomainConfiguration;
01114 typedef typename result_of::ncell<DomainConfiguration, 0>::type VertexType;
01115
01116 data_accessor_wrapper<VertexType> wrapper(new segment_scalar_data_accessor<VertexType, KeyType, DataType>(key));
01117 reader.add_scalar_data_on_vertices(wrapper, quantity_name);
01118 return reader;
01119 }
01120
01121
01131 template <typename KeyType, typename DataType, typename DomainType>
01132 vtk_reader<DomainType> & add_vector_data_on_vertices(vtk_reader<DomainType> & reader,
01133 KeyType const & key,
01134 std::string quantity_name)
01135 {
01136 typedef typename DomainType::config_type DomainConfiguration;
01137 typedef typename result_of::ncell<DomainConfiguration, 0>::type VertexType;
01138
01139 data_accessor_wrapper<VertexType> wrapper(new global_vector_data_accessor<VertexType, KeyType, DataType>(key));
01140 reader.add_vector_data_on_vertices(wrapper, quantity_name);
01141 return reader;
01142 }
01143
01153 template <typename KeyType, typename DataType, typename DomainType>
01154 vtk_reader<DomainType> & add_vector_data_on_vertices_per_segment(vtk_reader<DomainType> & reader,
01155 KeyType const & key,
01156 std::string quantity_name)
01157 {
01158 typedef typename DomainType::config_type DomainConfiguration;
01159 typedef typename result_of::ncell<DomainConfiguration, 0>::type VertexType;
01160
01161 data_accessor_wrapper<VertexType> wrapper(new segment_vector_data_accessor<VertexType, KeyType, DataType>(key));
01162 reader.add_vector_data_on_vertices(wrapper, quantity_name);
01163 return reader;
01164 }
01165
01166
01167
01177 template <typename KeyType, typename DataType, typename DomainType>
01178 vtk_reader<DomainType> & add_normal_data_on_vertices(vtk_reader<DomainType> & reader,
01179 KeyType const & key,
01180 std::string quantity_name)
01181 {
01182 typedef typename DomainType::config_type DomainConfiguration;
01183 typedef typename result_of::ncell<DomainConfiguration, 0>::type VertexType;
01184
01185 data_accessor_wrapper<VertexType> wrapper(new global_vector_data_accessor<VertexType, KeyType, DataType>(key));
01186 reader.add_vector_data_on_vertices(wrapper, quantity_name);
01187 return reader;
01188 }
01189
01199 template <typename KeyType, typename DataType, typename DomainType>
01200 vtk_reader<DomainType> & add_normal_data_on_vertices_per_segment(vtk_reader<DomainType> & reader,
01201 KeyType const & key,
01202 std::string quantity_name)
01203 {
01204 typedef typename DomainType::config_type DomainConfiguration;
01205 typedef typename result_of::ncell<DomainConfiguration, 0>::type VertexType;
01206
01207 data_accessor_wrapper<VertexType> wrapper(new segment_vector_data_accessor<VertexType, KeyType, DataType>(key));
01208 reader.add_vector_data_on_vertices(wrapper, quantity_name);
01209 return reader;
01210 }
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01231 template <typename KeyType, typename DataType, typename DomainType>
01232 vtk_reader<DomainType> & add_scalar_data_on_cells(vtk_reader<DomainType> & reader,
01233 KeyType const & key,
01234 std::string quantity_name)
01235 {
01236 typedef typename DomainType::config_type DomainConfiguration;
01237 typedef typename DomainConfiguration::cell_tag CellTag;
01238 typedef typename result_of::ncell<DomainConfiguration, CellTag::dim>::type CellType;
01239
01240 data_accessor_wrapper<CellType> wrapper(new global_scalar_data_accessor<CellType, KeyType, DataType>(key));
01241 reader.add_scalar_data_on_cells(wrapper, quantity_name);
01242 return reader;
01243 }
01244
01254 template <typename KeyType, typename DataType, typename DomainType>
01255 vtk_reader<DomainType> & add_scalar_data_on_cells_per_segment(vtk_reader<DomainType> & reader,
01256 KeyType const & key,
01257 std::string quantity_name)
01258 {
01259 typedef typename DomainType::config_type DomainConfiguration;
01260 typedef typename DomainConfiguration::cell_tag CellTag;
01261 typedef typename result_of::ncell<DomainConfiguration, CellTag::dim>::type CellType;
01262
01263 data_accessor_wrapper<CellType> wrapper(new segment_scalar_data_accessor<CellType, KeyType, DataType>(key));
01264 reader.add_scalar_data_on_cells(wrapper, quantity_name);
01265 return reader;
01266 }
01267
01268
01278 template <typename KeyType, typename DataType, typename DomainType>
01279 vtk_reader<DomainType> & add_vector_data_on_cells(vtk_reader<DomainType> & reader,
01280 KeyType const & key,
01281 std::string quantity_name)
01282 {
01283 typedef typename DomainType::config_type DomainConfiguration;
01284 typedef typename DomainConfiguration::cell_tag CellTag;
01285 typedef typename result_of::ncell<DomainConfiguration, CellTag::dim>::type CellType;
01286
01287 data_accessor_wrapper<CellType> wrapper(new global_vector_data_accessor<CellType, KeyType, DataType>(key));
01288 reader.add_vector_data_on_cells(wrapper, quantity_name);
01289 return reader;
01290 }
01291
01301 template <typename KeyType, typename DataType, typename DomainType>
01302 vtk_reader<DomainType> & add_vector_data_on_cells_per_segment(vtk_reader<DomainType> & reader,
01303 KeyType const & key,
01304 std::string quantity_name)
01305 {
01306 typedef typename DomainType::config_type DomainConfiguration;
01307 typedef typename DomainConfiguration::cell_tag CellTag;
01308 typedef typename result_of::ncell<DomainConfiguration, CellTag::dim>::type CellType;
01309
01310 data_accessor_wrapper<CellType> wrapper(new segment_vector_data_accessor<CellType, KeyType, DataType>(key));
01311 reader.add_vector_data_on_cells(wrapper, quantity_name);
01312 return reader;
01313 }
01314
01315
01325 template <typename KeyType, typename DataType, typename DomainType>
01326 vtk_reader<DomainType> & add_normal_data_on_cells(vtk_reader<DomainType> & reader,
01327 KeyType const & key,
01328 std::string quantity_name)
01329 {
01330 typedef typename DomainType::config_type DomainConfiguration;
01331 typedef typename DomainConfiguration::cell_tag CellTag;
01332 typedef typename result_of::ncell<DomainConfiguration, CellTag::dim>::type CellType;
01333
01334 data_accessor_wrapper<CellType> wrapper(new global_vector_data_accessor<CellType, KeyType, DataType>(key));
01335 reader.add_vector_data_on_cells(wrapper, quantity_name);
01336 return reader;
01337 }
01338
01348 template <typename KeyType, typename DataType, typename DomainType>
01349 vtk_reader<DomainType> & add_normal_data_on_cells_per_segment(vtk_reader<DomainType> & reader,
01350 KeyType const & key,
01351 std::string quantity_name)
01352 {
01353 typedef typename DomainType::config_type DomainConfiguration;
01354 typedef typename DomainConfiguration::cell_tag CellTag;
01355 typedef typename result_of::ncell<DomainConfiguration, CellTag::dim>::type CellType;
01356
01357 data_accessor_wrapper<CellType> wrapper(new segment_vector_data_accessor<CellType, KeyType, DataType>(key));
01358 reader.add_vector_data_on_cells(wrapper, quantity_name);
01359 return reader;
01360 }
01361
01363 template <typename ReaderType>
01364 std::vector<std::pair<std::size_t, std::string> > const &
01365 get_scalar_data_on_vertices(ReaderType const & reader) { return reader.get_scalar_data_on_vertices(); }
01366
01368 template <typename ReaderType>
01369 std::vector<std::pair<std::size_t, std::string> > const &
01370 get_vector_data_on_vertices(ReaderType const & reader) { return reader.get_vector_data_on_vertices(); }
01371
01373 template <typename ReaderType>
01374 std::vector<std::pair<std::size_t, std::string> > const &
01375 get_scalar_data_on_cells(ReaderType const & reader) { return reader.get_scalar_data_on_cells(); }
01376
01378 template <typename ReaderType>
01379 std::vector<std::pair<std::size_t, std::string> > const &
01380 get_vector_data_on_cells(ReaderType const & reader) { return reader.get_vector_data_on_cells(); }
01381
01382
01383 }
01384 }
01385
01386 #endif