dune-grid  2.2.0
vtkwriter.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #ifndef DUNE_VTKWRITER_HH
5 #define DUNE_VTKWRITER_HH
6 
7 #include <cstring>
8 #include <iostream>
9 #include <string>
10 #include <fstream>
11 #include <sstream>
12 #include <iomanip>
13 
14 #include <vector>
15 #include <list>
16 
17 #include <dune/common/deprecated.hh>
18 #include <dune/common/exceptions.hh>
19 #include <dune/common/indent.hh>
20 #include <dune/common/iteratorfacades.hh>
21 #include <dune/common/path.hh>
22 #include <dune/common/shared_ptr.hh>
23 #include <dune/geometry/referenceelements.hh>
32 
50 namespace Dune
51 {
60  template< class GridView >
61  class VTKWriter {
62  // extract types
63  typedef typename GridView::Grid Grid;
64  typedef typename GridView::ctype DT;
65  enum { n = GridView::dimension };
66  enum { w = GridView::dimensionworld };
67 
68  typedef typename GridView::template Codim< 0 >::Entity Cell;
69  typedef typename GridView::template Codim< n >::Entity Vertex;
70  typedef Cell Entity;
71 
72  typedef typename GridView::IndexSet IndexSet;
73 
74  static const PartitionIteratorType VTK_Partition = InteriorBorder_Partition;
75 
76  typedef typename GridView::template Codim< 0 >
77  ::template Partition< VTK_Partition >::Iterator
78  GridCellIterator;
79  typedef typename GridView::template Codim< n >
80  ::template Partition< VTK_Partition >::Iterator
81  GridVertexIterator;
82 
84 
85  public:
87  typedef shared_ptr< const VTKFunction > VTKFunctionPtr;
88 
89  protected:
90  typedef typename std::list<VTKFunctionPtr>::const_iterator FunctionIterator;
91 
93 
98  class CellIterator : public GridCellIterator
99  {
100  public:
102  CellIterator(const GridCellIterator & x) : GridCellIterator(x) {};
105  const FieldVector<DT,n> position() const
106  {
107  return GenericReferenceElements<DT,n>::general((*this)->type()).position(0,0);
108  }
109  };
110 
112  {
113  return gridView_.template begin< 0, VTK_Partition >();
114  }
115 
117  {
118  return gridView_.template end< 0, VTK_Partition >();
119  }
120 
122 
137  public ForwardIteratorFacade<VertexIterator, const Entity, const Entity&, int>
138  {
139  GridCellIterator git;
140  GridCellIterator gend;
141  VTK::DataMode datamode;
142  // Index of the currently visited corner within the current element.
143  // NOTE: this is in Dune-numbering, in contrast to CornerIterator.
144  int cornerIndexDune;
145  const VertexMapper & vertexmapper;
146  std::vector<bool> visited;
147  // in conforming mode, for each vertex id (as obtained by vertexmapper)
148  // hold its number in the iteration order (VertexIterator)
149  int offset;
150  protected:
152  {
153  if( git == gend )
154  return;
155  ++cornerIndexDune;
156  const int numCorners = git->template count< n >();
157  if( cornerIndexDune == numCorners )
158  {
159  offset += numCorners;
160  cornerIndexDune = 0;
161 
162  ++git;
163  while( (git != gend) && (git->partitionType() != InteriorEntity) )
164  ++git;
165  }
166  }
167  public:
168  VertexIterator(const GridCellIterator & x,
169  const GridCellIterator & end,
170  const VTK::DataMode & dm,
171  const VertexMapper & vm) :
172  git(x), gend(end), datamode(dm), cornerIndexDune(0),
173  vertexmapper(vm), visited(vm.size(), false),
174  offset(0)
175  {
176  if (datamode == VTK::conforming && git != gend)
177  visited[vertexmapper.map(*git,cornerIndexDune,n)] = true;
178  };
179  void increment ()
180  {
181  switch (datamode)
182  {
183  case VTK::conforming:
184  while(visited[vertexmapper.map(*git,cornerIndexDune,n)])
185  {
186  basicIncrement();
187  if (git == gend) return;
188  }
189  visited[vertexmapper.map(*git,cornerIndexDune,n)] = true;
190  break;
191  case VTK::nonconforming:
192  basicIncrement();
193  break;
194  }
195  }
196  bool equals (const VertexIterator & cit) const
197  {
198  return git == cit.git
199  && cornerIndexDune == cit.cornerIndexDune
200  && datamode == cit.datamode;
201  }
202  const Entity& dereference() const
203  {
204  return *git;
205  }
207  int localindex () const
208  {
209  return cornerIndexDune;
210  }
212  const FieldVector<DT,n> & position () const
213  {
214  return GenericReferenceElements<DT,n>::general(git->type())
215  .position(cornerIndexDune,n);
216  }
217  };
218 
220  {
221  return VertexIterator( gridView_.template begin< 0, VTK_Partition >(),
222  gridView_.template end< 0, VTK_Partition >(),
223  datamode, *vertexmapper );
224  }
225 
227  {
228  return VertexIterator( gridView_.template end< 0, VTK_Partition >(),
229  gridView_.template end< 0, VTK_Partition >(),
230  datamode, *vertexmapper );
231  }
232 
234 
249  public ForwardIteratorFacade<CornerIterator, const Entity, const Entity&, int>
250  {
251  GridCellIterator git;
252  GridCellIterator gend;
253  VTK::DataMode datamode;
254  // Index of the currently visited corner within the current element.
255  // NOTE: this is in VTK-numbering, in contrast to VertexIterator.
256  int cornerIndexVTK;
257  const VertexMapper & vertexmapper;
258  // in conforming mode, for each vertex id (as obtained by vertexmapper)
259  // hold its number in the iteration order of VertexIterator (*not*
260  // CornerIterator)
261  const std::vector<int> & number;
262  // holds the number of corners of all the elements we have seen so far,
263  // excluding the current element
264  int offset;
265 
266  public:
267  CornerIterator(const GridCellIterator & x,
268  const GridCellIterator & end,
269  const VTK::DataMode & dm,
270  const VertexMapper & vm,
271  const std::vector<int> & num) :
272  git(x), gend(end), datamode(dm), cornerIndexVTK(0),
273  vertexmapper(vm),
274  number(num), offset(0) {};
275  void increment ()
276  {
277  if( git == gend )
278  return;
279  ++cornerIndexVTK;
280  const int numCorners = git->template count< n >();
281  if( cornerIndexVTK == numCorners )
282  {
283  offset += numCorners;
284  cornerIndexVTK = 0;
285 
286  ++git;
287  while( (git != gend) && (git->partitionType() != InteriorEntity) )
288  ++git;
289  }
290  }
291  bool equals (const CornerIterator & cit) const
292  {
293  return git == cit.git
294  && cornerIndexVTK == cit.cornerIndexVTK
295  && datamode == cit.datamode;
296  }
297  const Entity& dereference() const
298  {
299  return *git;
300  }
302 
306  int id () const
307  {
308  switch (datamode)
309  {
310  case VTK::conforming:
311  return
312  number[vertexmapper.map(*git,VTK::renumber(*git,cornerIndexVTK),
313  n)];
314  case VTK::nonconforming:
315  return offset + VTK::renumber(*git,cornerIndexVTK);
316  default:
317  DUNE_THROW(IOError,"VTKWriter: unsupported DataMode" << datamode);
318  }
319  }
320  };
321 
323  {
324  return CornerIterator( gridView_.template begin< 0, VTK_Partition >(),
325  gridView_.template end< 0, VTK_Partition >(),
326  datamode, *vertexmapper, number );
327  }
328 
330  {
331  return CornerIterator( gridView_.template end< 0, VTK_Partition >(),
332  gridView_.template end< 0, VTK_Partition >(),
333  datamode, *vertexmapper, number );
334  }
335 
336  public:
344  explicit VTKWriter ( const GridView &gridView,
346  : gridView_( gridView ),
347  datamode( dm )
348  { }
349 
354  void addCellData (const VTKFunctionPtr & p)
355  {
356  celldata.push_back(p);
357  }
358 
364  void addCellData (VTKFunction* p) // DUNE_DEPRECATED
365  {
366  celldata.push_back(VTKFunctionPtr(p));
367  }
368 
384  template<class V>
385  void addCellData (const V& v, const std::string &name, int ncomps = 1)
386  {
387  typedef P0VTKFunction<GridView, V> Function;
388  for (int c=0;c<ncomps;++c) {
389  std::stringstream compName;
390  compName << name;
391  if (ncomps>1)
392  compName << "[" << c << "]";
393  VTKFunction* p = new Function(gridView_, v, compName.str(), ncomps, c);
394  celldata.push_back(VTKFunctionPtr(p));
395  }
396  }
397 
403  void addVertexData (VTKFunction* p) // DUNE_DEPRECATED
404  {
405  vertexdata.push_back(VTKFunctionPtr(p));
406  }
407 
412  void addVertexData (const VTKFunctionPtr & p)
413  {
414  vertexdata.push_back(p);
415  }
416 
432  template<class V>
433  void addVertexData (const V& v, const std::string &name, int ncomps=1)
434  {
435  typedef P1VTKFunction<GridView, V> Function;
436  for (int c=0;c<ncomps;++c) {
437  std::stringstream compName;
438  compName << name;
439  if (ncomps>1)
440  compName << "[" << c << "]";
441  VTKFunction* p = new Function(gridView_, v, compName.str(), ncomps, c);
442  vertexdata.push_back(VTKFunctionPtr(p));
443  }
444  }
445 
447  void clear ()
448  {
449  celldata.clear();
450  vertexdata.clear();
451  }
452 
454  virtual ~VTKWriter ()
455  {
456  this->clear();
457  }
458 
470  std::string write ( const std::string &name,
471  VTK::OutputType type = VTK::ascii )
472  {
473  return write( name, type, gridView_.comm().rank(), gridView_.comm().size() );
474  }
475 
502  std::string pwrite ( const std::string & name, const std::string & path, const std::string & extendpath,
503  VTK::OutputType type = VTK::ascii )
504  {
505  return pwrite( name, path, extendpath, type, gridView_.comm().rank(), gridView_.comm().size() );
506  }
507 
508  protected:
510 
521  std::string getParallelPieceName(const std::string& name,
522  const std::string& path,
523  int commRank, int commSize) const
524  {
525  std::ostringstream s;
526  if(path.size() > 0) {
527  s << path;
528  if(path[path.size()-1] != '/')
529  s << '/';
530  }
531  s << 's' << std::setw(4) << std::setfill('0') << commSize << '-';
532  s << 'p' << std::setw(4) << std::setfill('0') << commRank << '-';
533  s << name;
534  if(GridView::dimension > 1)
535  s << ".vtu";
536  else
537  s << ".vtp";
538  return s.str();
539  }
540 
542 
552  std::string getParallelHeaderName(const std::string& name,
553  const std::string& path,
554  int commSize) const
555  {
556  std::ostringstream s;
557  if(path.size() > 0) {
558  s << path;
559  if(path[path.size()-1] != '/')
560  s << '/';
561  }
562  s << 's' << std::setw(4) << std::setfill('0') << commSize << '-';
563  s << name;
564  if(GridView::dimension > 1)
565  s << ".pvtu";
566  else
567  s << ".pvtp";
568  return s.str();
569  }
570 
572 
584  std::string getSerialPieceName(const std::string& name,
585  const std::string& path) const
586  {
587  static const std::string extension =
588  GridView::dimension == 1 ? ".vtp" : ".vtu";
589 
590  return concatPaths(path, name+extension);
591  }
592 
608  std::string write ( const std::string &name,
609  VTK::OutputType type,
610  const int commRank,
611  const int commSize )
612  {
613  // in the parallel case, just use pwrite, it has all the necessary
614  // stuff, so we don't need to reimplement it here.
615  if(commSize > 1)
616  return pwrite(name, "", "", type, commRank, commSize);
617 
618  // make data mode visible to private functions
619  outputtype = type;
620 
621  // generate filename for process data
622  std::string pieceName = getSerialPieceName(name, "");
623 
624  // write process data
625  std::ofstream file;
626  file.open( pieceName.c_str(), std::ios::binary );
627  if (! file.is_open())
628  DUNE_THROW(IOError, "Could not write to piece file " << pieceName);
629  writeDataFile( file );
630  file.close();
631 
632  return pieceName;
633  }
634 
636 
659  std::string pwrite(const std::string& name, const std::string& path,
660  const std::string& extendpath,
661  VTK::OutputType ot, const int commRank,
662  const int commSize )
663  {
664  // make data mode visible to private functions
665  outputtype=ot;
666 
667  // do some magic because paraview can only cope with relative pathes to piece files
668  std::ofstream file;
669  std::string piecepath = concatPaths(path, extendpath);
670  std::string relpiecepath = relativePath(path, piecepath);
671 
672  // write this processes .vtu/.vtp piece file
673  std::string fullname = getParallelPieceName(name, piecepath, commRank,
674  commSize);
675  file.open(fullname.c_str(),std::ios::binary);
676  if (! file.is_open())
677  DUNE_THROW(IOError, "Could not write to piecefile file " << fullname);
678  writeDataFile(file);
679  file.close();
680  gridView_.comm().barrier();
681 
682  // if we are rank 0, write .pvtu/.pvtp parallel header
683  fullname = getParallelHeaderName(name, path, commSize);
684  if( commRank ==0 )
685  {
686  file.open(fullname.c_str());
687  if (! file.is_open())
688  DUNE_THROW(IOError, "Could not write to parallel file " << fullname);
689  writeParallelHeader(file,name,relpiecepath, commSize );
690  file.close();
691  }
692  gridView_.comm().barrier();
693  return fullname;
694  }
695 
696  private:
698 
715  void writeParallelHeader(std::ostream& s, const std::string& piecename,
716  const std::string& piecepath, const int commSize)
717  {
718  VTK::FileType fileType =
720 
721  VTK::PVTUWriter writer(s, fileType);
722 
723  writer.beginMain();
724 
725  // PPointData
726  {
727  std::string scalars;
728  for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end();
729  ++it)
730  if ((*it)->ncomps()==1)
731  {
732  scalars = (*it)->name();
733  break;
734  }
735  std::string vectors;
736  for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end();
737  ++it)
738  if ((*it)->ncomps()>1)
739  {
740  vectors = (*it)->name();
741  break;
742  }
743  writer.beginPointData(scalars, vectors);
744  }
745  for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end();
746  ++it)
747  {
748  unsigned writecomps = (*it)->ncomps();
749  if(writecomps == 2) writecomps = 3;
750  writer.addArray<float>((*it)->name(), writecomps);
751  }
752  writer.endPointData();
753 
754  // PCellData
755  {
756  std::string scalars;
757  for (FunctionIterator it=celldata.begin(); it!=celldata.end();
758  ++it)
759  if ((*it)->ncomps()==1)
760  {
761  scalars = (*it)->name();
762  break;
763  }
764  std::string vectors;
765  for (FunctionIterator it=celldata.begin(); it!=celldata.end();
766  ++it)
767  if ((*it)->ncomps()>1)
768  {
769  vectors = (*it)->name();
770  break;
771  }
772  writer.beginCellData(scalars, vectors);
773  }
774  for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it) {
775  unsigned writecomps = (*it)->ncomps();
776  if(writecomps == 2) writecomps = 3;
777  writer.addArray<float>((*it)->name(), writecomps);
778  }
779  writer.endCellData();
780 
781  // PPoints
782  writer.beginPoints();
783  writer.addArray<float>("Coordinates", 3);
784  writer.endPoints();
785 
786  // Pieces
787  for( int i = 0; i < commSize; ++i )
788  {
789  const std::string& fullname = getParallelPieceName(piecename,
790  piecepath, i,
791  commSize);
792  writer.addPiece(fullname);
793  }
794 
795  writer.endMain();
796  }
797 
799  void writeDataFile (std::ostream& s)
800  {
801  VTK::FileType fileType =
803 
804  VTK::VTUWriter writer(s, outputtype, fileType);
805 
806  // Grid characteristics
807  vertexmapper = new VertexMapper( gridView_ );
808  if (datamode == VTK::conforming)
809  {
810  number.resize(vertexmapper->size());
811  for (std::vector<int>::size_type i=0; i<number.size(); i++) number[i] = -1;
812  }
814 
815  writer.beginMain(ncells, nvertices);
816  writeAllData(writer);
817  writer.endMain();
818 
819  // write appended binary data section
820  if(writer.beginAppended())
821  writeAllData(writer);
822  writer.endAppended();
823 
824  delete vertexmapper; number.clear();
825  }
826 
827  void writeAllData(VTK::VTUWriter& writer) {
828  // PointData
829  writeVertexData(writer);
830 
831  // CellData
832  writeCellData(writer);
833 
834  // Points
835  writeGridPoints(writer);
836 
837  // Cells
838  writeGridCells(writer);
839  }
840 
841  protected:
842  std::string getFormatString() const
843  {
844  if (outputtype==VTK::ascii)
845  return "ascii";
846  if (outputtype==VTK::base64)
847  return "binary";
849  return "appended";
851  return "appended";
852  DUNE_THROW(IOError, "VTKWriter: unsupported OutputType" << outputtype);
853  }
854 
855  std::string getTypeString() const
856  {
857  if (n==1)
858  return "PolyData";
859  else
860  return "UnstructuredGrid";
861  }
862 
864  virtual void countEntities(int &nvertices, int &ncells, int &ncorners)
865  {
866  nvertices = 0;
867  ncells = 0;
868  ncorners = 0;
869  for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
870  {
871  ncells++;
872  // because of the use of vertexmapper->map(), this iteration must be
873  // in the order of Dune's numbering.
874  for (int i=0; i<it->template count<n>(); ++i)
875  {
876  ncorners++;
877  if (datamode == VTK::conforming)
878  {
879  int alpha = vertexmapper->map(*it,i,n);
880  if (number[alpha]<0)
881  number[alpha] = nvertices++;
882  }
883  else
884  {
885  nvertices++;
886  }
887  }
888  }
889  }
890 
892  virtual void writeCellData(VTK::VTUWriter& writer)
893  {
894  if(celldata.size() == 0)
895  return;
896 
897  std::string scalars = "";
898  for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it)
899  if ((*it)->ncomps()==1)
900  {
901  scalars = (*it)->name();
902  break;
903  }
904  std::string vectors = "";
905  for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it)
906  if ((*it)->ncomps()>1)
907  {
908  vectors = (*it)->name();
909  break;
910  }
911 
912  writer.beginCellData(scalars, vectors);
913  for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it)
914  {
915  // vtk file format: a vector data always should have 3 comps (with
916  // 3rd comp = 0 in 2D case)
917  unsigned writecomps = (*it)->ncomps();
918  if(writecomps == 2) writecomps = 3;
919  shared_ptr<VTK::DataArrayWriter<float> > p
920  (writer.makeArrayWriter<float>((*it)->name(), writecomps,
921  ncells));
922  if(!p->writeIsNoop())
923  for (CellIterator i=cellBegin(); i!=cellEnd(); ++i)
924  {
925  for (int j=0; j<(*it)->ncomps(); j++)
926  p->write((*it)->evaluate(j,*i,i.position()));
927  // vtk file format: a vector data always should have 3 comps
928  // (with 3rd comp = 0 in 2D case)
929  for (unsigned j=(*it)->ncomps(); j < writecomps; ++j)
930  p->write(0.0);
931  }
932  }
933  writer.endCellData();
934  }
935 
937  virtual void writeVertexData(VTK::VTUWriter& writer)
938  {
939  if(vertexdata.size() == 0)
940  return;
941 
942  std::string scalars = "";
943  for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it)
944  if ((*it)->ncomps()==1)
945  {
946  scalars = (*it)->name();
947  break;
948  }
949  std::string vectors = "";
950  for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it)
951  if ((*it)->ncomps()>1)
952  {
953  vectors = (*it)->name();
954  break;
955  }
956 
957  writer.beginPointData(scalars, vectors);
958  for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it)
959  {
960  // vtk file format: a vector data always should have 3 comps (with
961  // 3rd comp = 0 in 2D case)
962  unsigned writecomps = (*it)->ncomps();
963  if(writecomps == 2) writecomps = 3;
964  shared_ptr<VTK::DataArrayWriter<float> > p
965  (writer.makeArrayWriter<float>((*it)->name(), writecomps,
966  nvertices));
967  if(!p->writeIsNoop())
968  for (VertexIterator vit=vertexBegin(); vit!=vertexEnd(); ++vit)
969  {
970  for (int j=0; j<(*it)->ncomps(); j++)
971  p->write((*it)->evaluate(j,*vit,vit.position()));
972  // vtk file format: a vector data always should have 3 comps
973  // (with 3rd comp = 0 in 2D case)
974  for (unsigned j=(*it)->ncomps(); j < writecomps; ++j)
975  p->write(0.0);
976  }
977  }
978  writer.endPointData();
979  }
980 
982  virtual void writeGridPoints(VTK::VTUWriter& writer)
983  {
984  writer.beginPoints();
985 
986  shared_ptr<VTK::DataArrayWriter<float> > p
987  (writer.makeArrayWriter<float>("Coordinates", 3, nvertices));
988  if(!p->writeIsNoop()) {
989  VertexIterator vEnd = vertexEnd();
990  for (VertexIterator vit=vertexBegin(); vit!=vEnd; ++vit)
991  {
992  int dimw=w;
993  for (int j=0; j<std::min(dimw,3); j++)
994  p->write(vit->geometry().corner(vit.localindex())[j]);
995  for (int j=std::min(dimw,3); j<3; j++)
996  p->write(0.0);
997  }
998  }
999  // free the VTK::DataArrayWriter before touching the stream
1000  p.reset();
1001 
1002  writer.endPoints();
1003  }
1004 
1006  virtual void writeGridCells(VTK::VTUWriter& writer)
1007  {
1008  writer.beginCells();
1009 
1010  // connectivity
1011  {
1012  shared_ptr<VTK::DataArrayWriter<int> > p1
1013  (writer.makeArrayWriter<int>("connectivity", 1, ncorners));
1014  if(!p1->writeIsNoop())
1015  for (CornerIterator it=cornerBegin(); it!=cornerEnd(); ++it)
1016  p1->write(it.id());
1017  }
1018 
1019  // offsets
1020  {
1021  shared_ptr<VTK::DataArrayWriter<int> > p2
1022  (writer.makeArrayWriter<int>("offsets", 1, ncells));
1023  if(!p2->writeIsNoop()) {
1024  int offset = 0;
1025  for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
1026  {
1027  offset += it->template count<n>();
1028  p2->write(offset);
1029  }
1030  }
1031  }
1032 
1033  // types
1034  if (n>1)
1035  {
1036  shared_ptr<VTK::DataArrayWriter<unsigned char> > p3
1037  (writer.makeArrayWriter<unsigned char>("types", 1, ncells));
1038  if(!p3->writeIsNoop())
1039  for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
1040  {
1041  int vtktype = VTK::geometryType(it->type());
1042  p3->write(vtktype);
1043  }
1044  }
1045 
1046  writer.endCells();
1047  }
1048 
1049  protected:
1050  // the list of registered functions
1051  std::list<VTKFunctionPtr> celldata;
1052  std::list<VTKFunctionPtr> vertexdata;
1053 
1054  private:
1055  // the grid
1056  GridView gridView_;
1057 
1058  // temporary grid information
1059  protected:
1060  int ncells;
1063  private:
1064  VertexMapper* vertexmapper;
1065  // in conforming mode, for each vertex id (as obtained by vertexmapper)
1066  // hold its number in the iteration order (VertexIterator)
1067  std::vector<int> number;
1068  VTK::DataMode datamode;
1069  protected:
1071  };
1072 
1073 }
1074 
1075 #endif