1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3
4 #include "config.h"
5
6 #include <memory>
7 #include <vector>
8 #include <fstream>
9 #include <string>
10 #include <unistd.h>
11
12 // dune headers
13 #include <dune/grid/yaspgrid.hh>
14 #include <dune/grid/io/file/vtk/vtksequencewriter.hh>
15
16 // check that the files have equal amount of lines
checkFileEqualNumLines(const std::string & name1,const std::string & name2)17 void checkFileEqualNumLines(const std::string& name1, const std::string& name2)
18 {
19 std::ifstream pvdFile1(name1);
20 std::ifstream pvdFile2(name2);
21
22 if (pvdFile1.fail())
23 DUNE_THROW(Dune::Exception, "File " << name1 << " could not be opened!");
24
25 if (pvdFile2.fail())
26 DUNE_THROW(Dune::Exception, "File " << name2 << " could not be opened!");
27
28 std::string line;
29 unsigned int lines1 = 0, lines2 = 0;
30 while (std::getline(pvdFile1, line)) ++lines1;
31 while (std::getline(pvdFile2, line)) ++lines2;
32
33 if (lines1 != lines2)
34 DUNE_THROW(Dune::Exception, "Not the same number of lines (comparing " << name1 << " and " << name2 << ")");
35 }
36
VTKDataMode(Dune::VTK::DataMode dm)37 std::string VTKDataMode(Dune::VTK::DataMode dm)
38 {
39 switch(dm)
40 {
41 case Dune::VTK::conforming :
42 return "conforming";
43 case Dune::VTK::nonconforming :
44 return "nonconforming";
45 }
46 return "";
47 }
48
49 template< class GridView >
50 class VTKVectorFunction
51 : public Dune :: VTKWriter< GridView > :: VTKFunction
52 {
53 // extract types
54 enum { n = GridView :: dimension };
55 enum { w = GridView :: dimensionworld };
56 typedef typename GridView :: Grid :: ctype DT;
57 typedef typename GridView :: template Codim< 0 > :: Entity Entity;
58 double time_;
59 public:
VTKVectorFunction()60 VTKVectorFunction() : time_(0) {}
setTime(double time)61 void setTime(double time) {
62 time_ = time;
63 }
64 //! return number of components
ncomps() const65 virtual int ncomps () const { return n; }
66
67 //! evaluate single component comp in the entity e at local coordinates xi
68 /*! Evaluate the function in an entity at local coordinates.
69 @param[in] comp number of component to be evaluated
70 @param[in] e reference to grid entity of codimension 0
71 @param[in] xi point in local coordinates of the reference element of e
72 \return value of the component
73 */
evaluate(int comp,const Entity & e,const Dune::FieldVector<DT,n> & xi) const74 virtual double evaluate (int comp, [[maybe_unused]] const Entity& e,
75 [[maybe_unused]] const Dune::FieldVector<DT,n>& xi) const
76 {
77 return comp*0.1*sin(time_*2.*M_PI);
78 }
79
80 // get name
name() const81 virtual std::string name () const
82 {
83 char _name[256];
84 snprintf(_name, 256, "vector-%iD", ncomps());
85 return std::string(_name);
86 }
87 };
88
89 template< class GridView >
doWrite(const GridView & gridView,Dune::VTK::DataMode dm,bool testRestart=false)90 std::string doWrite( const GridView &gridView, Dune::VTK::DataMode dm,
91 bool testRestart = false)
92 {
93 enum { dim = GridView :: dimension };
94
95 const typename GridView :: IndexSet &is = gridView.indexSet();
96 std::vector<int> vertexdata(is.size(dim),dim);
97 std::vector<int> celldata(is.size(0),0);
98
99 std::stringstream name;
100 name << "vtktest-" << dim << "D-" << VTKDataMode(dm);
101
102 if (testRestart)
103 name << "-restart";
104
105 auto vtkWriter = std::make_shared<Dune::VTKWriter<GridView> >(gridView, dm);
106 Dune :: VTKSequenceWriter< GridView > vtk( vtkWriter, name.str(), ".", "" );
107
108 vtk.addVertexData(vertexdata,"vertexData");
109 vtk.addCellData(celldata,"cellData");
110 auto vectordata = std::make_shared<VTKVectorFunction<GridView> >();
111 vtk.addVertexData(vectordata);
112 double time = 0;
113 double tEnd = testRestart ? 0.5 : 1.0;
114 while (time < tEnd) {
115 vectordata->setTime(time);
116 vtk.write(time);
117 time += 0.1;
118 }
119
120 if (testRestart)
121 {
122 const auto& timeSteps = vtk.getTimeSteps();
123 auto vtkWriter2 = std::make_shared<Dune::VTKWriter<GridView> >(gridView, dm);
124 Dune :: VTKSequenceWriter< GridView > vtk2( vtkWriter2, name.str(), ".", "" );
125 vtk2.setTimeSteps(timeSteps);
126
127 vtk2.addVertexData(vertexdata,"vertexData");
128 vtk2.addCellData(celldata,"cellData");
129 vtk2.addVertexData(vectordata);
130
131 while (time < 1) {
132 vectordata->setTime(time);
133 vtk2.write(time);
134 time += 0.1;
135 }
136 }
137
138 return name.str();
139 }
140
141 template<int dim>
vtkCheck(const std::array<int,dim> & n,const Dune::FieldVector<double,dim> & h,bool testRestart=false)142 void vtkCheck(const std::array<int,dim>& n,
143 const Dune::FieldVector<double,dim>& h,
144 bool testRestart = false)
145 {
146 std::cout << std::endl << "vtkSequenceCheck dim=" << dim << std::endl << std::endl;
147 Dune::YaspGrid<dim> g(h, n);
148 g.globalRefine(1);
149
150 const auto fileName = doWrite( g.leafGridView(), Dune::VTK::conforming ) + ".pvd";
151 const auto fileNameRestart = doWrite( g.leafGridView(), Dune::VTK::conforming, testRestart ) + ".pvd";
152 checkFileEqualNumLines(fileName, fileNameRestart);
153
154 doWrite( g.leafGridView(), Dune::VTK::nonconforming );
155 doWrite( g.levelGridView( 0 ), Dune::VTK::conforming );
156 doWrite( g.levelGridView( 0 ), Dune::VTK::nonconforming );
157 doWrite( g.levelGridView( g.maxLevel() ), Dune::VTK::conforming );
158 doWrite( g.levelGridView( g.maxLevel() ), Dune::VTK::nonconforming );
159 }
160
main(int argc,char ** argv)161 int main(int argc, char **argv)
162 {
163 try {
164
165 const Dune::MPIHelper &mpiHelper = Dune::MPIHelper::instance(argc, argv);
166
167 if(mpiHelper.rank() == 0)
168 std::cout << "subsamplingvtktest: MPI_Comm_size == " << mpiHelper.size()
169 << std::endl;
170
171 {
172 std::array<int,1> n = { { 5 } };
173 Dune::FieldVector<double,1> h = { 1.0 };
174 vtkCheck<1>(n,h);
175 }
176 {
177 std::array<int,2> n = { { 5, 5 } };
178 Dune::FieldVector<double,2> h = { 1.0, 2.0 };
179 vtkCheck<2>(n,h);
180 }
181 {
182 std::array<int,3> n = { { 5, 5, 5 } };
183 Dune::FieldVector<double,3> h = { 1.0, 2.0, 3.0 };
184 vtkCheck<3>(n,h, /*testRestart=*/true);
185 }
186
187 } catch (Dune::Exception &e) {
188 std::cerr << e << std::endl;
189 return 1;
190 } catch (std::exception &e) {
191 std::cerr << e.what() << std::endl;
192 return 1;
193 } catch (...) {
194 std::cerr << "Generic exception!" << std::endl;
195 return 2;
196 }
197
198 return 0;
199 }
200