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