ViennaGrid - The Vienna Grid Library
2.1.0
|
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