1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3
4 #if HAVE_CONFIG_H
5 #include "config.h" // autoconf defines, needed by the dune headers
6 #endif
7
8 #include <algorithm>
9 #include <iostream>
10 #include <ostream>
11 #include <sstream>
12 #include <string>
13 #include <vector>
14
15 #include <dune/common/std/make_array.hh>
16 #include <dune/common/fvector.hh>
17 #include <dune/common/parallel/mpihelper.hh>
18
19 #include <dune/grid/io/file/test/checkvtkfile.hh>
20 #include <dune/grid/io/file/vtk/vtkwriter.hh>
21 #include <dune/grid/yaspgrid.hh>
22
VTKDataMode(Dune::VTK::DataMode dm)23 const char* VTKDataMode(Dune::VTK::DataMode dm)
24 {
25 switch(dm)
26 {
27 case Dune::VTK::conforming :
28 return "conforming";
29 case Dune::VTK::nonconforming :
30 return "nonconforming";
31 }
32 return "";
33 }
34
35 template< class GridView >
36 class VTKVectorFunction
37 : public Dune :: VTKWriter< GridView > :: VTKFunction
38 {
39 // extract types
40 enum { n = GridView :: dimension };
41 typedef typename GridView :: Grid :: ctype DT;
42 typedef typename GridView :: template Codim< 0 > :: Entity Entity;
43 const std::string type_;
44 public:
45 /** @brief Make a new VTKVectorFunction
46 *
47 * @param type_ Type of the function for use in its name (hint: "cell" or
48 * "vertex")
49 */
VTKVectorFunction(std::string type)50 VTKVectorFunction(std::string type)
51 : type_(type)
52 { }
53
54 //! return number of components
ncomps() const55 virtual int ncomps () const { return n; }
56
57 //! evaluate single component comp in the entity e at local coordinates xi
58 /*! Evaluate the function in an entity at local coordinates.
59 @param[in] comp number of component to be evaluated
60 @param[in] e reference to grid entity of codimension 0
61 @param[in] xi point in local coordinates of the reference element of e
62 \return value of the component
63 */
evaluate(int comp,const Entity & e,const Dune::FieldVector<DT,n> & xi) const64 virtual double evaluate (int comp, [[maybe_unused]] const Entity& e,
65 [[maybe_unused]] const Dune::FieldVector<DT,n>& xi) const
66 {
67 return comp*0.1;
68 }
69
70 // get name
name() const71 virtual std::string name () const
72 {
73 std::ostringstream os;
74 os << type_ << "-vector-" << ncomps() << "D";
75 return os.str();
76 }
77
78 };
79
80 // accumulate exit status
acc(int & accresult,int result)81 void acc(int &accresult, int result)
82 {
83 if(accresult == 0 || (accresult == 77 && result != 0))
84 accresult = result;
85 }
86
87 struct Acc
88 {
operator ()Acc89 int operator()(int v1, int v2) const
90 {
91 acc(v1, v2);
92 return v1;
93 }
94 };
95
96 template< class GridView >
doWrite(Dune::VTKChecker & vtkChecker,const std::string & gridViewName,const GridView & gridView,Dune::VTK::DataMode dm,const std::string & path="./")97 int doWrite( Dune::VTKChecker& vtkChecker, const std::string& gridViewName, const GridView &gridView, Dune :: VTK :: DataMode dm, const std::string& path = "./" )
98 {
99 enum { dim = GridView :: dimension };
100
101 Dune :: VTKWriter< GridView > vtk( gridView, dm );
102
103 const typename GridView :: IndexSet &is = gridView.indexSet();
104 std::vector<int> vertexdata(is.size(dim),dim);
105 std::vector<int> celldata(is.size(0),0);
106 std::vector<double> celldoubledata(is.size(0), 1.00000000000346e-12);
107
108 // Write algebraic data
109 vtk.addVertexData(vertexdata,"vertexData");
110 vtk.addCellData(celldata,"cellData");
111 vtk.addCellData(celldoubledata, "cellData_double", 1, Dune::VTK::Precision::float64);
112
113 // Write functions given as VTKFunction objects
114 vtk.addVertexData(std::make_shared< VTKVectorFunction<GridView> >("vertex"));
115 vtk.addCellData(std::make_shared< VTKVectorFunction<GridView> >("cell"));
116
117 // Write globally defined C++ functions
118 auto f1 = [](const Dune :: FieldVector<double,dim>& x)
119 {
120 return std::sin(x.two_norm());
121 };
122
123 vtk.addCellData(f1,
124 Dune :: VTK :: FieldInfo("scalar-valued lambda",
125 Dune :: VTK :: FieldInfo :: Type :: scalar, 1));
126
127 auto f2 = [](const Dune :: FieldVector<double,dim>& x)
128 {
129 return x;
130 };
131
132 vtk.addVertexData(f2,
133 Dune :: VTK :: FieldInfo("vector-valued lambda",
134 Dune :: VTK :: FieldInfo :: Type :: vector, 2,
135 Dune :: VTK :: Precision::float64));
136
137 int result = 0;
138 std::string name;
139 std::ostringstream prefix;
140 prefix << path << "/";
141
142 prefix << "vtktest-" << dim << "D-" << gridViewName << "-" << VTKDataMode(dm);
143 int rank = gridView.comm().rank();
144
145 name = vtk.write(prefix.str() + "-ascii");
146 if(rank == 0) vtkChecker.push(name);
147
148 name = vtk.write(prefix.str() + "-base64", Dune::VTK::base64);
149 if(rank == 0) vtkChecker.push(name);
150
151 name = vtk.write(prefix.str() + "-appendedraw", Dune::VTK::appendedraw);
152 if(rank == 0) vtkChecker.push(name);
153
154 name = vtk.write(prefix.str() + "-appendedbase64",
155 Dune::VTK::appendedbase64);
156 if(rank == 0) vtkChecker.push(name);
157
158 return result;
159 }
160
161 template<int dim>
vtkCheck(Dune::VTKChecker & vtkChecker,const std::array<int,dim> & elements,const Dune::FieldVector<double,dim> & upperRight)162 int vtkCheck(Dune::VTKChecker& vtkChecker, const std::array<int, dim>& elements,
163 const Dune::FieldVector<double, dim>& upperRight)
164 {
165 Dune::YaspGrid<dim> g(upperRight, elements);
166
167 if(g.comm().rank() == 0)
168 std::cout << std::endl
169 << "vtkCheck dim=" << dim << std::endl
170 << std::endl;
171
172 g.globalRefine(1);
173
174 int result = 0;
175
176 acc(result, doWrite( vtkChecker, "leafview", g.leafGridView(), Dune::VTK::conforming ));
177 acc(result, doWrite( vtkChecker, "leafview", g.leafGridView(), Dune::VTK::nonconforming ));
178 acc(result, doWrite( vtkChecker, "coarselevelview", g.levelGridView( 0 ),
179 Dune::VTK::conforming ));
180 acc(result, doWrite( vtkChecker, "coarselevelview", g.levelGridView( 0 ),
181 Dune::VTK::nonconforming ));
182 acc(result, doWrite( vtkChecker, "finelevelview", g.levelGridView( g.maxLevel() ),
183 Dune::VTK::conforming ));
184 acc(result, doWrite( vtkChecker, "finelevelview", g.levelGridView( g.maxLevel() ),
185 Dune::VTK::nonconforming ));
186
187 // write with filename including a path which is the current directory
188 // this is for testing the parallel version of 'write'
189 acc(result, doWrite( vtkChecker, "leafview", g.leafGridView(), Dune::VTK::conforming, "../test" ));
190 acc(result, doWrite( vtkChecker, "leafview", g.leafGridView(), Dune::VTK::nonconforming, "../test" ));
191
192 return result;
193 }
194
main(int argc,char ** argv)195 int main(int argc, char **argv)
196 {
197 try {
198
199 const Dune::MPIHelper &mpiHelper = Dune::MPIHelper::instance(argc, argv);
200
201 if(mpiHelper.rank() == 0)
202 std::cout << "vtktest: MPI_Comm_size == " << mpiHelper.size()
203 << std::endl;
204
205 int result = 0; // pass by default
206 using Dune::Std::make_array;
207
208 Dune::VTKChecker vtkChecker;
209
210 acc(result, vtkCheck<1>(vtkChecker, make_array(8), {1.0}));
211 acc(result, vtkCheck<2>(vtkChecker, make_array(8,4), {1.0, 2.0}));
212 acc(result, vtkCheck<3>(vtkChecker, make_array(8,4,4), {1.0, 2.0, 3.0}));
213
214 acc(result, vtkChecker.check());
215
216 mpiHelper.getCollectiveCommunication().allreduce<Acc>(&result, 1);
217
218 return result;
219
220 } catch (Dune::Exception &e) {
221 std::cerr << e << std::endl;
222 return 1;
223 } catch (std::exception &e) {
224 std::cerr << e.what() << std::endl;
225 return 1;
226 } catch (...) {
227 std::cerr << "Generic exception!" << std::endl;
228 return 2;
229 }
230
231 }
232