1 #ifndef DUNE_FEM_VTKIO_HH
2 #define DUNE_FEM_VTKIO_HH
3 
4 #include <type_traits>
5 
6 #include <dune/common/deprecated.hh>
7 
8 #include <dune/grid/io/file/vtk/vtkwriter.hh>
9 #include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
10 
11 #include <dune/fem/version.hh>
12 #include <dune/fem/io/parameter.hh>
13 
14 #include <dune/fem/misc/threads/threadmanager.hh>
15 #include <dune/fem/misc/threads/domainthreaditerator.hh>
16 
17 namespace Dune
18 {
19 
20   namespace Fem
21   {
22 
23     // Internal Forward Declarations
24     // -----------------------------
25 
26     template< class GridPart, bool subsampling = false >
27     class VTKIO;
28 
29 
30 
31     // VTKFunctionWrapper
32     // ------------------
33 
34     template< class DF >
35     class VTKFunctionWrapper
36     : public VTKFunction< typename DF::GridPartType::GridViewType >
37     {
38       typedef VTKFunctionWrapper< DF > ThisType;
39 
40     public:
41       enum TypeOfField { real, complex_real, complex_imag };
42       typedef DF DiscreteFunctionType;
43 
44       typedef ConstLocalFunction< DF > LocalFunctionType;
45       typedef typename DiscreteFunctionType::FunctionSpaceType FunctionSpaceType;
46 
47       static const int dimRange = FunctionSpaceType::dimRange;
48       static const int dimDomain = FunctionSpaceType::dimDomain;
49 
50       typedef typename FunctionSpaceType::DomainType DomainType;
51       typedef typename FunctionSpaceType::RangeType RangeType;
52 
53       typedef typename DiscreteFunctionType::GridPartType GridPartType;
54       typedef typename GridPartType::template Codim< 0 >::EntityType EntityType;
55 
56       typedef typename EntityType::Geometry::LocalCoordinate LocalCoordinateType;
57 
58       //! constructor taking discrete function
VTKFunctionWrapper(const DiscreteFunctionType & df,const std::string & dataName,int component,bool vector,TypeOfField typeOfField)59       VTKFunctionWrapper ( const DiscreteFunctionType& df,
60                            const std::string& dataName,
61                            int component, bool vector, TypeOfField typeOfField )
62       : localFunction_( df ),
63         name_( ( dataName.size() > 0 ) ? dataName : df.name() ),
64         vector_( vector ),
65         typeOfField_( typeOfField ),
66         component_( component )
67       {}
68 
69       //! virtual destructor
~VTKFunctionWrapper()70       virtual ~VTKFunctionWrapper ()
71       {}
72 
73       //! return number of components
ncomps() const74       virtual int ncomps () const
75       {
76         return (!vector_) ? 1 : 3; // dimDomain;
77       }
78 
79       //! evaluate single component comp in
80       //! the entity
evaluate(int comp,const EntityType & e,const LocalCoordinateType & xi) const81       virtual double evaluate ( int comp, const EntityType &e, const LocalCoordinateType &xi ) const
82       {
83         localFunction_.bind( e );
84         typedef typename LocalFunctionType::RangeFieldType RangeFieldType;
85         RangeType val;
86         localFunction_.evaluate(xi,val);
87         localFunction_.unbind();
88 
89         RangeFieldType outVal( 0 );
90         if (vector_)
91         {
92           if( comp <= dimDomain )
93             outVal = val[ comp + component_ ] ;
94         }
95         else
96           outVal = val[ component_ ] ;
97         if (typeOfField_ == TypeOfField::real || typeOfField_ == TypeOfField::complex_real )
98           return std::real( outVal );
99         else
100           return std::imag( outVal );
101       }
102 
103       //! get name
name() const104       virtual std::string name () const
105       {
106         std::stringstream ret;
107         ret << name_;
108         if (typeOfField_ == TypeOfField::complex_real)
109           ret << "_real_";
110         if (typeOfField_ == TypeOfField::complex_imag)
111           ret << "_imag_";
112         if (vector_)
113           ret << "_vec" << component_;
114         else
115           ret << component_;
116         return ret.str();
117       }
118 
119     private:
120       mutable LocalFunctionType localFunction_;
121       const std::string name_ ;
122       const bool vector_;
123       const TypeOfField typeOfField_;
124       const int component_;
125     };
126 
127 
128 
129     // VTKIOBase
130     // ---------
131 
132     //! \brief Output using VTK
133     template< class GridPart, bool subsampling >
134     class VTKIOBase
135     {
136       typedef VTKIOBase< GridPart, subsampling > ThisType;
137 
138     protected:
139       typedef typename GridPart::GridViewType GridViewType;
140 
141       class VTKWriter;
142       class SubsamplingVTKWriter;
143 
144       typedef typename std::conditional< subsampling, SubsamplingVTKWriter, VTKWriter >::type
145         VTKWriterType;
146 
147     public:
148       typedef GridPart GridPartType;
149 
150       typedef typename GridPartType::GridType GridType;
151       typedef typename GridPartType::IndexSetType IndexSetType;
152 
153     protected:
154       class PartitioningData
155         : public VTKFunction< GridViewType >
156       {
157         typedef PartitioningData   ThisType;
158 
159       public:
160         typedef typename GridViewType :: template Codim< 0 >::Entity EntityType;
161         typedef typename EntityType::Geometry::LocalCoordinate LocalCoordinateType;
162 
163         typedef DomainDecomposedIteratorStorage< GridPartType >  ThreadIteratorType;
164 
165         //! constructor taking discrete function
PartitioningData(const GridPartType & gridPart,const std::string & name,const int rank,const int nThreads)166         PartitioningData( const GridPartType& gridPart,
167                           const std::string& name,
168                           const int rank, const int nThreads )
169           : iterators_( gridPart ), name_( name ), rank_( rank ), nThreads_( nThreads ) {}
170 
171         //! virtual destructor
~PartitioningData()172         virtual ~PartitioningData () {}
173 
174         //! return number of components
ncomps() const175         virtual int ncomps () const { return 1; }
176 
177         //! evaluate single component comp in
178         //! the entity
evaluate(int comp,const EntityType & e,const LocalCoordinateType & xi) const179         virtual double evaluate ( int comp, const EntityType &e, const LocalCoordinateType &xi ) const
180         {
181           const int thread = iterators_.thread( e );
182           return (nThreads_ < 0) ? double( rank_ ) : double( rank_ * nThreads_ + thread );
183         }
184 
185         //! get name
name() const186         virtual std::string name () const
187         {
188           return name_;
189         }
190 
191       private:
192         ThreadIteratorType iterators_;
193         const std::string name_;
194         const int rank_;
195         const int nThreads_;
196 
197       };
198 
199       class VolumeData
200         : public VTKFunction< GridViewType >
201       {
202         typedef PartitioningData   ThisType;
203 
204       public:
205         typedef typename GridViewType :: template Codim< 0 >::Entity EntityType;
206         typedef typename EntityType::Geometry::LocalCoordinate LocalCoordinateType;
207 
208         typedef DomainDecomposedIteratorStorage< GridPartType >  ThreadIteratorType;
209 
210         //! constructor taking discrete function
VolumeData()211         VolumeData() {}
212 
213         //! virtual destructor
~VolumeData()214         virtual ~VolumeData () {}
215 
216         //! return number of components
ncomps() const217         virtual int ncomps () const { return 1; }
218 
219         //! evaluate single component comp in
220         //! the entity
evaluate(int comp,const EntityType & e,const LocalCoordinateType & xi) const221         virtual double evaluate ( int comp, const EntityType &e, const LocalCoordinateType &xi ) const
222         {
223           return e.geometry().volume();
224         }
225 
226         //! get name
name() const227         virtual std::string name () const
228         {
229           return std::string("volume");
230         }
231       };
232 
getPartitionParameter(const ParameterReader & parameter=Parameter::container ()) const233       int getPartitionParameter ( const ParameterReader &parameter = Parameter::container() ) const
234       {
235         // 0 = none, 1 = MPI ranks only, 2 = ranks + threads, 3 = like 1 and also threads only
236         const std::string names[] = { "none", "rank", "rank+thread", "rank/thread" };
237         return parameter.getEnum( "fem.io.partitioning", names, 0 );
238       }
239 
240     protected :
VTKIOBase(const GridPartType & gridPart,VTKWriterType * vtkWriter,const ParameterReader & parameter=Parameter::container ())241       VTKIOBase ( const GridPartType &gridPart, VTKWriterType *vtkWriter, const ParameterReader &parameter = Parameter::container() )
242       : gridPart_( gridPart ),
243         vtkWriter_( vtkWriter ),
244         addPartition_( getPartitionParameter( parameter ) )
245       {
246         static const std::string typeTable[] = { "ascii", "base64", "appended-raw", "appended-base64" };
247         static const VTK::OutputType typeValue[] = { VTK::ascii, VTK::base64, VTK::appendedraw, VTK::appendedbase64 };
248         type_ = typeValue[ parameter.getEnum( "fem.io.vtk.type", typeTable, 2 ) ];
249       }
250 
addPartitionData(const int myRank=-1)251       void addPartitionData( const int myRank = -1 )
252       {
253         if( addPartition_ > 0 )
254         {
255           std::shared_ptr<VolumeData> volumePtr( std::make_shared<VolumeData>() );
256           vtkWriter_->addCellData( volumePtr );
257 
258           const int rank = ( myRank < 0 ) ? gridPart_.comm().rank() : myRank ;
259           const int nThreads = ( addPartition_ > 1 ) ? ThreadManager::maxThreads() : 1 ;
260           if( addPartition_ <= 2 )
261           {
262             std::shared_ptr<PartitioningData> dataRankPtr( std::make_shared<PartitioningData>(gridPart_, "rank", rank, nThreads) );
263             vtkWriter_->addCellData( dataRankPtr );
264           }
265           else
266           {
267             // rank only visualization
268             std::shared_ptr<PartitioningData> dataRankPtr( std::make_shared<PartitioningData>(gridPart_, "rank", rank, -1) );
269             vtkWriter_->addCellData( dataRankPtr );
270             // thread only visualization
271             std::shared_ptr<PartitioningData> dataThreadPtr( std::make_shared<PartitioningData>(gridPart_, "thread", 0, nThreads) );
272             vtkWriter_->addCellData( dataThreadPtr );
273           }
274           addPartition_ = 0 ;
275         }
276       }
277 
278       template < class DF >
notComplex()279       static bool notComplex()
280       {
281         typedef typename DF::RangeFieldType RangeFieldType;
282         typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
283         return ! std::is_same< typename std::remove_cv<RangeFieldType>::type, std::complex<RealType> >::value;
284       }
285 
286     public:
~VTKIOBase()287       ~VTKIOBase ()
288       {
289         delete vtkWriter_;
290       }
291 
292       //! return grid part
gridPart() const293       const GridPartType &gridPart () const
294       {
295         return gridPart_;
296       }
297 
298       template< class DF >
addCellData(DF & df,const std::string & dataName="")299       void addCellData( DF &df , const std::string& dataName = "" )
300       {
301         static const int dimRange = DF::FunctionSpaceType::dimRange;
302         for( int i = 0;i < dimRange; ++i )
303         {
304           if ( notComplex<DF>() )
305           {
306             std::shared_ptr<VTKFunctionWrapper< DF > > ptr( std::make_shared<VTKFunctionWrapper< DF > >( df, dataName, i,
307                   false, VTKFunctionWrapper<DF>::TypeOfField::real) );
308             vtkWriter_->addCellData( ptr );
309           }
310           else
311           {
312             std::shared_ptr<VTKFunctionWrapper< DF > > ptrR( std::make_shared<VTKFunctionWrapper< DF > >( df, dataName, i,
313                   false, VTKFunctionWrapper<DF>::TypeOfField::complex_real) );
314             vtkWriter_->addCellData( ptrR );
315             std::shared_ptr<VTKFunctionWrapper< DF > > ptrI( std::make_shared<VTKFunctionWrapper< DF > >( df, dataName, i,
316                   false, VTKFunctionWrapper<DF>::TypeOfField::complex_imag) );
317             vtkWriter_->addCellData( ptrI );
318           }
319         }
320       }
321 
322       template< class DF >
addVectorCellData(DF & df,const std::string & dataName="",int startPoint=0)323       void addVectorCellData( DF &df,
324                               const std::string& dataName = "" ,
325                               int startPoint = 0 )
326       {
327         if ( notComplex<DF>() )
328         {
329           std::shared_ptr<VTKFunctionWrapper< DF > > ptr( std::make_shared<VTKFunctionWrapper< DF > >( df, dataName, startPoint,
330                   true, VTKFunctionWrapper<DF>::TypeOfField::real) );
331           vtkWriter_->addCellData( ptr );
332         }
333         else
334         {
335           std::shared_ptr<VTKFunctionWrapper< DF > > ptrR( std::make_shared<VTKFunctionWrapper< DF > >( df, dataName, startPoint,
336                   true, VTKFunctionWrapper<DF>::TypeOfField::complex_real) );
337           vtkWriter_->addCellData( ptrR );
338           std::shared_ptr<VTKFunctionWrapper< DF > > ptrI( std::make_shared<VTKFunctionWrapper< DF > >( df, dataName, startPoint,
339                   true, VTKFunctionWrapper<DF>::TypeOfField::complex_imag) );
340           vtkWriter_->addCellData( ptrI );
341         }
342       }
343 
344       template< class DF >
addVertexData(DF & df,const std::string & dataName="")345       void addVertexData( DF &df, const std::string& dataName = "" )
346       {
347         static const int dimRange = DF::FunctionSpaceType::dimRange;
348         std::string name = ( dataName.size() > 0 ) ? dataName : df.name() ;
349         for( int i = 0;i < dimRange; ++i )
350         {
351           if ( notComplex<DF>() )
352           {
353             std::shared_ptr<VTKFunctionWrapper< DF > > ptr( std::make_shared<VTKFunctionWrapper< DF > >( df, dataName, i,
354                   false, VTKFunctionWrapper<DF>::TypeOfField::real) );
355             vtkWriter_->addVertexData( ptr );
356           }
357           else
358           {
359             std::shared_ptr<VTKFunctionWrapper< DF > > ptrR( std::make_shared<VTKFunctionWrapper< DF > >( df, dataName, i,
360                   false, VTKFunctionWrapper<DF>::TypeOfField::complex_real) );
361             vtkWriter_->addVertexData( ptrR );
362             std::shared_ptr<VTKFunctionWrapper< DF > > ptrI( std::make_shared<VTKFunctionWrapper< DF > >( df, dataName, i,
363                   false, VTKFunctionWrapper<DF>::TypeOfField::complex_imag) );
364             vtkWriter_->addVertexData( ptrI );
365           }
366         }
367       }
368 
369       template< class DF >
addVectorVertexData(DF & df,const std::string & dataName="",int startPoint=0)370       void addVectorVertexData( DF &df,
371                                 const std::string& dataName = "" ,
372                                 int startPoint = 0 )
373       {
374         if ( notComplex<DF>() )
375         {
376           std::shared_ptr<VTKFunctionWrapper< DF > > ptr( std::make_shared<VTKFunctionWrapper< DF > >( df, dataName, startPoint,
377                   true, VTKFunctionWrapper<DF>::TypeOfField::real) );
378           vtkWriter_->addVertexData( ptr );
379         }
380         else
381         {
382           std::shared_ptr<VTKFunctionWrapper< DF > > ptrR( std::make_shared<VTKFunctionWrapper< DF > >( df, dataName, startPoint,
383                   true, VTKFunctionWrapper<DF>::TypeOfField::complex_real) );
384           vtkWriter_->addVertexData( ptrR );
385           std::shared_ptr<VTKFunctionWrapper< DF > > ptrI( std::make_shared<VTKFunctionWrapper< DF > >( df, dataName, startPoint,
386                   true, VTKFunctionWrapper<DF>::TypeOfField::complex_imag) );
387           vtkWriter_->addVertexData( ptrI );
388         }
389       }
390 
clear()391       void clear ()
392       {
393         vtkWriter_->clear();
394       }
395 
write(const std::string & name,VTK::OutputType type)396       std::string write ( const std::string &name, VTK::OutputType type )
397       {
398         addPartitionData();
399         size_t pos = name.find_last_of( '/' );
400         if( pos != name.npos )
401           return vtkWriter_->pwrite( name.substr( pos+1, name.npos ), name.substr( 0, pos ), "", type );
402         else
403           return vtkWriter_->write( name, type );
404       }
405 
write(const std::string & name)406       std::string write ( const std::string &name )
407       {
408         return write( name, type_ );
409       }
410 
pwrite(const std::string & name,const std::string & path,const std::string & extendpath,VTK::OutputType type)411       std::string pwrite ( const std::string &name,
412                            const std::string &path,
413                            const std::string &extendpath,
414                            VTK::OutputType type )
415       {
416         addPartitionData();
417         return vtkWriter_->pwrite( name, path, extendpath, type );
418       }
419 
pwrite(const std::string & name,const std::string & path,const std::string & extendpath)420       std::string pwrite ( const std::string &name,
421                            const std::string &path,
422                            const std::string &extendpath )
423       {
424         return pwrite( name, path, extendpath, type_ );
425       }
426 
write(const std::string & name,VTK::OutputType type,const int rank,const int size)427       std::string write ( const std::string &name,
428                           VTK::OutputType type,
429                           const int rank,
430                           const int size )
431       {
432         addPartitionData( rank );
433         return vtkWriter_->write( name, type, rank, size );
434       }
435 
write(const std::string & name,const int rank,const int size)436       std::string write ( const std::string &name,
437                           const int rank,
438                           const int size )
439       {
440         return write( name, type_, rank, size );
441       }
442 
443     protected:
444       const GridPartType &gridPart_;
445       VTKWriterType *vtkWriter_;
446       int addPartition_;
447       VTK::OutputType type_;
448     };
449 
450 
451 
452     // VTKIOBase::VTKWriter
453     // --------------------
454 
455     template< class GridPart, bool subsampling >
456     class VTKIOBase< GridPart, subsampling >::VTKWriter
457     : public Dune::VTKWriter< GridViewType >
458     {
459       typedef Dune::VTKWriter< GridViewType > BaseType;
460 
461     public:
462       // make all write methods public for data convert
463       using BaseType::write;
464       using BaseType::pwrite;
465 
466       //! constructor
VTKWriter(const GridPartType & gridPart,VTK::DataMode dm=VTK::conforming)467       VTKWriter( const GridPartType &gridPart,
468                  VTK::DataMode dm = VTK::conforming )
469       : BaseType( gridPart.gridView(), dm )
470       {}
471     };
472 
473 
474 
475     // VTKIOBase::SubSamplingVTKWriter
476     // -------------------------------
477 
478     template< class GridPart, bool subsampling >
479     class VTKIOBase< GridPart, subsampling >::SubsamplingVTKWriter
480       : public Dune::SubsamplingVTKWriter< GridViewType >
481     {
482       typedef Dune::SubsamplingVTKWriter< GridViewType > BaseType;
483 
484     public:
485       // make all write methods public for data convert
486       using BaseType::write;
487       using BaseType::pwrite;
488 
489       //! constructor
SubsamplingVTKWriter(const GridPartType & gridPart,Dune::RefinementIntervals intervals,bool coerceToSimplex=false)490       SubsamplingVTKWriter( const GridPartType &gridPart,
491                             Dune::RefinementIntervals intervals,
492                             bool coerceToSimplex = false )
493         : BaseType( gridPart.gridView(), intervals, coerceToSimplex )
494       {}
495     };
496 
497 
498 
499     // VTKIO (without subsampling)
500     // ---------------------------
501 
502     template< class GridPart >
503     class VTKIO< GridPart, false >
504       : public VTKIOBase< GridPart, false >
505     {
506       typedef VTKIO< GridPart > ThisType;
507       typedef VTKIOBase< GridPart, false > BaseType;
508 
509       typedef typename BaseType::VTKWriterType VTKWriterType;
510 
511     public:
512       typedef typename BaseType::GridPartType GridPartType;
513 
VTKIO(const GridPartType & gridPart,VTK::DataMode dm,const ParameterReader & parameter=Parameter::container ())514       VTKIO ( const GridPartType &gridPart, VTK::DataMode dm, const ParameterReader &parameter = Parameter::container() )
515         : BaseType( gridPart, new VTKWriterType( gridPart, dm ), parameter )
516       {}
517 
VTKIO(const GridPartType & gridPart,const ParameterReader & parameter=Parameter::container ())518       explicit VTKIO ( const GridPartType &gridPart, const ParameterReader &parameter = Parameter::container() )
519         : VTKIO( gridPart, VTK::conforming, parameter )
520       {}
521     };
522 
523 
524 
525     // VTKIO (with subsampling)
526     // ------------------------
527 
528     template< class GridPart >
529     class VTKIO< GridPart, true >
530       : public VTKIOBase< GridPart , true >
531     {
532       typedef VTKIO< GridPart, true > ThisType;
533       typedef VTKIOBase< GridPart, true > BaseType;
534 
535        typedef typename BaseType::VTKWriterType VTKWriterType;
536 
537     public:
538       typedef typename BaseType::GridPartType GridPartType;
539 
VTKIO(const GridPartType & gridPart,Dune::RefinementIntervals intervals,bool coerceToSimplex,const ParameterReader & parameter=Parameter::container ())540       explicit VTKIO ( const GridPartType &gridPart, Dune::RefinementIntervals intervals, bool coerceToSimplex, const ParameterReader &parameter = Parameter::container() )
541         : BaseType( gridPart, new VTKWriterType( gridPart, intervals, coerceToSimplex ), parameter )
542       {}
543 
VTKIO(const GridPartType & gridPart,unsigned int level,bool coerceToSimplex,const ParameterReader & parameter=Parameter::container ())544       explicit VTKIO ( const GridPartType &gridPart, unsigned int level, bool coerceToSimplex, const ParameterReader &parameter = Parameter::container() )
545         : VTKIO( gridPart, Dune::refinementLevels(level), coerceToSimplex, parameter )
546       {}
547 
VTKIO(const GridPartType & gridPart,unsigned int level,const ParameterReader & parameter=Parameter::container ())548       VTKIO ( const GridPartType &gridPart, unsigned int level, const ParameterReader &parameter = Parameter::container() )
549         : VTKIO( gridPart, level, false, parameter )
550       {}
551 
VTKIO(const GridPartType & gridPart,int level,const ParameterReader & parameter=Parameter::container ())552       VTKIO ( const GridPartType &gridPart, int level, const ParameterReader &parameter = Parameter::container() ) DUNE_DEPRECATED_MSG( "pass level as unsigned int" )
553         : VTKIO( gridPart, level, false, parameter )
554       {}
555 
VTKIO(const GridPartType & gridPart,const ParameterReader & parameter=Parameter::container ())556       VTKIO ( const GridPartType &gridPart, const ParameterReader &parameter = Parameter::container() )
557         : VTKIO( gridPart, 0, false, parameter )
558       {}
559 
VTKIO(const GridPartType & gridPart,bool coerceToSimplex,const ParameterReader & parameter=Parameter::container ())560       VTKIO ( const GridPartType &gridPart, bool coerceToSimplex, const ParameterReader &parameter = Parameter::container() )
561         : VTKIO( gridPart, 0, coerceToSimplex, parameter )
562       {}
563     };
564 
565 
566 
567     // SubsamplingVTKIO
568     // ----------------
569 
570     template< class GridPart >
571     using SubsamplingVTKIO = VTKIO< GridPart, true >;
572 
573   } // namespace Fem
574 
575 } // namespace Dune
576 
577 #endif // #ifndef DUNE_FEM_VTKIO_HH
578