pointcloudvtkwriter.hh 12.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
21
 * \ingroup InputOutput
22
 * \brief A VTK writer specialized for staggered grid implementations with dofs on the faces
23
 */
Timo Koch's avatar
Timo Koch committed
24
25
#ifndef DUMUX_POINTCLOUD_VTK_WRITER_HH
#define DUMUX_POINTCLOUD_VTK_WRITER_HH
26
27
28

#include <dune/common/fvector.hh>

29
#include <dumux/io/vtkoutputmodule.hh>
30
#include <dumux/io/staggeredvtkoutputmodule.hh>
31
#include <dune/grid/io/file/vtk/common.hh>
32
#include <dune/common/path.hh>
33

Timo Koch's avatar
Timo Koch committed
34
namespace Dumux {
35
36
37
38
39
40
41
42
43
44

/*!
 * \ingroup InputOutput
 * \brief A VTK output module to simplify writing dumux simulation data to VTK format
 *
 * Handles the output of scalar and vector fields to VTK formatted file for multiple
 * variables and timesteps. Certain predefined fields can be registered on problem / model
 * initialization and/or be turned on/off using the designated properties. Additionally
 * non-standardized scalar and vector fields can be added to the writer manually.
 */
45
46
template<class Scalar, int dim>
class PointCloudVtkWriter
47
{
48
    using GlobalPosition = Dune::FieldVector<Scalar, dim>;
49
50
51
52

    static constexpr unsigned int precision = 6;
    static constexpr unsigned int numBeforeLineBreak = 15;

53
54
55
     /*!
     * \brief A class holding a data container and additional information
     */
56
    template<class ContainerType>
57
    class VTKFunction
58
59
    {
    public:
60
61
62
63
64
65
66
        /*!
        * \brief A class holding a data container and additional information
        *
        * \param data The data container
        * \param name The name of the data set
        * \param numComponents The number of components of the data set
        */
67
        VTKFunction(const ContainerType& data, const std::string& name, const int numComponents) : data_(data), name_(name), numComponents_(numComponents)
68
69
        {}

70
71
72
        /*!
        * \brief Returns the name of the data set
        */
73
74
75
76
77
        const std::string& name() const
        {
            return name_;
        }

78
79
80
        /*!
        * \brief Returns the number of components of the data set
        */
81
82
83
84
85
        int numComponents() const
        {
            return numComponents_;
        }

86
87
88
89
90
91
        /*!
        * \brief Allows random acces to data
        *
        * \param dofIdx The index
        */
        auto& operator() (const int idx) const  { return data_[idx]; }
92
93
94
95
96
97
98
99
100
101
102

        decltype(auto) begin() const
        {
            return data_.begin();
        }

        decltype(auto) end() const
        {
            return data_.end();
        }

103
104
105
        /*!
        * \brief Returns the size of the data container
        */
106
107
108
109
110
111
112
113
114
115
116
        int size() const
        {
            return data_.size();
        }

    private:
        const ContainerType& data_;
        const std::string name_;
        const int numComponents_;
    };

117

118
public:
119
120
    using ScalarFunction = VTKFunction<std::vector<Scalar>>;
    using VectorFunction = VTKFunction<std::vector<GlobalPosition>>;
121
122


123
    PointCloudVtkWriter(const std::vector<GlobalPosition>& coordinates) : coordinates_(coordinates)
124
125
    {}

126
127
128
129
130
131
132
     /*!
     * \brief Create an output file
     *
     * \param name The base name
     * \param type The output type
     */
    void write(const std::string& name, Dune::VTK::OutputType type = Dune::VTK::ascii)
133
    {
134
        auto filename = getSerialPieceName(name, "");
135

136
137
138
139
140
        std::ofstream file;
        file.open(filename);
        writeHeader_(file);
        writeCoordinates_(file, coordinates_);
        writeDataInfo_(file);
141
142
143

        for(auto&& data : scalarPointData_)
        {
144
            writeData_(file, data);
145
146
147
        }
        for(auto&& data :vectorPointData_)
        {
148
            writeData_(file, data);
149
150
        }

151
152
153
154
        file << "</PointData>\n";
        file << "</Piece>\n";
        file << "</PolyData>\n";
        file << "</VTKFile>";
155
156
157

        clear();

158
        file.close();
159
160
    }

161
162
163
164
165
166
167
168
     /*!
     * \brief Create an output file in parallel
     *
     * \param name The base name
     * \param type The output type
     */
    void pwrite(const std::string & name,  const std::string & path, const std::string & extendpath,
                Dune::VTK::OutputType type = Dune::VTK::ascii)
169
    {
170
        DUNE_THROW(Dune::NotImplemented, "parallel point cloud vtk output not supported yet");
171
172
    }

173
174
175
176
177
178
179
     /*!
     * \brief Add a vector of scalar data that live on arbitrary points to the visualization.
     *
     * \param v The vector containing the data
     * \param name The name of the data set
     * \param ncomps The number of components of the data set
     */
180
    void addPointData(const std::vector<Scalar>& v, const std::string &name)
181
    {
182
183
        assert(v.size() == coordinates_.size());
        scalarPointData_.push_back(ScalarFunction(v, name, 1));
184
185
    }

186
187
188
189
190
191
192
     /*!
     * \brief Add a vector of vector data that live on arbitrary points to the visualization.
     *
     * \param v The vector containing the data
     * \param name The name of the data set
     * \param ncomps The number of components of the data set
     */
193
    void addPointData(const std::vector<GlobalPosition>& v, const std::string &name)
194
195
    {
        assert(v.size() == coordinates_.size());
196
        vectorPointData_.push_back(VectorFunction(v, name, 3));
197
198
    }

199
200
201
     /*!
     * \brief Clears the data lists
     */
202
203
204
205
206
207
    void clear()
    {
        scalarPointData_.clear();
        vectorPointData_.clear();
    }

208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
     /*!
     * \brief Return name of a serial header file
     *
     * \param name     Base name of the VTK output.  This should be without
     *                 any directory parts and without a filename extension.
     * \param path     Directory part of the resulting header name.  May be
     *                 empty, in which case the resulting name will not have a
     *                 directory part.  If non-empty, may or may not have a
     *                 trailing '/'.  If a trailing slash is missing, one is
     *                 appended implicitly.
     */
    std::string getSerialPieceName(const std::string& name,
                                   const std::string& path) const
    {
      static const std::string extension = ".vtp";

      return Dune::concatPaths(path, name+extension);
    }

     /*!
     * \brief Return name of a parallel header file
     *
     * \param name     Base name of the VTK output.  This should be without
     *                 any directory parts and without a filename extension.
     * \param path     Directory part of the resulting header name.  May be
     *                 empty, in which case the resulting name will not have a
     *                 directory part.  If non-empty, may or may not have a
     *                 trailing '/'.  If a trailing slash is missing, one is
     *                 appended implicitly.
     * \param commSize Number of processes writing a parallel vtk output.
     */
    std::string getParallelHeaderName(const std::string& name,
                                      const std::string& path,
                                      int commSize) const
    {
      std::ostringstream s;
      if(path.size() > 0) {
        s << path;
        if(path[path.size()-1] != '/')
          s << '/';
      }
      s << 's' << std::setw(4) << std::setfill('0') << commSize << '-';
      s << name;
      s << ".pvtp";
      return s.str();
    }

255

256
private:
257
258
259
     /*!
     * \brief Writes the header to the file
     */
260
    void writeHeader_(std::ostream& file)
261
262
263
264
    {
        std::string header = "<?xml version=\"1.0\"?>\n";
                    header += "<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
                    header += "<PolyData>\n";
265
                    header += "<Piece NumberOfLines=\"0\" NumberOfPoints=\"" + std::to_string(coordinates_.size()) + "\">\n";
266
        file << header;
267
268
    }

269
     /*!
270
     * \brief Writes information about the data to the file
271
     */
272
    void writeDataInfo_(std::ostream& file)
273
    {
274
275
        std::string scalarName;
        std::string vectorName;
276
277
        bool foundScalar = false;
        bool foundVector = false;
278

279
        for(auto&& data : scalarPointData_)
280
        {
281
282
283
284
285
286
287
288
289
290
291
292
            if(data.numComponents() == 1 && !foundScalar)
            {
                scalarName = data.name();
                foundScalar = true;
                continue;
            }

            if(data.numComponents() > 1 && !foundVector)
            {
                vectorName = data.name();
                foundVector = true;
            }
293
        }
294
295

        for(auto&& data : vectorPointData_)
296
        {
297
298
299
300
301
            if(data.numComponents() > 1 && !foundVector)
            {
                vectorName = data.name();
                foundVector = true;
            }
302
303
        }

304
305
        if(foundScalar)
            if(foundVector)
306
                file << "<PointData Scalars=\"" << scalarName << "\" Vectors=\"" << vectorName <<"\">\n";
307
            else
308
                file << "<PointData Scalars=\"" << scalarName << "\">\n";
309
        else if(foundVector)
310
            file << "<PointData Vectors=\"" << vectorName << "\">\n";
311
312
        else
            return;
313
314
    }

315
316
317
318
319
     /*!
     * \brief Writes the coordinates to the file
     *
     * \param positions Container to store the positions
     */
320
    void writeCoordinates_(std::ostream& file, const std::vector<GlobalPosition>& positions)
321
322
    {
        // write the positions to the file
323
324
        file << "<Points>\n";
        file << "<DataArray type=\"Float32\" Name=\"Coordinates\" NumberOfComponents=\"3\" format=\"ascii\">\n";
325
326
327
        int counter = 0;
        for(auto&& x : positions)
        {
328
            file << x ;
329
330

            if(x.size() == 1)
331
                file << " 0 0 ";
332
            if(x.size() == 2)
333
                file << " 0 ";
334
335
            if(x.size() == 3)
                file << " ";
336

337
            // introduce a line break after a certain time
338
            if((++counter)  > numBeforeLineBreak)
339
            {
340
                file << std::endl;
341
342
343
                counter = 0;
            }
        }
344
345
        file << "\n</DataArray>\n";
        file << "</Points>\n";
346
347
    }

348
     /*!
349
     * \brief Writes data to the file
350
     *
351
     * \param data The data container which hold the data itself, as well as the name of the data set and the number of its components
352
     */
353
    template<class T>
354
    void writeData_(std::ostream& file, const T& data)
355
    {
356
        file << "<DataArray type=\"Float32\" Name=\"" << data.name() << "\" NumberOfComponents=\"" << data.numComponents() << "\" format=\"ascii\">\n";
357
358
        int counter = 0;
        for(auto&& value : data)
359
        {
360
            // forward to specialized function
361
            writeToFile_(file, value);
362
363
364

            // introduce a line break after a certain time
            if((++counter)  > numBeforeLineBreak)
365
            {
366
                file << std::endl;
367
                counter = 0;
368
369
            }
        }
370
        file << "\n</DataArray>\n";
371
    }
372

373
     /*!
374
     * \brief Writes a scalar to the file
375
     *
376
     * \param s The scalar
377
     */
378
    void writeToFile_(std::ostream& file, const Scalar& s)
379
    {
380
        file << s << " ";
381
382
383
384
385
386
387
    }

     /*!
     * \brief Writes a vector to the file
     *
     * \param g The vector
     */
388
    void writeToFile_(std::ostream& file, const GlobalPosition& g)
389
390
391
    {
        assert(g.size() > 1 && g.size() < 4);
        if(g.size() < 3)
392
            file << g << " 0 ";
393
        else
394
            file << g;
395
396
    }

397
    const std::vector<GlobalPosition>& coordinates_;
398
399
    std::list<ScalarFunction> scalarPointData_;
    std::list<VectorFunction> vectorPointData_;
400
401
402
403
};
} // end namespace Dumux

#endif