ViennaGrid - The Vienna Grid Library  2.1.0
viennagrid/io/vtk_reader.hpp
Go to the documentation of this file.
00001 #ifndef VIENNAGRID_IO_VTK_READER_HPP
00002 #define VIENNAGRID_IO_VTK_READER_HPP
00003 
00004 /* =======================================================================
00005    Copyright (c) 2011-2014, Institute for Microelectronics,
00006                             Institute for Analysis and Scientific Computing,
00007                             TU Wien.
00008 
00009                             -----------------
00010                      ViennaGrid - The Vienna Grid Library
00011                             -----------------
00012 
00013    License:      MIT (X11), see file LICENSE in the base directory
00014 ======================================================================= */
00015 
00016 
00021 #include <fstream>
00022 #include <iostream>
00023 #include <sstream>
00024 #include <vector>
00025 #include <deque>
00026 #include <string>
00027 #include <algorithm>
00028 
00029 #include "viennagrid/forwards.hpp"
00030 #include "viennagrid/point.hpp"
00031 #include "viennagrid/io/vtk_common.hpp"
00032 #include "viennagrid/io/helper.hpp"
00033 #include "viennagrid/io/xml_tag.hpp"
00034 #include "viennagrid/mesh/element_creation.hpp"
00035 
00036 namespace viennagrid
00037 {
00038   namespace io
00039   {
00040 
00046     template <typename MeshType, typename SegmentationType = typename viennagrid::result_of::segmentation<MeshType>::type >
00047     class vtk_reader
00048     {
00049     protected:
00050 
00051       typedef MeshType mesh_type;
00052       typedef SegmentationType segmentation_type;
00053 
00054       typedef typename SegmentationType::segment_handle_type SegmentHandleType;
00055       typedef typename SegmentationType::segment_id_type segment_id_type;
00056 
00057 
00058       typedef typename viennagrid::result_of::point<MeshType>::type PointType;
00059       typedef typename viennagrid::result_of::coord<PointType>::type CoordType;
00060       static const int geometric_dim = viennagrid::result_of::static_size<PointType>::value;
00061 
00062       typedef typename viennagrid::result_of::cell_tag<MeshType>::type CellTag;
00063 
00064 
00065       typedef typename result_of::vertex<MeshType>::type                          VertexType;
00066       typedef typename result_of::vertex_handle<MeshType>::type                           VertexHandleType;
00067       typedef typename result_of::vertex_id<MeshType>::type                           VertexIDType;
00068       typedef typename result_of::element<MeshType, CellTag>::type     CellType;
00069       typedef typename result_of::handle<MeshType, CellTag>::type     CellHandleType;
00070 
00071       typedef typename viennagrid::result_of::vertex_range<MeshType>::type   VertexRange;
00072       typedef typename viennagrid::result_of::iterator<VertexRange>::type                             VertexIterator;
00073 
00074       typedef typename viennagrid::result_of::line_range<MeshType>::type     EdgeRange;
00075       typedef typename viennagrid::result_of::iterator<EdgeRange>::type                               EdgeIterator;
00076 
00077       typedef typename viennagrid::result_of::facet_range<MeshType>::type   FacetRange;
00078       typedef typename viennagrid::result_of::iterator<FacetRange>::type                                   FacetIterator;
00079 
00080       typedef typename viennagrid::result_of::cell_range<MeshType>::type     CellRange;
00081       typedef typename viennagrid::result_of::iterator<CellRange>::type                  CellIterator;
00082 
00083 
00084       typedef std::vector<double> vector_data_type;
00085 
00086       typedef std::map< std::string, base_dynamic_field<double, VertexType> * >             VertexScalarOutputFieldContainer;
00087       typedef std::map< std::string, base_dynamic_field<vector_data_type, VertexType> * >   VertexVectorOutputFieldContainer;
00088 
00089       typedef std::map< std::string, base_dynamic_field<double, CellType> * >               CellScalarOutputFieldContainer;
00090       typedef std::map< std::string, base_dynamic_field<vector_data_type, CellType> * >     CellVectorOutputFieldContainer;
00091 
00092 
00093 
00094 
00095 
00096       std::ifstream                                        reader;
00097       std::map<PointType, std::size_t, point_less>         global_points;
00098       std::map<std::size_t, PointType>                     global_points_2;
00099       std::map<int, std::deque<std::size_t> >              local_to_global_map;
00100       std::map<int, std::deque<std::size_t> >              local_cell_vertices;
00101       std::map<int, std::deque<std::size_t> >              local_cell_offsets;
00102       std::map<int, std::size_t>                           local_cell_num;
00103       std::map<int, std::deque<CellHandleType> >           local_cell_handle;
00104 
00105       typedef viennagrid::element_key<CellType> CellElementKeyType;
00106       std::map<CellElementKeyType, CellHandleType>         global_cells;
00107 
00108       //data containers:
00109       std::map<int, std::deque<std::pair<std::string, std::deque<double> > > >  local_scalar_vertex_data;
00110       std::map<int, std::deque<std::pair<std::string, std::deque<double> > > >  local_vector_vertex_data;
00111       std::map<int, std::deque<std::pair<std::string, std::deque<double> > > >  local_scalar_cell_data;
00112       std::map<int, std::deque<std::pair<std::string, std::deque<double> > > >  local_vector_cell_data;
00113 
00114 
00115       template<typename map_type>
00116       void clear_map( map_type & map )
00117       {
00118         for (typename map_type::iterator it = map.begin(); it != map.end(); ++it)
00119           delete it->second;
00120 
00121         map.clear();
00122       }
00123 
00124       void post_clear()
00125       {
00126         clear_map(registered_vertex_scalar_data);
00127         clear_map(registered_vertex_vector_data);
00128 
00129         clear_map(registered_cell_scalar_data);
00130         clear_map(registered_cell_vector_data);
00131 
00132 
00133         for (typename std::map< segment_id_type, VertexScalarOutputFieldContainer>::iterator it = registered_segment_vertex_scalar_data.begin(); it != registered_segment_vertex_scalar_data.end(); ++it)
00134           clear_map(it->second);
00135 
00136         for (typename std::map< segment_id_type, VertexVectorOutputFieldContainer>::iterator it = registered_segment_vertex_vector_data.begin(); it != registered_segment_vertex_vector_data.end(); ++it)
00137           clear_map(it->second);
00138 
00139 
00140         for (typename std::map< segment_id_type, CellScalarOutputFieldContainer>::iterator it = registered_segment_cell_scalar_data.begin(); it != registered_segment_cell_scalar_data.end(); ++it)
00141           clear_map(it->second);
00142 
00143         for (typename std::map< segment_id_type, CellVectorOutputFieldContainer>::iterator it = registered_segment_cell_vector_data.begin(); it != registered_segment_cell_vector_data.end(); ++it)
00144           clear_map(it->second);
00145 
00146 
00147         registered_segment_vertex_scalar_data.clear();
00148         registered_segment_vertex_vector_data.clear();
00149 
00150         registered_segment_cell_scalar_data.clear();
00151         registered_segment_cell_vector_data.clear();
00152       }
00153 
00154 
00155       void pre_clear()
00156       {
00157         vertex_data_scalar_read.clear();
00158         vertex_data_vector_read.clear();
00159 
00160         cell_data_scalar_read.clear();
00161         cell_data_vector_read.clear();
00162 
00163         vertex_scalar_data.clear();
00164         vertex_vector_data.clear();
00165 
00166         cell_scalar_data.clear();
00167         cell_vector_data.clear();
00168 
00169         global_points.clear();
00170         global_points_2.clear();
00171         local_to_global_map.clear();
00172         local_cell_vertices.clear();
00173         local_cell_offsets.clear();
00174         local_cell_num.clear();
00175         local_cell_handle.clear();
00176 
00177         global_cells.clear();
00178 
00179         local_scalar_vertex_data.clear();
00180         local_vector_vertex_data.clear();
00181         local_scalar_cell_data.clear();
00182         local_vector_cell_data.clear();
00183       }
00184 
00185 
00186 
00187 
00188 
00190       void openFile(std::string const & filename)
00191       {
00192         reader.open(filename.c_str());
00193         if(!reader)
00194         {
00195           throw cannot_open_file_exception("* ViennaGrid: vtk_reader::openFile(): File " + filename + ": Cannot open file!");
00196         }
00197       }
00198 
00200       void closeFile()
00201       {
00202         reader.close();
00203       }
00204 
00206       bool lowercase_compare(std::string const & s1, std::string const & s2)
00207       {
00208         std::string s1_lower = s1;
00209         std::transform(s1.begin(), s1.end(), s1_lower.begin(), char_to_lower<>());
00210 
00211         std::string s2_lower = s2;
00212         std::transform(s2.begin(), s2.end(), s2_lower.begin(), char_to_lower<>());
00213 
00214         return s1_lower == s2_lower;
00215       }
00216 
00218       void checkNextToken(std::string const & expectedToken)
00219       {
00220         std::string token;
00221         reader >> token;
00222 
00223         if ( !lowercase_compare(token, expectedToken) )
00224         {
00225           std::cerr << "* vtk_reader::operator(): Expected \"" << expectedToken << "\", but got \"" << token << "\"" << std::endl;
00226           std::stringstream ss;
00227           ss << "* ViennaGrid: vtk_reader::operator(): Expected \"" << expectedToken << "\", but got \"" << token << "\"";
00228           throw bad_file_format_exception(ss.str());
00229         }
00230       }
00231 
00233       void readNodeCoordinates(std::size_t nodeNum, std::size_t numberOfComponents, segment_id_type seg_id)
00234       {
00235         double nodeCoord;
00236         local_to_global_map[seg_id].resize(nodeNum);
00237 
00238         for(std::size_t i = 0; i < nodeNum; i++)
00239         {
00240           PointType p;
00241 
00242           for(std::size_t j = 0; j < numberOfComponents; j++)
00243           {
00244             reader >> nodeCoord;
00245             if (j < static_cast<std::size_t>(geometric_dim))
00246               p[j] = nodeCoord;
00247           }
00248 
00249           //add point to global list if not already there
00250           if (global_points.find(p) == global_points.end())
00251           {
00252             std::size_t new_global_id = global_points.size();
00253             global_points.insert( std::make_pair(p, new_global_id) );
00254             global_points_2.insert( std::make_pair(new_global_id, p) );
00255             local_to_global_map[seg_id][i] = new_global_id;
00256           }
00257           else
00258           {
00259             local_to_global_map[seg_id][i] = global_points[p];
00260           }
00261         }
00262       }
00263 
00265       void readCellIndices(segment_id_type seg_id)
00266       {
00267         std::size_t cellNode = 0;
00268         std::string token;
00269         reader >> token;
00270 
00271         while (token != "</DataArray>")
00272         {
00273           assert( strChecker::myIsNumber(token) && "Cell vertex index is not a number!" );
00274 
00275           cellNode = static_cast<std::size_t>(atoi(token.c_str()));
00276           local_cell_vertices[seg_id].push_back(cellNode);
00277 
00278           reader >> token;
00279         }
00280       }
00281 
00283       void readOffsets(segment_id_type seg_id)
00284       {
00285           //****************************************************************************
00286           // read in the offsets: describe the affiliation of the nodes to the cells
00287           // (see: http://www.vtk.org/pdf/file-formats.pdf , page 9)
00288           //****************************************************************************
00289 
00290           std::string token;
00291           std::size_t offset = 0;
00292           reader >> token;
00293 
00294           while(token != "</DataArray>")
00295           {
00296             assert( strChecker::myIsNumber(token) && "Cell offset is not a number!" );
00297 
00298             offset = static_cast<std::size_t>(atoi(token.c_str()));
00299             local_cell_offsets[seg_id].push_back(offset);
00300 
00301             //std::cout << "Vertex#: " << offset << std::endl;
00302             reader >> token;
00303           }
00304       }
00305 
00307       void readTypes()
00308       {
00309           //****************************************************************************
00310           // read in the offsets: describe the affiliation of the nodes to the cells
00311           // (see: http://www.vtk.org/pdf/file-formats.pdf , page 9)
00312           //****************************************************************************
00313 
00314           std::string token;
00315 #ifndef NDEBUG
00316           long type = 0;
00317 #endif
00318           reader >> token;
00319 
00320           while(token != "</DataArray>")
00321           {
00322             assert( strChecker::myIsNumber(token) && "Cell type is not a number!" );
00323 #ifndef NDEBUG
00324             type = atoi(token.c_str());
00325             assert(type == detail::ELEMENT_TAG_TO_VTK_TYPE<CellTag>::value && "Error in VTK reader: Type mismatch!");
00326 #endif
00327             //std::cout << "Vertex#: " << offset << std::endl;
00328             reader >> token;
00329           }
00330       }
00331 
00333       void readData(std::deque<double> & container)
00334       {
00335         std::string token;
00336         reader >> token;
00337 
00338         while (string_to_lower(token) != "</dataarray>")
00339         {
00340           container.push_back( atof(token.c_str()) );
00341           reader >> token;
00342         }
00343       }
00344 
00346       template <typename ContainerType, typename NameContainerType>
00347       void readPointCellData(segment_id_type seg_id,
00348                              ContainerType & scalar_data,
00349                              ContainerType & vector_data,
00350                              NameContainerType & data_names_scalar,
00351                              NameContainerType & data_names_vector)
00352       {
00353         std::string name;
00354         std::size_t components = 1;
00355 
00356         xml_tag<> tag;
00357 
00358         tag.parse(reader);
00359 
00360         while (tag.name() == "dataarray")
00361         {
00362           tag.check_attribute("name", "");
00363           name = tag.get_value("name");
00364 
00365           if (tag.has_attribute("numberofcomponents"))
00366             components = static_cast<std::size_t>(atoi(tag.get_value("numberofcomponents").c_str()));
00367 
00368           //now read data:
00369           if (components == 1)
00370           {
00371             data_names_scalar.push_back(std::make_pair(seg_id, name));
00372             scalar_data[seg_id].push_back( std::make_pair(name, std::deque<double>()) );
00373             readData(scalar_data[seg_id].back().second);
00374           }
00375           else if (components == 3)
00376           {
00377             data_names_vector.push_back(std::make_pair(seg_id, name));
00378             vector_data[seg_id].push_back( std::make_pair(name, std::deque<double>()) );
00379             readData(vector_data[seg_id].back().second);
00380           }
00381           else
00382             throw bad_file_format_exception("* ViennaGrid: vtk_reader::readPointCellData(): Number of components for data invalid!");
00383 
00384           tag.parse(reader);
00385         }
00386 
00387 
00388         if (tag.name() != "/pointdata" && tag.name() != "/celldata")
00389             throw bad_file_format_exception("* ViennaGrid: vtk_reader::readPointCellData(): XML Parse error: Expected </PointData> or </CellData>!");
00390 
00391       }
00392 
00394 
00396       void setupVertices(MeshType & mesh_obj)
00397       {
00398         for (std::size_t i=0; i<global_points_2.size(); ++i)
00399           viennagrid::make_vertex_with_id( mesh_obj, typename VertexType::id_type(typename VertexType::id_type::base_id_type(i)), global_points_2[i] );
00400       }
00401 
00403       void setupCells(MeshType & mesh_obj, SegmentationType & segmentation, segment_id_type seg_id)
00404       {
00405         //***************************************************
00406         // building up the cells in ViennaGrid
00407         // -------------------------------------------------
00408         // "cells"   ... contain the indices of the nodes
00409         // "offsets" ... contain the information about
00410         //               the affiliation of the nodes
00411         //               to the cells
00412         //***************************************************
00413         //CellType cell;
00414         std::size_t numVertices = 0;
00415         std::size_t offsetIdx = 0;
00416 
00417         std::deque<std::size_t> const & offsets = local_cell_offsets[seg_id];
00418 
00419         for (std::size_t i = 0; i < local_cell_num[seg_id]; i++)
00420         {
00421           //*************************************************
00422           // choose the offset index for the i-th cell
00423           // in the "cells"-vector and calculate the
00424           // number of vertices belonging to the i-th cell
00425           //*************************************************
00426           if(i == 0)
00427           {
00428             offsetIdx = 0;
00429             numVertices = offsets[i];
00430           }
00431           else
00432           {
00433             offsetIdx = offsets[i-1];
00434             numVertices = offsets[i]-offsets[i-1];
00435           }
00436 
00437 
00438           //****************************************************
00439           // read out the node indices form the "cells"-vector
00440           // and add the cells to the "vertices"-array
00441           //****************************************************
00442 
00443           viennagrid::static_array<VertexHandleType, boundary_elements<CellTag, vertex_tag>::num> cell_vertex_handles;
00444           std::vector<VertexIDType> cell_vertex_ids(numVertices);
00445 
00446           detail::vtk_to_viennagrid_orientations<CellTag> reorderer;
00447           for (std::size_t j = 0; j < numVertices; j++)
00448           {
00449             std::size_t reordered_j = reorderer(j);
00450             std::size_t local_index = local_cell_vertices[seg_id][reordered_j + offsetIdx];
00451             std::size_t global_vertex_index = local_to_global_map[seg_id][local_index];
00452 
00453             cell_vertex_handles[j] = viennagrid::elements<viennagrid::vertex_tag>(mesh_obj).handle_at(global_vertex_index);
00454             viennagrid::add( segmentation[seg_id], viennagrid::dereference_handle(segmentation, cell_vertex_handles[j]) );
00455 
00456             cell_vertex_ids[j] = viennagrid::dereference_handle(mesh_obj, cell_vertex_handles[j]).id();
00457           }
00458 
00459 
00460           CellElementKeyType cell_key(cell_vertex_ids);
00461           typename std::map<CellElementKeyType, CellHandleType>::iterator chit = global_cells.find( cell_key );
00462           if (chit != global_cells.end())
00463           {
00464             local_cell_handle[seg_id].push_back( chit->second );
00465             viennagrid::add( segmentation[seg_id], viennagrid::dereference_handle(mesh_obj, chit->second) );
00466           }
00467           else
00468           {
00469             CellHandleType cell_handle = viennagrid::make_element<CellType>(segmentation[seg_id], cell_vertex_handles.begin(), cell_vertex_handles.end());
00470             global_cells[cell_key] = cell_handle;
00471 
00472             local_cell_handle[seg_id].push_back(cell_handle);
00473           }
00474         }
00475       }
00476 
00478       template <typename ContainerType>
00479       void setupDataVertex(MeshType & mesh_obj, SegmentHandleType &
00480 #ifndef NDEBUG
00481         segment
00482 #endif
00483         , segment_id_type seg_id, ContainerType const & container, std::size_t num_components)
00484       {
00485         std::string const & name = container.first;
00486 
00487         if (num_components == 1)
00488         {
00489           VertexScalarOutputFieldContainer & current_registered_segment_vertex_scalar_data = registered_segment_vertex_scalar_data[seg_id];
00490           if (registered_vertex_scalar_data.find(name) != registered_vertex_scalar_data.end())
00491           {
00492             for (std::size_t i=0; i<container.second.size(); ++i)
00493             {
00494               std::size_t global_vertex_id = local_to_global_map[seg_id][i];
00495               VertexType const & vertex = viennagrid::elements<viennagrid::vertex_tag>(mesh_obj)[global_vertex_id];
00496 
00497               (*registered_vertex_scalar_data[name])(vertex) = (container.second)[i];
00498             }
00499           }
00500           else if (current_registered_segment_vertex_scalar_data.find(name) != current_registered_segment_vertex_scalar_data.end())
00501           {
00502             for (std::size_t i=0; i<container.second.size(); ++i)
00503             {
00504               std::size_t global_vertex_id = local_to_global_map[seg_id][i];
00505               VertexType const & vertex = viennagrid::elements<viennagrid::vertex_tag>(mesh_obj)[global_vertex_id];
00506 
00507               (*current_registered_segment_vertex_scalar_data[name])(vertex) = (container.second)[i];
00508             }
00509           }
00510           else
00511           {
00512             #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00513             std::cout << "* vtk_reader::operator(): Reading scalar quantity "
00514                       << container.first << " to vertices." << std::endl;
00515             #endif
00516             for (std::size_t i=0; i<container.second.size(); ++i)
00517             {
00518               std::size_t global_vertex_id = local_to_global_map[seg_id][i];
00519               VertexType const & vertex = viennagrid::elements<viennagrid::vertex_tag>(mesh_obj)[global_vertex_id];
00520 
00521               if ( static_cast<typename VertexType::id_type::base_id_type>(vertex_scalar_data[container.first][seg_id].size()) <= vertex.id().get())
00522                 vertex_scalar_data[container.first][seg_id].resize(static_cast<std::size_t>(vertex.id().get()+1));
00523               vertex_scalar_data[container.first][seg_id][static_cast<std::size_t>(vertex.id().get())] = (container.second)[i];
00524             }
00525           }
00526         }
00527         else
00528         {
00529           VertexVectorOutputFieldContainer & current_registered_segment_vertex_vector_data = registered_segment_vertex_vector_data[seg_id];
00530           if (registered_vertex_vector_data.find(name) != registered_vertex_vector_data.end())
00531           {
00532             for (std::size_t i=0; i<container.second.size()/3; ++i)
00533             {
00534               std::size_t global_vertex_id = local_to_global_map[seg_id][i];
00535               VertexType const & vertex = viennagrid::elements<viennagrid::vertex_tag>(mesh_obj)[global_vertex_id];
00536 
00537               (*registered_vertex_vector_data[name])(vertex).resize(3);
00538               (*registered_vertex_vector_data[name])(vertex)[0] = (container.second)[3*i+0];
00539               (*registered_vertex_vector_data[name])(vertex)[1] = (container.second)[3*i+1];
00540               (*registered_vertex_vector_data[name])(vertex)[2] = (container.second)[3*i+2];
00541             }
00542           }
00543           else if (current_registered_segment_vertex_vector_data.find(name) != current_registered_segment_vertex_vector_data.end())
00544           {
00545             for (std::size_t i=0; i<container.second.size(); ++i)
00546             {
00547               std::size_t global_vertex_id = local_to_global_map[seg_id][i];
00548               VertexType const & vertex = viennagrid::elements<viennagrid::vertex_tag>(mesh_obj)[global_vertex_id];
00549 
00550               (*current_registered_segment_vertex_vector_data[name])(vertex).resize(3);
00551               (*current_registered_segment_vertex_vector_data[name])(vertex)[0] = (container.second)[3*i+0];
00552               (*current_registered_segment_vertex_vector_data[name])(vertex)[1] = (container.second)[3*i+1];
00553               (*current_registered_segment_vertex_vector_data[name])(vertex)[2] = (container.second)[3*i+2];
00554             }
00555           }
00556           else
00557           {
00558             #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00559             std::cout << "* vtk_reader::operator(): Reading vector quantity "
00560                       << container.first << " to vertices." << std::endl;
00561             #endif
00562             assert( 3 * viennagrid::elements<viennagrid::vertex_tag>(segment).size() == container.second.size());
00563             for (std::size_t i=0; i<container.second.size() / 3; ++i)
00564             {
00565               std::size_t global_vertex_id = local_to_global_map[seg_id][i];
00566               VertexType const & vertex = viennagrid::elements<viennagrid::vertex_tag>(mesh_obj)[global_vertex_id];
00567 
00568               if ( static_cast<typename VertexType::id_type::base_id_type>(vertex_vector_data[container.first][seg_id].size()) <= vertex.id().get())
00569                 vertex_vector_data[container.first][seg_id].resize(static_cast<std::size_t>(vertex.id().get()+1));
00570               vertex_vector_data[container.first][seg_id][static_cast<std::size_t>(vertex.id().get())].resize(3);
00571               vertex_vector_data[container.first][seg_id][static_cast<std::size_t>(vertex.id().get())][0] = (container.second)[3*i+0];
00572               vertex_vector_data[container.first][seg_id][static_cast<std::size_t>(vertex.id().get())][1] = (container.second)[3*i+1];
00573               vertex_vector_data[container.first][seg_id][static_cast<std::size_t>(vertex.id().get())][2] = (container.second)[3*i+2];
00574             }
00575           }
00576         }
00577       }
00578 
00580       template <typename ContainerType>
00581       void setupDataCell(MeshType &, SegmentHandleType & segment, segment_id_type seg_id, ContainerType const & container, std::size_t num_components)
00582       {
00583         std::string const & name = container.first;
00584 
00585         if (num_components == 1)
00586         {
00587           CellScalarOutputFieldContainer & current_registered_segment_cell_scalar_data = registered_segment_cell_scalar_data[seg_id];
00588           if (registered_cell_scalar_data.find(name) != registered_cell_scalar_data.end())
00589           {
00590             for (std::size_t i=0; i<container.second.size(); ++i)
00591             {
00592               CellType const & cell = viennagrid::dereference_handle( segment, local_cell_handle[seg_id][i] );
00593 
00594               (*registered_cell_scalar_data[name])(cell) = (container.second)[i];
00595             }
00596           }
00597           else if (current_registered_segment_cell_scalar_data.find(name) != current_registered_segment_cell_scalar_data.end())
00598           {
00599             for (std::size_t i=0; i<container.second.size(); ++i)
00600             {
00601               CellType const & cell = viennagrid::dereference_handle( segment, local_cell_handle[seg_id][i] );
00602 
00603               (*current_registered_segment_cell_scalar_data[name])(cell) = (container.second)[i];
00604             }
00605           }
00606           else
00607           {
00608             #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00609             std::cout << "* vtk_reader::operator(): Reading scalar quantity "
00610                       << container.first << " to vertices." << std::endl;
00611             #endif
00612             for (std::size_t i=0; i<container.second.size(); ++i)
00613             {
00614               CellType const & cell = viennagrid::dereference_handle( segment, local_cell_handle[seg_id][i] );
00615 
00616               if ( static_cast<typename CellType::id_type::base_id_type>(cell_scalar_data[container.first][seg_id].size()) <= cell.id().get())
00617                 cell_scalar_data[container.first][seg_id].resize(static_cast<std::size_t>(cell.id().get()+1));
00618               cell_scalar_data[container.first][seg_id][static_cast<std::size_t>(cell.id().get())] = (container.second)[i];
00619             }
00620           }
00621         }
00622         else
00623         {
00624           CellVectorOutputFieldContainer & current_registered_segment_cell_vector_data = registered_segment_cell_vector_data[seg_id];
00625           if (registered_cell_vector_data.find(name) != registered_cell_vector_data.end())
00626           {
00627             for (std::size_t i=0; i<container.second.size()/3; ++i)
00628             {
00629               CellType const & cell = viennagrid::dereference_handle( segment, local_cell_handle[seg_id][i] );
00630 
00631               (*registered_cell_vector_data[name])(cell).resize(3);
00632               (*registered_cell_vector_data[name])(cell)[0] = (container.second)[3*i+0];
00633               (*registered_cell_vector_data[name])(cell)[1] = (container.second)[3*i+1];
00634               (*registered_cell_vector_data[name])(cell)[2] = (container.second)[3*i+2];
00635             }
00636           }
00637           else if (current_registered_segment_cell_vector_data.find(name) != current_registered_segment_cell_vector_data.end())
00638           {
00639             for (std::size_t i=0; i<container.second.size(); ++i)
00640             {
00641               CellType const & cell = viennagrid::dereference_handle( segment, local_cell_handle[seg_id][i] );
00642 
00643               (*current_registered_segment_cell_vector_data[name])(cell).resize(3);
00644               (*current_registered_segment_cell_vector_data[name])(cell)[0] = (container.second)[3*i+0];
00645               (*current_registered_segment_cell_vector_data[name])(cell)[1] = (container.second)[3*i+1];
00646               (*current_registered_segment_cell_vector_data[name])(cell)[2] = (container.second)[3*i+2];
00647             }
00648           }
00649           else
00650           {
00651             #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00652             std::cout << "* vtk_reader::operator(): Reading vector quantity "
00653                       << container.first << " to vertices." << std::endl;
00654             #endif
00655             assert( 3 * viennagrid::elements<CellTag>(segment).size() == container.second.size());
00656             for (std::size_t i=0; i<container.second.size() / 3; ++i)
00657             {
00658               CellType const & cell = viennagrid::dereference_handle( segment, local_cell_handle[seg_id][i] );
00659 
00660               if ( static_cast<typename CellType::id_type::base_id_type>(cell_vector_data[container.first][seg_id].size()) <= cell.id().get())
00661                 cell_vector_data[container.first][seg_id].resize(static_cast<std::size_t>(cell.id().get()+1));
00662               cell_vector_data[container.first][seg_id][static_cast<std::size_t>(cell.id().get())].resize(3);
00663               cell_vector_data[container.first][seg_id][static_cast<std::size_t>(cell.id().get())][0] = (container.second)[3*i+0];
00664               cell_vector_data[container.first][seg_id][static_cast<std::size_t>(cell.id().get())][1] = (container.second)[3*i+1];
00665               cell_vector_data[container.first][seg_id][static_cast<std::size_t>(cell.id().get())][2] = (container.second)[3*i+2];
00666             }
00667           }
00668         }
00669       }
00670 
00672       void setupData(MeshType & mesh_obj, SegmentationType & segmentation, segment_id_type seg_id)
00673       {
00674         for (size_t i=0; i<local_scalar_vertex_data[seg_id].size(); ++i)
00675         {
00676           setupDataVertex(mesh_obj, segmentation[seg_id], seg_id, local_scalar_vertex_data[seg_id][i], 1);
00677         }
00678 
00679         for (size_t i=0; i<local_vector_vertex_data[seg_id].size(); ++i)
00680         {
00681           setupDataVertex(mesh_obj, segmentation[seg_id], seg_id, local_vector_vertex_data[seg_id][i], 3);
00682         }
00683 
00684 
00685         for (size_t i=0; i<local_scalar_cell_data[seg_id].size(); ++i)
00686         {
00687           setupDataCell(mesh_obj, segmentation[seg_id], seg_id, local_scalar_cell_data[seg_id][i], 1);
00688         }
00689 
00690         for (size_t i=0; i<local_vector_cell_data[seg_id].size(); ++i)
00691         {
00692           setupDataCell(mesh_obj, segmentation[seg_id], seg_id, local_vector_cell_data[seg_id][i], 3);
00693         }
00694 
00695       }
00696 
00698       void parse_vtu_segment(std::string filename, segment_id_type seg_id)
00699       {
00700 
00701         try
00702         {
00703           openFile(filename);
00704 
00705           std::size_t nodeNum = 0;
00706           std::size_t numberOfComponents = 0;
00707 
00708           xml_tag<> tag;
00709 
00710           tag.parse(reader);
00711           if (tag.name() != "?xml" && tag.name() != "?xml?")
00712             throw bad_file_format_exception("* ViennaGrid: vtk_reader::parse_vtu_segment(): Parse error: No opening ?xml tag!");
00713 
00714           tag.parse_and_check_name(reader, "vtkfile", filename);
00715           tag.parse_and_check_name(reader, "unstructuredgrid", filename);
00716 
00717           tag.parse_and_check_name(reader, "piece", filename);
00718 
00719           tag.check_attribute("numberofpoints", filename);
00720 
00721           nodeNum = static_cast<std::size_t>(atoi(tag.get_value("numberofpoints").c_str()));
00722           #ifdef VIENNAGRID_DEBUG_IO
00723           std::cout << "#Nodes: " << nodeNum << std::endl;
00724           #endif
00725 
00726           tag.check_attribute("numberofcells", filename);
00727 
00728           local_cell_num[seg_id] = static_cast<std::size_t>(atoi(tag.get_value("numberofcells").c_str()));
00729           #ifdef VIENNAGRID_DEBUG_IO
00730           std::cout << "#Cells: " << local_cell_num[seg_id] << std::endl;
00731           #endif
00732 
00733           tag.parse_and_check_name(reader, "points", filename);
00734 
00735           tag.parse_and_check_name(reader, "dataarray", filename);
00736           tag.check_attribute("numberofcomponents", filename);
00737 
00738           numberOfComponents = static_cast<std::size_t>(atoi(tag.get_value("numberofcomponents").c_str()));
00739           readNodeCoordinates(nodeNum, numberOfComponents, seg_id);
00740 
00741           tag.parse_and_check_name(reader, "/dataarray", filename);
00742           tag.parse_and_check_name(reader, "/points", filename);
00743 
00744           tag.parse(reader);
00745           if (tag.name() == "pointdata")
00746           {
00747              readPointCellData(seg_id, local_scalar_vertex_data, local_vector_vertex_data,
00748                                       vertex_data_scalar_read, vertex_data_vector_read);
00749             tag.parse(reader);
00750           }
00751 
00752           if (tag.name() != "cells")
00753             throw bad_file_format_exception("* ViennaGrid: vtk_reader::parse_vtu_segment(): Parse error: Expected Cells tag!");
00754 
00755           for (std::size_t i=0; i<3; ++i)
00756           {
00757             tag.parse_and_check_name(reader, "dataarray", filename);
00758             tag.check_attribute("name", filename);
00759 
00760             if (tag.get_value("name") == "connectivity")
00761               readCellIndices(seg_id);
00762             else if (tag.get_value("name") == "offsets")
00763               readOffsets(seg_id);
00764             else if (tag.get_value("name") == "types")
00765               readTypes();
00766             else
00767               throw bad_file_format_exception("* ViennaGrid: vtk_reader::parse_vtu_segment(): Parse error: <DataArray> is not named 'connectivity', 'offsets' or 'types'!");
00768           }
00769 
00770           tag.parse_and_check_name(reader, "/cells", filename);
00771 
00772           tag.parse(reader);
00773           if (tag.name() == "celldata")
00774           {
00775             readPointCellData(seg_id, local_scalar_cell_data, local_vector_cell_data,
00776                                       cell_data_scalar_read, cell_data_vector_read);
00777             tag.parse(reader);
00778           }
00779 
00780 
00781           if (tag.name() != "/piece")
00782             throw bad_file_format_exception("* ViennaGrid: vtk_reader::parse_vtu_segment(): Parse error: Expected </Piece> tag!");
00783 
00784           tag.parse_and_check_name(reader, "/unstructuredgrid", filename);
00785           tag.parse_and_check_name(reader, "/vtkfile", filename);
00786 
00787           closeFile();
00788         }
00789         catch (std::exception const & ex) {
00790           std::cerr << "Problems while reading file " << filename << std::endl;
00791           std::cerr << "what(): " << ex.what() << std::endl;
00792         }
00793 
00794       }
00795 
00797       void process_vtu(std::string const & filename)
00798       {
00799         parse_vtu_segment(filename, 0);
00800       }
00801 
00803       void process_pvd(std::string const & filename)
00804       {
00805         std::map<int, std::string> filenames;
00806 
00807         //extract path from .pvd file:
00808         std::string::size_type pos = filename.rfind("/");
00809         std::string path_to_pvd;
00810         if (pos == std::string::npos)
00811           pos = filename.rfind("\\"); //a tribute to Windows... ;-)
00812 
00813         if (pos != std::string::npos)
00814           path_to_pvd = filename.substr(0, pos + 1);
00815 
00816         openFile(filename);
00817 
00818         //
00819         // Step 1: Get segments from pvd file:
00820         //
00821         xml_tag<> tag;
00822 
00823         tag.parse(reader);
00824         if (tag.name() != "?xml" && tag.name() != "?xml?")
00825           throw bad_file_format_exception("* ViennaGrid: vtk_reader::process_pvd(): Parse error: No opening <?xml?> tag!");
00826 
00827         tag.parse(reader);
00828         if (tag.name() != "vtkfile")
00829           throw bad_file_format_exception("* ViennaGrid: vtk_reader::process_pvd(): Parse error: VTKFile tag expected!");
00830 
00831         if (!tag.has_attribute("type"))
00832           throw bad_file_format_exception("* ViennaGrid: vtk_reader::process_pvd(): Parse error: VTKFile tag has no attribute 'type'!");
00833 
00834         if (string_to_lower(tag.get_value("type")) != "collection")
00835           throw bad_file_format_exception("* ViennaGrid: vtk_reader::process_pvd(): Parse error: Type-attribute of VTKFile tag is not 'Collection'!");
00836 
00837         tag.parse(reader);
00838         if (tag.name() != "collection")
00839           throw bad_file_format_exception("* ViennaGrid: vtk_reader::process_pvd(): Parse error: Collection tag expected!");
00840 
00841         while (reader.good())
00842         {
00843           tag.parse(reader);
00844 
00845           if (tag.name() == "/collection")
00846             break;
00847 
00848           if (tag.name() != "dataset")
00849             throw bad_file_format_exception("* ViennaGrid: vtk_reader::process_pvd(): Parse error: DataSet tag expected!");
00850 
00851           if (!tag.has_attribute("file"))
00852             throw bad_file_format_exception("* ViennaGrid: vtk_reader::process_pvd(): Parse error: DataSet tag has no file attribute!");
00853 
00854           filenames[ atoi(tag.get_value("part").c_str()) ] = tag.get_value("file");
00855 
00856         }
00857 
00858         tag.parse(reader);
00859         if (tag.name() != "/vtkfile")
00860           throw bad_file_format_exception("* ViennaGrid: vtk_reader::process_pvd(): Parse error: Closing VTKFile tag expected!");
00861 
00862         closeFile();
00863 
00864         assert(filenames.size() > 0 && "No segments in pvd-file specified!");
00865 
00866         //
00867         // Step 2: Parse .vtu files:
00868         //
00869         for (std::map<int, std::string>::iterator it = filenames.begin(); it != filenames.end(); ++it)
00870         {
00871           #if defined VIENNAGRID_DEBUG_ALL || defined VIENNAGRID_DEBUG_IO
00872           std::cout << "Parsing file " << path_to_pvd + it->second << std::endl;
00873           #endif
00874           parse_vtu_segment(path_to_pvd + it->second, it->first);
00875         }
00876 
00877       }
00878 
00879 
00880     public:
00881 
00882 
00883       ~vtk_reader() { pre_clear(); post_clear(); }
00884 
00885 
00892       void operator()(MeshType & mesh_obj, SegmentationType & segmentation, std::string const & filename)
00893       {
00894         pre_clear();
00895 
00896         std::string::size_type pos  = filename.rfind(".")+1;
00897         std::string extension = filename.substr(pos, filename.size());
00898 
00899         if(extension == "vtu")
00900         {
00901           process_vtu(filename);
00902         }
00903         else if(extension == "pvd")
00904         {
00905           process_pvd(filename);
00906         }
00907         else
00908         {
00909           std::stringstream ss;
00910           ss << "* ViennaGrid: " << filename << " - ";
00911           ss << "Error: vtk-reader does not support file extension: " << extension << "\n";
00912           ss << "       the vtk-reader supports: \n";
00913           ss << "       *.vtu -- single-segment meshs\n";
00914           ss << "       *.pvd -- multi-segment meshs\n";
00915 
00916           throw bad_file_format_exception(ss.str());
00917         }
00918 
00919         //
00920         // push everything to the ViennaGrid mesh:
00921         //
00922         setupVertices(mesh_obj);
00923 //         for (size_t seg_id = 0; seg_id < local_cell_num.size(); ++seg_id)
00924         for (std::map<int, std::size_t>::iterator it = local_cell_num.begin(); it != local_cell_num.end(); ++it)
00925         {
00926           segmentation.get_make_segment( it->first );
00927         }
00928 
00929 //         for (size_t seg_id = 0; seg_id < local_cell_num.size(); ++seg_id)
00930         for (std::map<int, std::size_t>::iterator it = local_cell_num.begin(); it != local_cell_num.end(); ++it)
00931         {
00932           setupCells(mesh_obj, segmentation, it->first);
00933           setupData(mesh_obj, segmentation, it->first);
00934         }
00935 
00936         post_clear();
00937       } //operator()
00938 
00939 
00945       void operator()(MeshType & mesh_obj, std::string const & filename)
00946       {
00947         SegmentationType tmp(mesh_obj);
00948         (*this)(mesh_obj, tmp, filename);
00949       }
00950 
00951 
00952 
00953 
00955       std::vector<std::string> scalar_vertex_data_names(segment_id_type segment_id) const
00956       {
00957         std::vector<std::string> ret;
00958 
00959         std::map<int, std::deque<std::pair<std::string, std::deque<double> > > >::const_iterator it = local_scalar_vertex_data.find(segment_id);
00960         if (it == local_scalar_vertex_data.end())
00961           return ret;
00962 
00963         for (std::size_t i=0; i<it->second.size(); ++i)
00964           ret.push_back(it->second[i].first);
00965 
00966         return ret;
00967       }
00968 
00970       std::vector<std::string> vector_vertex_data_names(segment_id_type segment_id) const
00971       {
00972         std::vector<std::string> ret;
00973 
00974         std::map<int, std::deque<std::pair<std::string, std::deque<double> > > >::const_iterator it = local_vector_vertex_data.find(segment_id);
00975         if (it == local_vector_vertex_data.end())
00976           return ret;
00977 
00978         for (std::size_t i=0; i<it->second.size(); ++i)
00979           ret.push_back(it->second[i].first);
00980 
00981         return ret;
00982       }
00983 
00985       std::vector<std::string> scalar_cell_data_names(segment_id_type segment_id) const
00986       {
00987         std::vector<std::string> ret;
00988 
00989         std::map<int, std::deque<std::pair<std::string, std::deque<double> > > >::const_iterator it = local_scalar_cell_data.find(segment_id);
00990         if (it == local_scalar_cell_data.end())
00991           return ret;
00992 
00993         for (std::size_t i=0; i<it->second.size(); ++i)
00994           ret.push_back(it->second[i].first);
00995 
00996         return ret;
00997       }
00998 
01000       std::vector<std::string> vector_cell_data_names(segment_id_type segment_id) const
01001       {
01002         std::vector<std::string> ret;
01003 
01004         std::map<int, std::deque<std::pair<std::string, std::deque<double> > > >::const_iterator it = local_vector_cell_data.find(segment_id);
01005         if (it == local_vector_cell_data.end())
01006           return ret;
01007 
01008         for (std::size_t i=0; i<it->second.size(); ++i)
01009           ret.push_back(it->second[i].first);
01010 
01011         return ret;
01012       }
01013 
01014         // Extract data read from file:
01015 
01017         std::vector<std::pair<std::size_t, std::string> > const & get_scalar_data_on_vertices() const
01018         {
01019           return vertex_data_scalar_read;
01020         }
01021 
01023         std::vector<std::pair<std::size_t, std::string> > const & get_vector_data_on_vertices() const
01024         {
01025           return vertex_data_vector_read;
01026         }
01027 
01029         std::vector<std::pair<std::size_t, std::string> > const & get_scalar_data_on_cells() const
01030         {
01031           return cell_data_scalar_read;
01032         }
01033 
01035         std::vector<std::pair<std::size_t, std::string> > const & get_vector_data_on_cells() const
01036         {
01037           return cell_data_vector_read;
01038         }
01039 
01040 
01041 
01042     private:
01043 
01044       template<typename MapType, typename AccessorOrFieldType>
01045       void register_to_map(MapType & map, AccessorOrFieldType accessor_or_field, std::string const & name)
01046       {
01047           typename MapType::iterator it = map.find(name);
01048           if (it != map.end())
01049           {
01050             delete it->second;
01051             it->second = new dynamic_field_wrapper<AccessorOrFieldType>( accessor_or_field );
01052           }
01053           else
01054             map[name] = new dynamic_field_wrapper<AccessorOrFieldType>( accessor_or_field );
01055       }
01056 
01057 
01058     public:
01059 
01061       template<typename AccessorOrFieldType>
01062       void register_vertex_scalar(AccessorOrFieldType accessor_or_field, std::string const & quantity_name)
01063       { register_to_map(registered_vertex_scalar_data, accessor_or_field, quantity_name); }
01064 
01066       template<typename AccessorOrFieldType>
01067       void register_vertex_scalar(segment_id_type seg_id, AccessorOrFieldType accessor_or_field, std::string const & quantity_name)
01068       { register_to_map(registered_segment_vertex_scalar_data[seg_id], accessor_or_field, quantity_name); }
01069 
01071       template<typename AccessorOrFieldType>
01072       void register_vertex_scalar(SegmentHandleType const & segment, AccessorOrFieldType accessor_or_field, std::string const & quantity_name)
01073       { register_vertex_scalar(segment.id(), accessor_or_field, quantity_name); }
01074 
01075 
01077       template<typename AccessorOrFieldType>
01078       void register_vertex_vector(AccessorOrFieldType accessor_or_field, std::string const & quantity_name)
01079       { register_to_map(registered_vertex_vector_data, accessor_or_field, quantity_name); }
01080 
01082       template<typename AccessorOrFieldType>
01083       void register_vertex_vector(segment_id_type seg_id, AccessorOrFieldType accessor_or_field, std::string const & quantity_name)
01084       { register_to_map(registered_segment_vertex_vector_data[seg_id], accessor_or_field, quantity_name); }
01085 
01087       template<typename AccessorOrFieldType>
01088       void register_vertex_vector(SegmentHandleType const & segment, AccessorOrFieldType accessor_or_field, std::string const & quantity_name)
01089       { register_vertex_vector(segment.id(), accessor_or_field, quantity_name); }
01090 
01091 
01092 
01093 
01095       template<typename AccessorOrFieldType>
01096       void register_cell_scalar(AccessorOrFieldType accessor_or_field, std::string const & quantity_name)
01097       { register_to_map(registered_cell_scalar_data, accessor_or_field, quantity_name); }
01098 
01100       template<typename AccessorOrFieldType>
01101       void register_cell_scalar(segment_id_type seg_id, AccessorOrFieldType accessor_or_field, std::string const & quantity_name)
01102       { register_to_map(registered_segment_cell_scalar_data[seg_id], accessor_or_field, quantity_name); }
01103 
01105       template<typename AccessorOrFieldType>
01106       void register_cell_scalar(SegmentHandleType const & segment, AccessorOrFieldType accessor_or_field, std::string const & quantity_name)
01107       { register_cell_scalar(segment.id(), accessor_or_field, quantity_name); }
01108 
01109 
01111       template<typename AccessorOrFieldType>
01112       void register_cell_vector(AccessorOrFieldType accessor_or_field, std::string const & quantity_name)
01113       { register_to_map(registered_cell_vector_data, accessor_or_field, quantity_name); }
01114 
01116       template<typename AccessorOrFieldType>
01117       void register_cell_vector(segment_id_type seg_id, AccessorOrFieldType accessor_or_field, std::string const & quantity_name)
01118       { register_to_map(registered_segment_cell_vector_data[seg_id], accessor_or_field, quantity_name); }
01119 
01121       template<typename AccessorOrFieldType>
01122       void register_cell_vector(SegmentHandleType const & segment, AccessorOrFieldType accessor_or_field, std::string const & quantity_name)
01123       { register_cell_vector(segment.id(), accessor_or_field, quantity_name); }
01124 
01125 
01126 
01127 
01128 
01129 
01131       typename viennagrid::result_of::field<const std::deque<double>, VertexType >::type vertex_scalar_field( std::string const & quantity_name, segment_id_type seg_id ) const
01132       {
01133         typename std::map< std::string, std::map<segment_id_type, std::deque<double> > >::const_iterator it = vertex_scalar_data.find(quantity_name);
01134         if (it == vertex_scalar_data.end()) return typename viennagrid::result_of::field<const std::deque<double>, VertexType >::type();
01135 
01136         typename std::map<segment_id_type, std::deque<double> >::const_iterator jt = it->second.find( seg_id );
01137         if (jt == it->second.end()) return typename viennagrid::result_of::field<const std::deque<double>, VertexType >::type();
01138 
01139         return viennagrid::make_field<VertexType>( jt->second );
01140       }
01141 
01143       typename viennagrid::result_of::field<const std::deque<double>, VertexType >::type vertex_scalar_field( std::string const & quantity_name, SegmentHandleType const & segment ) const
01144       { return vertex_scalar_field(quantity_name, segment.id()); }
01145 
01146 
01148       typename viennagrid::result_of::field<const std::deque<vector_data_type>, VertexType >::type vertex_vector_field( std::string const & quantity_name, segment_id_type seg_id ) const
01149       {
01150         typename std::map< std::string, std::map<segment_id_type, std::deque<vector_data_type> > >::const_iterator it = vertex_vector_data.find(quantity_name);
01151         if (it == vertex_vector_data.end()) return typename viennagrid::result_of::field<const std::deque<vector_data_type>, VertexType >::type();
01152 
01153         typename std::map<segment_id_type, std::deque<vector_data_type> >::const_iterator jt = it->second.find( seg_id );
01154         if (jt == it->second.end()) return typename viennagrid::result_of::field<const std::deque<vector_data_type>, VertexType >::type();
01155 
01156         return viennagrid::make_field<VertexType>( jt->second );
01157       }
01158 
01160       typename viennagrid::result_of::field<const std::deque<vector_data_type>, VertexType >::type vertex_vector_field( std::string const & quantity_name, SegmentHandleType const & segment )
01161       { return vertex_vector_field(quantity_name, segment.id()); }
01162 
01163 
01164 
01166       typename viennagrid::result_of::field<const std::deque<double>, CellType >::type cell_scalar_field( std::string const & quantity_name, segment_id_type seg_id ) const
01167       {
01168         typename std::map< std::string, std::map<segment_id_type, std::deque<double> > >::const_iterator it = cell_scalar_data.find(quantity_name);
01169         if (it == cell_scalar_data.end()) return typename viennagrid::result_of::field<const std::deque<double>, CellType >::type();
01170 
01171         typename std::map<segment_id_type, std::deque<double> >::const_iterator jt = it->second.find( seg_id );
01172         if (jt == it->second.end()) return typename viennagrid::result_of::field<const std::deque<double>, CellType >::type();
01173 
01174         return viennagrid::make_field<CellType>( jt->second );
01175       }
01176 
01178       typename viennagrid::result_of::field<const std::deque<double>, CellType >::type cell_scalar_field( std::string const & quantity_name, SegmentHandleType const & segment ) const
01179       { return cell_scalar_field(quantity_name, segment.id()); }
01180 
01181 
01183       typename viennagrid::result_of::field<const std::deque<vector_data_type>, CellType >::type cell_vector_field( std::string const & quantity_name, segment_id_type seg_id ) const
01184       {
01185         typename std::map< std::string, std::map<segment_id_type, std::deque<vector_data_type> > >::const_iterator it = cell_vector_data.find(quantity_name);
01186         if (it == cell_vector_data.end()) return typename viennagrid::result_of::field<const std::deque<vector_data_type>, CellType >::type();
01187 
01188         typename std::map<segment_id_type, std::deque<vector_data_type> >::const_iterator jt = it->second.find( seg_id );
01189         if (jt == it->second.end()) return typename viennagrid::result_of::field<const std::deque<vector_data_type>, CellType >::type();
01190 
01191         return viennagrid::make_field<CellType>( jt->second );
01192       }
01193 
01195       typename viennagrid::result_of::field<const std::deque<vector_data_type>, CellType >::type cell_vector_field( std::string const & quantity_name, SegmentHandleType const & segment ) const
01196       { return cell_vector_field(quantity_name, segment.id()); }
01197 
01198 
01199     private:
01200 
01201       // Quantities read:
01202       std::vector<std::pair<std::size_t, std::string> >         vertex_data_scalar_read;
01203       std::vector<std::pair<std::size_t, std::string> >         vertex_data_vector_read;
01204 
01205       std::vector<std::pair<std::size_t, std::string> >         cell_data_scalar_read;
01206       std::vector<std::pair<std::size_t, std::string> >         cell_data_vector_read;
01207 
01208 
01209       std::map< std::string, std::map<segment_id_type, std::deque<double> > >           vertex_scalar_data;
01210       std::map< std::string, std::map<segment_id_type, std::deque<vector_data_type> > > vertex_vector_data;
01211 
01212       std::map< std::string, std::map<segment_id_type, std::deque<double> > >           cell_scalar_data;
01213       std::map< std::string, std::map<segment_id_type, std::deque<vector_data_type> > > cell_vector_data;
01214 
01215 
01216       VertexScalarOutputFieldContainer          registered_vertex_scalar_data;
01217       VertexVectorOutputFieldContainer          registered_vertex_vector_data;
01218 
01219       CellScalarOutputFieldContainer          registered_cell_scalar_data;
01220       CellVectorOutputFieldContainer          registered_cell_vector_data;
01221 
01222       std::map< segment_id_type, VertexScalarOutputFieldContainer > registered_segment_vertex_scalar_data;
01223       std::map< segment_id_type, VertexVectorOutputFieldContainer > registered_segment_vertex_vector_data;
01224 
01225       std::map< segment_id_type, CellScalarOutputFieldContainer >   registered_segment_cell_scalar_data;
01226       std::map< segment_id_type, CellVectorOutputFieldContainer >   registered_segment_cell_vector_data;
01227 
01228     }; //class vtk_reader
01229 
01230 
01232     template <typename MeshType, typename SegmentationType >
01233     int import_vtk(MeshType & mesh_obj, SegmentationType & segmentation, std::string const & filename)
01234     {
01235       vtk_reader<MeshType, SegmentationType> vtk_reader;
01236       return vtk_reader(mesh_obj, segmentation, filename);
01237     }
01238 
01240     template <typename MeshType>
01241     int import_vtk(MeshType & mesh_obj, std::string const & filename)
01242     {
01243       vtk_reader<MeshType> vtk_reader;
01244       return vtk_reader(mesh_obj, filename);
01245     }
01246 
01247 
01248 
01249 
01250 
01260     template <typename MeshT, typename SegmentationT, typename AccessorOrFieldT>
01261     vtk_reader<MeshT, SegmentationT> & add_scalar_data_on_vertices(vtk_reader<MeshT, SegmentationT> & reader,
01262                                                                     AccessorOrFieldT const accessor_or_field,
01263                                                                     std::string const & quantity_name)
01264     {
01265       reader.register_vertex_scalar(accessor_or_field, quantity_name);
01266       return reader;
01267     }
01268 
01278     template <typename MeshT, typename SegmentationT, typename AccessorOrFieldT>
01279     vtk_reader<MeshT, SegmentationT> & add_vector_data_on_vertices(vtk_reader<MeshT, SegmentationT> & reader,
01280                                                                     AccessorOrFieldT const accessor_or_field,
01281                                                                     std::string const & quantity_name)
01282     {
01283       reader.register_vertex_vector(accessor_or_field, quantity_name);
01284       return reader;
01285     }
01286 
01287 
01288 
01289 
01290 
01300     template <typename MeshT, typename SegmentationT, typename AccessorOrFieldT>
01301     vtk_reader<MeshT, SegmentationT> & add_scalar_data_on_cells(vtk_reader<MeshT, SegmentationT> & reader,
01302                                                                     AccessorOrFieldT const accessor_or_field,
01303                                                                     std::string const & quantity_name)
01304     {
01305       reader.register_cell_scalar(accessor_or_field, quantity_name);
01306       return reader;
01307     }
01308 
01318     template <typename MeshT, typename SegmentationT, typename AccessorOrFieldT>
01319     vtk_reader<MeshT, SegmentationT> & add_vector_data_on_cells(vtk_reader<MeshT, SegmentationT> & reader,
01320                                                                     AccessorOrFieldT const accessor_or_field,
01321                                                                     std::string const & quantity_name)
01322     {
01323       reader.register_cell_vector(accessor_or_field, quantity_name);
01324       return reader;
01325     }
01326 
01327 
01328 
01329 
01331     template<typename VTKReaderT>
01332     std::vector<std::pair<std::size_t, std::string> > const & get_scalar_data_on_vertices(VTKReaderT const & reader)
01333     {
01334       return reader.get_scalar_data_on_vertices();
01335     }
01336 
01338     template<typename VTKReaderT>
01339     std::vector<std::pair<std::size_t, std::string> > const & get_vector_data_on_vertices(VTKReaderT const & reader)
01340     {
01341       return reader.get_vector_data_on_vertices();
01342     }
01343 
01345     template<typename VTKReaderT>
01346     std::vector<std::pair<std::size_t, std::string> > const & get_scalar_data_on_cells(VTKReaderT const & reader)
01347     {
01348       return reader.get_scalar_data_on_cells();
01349     }
01350 
01352     template<typename VTKReaderT>
01353     std::vector<std::pair<std::size_t, std::string> > const & get_vector_data_on_cells(VTKReaderT const & reader)
01354     {
01355       return reader.get_vector_data_on_cells();
01356     }
01357 
01358 
01359 
01360   } //namespace io
01361 } //namespace viennagrid
01362 
01363 #endif