Dune Core Modules (2.6.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#include <memory>
14
15#include <vector>
16#include <list>
17#include <map>
18
21#include <dune/common/std/memory.hh>
22#include <dune/common/indent.hh>
24#include <dune/common/path.hh>
25#include <dune/geometry/referenceelements.hh>
27#include <dune/grid/common/gridenums.hh>
31#include <dune/grid/io/file/vtk/pvtuwriter.hh>
32#include <dune/grid/io/file/vtk/streams.hh>
33#include <dune/grid/io/file/vtk/vtuwriter.hh>
34
48namespace Dune
49{
50
51 namespace Impl
52 {
53
54 template< class F, class E, class = void >
55 struct IsBindable
56 : std::false_type
57 {};
58
59 template< class F, class E >
60 struct IsBindable< F, E, void_t< decltype( std::declval< F & >().bind( std::declval< const E & >() ) ), decltype( std::declval< F & >().unbind() ) > >
61 : std::true_type
62 {};
63
64 } // namespace Impl
65
66 namespace VTKWriteTypeTraits {
67 template<typename T>
68 struct IsLocalFunction
69 {
70 };
71 }
72
73 // Forward-declaration here, so the class can be friend of VTKWriter
74 template <class GridView>
75 class VTKSequenceWriterBase;
76 template <class GridView>
77 class VTKSequenceWriter;
78
87 template< class GridView >
88 class VTKWriter {
89
90 // VTKSequenceWriterBase needs getSerialPieceName
91 // and getParallelHeaderName
92 friend class VTKSequenceWriterBase<GridView>;
93 // VTKSequenceWriter needs the grid view, to get the MPI size and rank
94 friend class VTKSequenceWriter<GridView>;
95
96 // extract types
97 typedef typename GridView::Grid Grid;
98 typedef typename GridView::ctype DT;
99 enum { n = GridView::dimension };
100 enum { w = GridView::dimensionworld };
101
102 typedef typename GridView::template Codim< 0 >::Entity Cell;
103 typedef typename GridView::template Codim< n >::Entity Vertex;
104 typedef Cell Entity;
105
106 typedef typename GridView::IndexSet IndexSet;
107
108 static const PartitionIteratorType VTK_Partition = InteriorBorder_Partition;
109 //static const PartitionIteratorType VTK_Partition = All_Partition;
110
111 typedef typename GridView::template Codim< 0 >
112 ::template Partition< VTK_Partition >::Iterator
113 GridCellIterator;
114 typedef typename GridView::template Codim< n >
115 ::template Partition< VTK_Partition >::Iterator
116 GridVertexIterator;
117
118 typedef typename GridCellIterator::Reference EntityReference;
119
120 typedef typename GridView::template Codim< 0 >
121 ::Entity::Geometry::LocalCoordinate Coordinate;
122
124
125 // return true if entity should be skipped in Vertex and Corner iterator
126 static bool skipEntity( const PartitionType entityType )
127 {
128 switch( VTK_Partition )
129 {
130 // for All_Partition no entity has to be skipped
131 case All_Partition: return false;
132 case InteriorBorder_Partition: return ( entityType != InteriorEntity );
133 default: DUNE_THROW(NotImplemented,"Add check for this partition type");
134 }
135 return false ;
136 }
137
138 public:
139
141
142 protected:
143
145
149 {
150
151 public:
152
154
157 {
158
160 virtual void bind(const Entity& e) = 0;
161
163 virtual void unbind() = 0;
164
166
169 virtual void write(const Coordinate& pos, Writer& w, std::size_t count) const = 0;
170
171 virtual ~FunctionWrapperBase()
172 {}
173
174 };
175
177 template<typename F>
179 : public FunctionWrapperBase
180 {
181 using Function = typename std::decay<F>::type;
182
183 template<typename F_>
184 FunctionWrapper(F_&& f)
185 : _f(std::forward<F_>(f))
186 {}
187
188 virtual void bind(const Entity& e)
189 {
190 _f.bind(e);
191 }
192
193 virtual void unbind()
194 {
195 _f.unbind();
196 }
197
198 virtual void write(const Coordinate& pos, Writer& w, std::size_t count) const
199 {
200 auto r = _f(pos);
201 // we need to do different things here depending on whether r supports indexing into it or not.
202 do_write(w,r,count,is_indexable<decltype(r)>());
203 }
204
205 private:
206
207 template<typename R>
208 void do_write(Writer& w, const R& r, std::size_t count, std::true_type) const
209 {
210 for (std::size_t i = 0; i < count; ++i)
211 w.write(r[i]);
212 }
213
214 template<typename R>
215 void do_write(Writer& w, const R& r, std::size_t count, std::false_type) const
216 {
217 assert(count == 1);
218 w.write(r);
219 }
220
221 Function _f;
222 };
223
226 : public FunctionWrapperBase
227 {
228 VTKFunctionWrapper(const std::shared_ptr< const VTKFunction >& f)
229 : _f(f)
230 , _entity(nullptr)
231 {}
232
233 virtual void bind(const Entity& e)
234 {
235 _entity = &e;
236 }
237
238 virtual void unbind()
239 {
240 _entity = nullptr;
241 }
242
243 virtual void write(const Coordinate& pos, Writer& w, std::size_t count) const
244 {
245 for (std::size_t i = 0; i < count; ++i)
246 w.write(_f->evaluate(i,*_entity,pos));
247 }
248
249 private:
250
251 std::shared_ptr< const VTKFunction > _f;
252 const Entity* _entity;
253
254 };
255
257 template<typename F, std::enable_if_t<Impl::IsBindable<F, Entity>::value, int> = 0>
259 : _f(Dune::Std::make_unique<FunctionWrapper<F> >(std::forward<F>(f)))
260 , _fieldInfo(fieldInfo)
261 {}
262
264 template<typename F, std::enable_if_t<not Impl::IsBindable<F, Entity>::value, int> = 0>
266 : _f(Dune::Std::make_unique< FunctionWrapper<
267 typename std::decay<decltype(localFunction(std::forward<F>(f)))>::type
268 > >(localFunction(std::forward<F>(f))))
269 , _fieldInfo(fieldInfo)
270 {}
271
273 explicit VTKLocalFunction (const std::shared_ptr< const VTKFunction >& vtkFunctionPtr)
274 : _f(Dune::Std::make_unique<VTKFunctionWrapper>(vtkFunctionPtr))
275 , _fieldInfo(
276 vtkFunctionPtr->name(),
277 vtkFunctionPtr->ncomps() > 1 ? VTK::FieldInfo::Type::vector : VTK::FieldInfo::Type::scalar,
278 vtkFunctionPtr->ncomps()
279 )
280 {}
281
283 std::string name() const
284 {
285 return fieldInfo().name();
286 }
287
290 {
291 return _fieldInfo;
292 }
293
295 void bind(const Entity& e) const
296 {
297 _f->bind(e);
298 }
299
301 void unbind() const
302 {
303 _f->unbind();
304 }
305
307 void write(const Coordinate& pos, Writer& w) const
308 {
309 _f->write(pos,w,fieldInfo().size());
310 }
311
312 std::shared_ptr<FunctionWrapperBase> _f;
313 VTK::FieldInfo _fieldInfo;
314
315 };
316
317 typedef typename std::list<VTKLocalFunction>::const_iterator FunctionIterator;
318
320
325 class CellIterator : public GridCellIterator
326 {
327 public:
329 CellIterator(const GridCellIterator & x) : GridCellIterator(x) {}
333 {
334 return ReferenceElements<DT,n>::general((*this)->type()).position(0,0);
335 }
336 };
337
338 CellIterator cellBegin() const
339 {
340 return gridView_.template begin< 0, VTK_Partition >();
341 }
342
343 CellIterator cellEnd() const
344 {
345 return gridView_.template end< 0, VTK_Partition >();
346 }
347
349
364 public ForwardIteratorFacade<VertexIterator, const Entity, EntityReference, int>
365 {
366 GridCellIterator git;
367 GridCellIterator gend;
368 VTK::DataMode datamode;
369 // Index of the currently visited corner within the current element.
370 // NOTE: this is in Dune-numbering, in contrast to CornerIterator.
371 int cornerIndexDune;
372 const VertexMapper & vertexmapper;
373 std::vector<bool> visited;
374 // in conforming mode, for each vertex id (as obtained by vertexmapper)
375 // hold its number in the iteration order (VertexIterator)
376 int offset;
377
378 // hide operator ->
379 void operator->();
380 protected:
381 void basicIncrement ()
382 {
383 if( git == gend )
384 return;
385 ++cornerIndexDune;
386 const int numCorners = git->subEntities(n);
387 if( cornerIndexDune == numCorners )
388 {
389 offset += numCorners;
390 cornerIndexDune = 0;
391
392 ++git;
393 while( (git != gend) && skipEntity( git->partitionType() ) )
394 ++git;
395 }
396 }
397 public:
398 VertexIterator(const GridCellIterator & x,
399 const GridCellIterator & end,
400 const VTK::DataMode & dm,
401 const VertexMapper & vm) :
402 git(x), gend(end), datamode(dm), cornerIndexDune(0),
403 vertexmapper(vm), visited(vm.size(), false),
404 offset(0)
405 {
406 if (datamode == VTK::conforming && git != gend)
407 visited[vertexmapper.subIndex(*git,cornerIndexDune,n)] = true;
408 }
409 void increment ()
410 {
411 switch (datamode)
412 {
413 case VTK::conforming :
414 while(visited[vertexmapper.subIndex(*git,cornerIndexDune,n)])
415 {
416 basicIncrement();
417 if (git == gend) return;
418 }
419 visited[vertexmapper.subIndex(*git,cornerIndexDune,n)] = true;
420 break;
421 case VTK::nonconforming :
422 basicIncrement();
423 break;
424 }
425 }
426 bool equals (const VertexIterator & cit) const
427 {
428 return git == cit.git
429 && cornerIndexDune == cit.cornerIndexDune
430 && datamode == cit.datamode;
431 }
432 EntityReference dereference() const
433 {
434 return *git;
435 }
437 int localindex () const
438 {
439 return cornerIndexDune;
440 }
443 {
444 return referenceElement<DT,n>(git->type())
445 .position(cornerIndexDune,n);
446 }
447 };
448
449 VertexIterator vertexBegin () const
450 {
451 return VertexIterator( gridView_.template begin< 0, VTK_Partition >(),
452 gridView_.template end< 0, VTK_Partition >(),
453 datamode, *vertexmapper );
454 }
455
456 VertexIterator vertexEnd () const
457 {
458 return VertexIterator( gridView_.template end< 0, VTK_Partition >(),
459 gridView_.template end< 0, VTK_Partition >(),
460 datamode, *vertexmapper );
461 }
462
464
479 public ForwardIteratorFacade<CornerIterator, const Entity, EntityReference, int>
480 {
481 GridCellIterator git;
482 GridCellIterator gend;
483 VTK::DataMode datamode;
484 // Index of the currently visited corner within the current element.
485 // NOTE: this is in VTK-numbering, in contrast to VertexIterator.
486 int cornerIndexVTK;
487 const VertexMapper & vertexmapper;
488 // in conforming mode, for each vertex id (as obtained by vertexmapper)
489 // hold its number in the iteration order of VertexIterator (*not*
490 // CornerIterator)
491 const std::vector<int> & number;
492 // holds the number of corners of all the elements we have seen so far,
493 // excluding the current element
494 int offset;
495
496 // hide operator ->
497 void operator->();
498 public:
499 CornerIterator(const GridCellIterator & x,
500 const GridCellIterator & end,
501 const VTK::DataMode & dm,
502 const VertexMapper & vm,
503 const std::vector<int> & num) :
504 git(x), gend(end), datamode(dm), cornerIndexVTK(0),
505 vertexmapper(vm),
506 number(num), offset(0) {}
507 void increment ()
508 {
509 if( git == gend )
510 return;
511 ++cornerIndexVTK;
512 const int numCorners = git->subEntities(n);
513 if( cornerIndexVTK == numCorners )
514 {
515 offset += numCorners;
516 cornerIndexVTK = 0;
517
518 ++git;
519 while( (git != gend) && skipEntity( git->partitionType() ) )
520 ++git;
521 }
522 }
523 bool equals (const CornerIterator & cit) const
524 {
525 return git == cit.git
526 && cornerIndexVTK == cit.cornerIndexVTK
527 && datamode == cit.datamode;
528 }
529 EntityReference dereference() const
530 {
531 return *git;
532 }
534
538 int id () const
539 {
540 switch (datamode)
541 {
542 case VTK::conforming :
543 return
544 number[vertexmapper.subIndex(*git,VTK::renumber(*git,cornerIndexVTK),
545 n)];
546 case VTK::nonconforming :
547 return offset + VTK::renumber(*git,cornerIndexVTK);
548 default :
549 DUNE_THROW(IOError,"VTKWriter: unsupported DataMode" << datamode);
550 }
551 }
552 };
553
554 CornerIterator cornerBegin () const
555 {
556 return CornerIterator( gridView_.template begin< 0, VTK_Partition >(),
557 gridView_.template end< 0, VTK_Partition >(),
558 datamode, *vertexmapper, number );
559 }
560
561 CornerIterator cornerEnd () const
562 {
563 return CornerIterator( gridView_.template end< 0, VTK_Partition >(),
564 gridView_.template end< 0, VTK_Partition >(),
565 datamode, *vertexmapper, number );
566 }
567
568 public:
576 explicit VTKWriter ( const GridView &gridView,
577 VTK::DataMode dm = VTK::conforming )
578 : gridView_( gridView ),
579 datamode( dm ),
580 polyhedralCellsPresent_( checkForPolyhedralCells() )
581 { }
582
587 void addCellData (const std::shared_ptr< const VTKFunction > & p)
588 {
589 celldata.push_back(VTKLocalFunction(p));
590 }
591
595 template<typename F>
596 void addCellData(F&& f, VTK::FieldInfo vtkFieldInfo)
597 {
598 celldata.push_back(VTKLocalFunction(std::forward<F>(f),vtkFieldInfo));
599 }
600
616 template<class Container>
617 void addCellData (const Container& v, const std::string &name, int ncomps = 1)
618 {
620 for (int c=0; c<ncomps; ++c) {
621 std::stringstream compName;
622 compName << name;
623 if (ncomps>1)
624 compName << "[" << c << "]";
625 VTKFunction* p = new Function(gridView_, v, compName.str(), ncomps, c);
626 addCellData(std::shared_ptr< const VTKFunction >(p));
627 }
628 }
629
634 void addVertexData (const std::shared_ptr< const VTKFunction > & p)
635 {
636 vertexdata.push_back(VTKLocalFunction(p));
637 }
638
642 template<typename F>
643 void addVertexData(F&& f, VTK::FieldInfo vtkFieldInfo)
644 {
645 vertexdata.push_back(VTKLocalFunction(std::forward<F>(f),vtkFieldInfo));
646 }
647
648
664 template<class Container>
665 void addVertexData (const Container& v, const std::string &name, int ncomps=1)
666 {
668 for (int c=0; c<ncomps; ++c) {
669 std::stringstream compName;
670 compName << name;
671 if (ncomps>1)
672 compName << "[" << c << "]";
673 VTKFunction* p = new Function(gridView_, v, compName.str(), ncomps, c);
674 addVertexData(std::shared_ptr< const VTKFunction >(p));
675 }
676 }
677
679 void clear ()
680 {
681 celldata.clear();
682 vertexdata.clear();
683 }
684
686 virtual ~VTKWriter ()
687 {
688 this->clear();
689 }
690
702 std::string write ( const std::string &name,
703 VTK::OutputType type = VTK::ascii )
704 {
705 return write( name, type, gridView_.comm().rank(), gridView_.comm().size() );
706 }
707
734 std::string pwrite ( const std::string & name, const std::string & path, const std::string & extendpath,
735 VTK::OutputType type = VTK::ascii )
736 {
737 return pwrite( name, path, extendpath, type, gridView_.comm().rank(), gridView_.comm().size() );
738 }
739
740 protected:
742
753 std::string getParallelPieceName(const std::string& name,
754 const std::string& path,
755 int commRank, int commSize) const
756 {
757 std::ostringstream s;
758 if(path.size() > 0) {
759 s << path;
760 if(path[path.size()-1] != '/')
761 s << '/';
762 }
763 s << 's' << std::setw(4) << std::setfill('0') << commSize << '-';
764 s << 'p' << std::setw(4) << std::setfill('0') << commRank << '-';
765 s << name;
766 if(GridView::dimension > 1)
767 s << ".vtu";
768 else
769 s << ".vtp";
770 return s.str();
771 }
772
774
784 std::string getParallelHeaderName(const std::string& name,
785 const std::string& path,
786 int commSize) const
787 {
788 std::ostringstream s;
789 if(path.size() > 0) {
790 s << path;
791 if(path[path.size()-1] != '/')
792 s << '/';
793 }
794 s << 's' << std::setw(4) << std::setfill('0') << commSize << '-';
795 s << name;
796 if(GridView::dimension > 1)
797 s << ".pvtu";
798 else
799 s << ".pvtp";
800 return s.str();
801 }
802
804
816 std::string getSerialPieceName(const std::string& name,
817 const std::string& path) const
818 {
819 static const std::string extension =
820 GridView::dimension == 1 ? ".vtp" : ".vtu";
821
822 return concatPaths(path, name+extension);
823 }
824
840 std::string write ( const std::string &name,
841 VTK::OutputType type,
842 const int commRank,
843 const int commSize )
844 {
845 // in the parallel case, just use pwrite, it has all the necessary
846 // stuff, so we don't need to reimplement it here.
847 if(commSize > 1)
848 return pwrite(name, "", "", type, commRank, commSize);
849
850 // make data mode visible to private functions
851 outputtype = type;
852
853 // generate filename for process data
854 std::string pieceName = getSerialPieceName(name, "");
855
856 // write process data
857 std::ofstream file;
858 file.exceptions(std::ios_base::badbit | std::ios_base::failbit |
859 std::ios_base::eofbit);
860 // check if file can be opened
861 try {
862 file.open( pieceName.c_str(), std::ios::binary );
863 }
864 catch(...) {
865 std::cerr << "Filename: " << pieceName << " could not be opened" << std::endl;
866 throw;
867 }
868 if (! file.is_open())
869 DUNE_THROW(IOError, "Could not write to piece file " << pieceName);
870 writeDataFile( file );
871 file.close();
872
873 return pieceName;
874 }
875
877
900 std::string pwrite(const std::string& name, const std::string& path,
901 const std::string& extendpath,
902 VTK::OutputType ot, const int commRank,
903 const int commSize )
904 {
905 // make data mode visible to private functions
906 outputtype=ot;
907
908 // do some magic because paraview can only cope with relative paths to piece files
909 std::ofstream file;
910 file.exceptions(std::ios_base::badbit | std::ios_base::failbit |
911 std::ios_base::eofbit);
912 std::string piecepath = concatPaths(path, extendpath);
913 std::string relpiecepath = relativePath(path, piecepath);
914
915 // write this processes .vtu/.vtp piece file
916 std::string fullname = getParallelPieceName(name, piecepath, commRank,
917 commSize);
918 // check if file can be opened
919 try {
920 file.open(fullname.c_str(),std::ios::binary);
921 }
922 catch(...) {
923 std::cerr << "Filename: " << fullname << " could not be opened" << std::endl;
924 throw;
925 }
926 if (! file.is_open())
927 DUNE_THROW(IOError, "Could not write to piecefile file " << fullname);
928 writeDataFile(file);
929 file.close();
930 gridView_.comm().barrier();
931
932 // if we are rank 0, write .pvtu/.pvtp parallel header
933 fullname = getParallelHeaderName(name, path, commSize);
934 if( commRank ==0 )
935 {
936 file.open(fullname.c_str());
937 if (! file.is_open())
938 DUNE_THROW(IOError, "Could not write to parallel file " << fullname);
939 writeParallelHeader(file,name,relpiecepath, commSize );
940 file.close();
941 }
942 gridView_.comm().barrier();
943 return fullname;
944 }
945
946 private:
948
965 void writeParallelHeader(std::ostream& s, const std::string& piecename,
966 const std::string& piecepath, const int commSize)
967 {
968 VTK::FileType fileType =
969 (n == 1) ? VTK::polyData : VTK::unstructuredGrid;
970
971 VTK::PVTUWriter writer(s, fileType);
972
973 writer.beginMain();
974
975 // PPointData
976 {
977 std::string scalars, vectors;
978 std::tie(scalars,vectors) = getDataNames(vertexdata);
979 writer.beginPointData(scalars, vectors);
980 }
981 for (auto it = vertexdata.begin(),
982 end = vertexdata.end();
983 it != end;
984 ++it)
985 {
986 unsigned writecomps = it->fieldInfo().size();
987 if(writecomps == 2) writecomps = 3;
988 writer.addArray<float>(it->name(), writecomps);
989 }
990 writer.endPointData();
991
992 // PCellData
993 {
994 std::string scalars, vectors;
995 std::tie(scalars,vectors) = getDataNames(celldata);
996 writer.beginCellData(scalars, vectors);
997 }
998 for (auto it = celldata.begin(),
999 end = celldata.end();
1000 it != end;
1001 ++it)
1002 {
1003 unsigned writecomps = it->fieldInfo().size();
1004 if(writecomps == 2) writecomps = 3;
1005 writer.addArray<float>(it->name(), writecomps);
1006 }
1007 writer.endCellData();
1008
1009 // PPoints
1010 writer.beginPoints();
1011 writer.addArray<float>("Coordinates", 3);
1012 writer.endPoints();
1013
1014 // Pieces
1015 for( int i = 0; i < commSize; ++i )
1016 {
1017 const std::string& fullname = getParallelPieceName(piecename,
1018 piecepath, i,
1019 commSize);
1020 writer.addPiece(fullname);
1021 }
1022
1023 writer.endMain();
1024 }
1025
1027 void writeDataFile (std::ostream& s)
1028 {
1029 VTK::FileType fileType =
1030 (n == 1) ? VTK::polyData : VTK::unstructuredGrid;
1031
1032 VTK::VTUWriter writer(s, outputtype, fileType);
1033
1034 // Grid characteristics
1035 vertexmapper = new VertexMapper( gridView_, mcmgVertexLayout() );
1036 if (datamode == VTK::conforming)
1037 {
1038 number.resize(vertexmapper->size());
1039 for (std::vector<int>::size_type i=0; i<number.size(); i++) number[i] = -1;
1040 }
1041 countEntities(nvertices, ncells, ncorners);
1042
1043 writer.beginMain(ncells, nvertices);
1044 writeAllData(writer);
1045 writer.endMain();
1046
1047 // write appended binary data section
1048 if(writer.beginAppended())
1049 writeAllData(writer);
1050 writer.endAppended();
1051
1052 delete vertexmapper; number.clear();
1053 }
1054
1055 void writeAllData(VTK::VTUWriter& writer) {
1056 // PointData
1057 writeVertexData(writer);
1058
1059 // CellData
1060 writeCellData(writer);
1061
1062 // Points
1063 writeGridPoints(writer);
1064
1065 // Cells
1066 writeGridCells(writer);
1067 }
1068
1069 protected:
1070 std::string getFormatString() const
1071 {
1072 if (outputtype==VTK::ascii)
1073 return "ascii";
1074 if (outputtype==VTK::base64)
1075 return "binary";
1076 if (outputtype==VTK::appendedraw)
1077 return "appended";
1078 if (outputtype==VTK::appendedbase64)
1079 return "appended";
1080 DUNE_THROW(IOError, "VTKWriter: unsupported OutputType" << outputtype);
1081 }
1082
1083 std::string getTypeString() const
1084 {
1085 if (n==1)
1086 return "PolyData";
1087 else
1088 return "UnstructuredGrid";
1089 }
1090
1092 virtual void countEntities(int &nvertices, int &ncells, int &ncorners)
1093 {
1094 nvertices = 0;
1095 ncells = 0;
1096 ncorners = 0;
1097 for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
1098 {
1099 ncells++;
1100 // because of the use of vertexmapper->map(), this iteration must be
1101 // in the order of Dune's numbering.
1102 const int subEntities = it->subEntities(n);
1103 for (int i=0; i<subEntities; ++i)
1104 {
1105 ncorners++;
1106 if (datamode == VTK::conforming)
1107 {
1108 int alpha = vertexmapper->subIndex(*it,i,n);
1109 if (number[alpha]<0)
1110 number[alpha] = nvertices++;
1111 }
1112 else
1113 {
1114 nvertices++;
1115 }
1116 }
1117 }
1118 }
1119
1120 template<typename T>
1121 std::tuple<std::string,std::string> getDataNames(const T& data) const
1122 {
1123 std::string scalars = "";
1124 for (auto it = data.begin(),
1125 end = data.end();
1126 it != end;
1127 ++it)
1128 if (it->fieldInfo().type() == VTK::FieldInfo::Type::scalar)
1129 {
1130 scalars = it->name();
1131 break;
1132 }
1133
1134 std::string vectors = "";
1135 for (auto it = data.begin(),
1136 end = data.end();
1137 it != end;
1138 ++it)
1139 if (it->fieldInfo().type() == VTK::FieldInfo::Type::vector)
1140 {
1141 vectors = it->name();
1142 break;
1143 }
1144 return std::make_tuple(scalars,vectors);
1145 }
1146
1147 template<typename Data, typename Iterator>
1148 void writeData(VTK::VTUWriter& writer, const Data& data, const Iterator begin, const Iterator end, int nentries)
1149 {
1150 for (auto it = data.begin(),
1151 iend = data.end();
1152 it != iend;
1153 ++it)
1154 {
1155 const auto& f = *it;
1156 VTK::FieldInfo fieldInfo = f.fieldInfo();
1157 std::size_t writecomps = fieldInfo.size();
1158 switch (fieldInfo.type())
1159 {
1161 break;
1163 // vtk file format: a vector data always should have 3 comps (with
1164 // 3rd comp = 0 in 2D case)
1165 if (writecomps > 3)
1166 DUNE_THROW(IOError,"Cannot write VTK vectors with more than 3 components (components was " << writecomps << ")");
1167 writecomps = 3;
1168 break;
1170 DUNE_THROW(NotImplemented,"VTK output for tensors not implemented yet");
1171 }
1172 std::shared_ptr<VTK::DataArrayWriter<float> > p
1173 (writer.makeArrayWriter<float>(f.name(), writecomps, nentries));
1174 if(!p->writeIsNoop())
1175 for (Iterator eit = begin; eit!=end; ++eit)
1176 {
1177 const Entity & e = *eit;
1178 f.bind(e);
1179 f.write(eit.position(),*p);
1180 f.unbind();
1181 // vtk file format: a vector data always should have 3 comps
1182 // (with 3rd comp = 0 in 2D case)
1183 for (std::size_t j=fieldInfo.size(); j < writecomps; ++j)
1184 p->write(0.0);
1185 }
1186 }
1187 }
1188
1190 virtual void writeCellData(VTK::VTUWriter& writer)
1191 {
1192 if(celldata.size() == 0)
1193 return;
1194
1195 std::string scalars, vectors;
1196 std::tie(scalars,vectors) = getDataNames(celldata);
1197
1198 writer.beginCellData(scalars, vectors);
1199 writeData(writer,celldata,cellBegin(),cellEnd(),ncells);
1200 writer.endCellData();
1201 }
1202
1204 virtual void writeVertexData(VTK::VTUWriter& writer)
1205 {
1206 if(vertexdata.size() == 0)
1207 return;
1208
1209 std::string scalars, vectors;
1210 std::tie(scalars,vectors) = getDataNames(vertexdata);
1211
1212 writer.beginPointData(scalars, vectors);
1213 writeData(writer,vertexdata,vertexBegin(),vertexEnd(),nvertices);
1214 writer.endPointData();
1215 }
1216
1218 virtual void writeGridPoints(VTK::VTUWriter& writer)
1219 {
1220 writer.beginPoints();
1221
1222 std::shared_ptr<VTK::DataArrayWriter<float> > p
1223 (writer.makeArrayWriter<float>("Coordinates", 3, nvertices));
1224 if(!p->writeIsNoop()) {
1225 VertexIterator vEnd = vertexEnd();
1226 for (VertexIterator vit=vertexBegin(); vit!=vEnd; ++vit)
1227 {
1228 int dimw=w;
1229 for (int j=0; j<std::min(dimw,3); j++)
1230 p->write((*vit).geometry().corner(vit.localindex())[j]);
1231 for (int j=std::min(dimw,3); j<3; j++)
1232 p->write(0.0);
1233 }
1234 }
1235 // free the VTK::DataArrayWriter before touching the stream
1236 p.reset();
1237
1238 writer.endPoints();
1239 }
1240
1242 virtual void writeGridCells(VTK::VTUWriter& writer)
1243 {
1244 writer.beginCells();
1245
1246 // connectivity
1247 {
1248 std::shared_ptr<VTK::DataArrayWriter<int> > p1
1249 (writer.makeArrayWriter<int>("connectivity", 1, ncorners));
1250 if(!p1->writeIsNoop())
1251 for (CornerIterator it=cornerBegin(); it!=cornerEnd(); ++it)
1252 p1->write(it.id());
1253 }
1254
1255 // offsets
1256 {
1257 std::shared_ptr<VTK::DataArrayWriter<int> > p2
1258 (writer.makeArrayWriter<int>("offsets", 1, ncells));
1259 if(!p2->writeIsNoop()) {
1260 int offset = 0;
1261 for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
1262 {
1263 offset += it->subEntities(n);
1264 p2->write(offset);
1265 }
1266 }
1267 }
1268
1269 // types
1270 if (n>1)
1271 {
1272 {
1273 std::shared_ptr<VTK::DataArrayWriter<unsigned char> > p3
1274 (writer.makeArrayWriter<unsigned char>("types", 1, ncells));
1275
1276 if(!p3->writeIsNoop())
1277 {
1278 for (CellIterator it=cellBegin(); it!=cellEnd(); ++it)
1279 {
1280 int vtktype = VTK::geometryType(it->type());
1281 p3->write(vtktype);
1282 }
1283 }
1284 }
1285
1286
1287 // if polyhedron cells found also cell faces need to be written
1288 if( polyhedralCellsPresent_ )
1289 {
1290 writeCellFaces( writer );
1291 }
1292 }
1293
1294 writer.endCells();
1295 }
1296
1297 protected:
1298 bool checkForPolyhedralCells() const
1299 {
1300 // check if polyhedron cells are present
1301 for( const auto& geomType : gridView_.indexSet().types( 0 ) )
1302 {
1303 if( VTK::geometryType( geomType ) == VTK::polyhedron )
1304 {
1305 return true;
1306 }
1307 }
1308 return false;
1309 }
1310
1312 virtual void writeCellFaces(VTK::VTUWriter& writer)
1313 {
1314 if( ! faceVertices_ )
1315 {
1316 faceVertices_.reset( new std::pair< std::vector<int>, std::vector<int> > () );
1317 // fill face vertex structure
1318 fillFaceVertices( cornerBegin(), cornerEnd(), gridView_.indexSet(),
1319 faceVertices_->first, faceVertices_->second );
1320 }
1321
1322 std::vector< int >& faces = faceVertices_->first;
1323 std::vector< int >& faceOffsets = faceVertices_->second;
1324 assert( int(faceOffsets.size()) == ncells );
1325
1326 {
1327 std::shared_ptr<VTK::DataArrayWriter< int > > p4
1328 (writer.makeArrayWriter< int >("faces", 1, faces.size()));
1329 if(!p4->writeIsNoop())
1330 {
1331 for( const auto& face : faces )
1332 p4->write( face );
1333 }
1334 }
1335
1336 {
1337 std::shared_ptr<VTK::DataArrayWriter< int > > p5
1338 (writer.makeArrayWriter< int >("faceoffsets", 1, ncells ));
1339 if(!p5->writeIsNoop())
1340 {
1341 for( const auto& offset : faceOffsets )
1342 p5->write( offset );
1343
1344 // clear face vertex structure
1345 faceVertices_.reset();
1346 }
1347 }
1348 }
1349
1350 template <class CornerIterator, class IndexSet, class T>
1351 inline void fillFaceVertices( CornerIterator it,
1352 const CornerIterator end,
1353 const IndexSet& indexSet,
1354 std::vector<T>& faces,
1355 std::vector<T>& faceOffsets )
1356 {
1357 if( n == 3 && it != end )
1358 {
1359 // clear output arrays
1360 faces.clear();
1361 faces.reserve( 15 * ncells );
1362 faceOffsets.clear();
1363 faceOffsets.reserve( ncells );
1364
1365 int offset = 0;
1366
1367 Cell element = *it;
1368 int elIndex = indexSet.index( element );
1369 std::vector< T > vertices;
1370 vertices.reserve( 30 );
1371 for( ; it != end; ++it )
1372 {
1373 const Cell& cell = *it ;
1374 const int cellIndex = indexSet.index( cell ) ;
1375 if( elIndex != cellIndex )
1376 {
1377 fillFacesForElement( element, indexSet, vertices, offset, faces, faceOffsets );
1378
1379 vertices.clear();
1380 element = cell ;
1381 elIndex = cellIndex ;
1382 }
1383 vertices.push_back( it.id() );
1384 }
1385
1386 // fill faces for last element
1387 fillFacesForElement( element, indexSet, vertices, offset, faces, faceOffsets );
1388 }
1389 }
1390
1391 template <class Entity, class IndexSet, class T>
1392 static void fillFacesForElement( const Entity& element,
1393 const IndexSet& indexSet,
1394 const std::vector<T>& vertices,
1395 T& offset,
1396 std::vector<T>& faces,
1397 std::vector<T>& faceOffsets )
1398 {
1399 const int dim = n;
1400
1401 std::map< T, T > vxMap;
1402
1403 // get number of local faces
1404 const int nVertices = element.subEntities( dim );
1405 for( int vx = 0; vx < nVertices; ++ vx )
1406 {
1407 const int vxIdx = indexSet.subIndex( element, vx, dim );
1408 vxMap[ vxIdx ] = vertices[ vx ];
1409 }
1410
1411 // get number of local faces
1412 const int nFaces = element.subEntities( 1 );
1413 // store number of faces for current element
1414 faces.push_back( nFaces );
1415 ++offset;
1416 // extract each face as a set of vertex indices
1417 for( int fce = 0; fce < nFaces; ++ fce )
1418 {
1419 // obtain face
1420 const auto face = element.template subEntity< 1 > ( fce );
1421
1422 // get all vertex indices from current face
1423 const int nVxFace = face.subEntities( dim );
1424 faces.push_back( nVxFace );
1425 ++offset ;
1426 for( int i=0; i<nVxFace; ++i )
1427 {
1428 const T vxIndex = indexSet.subIndex( face, i, dim );
1429 assert( vxMap.find( vxIndex ) != vxMap.end() );
1430 faces.push_back( vxMap[ vxIndex ] );
1431 ++offset ;
1432 }
1433 }
1434
1435 // store face offset for each element
1436 faceOffsets.push_back( offset );
1437 }
1438
1439 protected:
1440 // the list of registered functions
1441 std::list<VTKLocalFunction> celldata;
1442 std::list<VTKLocalFunction> vertexdata;
1443
1444 // the grid
1445 GridView gridView_;
1446
1447 // temporary grid information
1448 int ncells;
1449 int nvertices;
1450 int ncorners;
1451 private:
1452 VertexMapper* vertexmapper;
1453 // in conforming mode, for each vertex id (as obtained by vertexmapper)
1454 // hold its number in the iteration order (VertexIterator)
1455 std::vector<int> number;
1456 VTK::DataMode datamode;
1457
1458 // true if polyhedral cells are present in the grid
1459 const bool polyhedralCellsPresent_;
1460
1461 // pointer holding face vertex connectivity if needed
1462 std::shared_ptr< std::pair< std::vector<int>, std::vector<int> > > faceVertices_;
1463
1464 protected:
1465 VTK::OutputType outputtype;
1466 };
1467
1468}
1469
1470#endif
vector space out of a tensor product of fields.
Definition: fvector.hh:93
Base class for stl conformant forward iterators.
Definition: iteratorfacades.hh:144
Base class template for function classes.
Definition: function.hh:29
Grid view abstract base class.
Definition: gridview.hh:60
Default exception class for I/O errors.
Definition: exceptions.hh:229
Index Set Interface base class.
Definition: indexidset.hh:76
IndexType index(const typename Traits::template Codim< cc >::Entity &e) const
Map entity to index. The result of calling this method with an entity that is not in the index set is...
Definition: indexidset.hh:111
Implementation class for a multiple codim and multiple geometry type mapper.
Definition: mcmgmapper.hh:200
Index subIndex(const typename GV::template Codim< 0 >::Entity &e, int i, unsigned int codim) const
Map subentity of codim 0 entity to starting index in array for dof block.
Definition: mcmgmapper.hh:268
size_type size() const
Return total number of entities in the entity set managed by the mapper.
Definition: mcmgmapper.hh:287
Default exception for dummy implementations.
Definition: exceptions.hh:261
Take a vector and interpret it as cell data for the VTKWriter.
Definition: function.hh:90
Take a vector and interpret it as point data for the VTKWriter.
Definition: function.hh:188
A base class for grid functions with any return type and dimension.
Definition: function.hh:39
Base class to write pvd-files which contains a list of all collected vtk-files.
Definition: vtksequencewriterbase.hh:32
Writer for the ouput of grid functions in the vtk format.
Definition: vtksequencewriter.hh:28
Iterator over the grids elements.
Definition: vtkwriter.hh:326
CellIterator(const GridCellIterator &x)
construct a CellIterator from the gridview's Iterator.
Definition: vtkwriter.hh:329
const FieldVector< DT, n > position() const
Definition: vtkwriter.hh:332
Iterate over the elements' corners.
Definition: vtkwriter.hh:480
int id() const
Process-local consecutive zero-starting vertex id.
Definition: vtkwriter.hh:538
Type erasure wrapper for VTK data sets.
Definition: vtkwriter.hh:149
void unbind() const
Unbind the data set from the currently bound entity.
Definition: vtkwriter.hh:301
VTKLocalFunction(F &&f, VTK::FieldInfo fieldInfo)
Construct a VTKLocalFunction for a dune-functions style LocalFunction.
Definition: vtkwriter.hh:258
std::string name() const
Returns the name of the data set.
Definition: vtkwriter.hh:283
const VTK::FieldInfo & fieldInfo() const
Returns the VTK::FieldInfo for the data set.
Definition: vtkwriter.hh:289
void bind(const Entity &e) const
Bind the data set to grid entity e.
Definition: vtkwriter.hh:295
VTKLocalFunction(const std::shared_ptr< const VTKFunction > &vtkFunctionPtr)
Construct a VTKLocalFunction for a legacy VTKFunction.
Definition: vtkwriter.hh:273
void write(const Coordinate &pos, Writer &w) const
Write the value of the data set at local coordinate pos to the writer w.
Definition: vtkwriter.hh:307
Iterate over the grid's vertices.
Definition: vtkwriter.hh:365
FieldVector< DT, n > position() const
position of vertex inside the entity
Definition: vtkwriter.hh:442
int localindex() const
index of vertex within the entity, in Dune-numbering
Definition: vtkwriter.hh:437
Writer for the ouput of grid functions in the vtk format.
Definition: vtkwriter.hh:88
void clear()
clear list of registered functions
Definition: vtkwriter.hh:679
std::string write(const std::string &name, VTK::OutputType type=VTK::ascii)
write output (interface might change later)
Definition: vtkwriter.hh:702
std::string getParallelHeaderName(const std::string &name, const std::string &path, int commSize) const
return name of a parallel header file
Definition: vtkwriter.hh:784
void addVertexData(const std::shared_ptr< const VTKFunction > &p)
Add a grid function that lives on the vertices of the grid to the visualization.
Definition: vtkwriter.hh:634
void addVertexData(const Container &v, const std::string &name, int ncomps=1)
Add a grid function (represented by container) that lives on the vertices of the grid to the visualiz...
Definition: vtkwriter.hh:665
std::string getSerialPieceName(const std::string &name, const std::string &path) const
return name of a serial piece file
Definition: vtkwriter.hh:816
void addCellData(const std::shared_ptr< const VTKFunction > &p)
Add a grid function that lives on the cells of the grid to the visualization.
Definition: vtkwriter.hh:587
void addVertexData(F &&f, VTK::FieldInfo vtkFieldInfo)
Add a grid function that lives on the vertices of the grid to the visualization.
Definition: vtkwriter.hh:643
virtual void writeCellData(VTK::VTUWriter &writer)
write cell data
Definition: vtkwriter.hh:1190
VTKWriter(const GridView &gridView, VTK::DataMode dm=VTK::conforming)
Construct a VTKWriter working on a specific GridView.
Definition: vtkwriter.hh:576
std::string getParallelPieceName(const std::string &name, const std::string &path, int commRank, int commSize) const
return name of a parallel piece file
Definition: vtkwriter.hh:753
virtual void countEntities(int &nvertices, int &ncells, int &ncorners)
count the vertices, cells and corners
Definition: vtkwriter.hh:1092
void addCellData(const Container &v, const std::string &name, int ncomps=1)
Add a grid function (represented by container) that lives on the cells of the grid to the visualizati...
Definition: vtkwriter.hh:617
virtual void writeGridCells(VTK::VTUWriter &writer)
write the connectivity array
Definition: vtkwriter.hh:1242
virtual void writeCellFaces(VTK::VTUWriter &writer)
write the connectivity array
Definition: vtkwriter.hh:1312
std::string write(const std::string &name, VTK::OutputType type, const int commRank, const int commSize)
write output (interface might change later)
Definition: vtkwriter.hh:840
virtual void writeGridPoints(VTK::VTUWriter &writer)
write the positions of vertices
Definition: vtkwriter.hh:1218
virtual void writeVertexData(VTK::VTUWriter &writer)
write vertex data
Definition: vtkwriter.hh:1204
void addCellData(F &&f, VTK::FieldInfo vtkFieldInfo)
Add a grid function that lives on the cells of the grid to the visualization.
Definition: vtkwriter.hh:596
virtual ~VTKWriter()
destructor
Definition: vtkwriter.hh:686
std::string pwrite(const std::string &name, const std::string &path, const std::string &extendpath, VTK::OutputType ot, const int commRank, const int commSize)
write output; interface might change later
Definition: vtkwriter.hh:900
std::string pwrite(const std::string &name, const std::string &path, const std::string &extendpath, VTK::OutputType type=VTK::ascii)
write output (interface might change later)
Definition: vtkwriter.hh:734
virtual void write(T data)=0
write one data element
Descriptor struct for VTK fields.
Definition: common.hh:315
@ tensor
tensor field (always 3x3)
@ vector
vector-valued field (always 3D, will be padded if necessary)
std::string name() const
The name of the data field.
Definition: common.hh:338
Dump a .vtu/.vtp files contents to a stream.
Definition: pvtuwriter.hh:60
Dump a .vtu/.vtp files contents to a stream.
Definition: vtuwriter.hh:96
void endCellData()
finish CellData section
Definition: vtuwriter.hh:218
void beginCells()
start section for the grid cells/PolyData lines
Definition: vtuwriter.hh:272
DataArrayWriter< T > * makeArrayWriter(const std::string &name, unsigned ncomps, unsigned nitems)
acquire a DataArrayWriter
Definition: vtuwriter.hh:379
void endPointData()
finish PointData section
Definition: vtuwriter.hh:180
void beginCellData(const std::string &scalars="", const std::string &vectors="")
start CellData section
Definition: vtuwriter.hh:203
void beginPointData(const std::string &scalars="", const std::string &vectors="")
start PointData section
Definition: vtuwriter.hh:165
void endPoints()
finish section for the point coordinates
Definition: vtuwriter.hh:247
void endCells()
start section for the grid cells/PolyData lines
Definition: vtuwriter.hh:283
void beginPoints()
start section for the point coordinates
Definition: vtuwriter.hh:236
Data array writers for the VTKWriter.
A few common exception classes.
Common stuff for the VTKWriter.
OutputType
How the bulk data should be stored in the file.
Definition: common.hh:40
FileType
which type of VTK file to write
Definition: common.hh:298
DataMode
Whether to produce conforming or non-conforming output.
Definition: common.hh:64
Functions for VTK output.
typename Impl::voider< Types... >::type void_t
Is void for all valid input types (see N3911). The workhorse for C++11 SFINAE-techniques.
Definition: typetraits.hh:39
Impl::MakeUniqueHelper< T >::NonArrayUniquePtr make_unique(Args &&... args)
Implementation of std::make_unique to be introduced in C++14.
Definition: memory.hh:66
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
const IndexSet & indexSet() const
obtain the index set
Definition: gridview.hh:172
Traits::Grid Grid
type of the grid
Definition: gridview.hh:77
const CollectiveCommunication & comm() const
obtain collective communication object
Definition: gridview.hh:246
Traits::IndexSet IndexSet
type of the index set
Definition: gridview.hh:80
Grid::ctype ctype
type used for coordinates in grid
Definition: gridview.hh:124
@ dimensionworld
The dimension of the world the grid lives in.
Definition: gridview.hh:131
@ dimension
The dimension of the grid.
Definition: gridview.hh:127
PartitionIteratorType
Parameter to be used for the parallel level- and leaf iterators.
Definition: gridenums.hh:134
PartitionType
Attributes used in the generic overlap model.
Definition: gridenums.hh:28
@ All_Partition
all entities
Definition: gridenums.hh:139
@ InteriorBorder_Partition
interior and border entities
Definition: gridenums.hh:136
@ InteriorEntity
all interior entities
Definition: gridenums.hh:29
MCMGLayout mcmgVertexLayout()
layout for vertices (dim-0 entities)
Definition: mcmgmapper.hh:160
std::string relativePath(const std::string &newbase, const std::string &p)
compute a relative path between two paths
Definition: path.cc:151
std::string concatPaths(const std::string &base, const std::string &p)
concatenate two paths
Definition: path.cc:30
Utility class for handling nested indentation in output.
This file implements iterator facade classes for writing stl conformant iterators.
Mapper for multiple codim and multiple geometry types.
Dune namespace.
Definition: alignedallocator.hh:10
STL namespace.
Utilities for handling filesystem paths.
Static tag representing a codimension.
Definition: dimension.hh:22
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:196
Base class for polymorphic container of underlying data set.
Definition: vtkwriter.hh:157
virtual void write(const Coordinate &pos, Writer &w, std::size_t count) const =0
Evaluate data set at local position pos inside the current entity and write result to w.
virtual void unbind()=0
Unbind data set from current grid entity - mostly here for performance and symmetry reasons.
virtual void bind(const Entity &e)=0
Bind data set to grid entity - must be called before evaluating (i.e. calling write())
Type erasure implementation for functions conforming to the dune-functions LocalFunction interface.
Definition: vtkwriter.hh:180
virtual void write(const Coordinate &pos, Writer &w, std::size_t count) const
Evaluate data set at local position pos inside the current entity and write result to w.
Definition: vtkwriter.hh:198
virtual void unbind()
Unbind data set from current grid entity - mostly here for performance and symmetry reasons.
Definition: vtkwriter.hh:193
virtual void bind(const Entity &e)
Bind data set to grid entity - must be called before evaluating (i.e. calling write())
Definition: vtkwriter.hh:188
Type erasure implementation for legacy VTKFunctions.
Definition: vtkwriter.hh:227
virtual void unbind()
Unbind data set from current grid entity - mostly here for performance and symmetry reasons.
Definition: vtkwriter.hh:238
virtual void write(const Coordinate &pos, Writer &w, std::size_t count) const
Evaluate data set at local position pos inside the current entity and write result to w.
Definition: vtkwriter.hh:243
virtual void bind(const Entity &e)
Bind data set to grid entity - must be called before evaluating (i.e. calling write())
Definition: vtkwriter.hh:233
Type trait to determine whether an instance of T has an operator[](I), i.e. whether it can be indexed...
Definition: typetraits.hh:222
Traits for type conversions and type information.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)