← Back to branch summary

~ubuntu-branches/ubuntu/utopic/dune-grid/utopic-proposed

« back to all changes in this revision

Viewing changes to dune/grid/io/file/vtk/vtkwriter.hh

  • Committer: Package Import Robot
  • Author(s): Ansgar Burchardt
  • Date: 2014-02-14 10:49:16 UTC
  • mfrom: (1.3.1) (5.1.4 experimental)
  • Revision ID: package-import@ubuntu.com-20140214104916-2clchnny3eu7ks4w
Tags: 2.3.0-2
Install /usr/share/dune-grid.

Show diffs side-by-side

added added

removed removed

Lines of Context:
dune/grid/io/file/vtk/vtkwriter.hh
14
14
#include <vector>
15
15
#include <list>
16
16
 
17
 
#include <dune/common/deprecated.hh>
18
17
#include <dune/common/exceptions.hh>
19
18
#include <dune/common/indent.hh>
20
19
#include <dune/common/iteratorfacades.hh>
33
32
/** @file
34
33
    @author Peter Bastian, Christian Engwer
35
34
    @brief Provides file i/o for the visualization toolkit
36
 
*/
 
35
 */
37
36
 
38
37
/**
39
38
   @todo put vtk io intro here ...
45
44
 
46
45
   http://www.geophysik.uni-muenchen.de/~moder/Paraview/VTK_File_Formats.php
47
46
   (alternative)
48
 
*/
 
47
 */
49
48
 
50
49
namespace Dune
51
50
{
52
 
  /** 
 
51
  /**
53
52
   * @brief Writer for the ouput of grid functions in the vtk format.
54
53
   * @ingroup VTK
55
54
   *
56
55
   * Writes arbitrary grid functions (living on cells or vertices of a grid)
57
 
   * to a file suitable for easy visualization with 
 
56
   * to a file suitable for easy visualization with
58
57
   * <a href="http://public.kitware.com/VTK/">The Visualization Toolkit (VTK)</a>.
59
58
   */
60
59
  template< class GridView >
68
67
    typedef typename GridView::template Codim< 0 >::Entity Cell;
69
68
    typedef typename GridView::template Codim< n >::Entity Vertex;
70
69
    typedef Cell Entity;
71
 
    
 
70
 
72
71
    typedef typename GridView::IndexSet IndexSet;
73
 
    
 
72
 
74
73
    static const PartitionIteratorType VTK_Partition = InteriorBorder_Partition;
75
 
    
 
74
    //static const PartitionIteratorType VTK_Partition = All_Partition;
 
75
 
76
76
    typedef typename GridView::template Codim< 0 >
77
 
      ::template Partition< VTK_Partition >::Iterator
78
 
      GridCellIterator;
 
77
    ::template Partition< VTK_Partition >::Iterator
 
78
    GridCellIterator;
79
79
    typedef typename GridView::template Codim< n >
80
 
      ::template Partition< VTK_Partition >::Iterator
81
 
      GridVertexIterator;
82
 
    
 
80
    ::template Partition< VTK_Partition >::Iterator
 
81
    GridVertexIterator;
 
82
 
83
83
    typedef MultipleCodimMultipleGeomTypeMapper< GridView, MCMGVertexLayout > VertexMapper;
84
84
 
 
85
    // return true if entity should be skipped in Vertex and Corner iterator
 
86
    static bool skipEntity( const PartitionType entityType )
 
87
    {
 
88
      switch( VTK_Partition )
 
89
      {
 
90
        // for All_Partition no entity has to be skipped
 
91
        case All_Partition:             return false;
 
92
        case InteriorBorder_Partition:  return ( entityType != InteriorEntity );
 
93
        default: DUNE_THROW(NotImplemented,"Add check for this partition type");
 
94
      }
 
95
      return false ;
 
96
    }
 
97
 
85
98
  public:
86
99
    typedef Dune::VTKFunction< GridView > VTKFunction;
87
100
    typedef shared_ptr< const VTKFunction > VTKFunctionPtr;
88
101
 
89
102
  protected:
90
103
    typedef typename std::list<VTKFunctionPtr>::const_iterator FunctionIterator;
91
 
    
 
104
 
92
105
    //! Iterator over the grids elements
93
106
    /**
94
107
     * This class iterates over the gridview's elements.  It is the same as
103
116
      //! get the position of the center of the element, in element-local
104
117
      //! coordinates
105
118
      const FieldVector<DT,n> position() const
106
 
        {
107
 
            return GenericReferenceElements<DT,n>::general((*this)->type()).position(0,0);
108
 
        }
 
119
      {
 
120
        return ReferenceElements<DT,n>::general((*this)->type()).position(0,0);
 
121
      }
109
122
    };
110
123
 
111
124
    CellIterator cellBegin() const
117
130
    {
118
131
      return gridView_.template end< 0, VTK_Partition >();
119
132
    }
120
 
    
 
133
 
121
134
    //! Iterate over the grid's vertices
122
135
    /**
123
136
     * This class iterates over the elements, and within the elements over the
160
173
          cornerIndexDune = 0;
161
174
 
162
175
          ++git;
163
 
          while( (git != gend) && (git->partitionType() != InteriorEntity) )
 
176
          while( (git != gend) && skipEntity( git->partitionType() ) )
164
177
            ++git;
165
178
        }
166
179
      }
172
185
        git(x), gend(end), datamode(dm), cornerIndexDune(0),
173
186
        vertexmapper(vm), visited(vm.size(), false),
174
187
        offset(0)
175
 
        {
176
 
          if (datamode == VTK::conforming && git != gend)
177
 
            visited[vertexmapper.map(*git,cornerIndexDune,n)] = true;
178
 
        }
 
188
      {
 
189
        if (datamode == VTK::conforming && git != gend)
 
190
          visited[vertexmapper.map(*git,cornerIndexDune,n)] = true;
 
191
      }
179
192
      void increment ()
 
193
      {
 
194
        switch (datamode)
180
195
        {
181
 
          switch (datamode)
 
196
        case VTK::conforming :
 
197
          while(visited[vertexmapper.map(*git,cornerIndexDune,n)])
182
198
          {
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
199
            basicIncrement();
193
 
            break;
 
200
            if (git == gend) return;
194
201
          }
195
 
       }
 
202
          visited[vertexmapper.map(*git,cornerIndexDune,n)] = true;
 
203
          break;
 
204
        case VTK::nonconforming :
 
205
          basicIncrement();
 
206
          break;
 
207
        }
 
208
      }
196
209
      bool equals (const VertexIterator & cit) const
197
 
        {
198
 
          return git == cit.git
199
 
            && cornerIndexDune == cit.cornerIndexDune
200
 
            && datamode == cit.datamode;
201
 
        }
 
210
      {
 
211
        return git == cit.git
 
212
               && cornerIndexDune == cit.cornerIndexDune
 
213
               && datamode == cit.datamode;
 
214
      }
202
215
      const Entity& dereference() const
203
 
        {
204
 
          return *git;
205
 
        }
 
216
      {
 
217
        return *git;
 
218
      }
206
219
      //! index of vertex within the entity, in Dune-numbering
207
220
      int localindex () const
208
 
        {
209
 
          return cornerIndexDune;
210
 
        }
 
221
      {
 
222
        return cornerIndexDune;
 
223
      }
211
224
      //! position of vertex inside the entity
212
225
      const FieldVector<DT,n> & position () const
213
 
        {
214
 
          return GenericReferenceElements<DT,n>::general(git->type())
215
 
            .position(cornerIndexDune,n);
216
 
        }
217
 
    };    
218
 
    
 
226
      {
 
227
        return ReferenceElements<DT,n>::general(git->type())
 
228
               .position(cornerIndexDune,n);
 
229
      }
 
230
    };
 
231
 
219
232
    VertexIterator vertexBegin () const
220
233
    {
221
234
      return VertexIterator( gridView_.template begin< 0, VTK_Partition >(),
229
242
                             gridView_.template end< 0, VTK_Partition >(),
230
243
                             datamode, *vertexmapper );
231
244
    }
232
 
    
 
245
 
233
246
    //! Iterate over the elements' corners
234
247
    /**
235
248
     * This class iterates over the elements, and within the elements over the
284
297
          cornerIndexVTK = 0;
285
298
 
286
299
          ++git;
287
 
          while( (git != gend) && (git->partitionType() != InteriorEntity) )
 
300
          while( (git != gend) && skipEntity( git->partitionType() ) )
288
301
            ++git;
289
302
        }
290
303
      }
291
304
      bool equals (const CornerIterator & cit) const
292
 
        {
293
 
          return git == cit.git
294
 
            && cornerIndexVTK == cit.cornerIndexVTK
295
 
            && datamode == cit.datamode;
296
 
        }
 
305
      {
 
306
        return git == cit.git
 
307
               && cornerIndexVTK == cit.cornerIndexVTK
 
308
               && datamode == cit.datamode;
 
309
      }
297
310
      const Entity& dereference() const
298
 
        {
299
 
          return *git;
300
 
        }
 
311
      {
 
312
        return *git;
 
313
      }
301
314
      //! Process-local consecutive zero-starting vertex id
302
315
      /**
303
316
       * This method returns the number of this corners associated vertex, in
304
317
       * the numbering given by the iteration order of VertexIterator.
305
318
       */
306
319
      int id () const
 
320
      {
 
321
        switch (datamode)
307
322
        {
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
 
          }
 
323
        case VTK::conforming :
 
324
          return
 
325
            number[vertexmapper.map(*git,VTK::renumber(*git,cornerIndexVTK),
 
326
                                    n)];
 
327
        case VTK::nonconforming :
 
328
          return offset + VTK::renumber(*git,cornerIndexVTK);
 
329
        default :
 
330
          DUNE_THROW(IOError,"VTKWriter: unsupported DataMode" << datamode);
319
331
        }
 
332
      }
320
333
    };
321
 
    
 
334
 
322
335
    CornerIterator cornerBegin () const
323
336
    {
324
337
      return CornerIterator( gridView_.template begin< 0, VTK_Partition >(),
325
338
                             gridView_.template end< 0, VTK_Partition >(),
326
339
                             datamode, *vertexmapper, number );
327
340
    }
328
 
    
 
341
 
329
342
    CornerIterator cornerEnd () const
330
343
    {
331
344
      return CornerIterator( gridView_.template end< 0, VTK_Partition >(),
332
345
                             gridView_.template end< 0, VTK_Partition >(),
333
346
                             datamode, *vertexmapper, number );
334
 
    }    
 
347
    }
335
348
 
336
349
  public:
337
350
    /**
338
351
     * @brief Construct a VTKWriter working on a specific GridView.
339
 
     * 
340
 
     * 
 
352
     *
 
353
     *
341
354
     * @param gridView The gridView the grid functions live on. (E. g. a LevelGridView.)
342
355
     * @param dm The data mode.
343
356
     */
344
357
    explicit VTKWriter ( const GridView &gridView,
345
358
                         VTK::DataMode dm = VTK::conforming )
346
 
    : gridView_( gridView ),
347
 
      datamode( dm )
 
359
      : gridView_( gridView ),
 
360
        datamode( dm )
348
361
    { }
349
362
 
350
363
    /**
352
365
     * @param p Dune::shared_ptr to the function to visualize
353
366
     */
354
367
    void addCellData (const VTKFunctionPtr & p)
355
 
      {
356
 
        celldata.push_back(p);
357
 
      }
 
368
    {
 
369
      celldata.push_back(p);
 
370
    }
358
371
 
359
372
    /**
360
373
     * @brief Add a grid function that lives on the cells of the grid to the visualization.
361
374
     * @param p The function to visualize.  The VTKWriter object will take
362
375
     *          ownership of the VTKFunction *p and delete it when it's done.
363
376
     */
364
 
      void addCellData (VTKFunction* p) // DUNE_DEPRECATED
365
 
      {
366
 
        celldata.push_back(VTKFunctionPtr(p));
367
 
      }
 
377
    void addCellData (VTKFunction* p)
 
378
    {
 
379
      celldata.push_back(VTKFunctionPtr(p));
 
380
    }
368
381
 
369
382
    /**
370
 
     * @brief Add a grid function (represented by container) that lives on the cells of 
 
383
     * @brief Add a grid function (represented by container) that lives on the cells of
371
384
     * the grid to the visualization.
372
385
     *
373
 
     * The container has to have random access via operator[] (e. g. std::vector). The 
 
386
     * The container has to have random access via operator[] (e. g. std::vector). The
374
387
     * value of the grid function for an arbitrary element
375
388
     * will be accessed by calling operator[] with the index (corresponding
376
 
     * with the grid view) of the element. 
 
389
     * with the grid view) of the element.
377
390
     * For vector valued data all components for an element are assumed to
378
391
     * be consecutive.
379
392
     *
385
398
    void addCellData (const V& v, const std::string &name, int ncomps = 1)
386
399
    {
387
400
      typedef P0VTKFunction<GridView, V> Function;
388
 
      for (int c=0;c<ncomps;++c) {
 
401
      for (int c=0; c<ncomps; ++c) {
389
402
        std::stringstream compName;
390
403
        compName << name;
391
404
        if (ncomps>1)
400
413
     * @param p The function to visualize.  The VTKWriter object will take
401
414
     *          ownership of the VTKFunction *p and delete it when it's done.
402
415
     */
403
 
      void addVertexData (VTKFunction* p) // DUNE_DEPRECATED
404
 
      {
405
 
        vertexdata.push_back(VTKFunctionPtr(p));
406
 
      }
 
416
    void addVertexData (VTKFunction* p)
 
417
    {
 
418
      vertexdata.push_back(VTKFunctionPtr(p));
 
419
    }
407
420
 
408
421
    /**
409
422
     * @brief Add a grid function that lives on the vertices of the grid to the visualization.
410
423
     * @param p Dune::shared_ptr to the function to visualize
411
424
     */
412
425
    void addVertexData (const VTKFunctionPtr & p)
413
 
      {
414
 
        vertexdata.push_back(p);
415
 
      }
 
426
    {
 
427
      vertexdata.push_back(p);
 
428
    }
416
429
 
417
430
    /**
418
 
     * @brief Add a grid function (represented by container) that lives on the vertices of the 
 
431
     * @brief Add a grid function (represented by container) that lives on the vertices of the
419
432
     * grid to the visualization output.
420
433
     *
421
 
     * The container has to have random access via operator[] (e. g. std::vector). The value 
 
434
     * The container has to have random access via operator[] (e. g. std::vector). The value
422
435
     * of the grid function for an arbitrary element
423
436
     * will be accessed by calling operator[] with the index (corresponding
424
437
     * to the grid view) of the vertex.
433
446
    void addVertexData (const V& v, const std::string &name, int ncomps=1)
434
447
    {
435
448
      typedef P1VTKFunction<GridView, V> Function;
436
 
      for (int c=0;c<ncomps;++c) {
 
449
      for (int c=0; c<ncomps; ++c) {
437
450
        std::stringstream compName;
438
451
        compName << name;
439
452
        if (ncomps>1)
445
458
 
446
459
    //! clear list of registered functions
447
460
    void clear ()
448
 
      {
449
 
        celldata.clear();
450
 
        vertexdata.clear();
451
 
      }
 
461
    {
 
462
      celldata.clear();
 
463
      vertexdata.clear();
 
464
    }
452
465
 
453
466
    //! destructor
454
467
    virtual ~VTKWriter ()
455
 
      {
456
 
        this->clear();
457
 
      }
 
468
    {
 
469
      this->clear();
 
470
    }
458
471
 
459
472
    /** \brief write output (interface might change later)
460
473
     *
623
636
 
624
637
      // write process data
625
638
      std::ofstream file;
 
639
      file.exceptions(std::ios_base::badbit | std::ios_base::failbit |
 
640
                      std::ios_base::eofbit);
626
641
      file.open( pieceName.c_str(), std::ios::binary );
627
642
      if (! file.is_open())
628
643
        DUNE_THROW(IOError, "Could not write to piece file " << pieceName);
660
675
                       const std::string& extendpath,
661
676
                       VTK::OutputType ot, const int commRank,
662
677
                       const int commSize )
 
678
    {
 
679
      // make data mode visible to private functions
 
680
      outputtype=ot;
 
681
 
 
682
      // do some magic because paraview can only cope with relative pathes to piece files
 
683
      std::ofstream file;
 
684
      file.exceptions(std::ios_base::badbit | std::ios_base::failbit |
 
685
                      std::ios_base::eofbit);
 
686
      std::string piecepath = concatPaths(path, extendpath);
 
687
      std::string relpiecepath = relativePath(path, piecepath);
 
688
 
 
689
      // write this processes .vtu/.vtp piece file
 
690
      std::string fullname = getParallelPieceName(name, piecepath, commRank,
 
691
                                                  commSize);
 
692
      file.open(fullname.c_str(),std::ios::binary);
 
693
      if (! file.is_open())
 
694
        DUNE_THROW(IOError, "Could not write to piecefile file " << fullname);
 
695
      writeDataFile(file);
 
696
      file.close();
 
697
      gridView_.comm().barrier();
 
698
 
 
699
      // if we are rank 0, write .pvtu/.pvtp parallel header
 
700
      fullname = getParallelHeaderName(name, path, commSize);
 
701
      if( commRank  ==0 )
663
702
      {
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);
 
703
        file.open(fullname.c_str());
676
704
        if (! file.is_open())
677
 
          DUNE_THROW(IOError, "Could not write to piecefile file " << fullname);
678
 
        writeDataFile(file);
 
705
          DUNE_THROW(IOError, "Could not write to parallel file " << fullname);
 
706
        writeParallelHeader(file,name,relpiecepath, commSize );
679
707
        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
708
      }
 
709
      gridView_.comm().barrier();
 
710
      return fullname;
 
711
    }
695
712
 
696
713
  private:
697
714
    //! write header file in parallel case to stream
714
731
     */
715
732
    void writeParallelHeader(std::ostream& s, const std::string& piecename,
716
733
                             const std::string& piecepath, const int commSize)
717
 
      {
718
 
        VTK::FileType fileType =
719
 
          (n == 1) ? VTK::polyData : VTK::unstructuredGrid;
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
 
      }
 
734
    {
 
735
      VTK::FileType fileType =
 
736
        (n == 1) ? VTK::polyData : VTK::unstructuredGrid;
 
737
 
 
738
      VTK::PVTUWriter writer(s, fileType);
 
739
 
 
740
      writer.beginMain();
 
741
 
 
742
      // PPointData
 
743
      {
 
744
        std::string scalars;
 
745
        for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end();
 
746
             ++it)
 
747
          if ((*it)->ncomps()==1)
 
748
          {
 
749
            scalars = (*it)->name();
 
750
            break;
 
751
          }
 
752
        std::string vectors;
 
753
        for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end();
 
754
             ++it)
 
755
          if ((*it)->ncomps()>1)
 
756
          {
 
757
            vectors = (*it)->name();
 
758
            break;
 
759
          }
 
760
        writer.beginPointData(scalars, vectors);
 
761
      }
 
762
      for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end();
 
763
           ++it)
 
764
      {
 
765
        unsigned writecomps = (*it)->ncomps();
 
766
        if(writecomps == 2) writecomps = 3;
 
767
        writer.addArray<float>((*it)->name(), writecomps);
 
768
      }
 
769
      writer.endPointData();
 
770
 
 
771
      // PCellData
 
772
      {
 
773
        std::string scalars;
 
774
        for (FunctionIterator it=celldata.begin(); it!=celldata.end();
 
775
             ++it)
 
776
          if ((*it)->ncomps()==1)
 
777
          {
 
778
            scalars = (*it)->name();
 
779
            break;
 
780
          }
 
781
        std::string vectors;
 
782
        for (FunctionIterator it=celldata.begin(); it!=celldata.end();
 
783
             ++it)
 
784
          if ((*it)->ncomps()>1)
 
785
          {
 
786
            vectors = (*it)->name();
 
787
            break;
 
788
          }
 
789
        writer.beginCellData(scalars, vectors);
 
790
      }
 
791
      for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it) {
 
792
        unsigned writecomps = (*it)->ncomps();
 
793
        if(writecomps == 2) writecomps = 3;
 
794
        writer.addArray<float>((*it)->name(), writecomps);
 
795
      }
 
796
      writer.endCellData();
 
797
 
 
798
      // PPoints
 
799
      writer.beginPoints();
 
800
      writer.addArray<float>("Coordinates", 3);
 
801
      writer.endPoints();
 
802
 
 
803
      // Pieces
 
804
      for( int i = 0; i < commSize; ++i )
 
805
      {
 
806
        const std::string& fullname = getParallelPieceName(piecename,
 
807
                                                           piecepath, i,
 
808
                                                           commSize);
 
809
        writer.addPiece(fullname);
 
810
      }
 
811
 
 
812
      writer.endMain();
 
813
    }
797
814
 
798
815
    //! write data file to stream
799
816
    void writeDataFile (std::ostream& s)
 
817
    {
 
818
      VTK::FileType fileType =
 
819
        (n == 1) ? VTK::polyData : VTK::unstructuredGrid;
 
820
 
 
821
      VTK::VTUWriter writer(s, outputtype, fileType);
 
822
 
 
823
      // Grid characteristics
 
824
      vertexmapper = new VertexMapper( gridView_ );
 
825
      if (datamode == VTK::conforming)
800
826
      {
801
 
        VTK::FileType fileType =
802
 
          (n == 1) ? VTK::polyData : VTK::unstructuredGrid;
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
 
        }
813
 
        countEntities(nvertices, ncells, ncorners);
814
 
 
815
 
        writer.beginMain(ncells, nvertices);
 
827
        number.resize(vertexmapper->size());
 
828
        for (std::vector<int>::size_type i=0; i<number.size(); i++) number[i] = -1;
 
829
      }
 
830
      countEntities(nvertices, ncells, ncorners);
 
831
 
 
832
      writer.beginMain(ncells, nvertices);
 
833
      writeAllData(writer);
 
834
      writer.endMain();
 
835
 
 
836
      // write appended binary data section
 
837
      if(writer.beginAppended())
816
838
        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
 
      }
 
839
      writer.endAppended();
 
840
 
 
841
      delete vertexmapper; number.clear();
 
842
    }
826
843
 
827
844
    void writeAllData(VTK::VTUWriter& writer) {
828
845
      // PointData
840
857
 
841
858
  protected:
842
859
    std::string getFormatString() const
843
 
      {
844
 
          if (outputtype==VTK::ascii)
845
 
              return "ascii";
846
 
          if (outputtype==VTK::base64)
847
 
              return "binary";
848
 
          if (outputtype==VTK::appendedraw)
849
 
              return "appended";
850
 
          if (outputtype==VTK::appendedbase64)
851
 
              return "appended";
852
 
          DUNE_THROW(IOError, "VTKWriter: unsupported OutputType" << outputtype);
853
 
      }
854
 
      
 
860
    {
 
861
      if (outputtype==VTK::ascii)
 
862
        return "ascii";
 
863
      if (outputtype==VTK::base64)
 
864
        return "binary";
 
865
      if (outputtype==VTK::appendedraw)
 
866
        return "appended";
 
867
      if (outputtype==VTK::appendedbase64)
 
868
        return "appended";
 
869
      DUNE_THROW(IOError, "VTKWriter: unsupported OutputType" << outputtype);
 
870
    }
 
871
 
855
872
    std::string getTypeString() const
856
 
      {
857
 
          if (n==1)
858
 
              return "PolyData";
859
 
          else
860
 
              return "UnstructuredGrid";
861
 
      }
862
 
      
 
873
    {
 
874
      if (n==1)
 
875
        return "PolyData";
 
876
      else
 
877
        return "UnstructuredGrid";
 
878
    }
 
879
 
863
880
    //! count the vertices, cells and corners
864
881
    virtual void countEntities(int &nvertices, int &ncells, int &ncorners)
 
882
    {
 
883
      nvertices = 0;
 
884
      ncells = 0;
 
885
      ncorners = 0;
 
886
      for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
865
887
      {
866
 
        nvertices = 0;
867
 
        ncells = 0;
868
 
        ncorners = 0;
869
 
        for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
 
888
        ncells++;
 
889
        // because of the use of vertexmapper->map(), this iteration must be
 
890
        // in the order of Dune's numbering.
 
891
        for (int i=0; i<it->template count<n>(); ++i)
870
892
        {
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
 
            }
 
893
          ncorners++;
 
894
          if (datamode == VTK::conforming)
 
895
          {
 
896
            int alpha = vertexmapper->map(*it,i,n);
 
897
            if (number[alpha]<0)
 
898
              number[alpha] = nvertices++;
 
899
          }
 
900
          else
 
901
          {
 
902
            nvertices++;
887
903
          }
888
904
        }
889
905
      }
 
906
    }
890
907
 
891
908
    //! write cell data
892
909
    virtual void writeCellData(VTK::VTUWriter& writer)
 
910
    {
 
911
      if(celldata.size() == 0)
 
912
        return;
 
913
 
 
914
      std::string scalars = "";
 
915
      for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it)
 
916
        if ((*it)->ncomps()==1)
 
917
        {
 
918
          scalars = (*it)->name();
 
919
          break;
 
920
        }
 
921
      std::string vectors = "";
 
922
      for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it)
 
923
        if ((*it)->ncomps()>1)
 
924
        {
 
925
          vectors = (*it)->name();
 
926
          break;
 
927
        }
 
928
 
 
929
      writer.beginCellData(scalars, vectors);
 
930
      for (FunctionIterator it=celldata.begin(); it!=celldata.end(); ++it)
893
931
      {
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();
 
932
        // vtk file format: a vector data always should have 3 comps (with
 
933
        // 3rd comp = 0 in 2D case)
 
934
        unsigned writecomps = (*it)->ncomps();
 
935
        if(writecomps == 2) writecomps = 3;
 
936
        shared_ptr<VTK::DataArrayWriter<float> > p
 
937
          (writer.makeArrayWriter<float>((*it)->name(), writecomps,
 
938
                                         ncells));
 
939
        if(!p->writeIsNoop())
 
940
          for (CellIterator i=cellBegin(); i!=cellEnd(); ++i)
 
941
          {
 
942
            for (int j=0; j<(*it)->ncomps(); j++)
 
943
              p->write((*it)->evaluate(j,*i,i.position()));
 
944
            // vtk file format: a vector data always should have 3 comps
 
945
            // (with 3rd comp = 0 in 2D case)
 
946
            for (unsigned j=(*it)->ncomps(); j < writecomps; ++j)
 
947
              p->write(0.0);
 
948
          }
934
949
      }
 
950
      writer.endCellData();
 
951
    }
935
952
 
936
953
    //! write vertex data
937
954
    virtual void writeVertexData(VTK::VTUWriter& writer)
 
955
    {
 
956
      if(vertexdata.size() == 0)
 
957
        return;
 
958
 
 
959
      std::string scalars = "";
 
960
      for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it)
 
961
        if ((*it)->ncomps()==1)
 
962
        {
 
963
          scalars = (*it)->name();
 
964
          break;
 
965
        }
 
966
      std::string vectors = "";
 
967
      for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it)
 
968
        if ((*it)->ncomps()>1)
 
969
        {
 
970
          vectors = (*it)->name();
 
971
          break;
 
972
        }
 
973
 
 
974
      writer.beginPointData(scalars, vectors);
 
975
      for (FunctionIterator it=vertexdata.begin(); it!=vertexdata.end(); ++it)
938
976
      {
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();
 
977
        // vtk file format: a vector data always should have 3 comps (with
 
978
        // 3rd comp = 0 in 2D case)
 
979
        unsigned writecomps = (*it)->ncomps();
 
980
        if(writecomps == 2) writecomps = 3;
 
981
        shared_ptr<VTK::DataArrayWriter<float> > p
 
982
          (writer.makeArrayWriter<float>((*it)->name(), writecomps,
 
983
                                         nvertices));
 
984
        if(!p->writeIsNoop())
 
985
          for (VertexIterator vit=vertexBegin(); vit!=vertexEnd(); ++vit)
 
986
          {
 
987
            for (int j=0; j<(*it)->ncomps(); j++)
 
988
              p->write((*it)->evaluate(j,*vit,vit.position()));
 
989
            // vtk file format: a vector data always should have 3 comps
 
990
            // (with 3rd comp = 0 in 2D case)
 
991
            for (unsigned j=(*it)->ncomps(); j < writecomps; ++j)
 
992
              p->write(0.0);
 
993
          }
979
994
      }
 
995
      writer.endPointData();
 
996
    }
980
997
 
981
998
    //! write the positions of vertices
982
999
    virtual void writeGridPoints(VTK::VTUWriter& writer)
983
 
      {
984
 
        writer.beginPoints();
 
1000
    {
 
1001
      writer.beginPoints();
985
1002
 
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
 
          }
 
1003
      shared_ptr<VTK::DataArrayWriter<float> > p
 
1004
        (writer.makeArrayWriter<float>("Coordinates", 3, nvertices));
 
1005
      if(!p->writeIsNoop()) {
 
1006
        VertexIterator vEnd = vertexEnd();
 
1007
        for (VertexIterator vit=vertexBegin(); vit!=vEnd; ++vit)
 
1008
        {
 
1009
          int dimw=w;
 
1010
          for (int j=0; j<std::min(dimw,3); j++)
 
1011
            p->write(vit->geometry().corner(vit.localindex())[j]);
 
1012
          for (int j=std::min(dimw,3); j<3; j++)
 
1013
            p->write(0.0);
998
1014
        }
999
 
        // free the VTK::DataArrayWriter before touching the stream
1000
 
        p.reset();
1001
 
 
1002
 
        writer.endPoints();
1003
1015
      }
 
1016
      // free the VTK::DataArrayWriter before touching the stream
 
1017
      p.reset();
 
1018
 
 
1019
      writer.endPoints();
 
1020
    }
1004
1021
 
1005
1022
    //! write the connectivity array
1006
1023
    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
 
      }
 
1024
    {
 
1025
      writer.beginCells();
 
1026
 
 
1027
      // connectivity
 
1028
      {
 
1029
        shared_ptr<VTK::DataArrayWriter<int> > p1
 
1030
          (writer.makeArrayWriter<int>("connectivity", 1, ncorners));
 
1031
        if(!p1->writeIsNoop())
 
1032
          for (CornerIterator it=cornerBegin(); it!=cornerEnd(); ++it)
 
1033
            p1->write(it.id());
 
1034
      }
 
1035
 
 
1036
      // offsets
 
1037
      {
 
1038
        shared_ptr<VTK::DataArrayWriter<int> > p2
 
1039
          (writer.makeArrayWriter<int>("offsets", 1, ncells));
 
1040
        if(!p2->writeIsNoop()) {
 
1041
          int offset = 0;
 
1042
          for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
 
1043
          {
 
1044
            offset += it->template count<n>();
 
1045
            p2->write(offset);
 
1046
          }
 
1047
        }
 
1048
      }
 
1049
 
 
1050
      // types
 
1051
      if (n>1)
 
1052
      {
 
1053
        shared_ptr<VTK::DataArrayWriter<unsigned char> > p3
 
1054
          (writer.makeArrayWriter<unsigned char>("types", 1, ncells));
 
1055
        if(!p3->writeIsNoop())
 
1056
          for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
 
1057
          {
 
1058
            int vtktype = VTK::geometryType(it->type());
 
1059
            p3->write(vtktype);
 
1060
          }
 
1061
      }
 
1062
 
 
1063
      writer.endCells();
 
1064
    }
1048
1065
 
1049
1066
  protected:
1050
1067
    // the list of registered functions
1051
1068
    std::list<VTKFunctionPtr> celldata;
1052
1069
    std::list<VTKFunctionPtr> vertexdata;
1053
1070
 
1054
 
  private: 
1055
1071
    // the grid
1056
1072
    GridView gridView_;
1057
1073
 
1058
1074
    // temporary grid information
1059
 
  protected:
1060
1075
    int ncells;
1061
1076
    int nvertices;
1062
1077
    int ncorners;