dune-grid
2.2.0
|
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set et ts=4 sw=2 sts=2: 00003 00004 #ifndef DUNE_VTKWRITER_HH 00005 #define DUNE_VTKWRITER_HH 00006 00007 #include <cstring> 00008 #include <iostream> 00009 #include <string> 00010 #include <fstream> 00011 #include <sstream> 00012 #include <iomanip> 00013 00014 #include <vector> 00015 #include <list> 00016 00017 #include <dune/common/deprecated.hh> 00018 #include <dune/common/exceptions.hh> 00019 #include <dune/common/indent.hh> 00020 #include <dune/common/iteratorfacades.hh> 00021 #include <dune/common/path.hh> 00022 #include <dune/common/shared_ptr.hh> 00023 #include <dune/geometry/referenceelements.hh> 00024 #include <dune/grid/common/mcmgmapper.hh> 00025 #include <dune/grid/common/gridenums.hh> 00026 #include <dune/grid/io/file/vtk/common.hh> 00027 #include <dune/grid/io/file/vtk/dataarraywriter.hh> 00028 #include <dune/grid/io/file/vtk/function.hh> 00029 #include <dune/grid/io/file/vtk/pvtuwriter.hh> 00030 #include <dune/grid/io/file/vtk/streams.hh> 00031 #include <dune/grid/io/file/vtk/vtuwriter.hh> 00032 00050 namespace Dune 00051 { 00060 template< class GridView > 00061 class VTKWriter { 00062 // extract types 00063 typedef typename GridView::Grid Grid; 00064 typedef typename GridView::ctype DT; 00065 enum { n = GridView::dimension }; 00066 enum { w = GridView::dimensionworld }; 00067 00068 typedef typename GridView::template Codim< 0 >::Entity Cell; 00069 typedef typename GridView::template Codim< n >::Entity Vertex; 00070 typedef Cell Entity; 00071 00072 typedef typename GridView::IndexSet IndexSet; 00073 00074 static const PartitionIteratorType VTK_Partition = InteriorBorder_Partition; 00075 00076 typedef typename GridView::template Codim< 0 > 00077 ::template Partition< VTK_Partition >::Iterator 00078 GridCellIterator; 00079 typedef typename GridView::template Codim< n > 00080 ::template Partition< VTK_Partition >::Iterator 00081 GridVertexIterator; 00082 00083 typedef MultipleCodimMultipleGeomTypeMapper< GridView, MCMGVertexLayout > VertexMapper; 00084 00085 public: 00086 typedef Dune::VTKFunction< GridView > VTKFunction; 00087 typedef shared_ptr< const VTKFunction > VTKFunctionPtr; 00088 00089 protected: 00090 typedef typename std::list<VTKFunctionPtr>::const_iterator FunctionIterator; 00091 00093 00098 class CellIterator : public GridCellIterator 00099 { 00100 public: 00102 CellIterator(const GridCellIterator & x) : GridCellIterator(x) {}; 00105 const FieldVector<DT,n> position() const 00106 { 00107 return GenericReferenceElements<DT,n>::general((*this)->type()).position(0,0); 00108 } 00109 }; 00110 00111 CellIterator cellBegin() const 00112 { 00113 return gridView_.template begin< 0, VTK_Partition >(); 00114 } 00115 00116 CellIterator cellEnd() const 00117 { 00118 return gridView_.template end< 0, VTK_Partition >(); 00119 } 00120 00122 00136 class VertexIterator : 00137 public ForwardIteratorFacade<VertexIterator, const Entity, const Entity&, int> 00138 { 00139 GridCellIterator git; 00140 GridCellIterator gend; 00141 VTK::DataMode datamode; 00142 // Index of the currently visited corner within the current element. 00143 // NOTE: this is in Dune-numbering, in contrast to CornerIterator. 00144 int cornerIndexDune; 00145 const VertexMapper & vertexmapper; 00146 std::vector<bool> visited; 00147 // in conforming mode, for each vertex id (as obtained by vertexmapper) 00148 // hold its number in the iteration order (VertexIterator) 00149 int offset; 00150 protected: 00151 void basicIncrement () 00152 { 00153 if( git == gend ) 00154 return; 00155 ++cornerIndexDune; 00156 const int numCorners = git->template count< n >(); 00157 if( cornerIndexDune == numCorners ) 00158 { 00159 offset += numCorners; 00160 cornerIndexDune = 0; 00161 00162 ++git; 00163 while( (git != gend) && (git->partitionType() != InteriorEntity) ) 00164 ++git; 00165 } 00166 } 00167 public: 00168 VertexIterator(const GridCellIterator & x, 00169 const GridCellIterator & end, 00170 const VTK::DataMode & dm, 00171 const VertexMapper & vm) : 00172 git(x), gend(end), datamode(dm), cornerIndexDune(0), 00173 vertexmapper(vm), visited(vm.size(), false), 00174 offset(0) 00175 { 00176 if (datamode == VTK::conforming && git != gend) 00177 visited[vertexmapper.map(*git,cornerIndexDune,n)] = true; 00178 }; 00179 void increment () 00180 { 00181 switch (datamode) 00182 { 00183 case VTK::conforming: 00184 while(visited[vertexmapper.map(*git,cornerIndexDune,n)]) 00185 { 00186 basicIncrement(); 00187 if (git == gend) return; 00188 } 00189 visited[vertexmapper.map(*git,cornerIndexDune,n)] = true; 00190 break; 00191 case VTK::nonconforming: 00192 basicIncrement(); 00193 break; 00194 } 00195 } 00196 bool equals (const VertexIterator & cit) const 00197 { 00198 return git == cit.git 00199 && cornerIndexDune == cit.cornerIndexDune 00200 && datamode == cit.datamode; 00201 } 00202 const Entity& dereference() const 00203 { 00204 return *git; 00205 } 00207 int localindex () const 00208 { 00209 return cornerIndexDune; 00210 } 00212 const FieldVector<DT,n> & position () const 00213 { 00214 return GenericReferenceElements<DT,n>::general(git->type()) 00215 .position(cornerIndexDune,n); 00216 } 00217 }; 00218 00219 VertexIterator vertexBegin () const 00220 { 00221 return VertexIterator( gridView_.template begin< 0, VTK_Partition >(), 00222 gridView_.template end< 0, VTK_Partition >(), 00223 datamode, *vertexmapper ); 00224 } 00225 00226 VertexIterator vertexEnd () const 00227 { 00228 return VertexIterator( gridView_.template end< 0, VTK_Partition >(), 00229 gridView_.template end< 0, VTK_Partition >(), 00230 datamode, *vertexmapper ); 00231 } 00232 00234 00248 class CornerIterator : 00249 public ForwardIteratorFacade<CornerIterator, const Entity, const Entity&, int> 00250 { 00251 GridCellIterator git; 00252 GridCellIterator gend; 00253 VTK::DataMode datamode; 00254 // Index of the currently visited corner within the current element. 00255 // NOTE: this is in VTK-numbering, in contrast to VertexIterator. 00256 int cornerIndexVTK; 00257 const VertexMapper & vertexmapper; 00258 // in conforming mode, for each vertex id (as obtained by vertexmapper) 00259 // hold its number in the iteration order of VertexIterator (*not* 00260 // CornerIterator) 00261 const std::vector<int> & number; 00262 // holds the number of corners of all the elements we have seen so far, 00263 // excluding the current element 00264 int offset; 00265 00266 public: 00267 CornerIterator(const GridCellIterator & x, 00268 const GridCellIterator & end, 00269 const VTK::DataMode & dm, 00270 const VertexMapper & vm, 00271 const std::vector<int> & num) : 00272 git(x), gend(end), datamode(dm), cornerIndexVTK(0), 00273 vertexmapper(vm), 00274 number(num), offset(0) {}; 00275 void increment () 00276 { 00277 if( git == gend ) 00278 return; 00279 ++cornerIndexVTK; 00280 const int numCorners = git->template count< n >(); 00281 if( cornerIndexVTK == numCorners ) 00282 { 00283 offset += numCorners; 00284 cornerIndexVTK = 0; 00285 00286 ++git; 00287 while( (git != gend) && (git->partitionType() != InteriorEntity) ) 00288 ++git; 00289 } 00290 } 00291 bool equals (const CornerIterator & cit) const 00292 { 00293 return git == cit.git 00294 && cornerIndexVTK == cit.cornerIndexVTK 00295 && datamode == cit.datamode; 00296 } 00297 const Entity& dereference() const 00298 { 00299 return *git; 00300 } 00302 00306 int id () const 00307 { 00308 switch (datamode) 00309 { 00310 case VTK::conforming: 00311 return 00312 number[vertexmapper.map(*git,VTK::renumber(*git,cornerIndexVTK), 00313 n)]; 00314 case VTK::nonconforming: 00315 return offset + VTK::renumber(*git,cornerIndexVTK); 00316 default: 00317 DUNE_THROW(IOError,"VTKWriter: unsupported DataMode" << datamode); 00318 } 00319 } 00320 }; 00321 00322 CornerIterator cornerBegin () const 00323 { 00324 return CornerIterator( gridView_.template begin< 0, VTK_Partition >(), 00325 gridView_.template end< 0, VTK_Partition >(), 00326 datamode, *vertexmapper, number ); 00327 } 00328 00329 CornerIterator cornerEnd () const 00330 { 00331 return CornerIterator( gridView_.template end< 0, VTK_Partition >(), 00332 gridView_.template end< 0, VTK_Partition >(), 00333 datamode, *vertexmapper, number ); 00334 } 00335 00336 public: 00344 explicit VTKWriter ( const GridView &gridView, 00345 VTK::DataMode dm = VTK::conforming ) 00346 : gridView_( gridView ), 00347 datamode( dm ) 00348 { } 00349 00354 void addCellData (const VTKFunctionPtr & p) 00355 { 00356 celldata.push_back(p); 00357 } 00358 00364 void addCellData (VTKFunction* p) // DUNE_DEPRECATED 00365 { 00366 celldata.push_back(VTKFunctionPtr(p)); 00367 } 00368 00384 template<class V> 00385 void addCellData (const V& v, const std::string &name, int ncomps = 1) 00386 { 00387 typedef P0VTKFunction<GridView, V> Function; 00388 for (int c=0;c<ncomps;++c) { 00389 std::stringstream compName; 00390 compName << name; 00391 if (ncomps>1) 00392 compName << "[" << c << "]"; 00393 VTKFunction* p = new Function(gridView_, v, compName.str(), ncomps, c); 00394 celldata.push_back(VTKFunctionPtr(p)); 00395 } 00396 } 00397 00403 void addVertexData (VTKFunction* p) // DUNE_DEPRECATED 00404 { 00405 vertexdata.push_back(VTKFunctionPtr(p)); 00406 } 00407 00412 void addVertexData (const VTKFunctionPtr & p) 00413 { 00414 vertexdata.push_back(p); 00415 } 00416 00432 template<class V> 00433 void addVertexData (const V& v, const std::string &name, int ncomps=1) 00434 { 00435 typedef P1VTKFunction<GridView, V> Function; 00436 for (int c=0;c<ncomps;++c) { 00437 std::stringstream compName; 00438 compName << name; 00439 if (ncomps>1) 00440 compName << "[" << c << "]"; 00441 VTKFunction* p = new Function(gridView_, v, compName.str(), ncomps, c); 00442 vertexdata.push_back(VTKFunctionPtr(p)); 00443 } 00444 } 00445 00447 void clear () 00448 { 00449 celldata.clear(); 00450 vertexdata.clear(); 00451 } 00452 00454 virtual ~VTKWriter () 00455 { 00456 this->clear(); 00457 } 00458 00470 std::string write ( const std::string &name, 00471 VTK::OutputType type = VTK::ascii ) 00472 { 00473 return write( name, type, gridView_.comm().rank(), gridView_.comm().size() ); 00474 } 00475 00502 std::string pwrite ( const std::string & name, const std::string & path, const std::string & extendpath, 00503 VTK::OutputType type = VTK::ascii ) 00504 { 00505 return pwrite( name, path, extendpath, type, gridView_.comm().rank(), gridView_.comm().size() ); 00506 } 00507 00508 protected: 00510 00521 std::string getParallelPieceName(const std::string& name, 00522 const std::string& path, 00523 int commRank, int commSize) const 00524 { 00525 std::ostringstream s; 00526 if(path.size() > 0) { 00527 s << path; 00528 if(path[path.size()-1] != '/') 00529 s << '/'; 00530 } 00531 s << 's' << std::setw(4) << std::setfill('0') << commSize << '-'; 00532 s << 'p' << std::setw(4) << std::setfill('0') << commRank << '-'; 00533 s << name; 00534 if(GridView::dimension > 1) 00535 s << ".vtu"; 00536 else 00537 s << ".vtp"; 00538 return s.str(); 00539 } 00540 00542 00552 std::string getParallelHeaderName(const std::string& name, 00553 const std::string& path, 00554 int commSize) const 00555 { 00556 std::ostringstream s; 00557 if(path.size() > 0) { 00558 s << path; 00559 if(path[path.size()-1] != '/') 00560 s << '/'; 00561 } 00562 s << 's' << std::setw(4) << std::setfill('0') << commSize << '-'; 00563 s << name; 00564 if(GridView::dimension > 1) 00565 s << ".pvtu"; 00566 else 00567 s << ".pvtp"; 00568 return s.str(); 00569 } 00570 00572 00584 std::string getSerialPieceName(const std::string& name, 00585 const std::string& path) const 00586 { 00587 static const std::string extension = 00588 GridView::dimension == 1 ? ".vtp" : ".vtu"; 00589 00590 return concatPaths(path, name+extension); 00591 } 00592 00608 std::string write ( const std::string &name, 00609 VTK::OutputType type, 00610 const int commRank, 00611 const int commSize ) 00612 { 00613 // in the parallel case, just use pwrite, it has all the necessary 00614 // stuff, so we don't need to reimplement it here. 00615 if(commSize > 1) 00616 return pwrite(name, "", "", type, commRank, commSize); 00617 00618 // make data mode visible to private functions 00619 outputtype = type; 00620 00621 // generate filename for process data 00622 std::string pieceName = getSerialPieceName(name, ""); 00623 00624 // write process data 00625 std::ofstream file; 00626 file.open( pieceName.c_str(), std::ios::binary ); 00627 if (! file.is_open()) 00628 DUNE_THROW(IOError, "Could not write to piece file " << pieceName); 00629 writeDataFile( file ); 00630 file.close(); 00631 00632 return pieceName; 00633 } 00634 00636 00659 std::string pwrite(const std::string& name, const std::string& path, 00660 const std::string& extendpath, 00661 VTK::OutputType ot, const int commRank, 00662 const int commSize ) 00663 { 00664 // make data mode visible to private functions 00665 outputtype=ot; 00666 00667 // do some magic because paraview can only cope with relative pathes to piece files 00668 std::ofstream file; 00669 std::string piecepath = concatPaths(path, extendpath); 00670 std::string relpiecepath = relativePath(path, piecepath); 00671 00672 // write this processes .vtu/.vtp piece file 00673 std::string fullname = getParallelPieceName(name, piecepath, commRank, 00674 commSize); 00675 file.open(fullname.c_str(),std::ios::binary); 00676 if (! file.is_open()) 00677 DUNE_THROW(IOError, "Could not write to piecefile file " << fullname); 00678 writeDataFile(file); 00679 file.close(); 00680 gridView_.comm().barrier(); 00681 00682 // if we are rank 0, write .pvtu/.pvtp parallel header 00683 fullname = getParallelHeaderName(name, path, commSize); 00684 if( commRank ==0 ) 00685 { 00686 file.open(fullname.c_str()); 00687 if (! file.is_open()) 00688 DUNE_THROW(IOError, "Could not write to parallel file " << fullname); 00689 writeParallelHeader(file,name,relpiecepath, commSize ); 00690 file.close(); 00691 } 00692 gridView_.comm().barrier(); 00693 return fullname; 00694 } 00695 00696 private: 00698 00715 void writeParallelHeader(std::ostream& s, const std::string& piecename, 00716 const std::string& piecepath, const int commSize) 00717 { 00718 VTK::FileType fileType = 00719 (n == 1) ? VTK::polyData : VTK::unstructuredGrid; 00720 00721 VTK::PVTUWriter writer(s, fileType); 00722 00723 writer.beginMain(); 00724 00725 // PPointData 00726 { 00727 std::string scalars; 00728 for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); 00729 ++it) 00730 if ((*it)->ncomps()==1) 00731 { 00732 scalars = (*it)->name(); 00733 break; 00734 } 00735 std::string vectors; 00736 for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); 00737 ++it) 00738 if ((*it)->ncomps()>1) 00739 { 00740 vectors = (*it)->name(); 00741 break; 00742 } 00743 writer.beginPointData(scalars, vectors); 00744 } 00745 for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); 00746 ++it) 00747 { 00748 unsigned writecomps = (*it)->ncomps(); 00749 if(writecomps == 2) writecomps = 3; 00750 writer.addArray<float>((*it)->name(), writecomps); 00751 } 00752 writer.endPointData(); 00753 00754 // PCellData 00755 { 00756 std::string scalars; 00757 for (FunctionIterator it=celldata.begin(); it!=celldata.end(); 00758 ++it) 00759 if ((*it)->ncomps()==1) 00760 { 00761 scalars = (*it)->name(); 00762 break; 00763 } 00764 std::string vectors; 00765 for (FunctionIterator it=celldata.begin(); it!=celldata.end(); 00766 ++it) 00767 if ((*it)->ncomps()>1) 00768 { 00769 vectors = (*it)->name(); 00770 break; 00771 } 00772 writer.beginCellData(scalars, vectors); 00773 } 00774 for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it) { 00775 unsigned writecomps = (*it)->ncomps(); 00776 if(writecomps == 2) writecomps = 3; 00777 writer.addArray<float>((*it)->name(), writecomps); 00778 } 00779 writer.endCellData(); 00780 00781 // PPoints 00782 writer.beginPoints(); 00783 writer.addArray<float>("Coordinates", 3); 00784 writer.endPoints(); 00785 00786 // Pieces 00787 for( int i = 0; i < commSize; ++i ) 00788 { 00789 const std::string& fullname = getParallelPieceName(piecename, 00790 piecepath, i, 00791 commSize); 00792 writer.addPiece(fullname); 00793 } 00794 00795 writer.endMain(); 00796 } 00797 00799 void writeDataFile (std::ostream& s) 00800 { 00801 VTK::FileType fileType = 00802 (n == 1) ? VTK::polyData : VTK::unstructuredGrid; 00803 00804 VTK::VTUWriter writer(s, outputtype, fileType); 00805 00806 // Grid characteristics 00807 vertexmapper = new VertexMapper( gridView_ ); 00808 if (datamode == VTK::conforming) 00809 { 00810 number.resize(vertexmapper->size()); 00811 for (std::vector<int>::size_type i=0; i<number.size(); i++) number[i] = -1; 00812 } 00813 countEntities(nvertices, ncells, ncorners); 00814 00815 writer.beginMain(ncells, nvertices); 00816 writeAllData(writer); 00817 writer.endMain(); 00818 00819 // write appended binary data section 00820 if(writer.beginAppended()) 00821 writeAllData(writer); 00822 writer.endAppended(); 00823 00824 delete vertexmapper; number.clear(); 00825 } 00826 00827 void writeAllData(VTK::VTUWriter& writer) { 00828 // PointData 00829 writeVertexData(writer); 00830 00831 // CellData 00832 writeCellData(writer); 00833 00834 // Points 00835 writeGridPoints(writer); 00836 00837 // Cells 00838 writeGridCells(writer); 00839 } 00840 00841 protected: 00842 std::string getFormatString() const 00843 { 00844 if (outputtype==VTK::ascii) 00845 return "ascii"; 00846 if (outputtype==VTK::base64) 00847 return "binary"; 00848 if (outputtype==VTK::appendedraw) 00849 return "appended"; 00850 if (outputtype==VTK::appendedbase64) 00851 return "appended"; 00852 DUNE_THROW(IOError, "VTKWriter: unsupported OutputType" << outputtype); 00853 } 00854 00855 std::string getTypeString() const 00856 { 00857 if (n==1) 00858 return "PolyData"; 00859 else 00860 return "UnstructuredGrid"; 00861 } 00862 00864 virtual void countEntities(int &nvertices, int &ncells, int &ncorners) 00865 { 00866 nvertices = 0; 00867 ncells = 0; 00868 ncorners = 0; 00869 for (CellIterator it=cellBegin(); it!=cellEnd(); ++it) 00870 { 00871 ncells++; 00872 // because of the use of vertexmapper->map(), this iteration must be 00873 // in the order of Dune's numbering. 00874 for (int i=0; i<it->template count<n>(); ++i) 00875 { 00876 ncorners++; 00877 if (datamode == VTK::conforming) 00878 { 00879 int alpha = vertexmapper->map(*it,i,n); 00880 if (number[alpha]<0) 00881 number[alpha] = nvertices++; 00882 } 00883 else 00884 { 00885 nvertices++; 00886 } 00887 } 00888 } 00889 } 00890 00892 virtual void writeCellData(VTK::VTUWriter& writer) 00893 { 00894 if(celldata.size() == 0) 00895 return; 00896 00897 std::string scalars = ""; 00898 for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it) 00899 if ((*it)->ncomps()==1) 00900 { 00901 scalars = (*it)->name(); 00902 break; 00903 } 00904 std::string vectors = ""; 00905 for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it) 00906 if ((*it)->ncomps()>1) 00907 { 00908 vectors = (*it)->name(); 00909 break; 00910 } 00911 00912 writer.beginCellData(scalars, vectors); 00913 for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it) 00914 { 00915 // vtk file format: a vector data always should have 3 comps (with 00916 // 3rd comp = 0 in 2D case) 00917 unsigned writecomps = (*it)->ncomps(); 00918 if(writecomps == 2) writecomps = 3; 00919 shared_ptr<VTK::DataArrayWriter<float> > p 00920 (writer.makeArrayWriter<float>((*it)->name(), writecomps, 00921 ncells)); 00922 if(!p->writeIsNoop()) 00923 for (CellIterator i=cellBegin(); i!=cellEnd(); ++i) 00924 { 00925 for (int j=0; j<(*it)->ncomps(); j++) 00926 p->write((*it)->evaluate(j,*i,i.position())); 00927 // vtk file format: a vector data always should have 3 comps 00928 // (with 3rd comp = 0 in 2D case) 00929 for (unsigned j=(*it)->ncomps(); j < writecomps; ++j) 00930 p->write(0.0); 00931 } 00932 } 00933 writer.endCellData(); 00934 } 00935 00937 virtual void writeVertexData(VTK::VTUWriter& writer) 00938 { 00939 if(vertexdata.size() == 0) 00940 return; 00941 00942 std::string scalars = ""; 00943 for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it) 00944 if ((*it)->ncomps()==1) 00945 { 00946 scalars = (*it)->name(); 00947 break; 00948 } 00949 std::string vectors = ""; 00950 for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it) 00951 if ((*it)->ncomps()>1) 00952 { 00953 vectors = (*it)->name(); 00954 break; 00955 } 00956 00957 writer.beginPointData(scalars, vectors); 00958 for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it) 00959 { 00960 // vtk file format: a vector data always should have 3 comps (with 00961 // 3rd comp = 0 in 2D case) 00962 unsigned writecomps = (*it)->ncomps(); 00963 if(writecomps == 2) writecomps = 3; 00964 shared_ptr<VTK::DataArrayWriter<float> > p 00965 (writer.makeArrayWriter<float>((*it)->name(), writecomps, 00966 nvertices)); 00967 if(!p->writeIsNoop()) 00968 for (VertexIterator vit=vertexBegin(); vit!=vertexEnd(); ++vit) 00969 { 00970 for (int j=0; j<(*it)->ncomps(); j++) 00971 p->write((*it)->evaluate(j,*vit,vit.position())); 00972 // vtk file format: a vector data always should have 3 comps 00973 // (with 3rd comp = 0 in 2D case) 00974 for (unsigned j=(*it)->ncomps(); j < writecomps; ++j) 00975 p->write(0.0); 00976 } 00977 } 00978 writer.endPointData(); 00979 } 00980 00982 virtual void writeGridPoints(VTK::VTUWriter& writer) 00983 { 00984 writer.beginPoints(); 00985 00986 shared_ptr<VTK::DataArrayWriter<float> > p 00987 (writer.makeArrayWriter<float>("Coordinates", 3, nvertices)); 00988 if(!p->writeIsNoop()) { 00989 VertexIterator vEnd = vertexEnd(); 00990 for (VertexIterator vit=vertexBegin(); vit!=vEnd; ++vit) 00991 { 00992 int dimw=w; 00993 for (int j=0; j<std::min(dimw,3); j++) 00994 p->write(vit->geometry().corner(vit.localindex())[j]); 00995 for (int j=std::min(dimw,3); j<3; j++) 00996 p->write(0.0); 00997 } 00998 } 00999 // free the VTK::DataArrayWriter before touching the stream 01000 p.reset(); 01001 01002 writer.endPoints(); 01003 } 01004 01006 virtual void writeGridCells(VTK::VTUWriter& writer) 01007 { 01008 writer.beginCells(); 01009 01010 // connectivity 01011 { 01012 shared_ptr<VTK::DataArrayWriter<int> > p1 01013 (writer.makeArrayWriter<int>("connectivity", 1, ncorners)); 01014 if(!p1->writeIsNoop()) 01015 for (CornerIterator it=cornerBegin(); it!=cornerEnd(); ++it) 01016 p1->write(it.id()); 01017 } 01018 01019 // offsets 01020 { 01021 shared_ptr<VTK::DataArrayWriter<int> > p2 01022 (writer.makeArrayWriter<int>("offsets", 1, ncells)); 01023 if(!p2->writeIsNoop()) { 01024 int offset = 0; 01025 for (CellIterator it=cellBegin(); it!=cellEnd(); ++it) 01026 { 01027 offset += it->template count<n>(); 01028 p2->write(offset); 01029 } 01030 } 01031 } 01032 01033 // types 01034 if (n>1) 01035 { 01036 shared_ptr<VTK::DataArrayWriter<unsigned char> > p3 01037 (writer.makeArrayWriter<unsigned char>("types", 1, ncells)); 01038 if(!p3->writeIsNoop()) 01039 for (CellIterator it=cellBegin(); it!=cellEnd(); ++it) 01040 { 01041 int vtktype = VTK::geometryType(it->type()); 01042 p3->write(vtktype); 01043 } 01044 } 01045 01046 writer.endCells(); 01047 } 01048 01049 protected: 01050 // the list of registered functions 01051 std::list<VTKFunctionPtr> celldata; 01052 std::list<VTKFunctionPtr> vertexdata; 01053 01054 private: 01055 // the grid 01056 GridView gridView_; 01057 01058 // temporary grid information 01059 protected: 01060 int ncells; 01061 int nvertices; 01062 int ncorners; 01063 private: 01064 VertexMapper* vertexmapper; 01065 // in conforming mode, for each vertex id (as obtained by vertexmapper) 01066 // hold its number in the iteration order (VertexIterator) 01067 std::vector<int> number; 01068 VTK::DataMode datamode; 01069 protected: 01070 VTK::OutputType outputtype; 01071 }; 01072 01073 } 01074 01075 #endif