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 ¶meter = 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 ¶meter = 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 ¶meter = 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 ¶meter = 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 ¶meter = 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 ¶meter = 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 ¶meter = 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 ¶meter = 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 ¶meter = 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 ¶meter = 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