staggeredvtkoutputmodule.hh 11.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
// -*- 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
 * \brief A VTK output module to simplify writing dumux simulation data to VTK format. Specialization for staggered grids with dofs on faces.
 */
#ifndef STAGGERED_VTK_OUTPUT_MODULE_HH
#define STAGGERED_VTK_OUTPUT_MODULE_HH

#include <dune/common/fvector.hh>

#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/pointcloudvtkwriter.hh>
#include <dumux/io/vtksequencewriter.hh>


namespace Dumux
{

/*!
 * \ingroup InputOutput
 * \brief A VTK output module to simplify writing dumux simulation data to VTK format
 *        Specialization for staggered grids with dofs on faces.
 */
template<typename TypeTag>
class StaggeredVtkOutputModule : public VtkOutputModule<TypeTag>
{
    friend class VtkOutputModule<TypeTag>;
    using ParentType = VtkOutputModule<TypeTag>;
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
    using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
    using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);


    enum { dim = GridView::dimension };
    enum { dimWorld = GridView::dimensionworld };

    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;

    using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
    typename DofTypeIndices::CellCenterIdx cellCenterIdx;
    typename DofTypeIndices::FaceIdx faceIdx;

    struct FaceVarScalarDataInfo { std::function<Scalar(const FaceVariables&)> get; std::string name; };
    struct FaceVarVectorDataInfo { std::function<GlobalPosition(const SubControlVolumeFace& scvf, const FaceVariables&)> get; std::string name; };

68
69
70
71
72
73
74
75
76
77
78
79
80
    struct FaceFieldScalarDataInfo
    {
        FaceFieldScalarDataInfo(const std::vector<Scalar>& f, const std::string& n) : data(f), name(n) {}
        const std::vector<Scalar>& data;
        const std::string name;
    };

    struct FaceFieldVectorDataInfo
    {
        FaceFieldVectorDataInfo(const std::vector<GlobalPosition>& f, const std::string& n) : data(f), name(n) {}
        const std::vector<GlobalPosition>& data;
        const std::string name;
    };
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106

public:

    StaggeredVtkOutputModule(const Problem& problem,
                    const FVGridGeometry& fvGridGeometry,
                    const GridVariables& gridVariables,
                    const SolutionVector& sol,
                    const std::string& name,
                    bool verbose = true,
                    Dune::VTK::DataMode dm = Dune::VTK::conforming)
    : ParentType(problem, fvGridGeometry, gridVariables, sol, name, verbose, dm)
    , problem_(problem)
    , gridGeom_(fvGridGeometry)
    , gridVariables_(gridVariables)
    , sol_(sol)
    , faceWriter_(std::make_shared<PointCloudVtkWriter<Scalar, dim>>(coordinates_))
    , sequenceWriter_(faceWriter_, problem.name() + "-face", "","",
                      fvGridGeometry.gridView().comm().rank(),
                      fvGridGeometry.gridView().comm().size() )

    {
        writeFaceVars_ = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Vtk.WriteFaceData", false);
        coordinatesInitialized_ = false;
    }

    //////////////////////////////////////////////////////////////////////////////////////////////
107
    //! Methods to conveniently add face variables
108
109
110
    //! Do not call these methods after initialization
    //////////////////////////////////////////////////////////////////////////////////////////////

111
112
113
114
    //! Add a scalar valued field
    //! \param v The field to be added
    //! \param name The name of the vtk field
    void addFaceField(const std::vector<Scalar>& v, const std::string& name)
115
116
    {
        if (v.size() == this->gridGeom_.gridView().size(1))
117
118
119
120
            faceFieldScalarDataInfo_.emplace_back(v, name);
        else
            DUNE_THROW(Dune::RangeError, "Size mismatch of added field!");
    }
121

122
123
124
125
126
127
128
    //! Add a vector valued field
    //! \param v The field to be added
    //! \param name The name of the vtk field
    void addFaceField(const std::vector<GlobalPosition>& v, const std::string& name)
    {
        if (v.size() == this->gridGeom_.gridView().size(1))
                faceFieldVectorDataInfo_.emplace_back(v, name);
129
130
131
132
        else
            DUNE_THROW(Dune::RangeError, "Size mismatch of added field!");
    }

133
134
135
    //! Add a scalar-valued faceVarible
    //! \param f A function taking a FaceVariables object and returning the desired scalar
    //! \param name The name of the vtk field
136
137
138
139
140
    void addFaceVariable(std::function<Scalar(const FaceVariables&)>&& f, const std::string& name)
    {
        faceVarScalarDataInfo_.push_back(FaceVarScalarDataInfo{f, name});
    }

141
142
143
    //! Add a vector-valued faceVarible
    //! \param f A function taking a SubControlVolumeFace and FaceVariables object and returning the desired vector
    //! \param name The name of the vtk field
144
145
146
147
148
    void addFaceVariable(std::function<GlobalPosition(const SubControlVolumeFace& scvf, const FaceVariables&)>&& f, const std::string& name)
    {
        faceVarVectorDataInfo_.push_back(FaceVarVectorDataInfo{f, name});
    }

149
150
151
    //! Write the values to vtp files
    //! \param time The current time
    //! \param type The output type
152
153
154
155
156
157
158
159
160
161
    void write(double time, Dune::VTK::OutputType type = Dune::VTK::ascii)
    {
        ParentType::write(time, type);
        if(writeFaceVars_)
            getFaceDataAndWrite_(time);
    }


private:

162
    //! Update the coordinates (the face centers)
163
164
165
166
167
168
169
170
171
172
173
    void updateCoordinates_()
    {
        coordinates_.resize(gridGeom_.numFaceDofs());
        for(auto&& facet : facets(gridGeom_.gridView()))
        {
            const int dofIdxGlobal = gridGeom_.gridView().indexSet().index(facet);
            coordinates_[dofIdxGlobal] = facet.geometry().center();
        }
        coordinatesInitialized_ = true;
    }

174
175
     //! Gathers all face-related data and invokes the face vtk-writer using these data.
     //! \param time The current time
176
177
178
179
180
181
182
183
184
185
186
    void getFaceDataAndWrite_(const Scalar time)
    {
        const auto numPoints = gridGeom_.numFaceDofs();

        // make sure not to iterate over the same dofs twice
        std::vector<bool> dofVisited(numPoints, false);

        // get fields for all primary coordinates and variables
        if(!coordinatesInitialized_)
            updateCoordinates_();

187
        // prepare some containers to store the relevant data
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
        std::vector<std::vector<Scalar>> faceVarScalarData;
        std::vector<std::vector<GlobalPosition>> faceVarVectorData;

        if(!faceVarScalarDataInfo_.empty())
            faceVarScalarData.resize(faceVarScalarDataInfo_.size(), std::vector<Scalar>(numPoints));

        if(!faceVarVectorDataInfo_.empty())
            faceVarVectorData.resize(faceVarVectorDataInfo_.size(), std::vector<GlobalPosition>(numPoints));

        for (const auto& element : elements(gridGeom_.gridView(), Dune::Partitions::interior))
        {
            auto fvGeometry = localView(gridGeom_);
            auto elemFaceVars = localView(gridVariables_.curGridFaceVars());

            if (!faceVarScalarDataInfo_.empty() || !faceVarVectorDataInfo_.empty())
            {
                fvGeometry.bind(element);
                elemFaceVars.bindElement(element, fvGeometry, sol_);

                for (auto&& scvf : scvfs(fvGeometry))
                {
                    const auto dofIdxGlobal = scvf.dofIndex();
                    if(dofVisited[dofIdxGlobal])
                        continue;

                    dofVisited[dofIdxGlobal] = true;

                    const auto& faceVars = elemFaceVars[scvf];

217
                    // get the scalar-valued data
218
219
220
                    for (std::size_t i = 0; i < faceVarScalarDataInfo_.size(); ++i)
                        faceVarScalarData[i][dofIdxGlobal] = faceVarScalarDataInfo_[i].get(faceVars);

221
                    // get the vector-valued data
222
223
224
225
226
227
                    for (std::size_t i = 0; i < faceVarVectorDataInfo_.size(); ++i)
                            faceVarVectorData[i][dofIdxGlobal] = faceVarVectorDataInfo_[i].get(scvf, faceVars);
                }
            }
        }

228
        // transfer the data to the point writer
229
230
231
232
233
234
        if(!faceVarScalarDataInfo_.empty())
            for (std::size_t i = 0; i < faceVarScalarDataInfo_.size(); ++i)
                faceWriter_->addPointData(faceVarScalarData[i], faceVarScalarDataInfo_[i].name);

        if(!faceVarVectorDataInfo_.empty())
            for (std::size_t i = 0; i < faceVarVectorDataInfo_.size(); ++i)
235
                faceWriter_->addPointData(faceVarVectorData[i], faceVarVectorDataInfo_[i].name);
236

237
238
239
240
241
242
243
244
        // account for the custom fields
        for(const auto& field: faceFieldScalarDataInfo_)
            faceWriter_->addPointData(field.data, field.name);

        for(const auto& field: faceFieldVectorDataInfo_)
            faceWriter_->addPointData(field.data, field.name);

        // write for the current time step
245
        sequenceWriter_.write(time);
246
247

        // clear coordinates to save some memory
248
249
250
        coordinates_.clear();
        coordinates_.shrink_to_fit();
        coordinatesInitialized_ = false;
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
    }


    const Problem& problem_;
    const FVGridGeometry& gridGeom_;
    const GridVariables& gridVariables_;
    const SolutionVector& sol_;

    std::shared_ptr<PointCloudVtkWriter<Scalar, dim>> faceWriter_;

    VTKSequenceWriter<PointCloudVtkWriter<Scalar, dim>> sequenceWriter_;

    bool writeFaceVars_;

    std::vector<GlobalPosition> coordinates_;
    bool coordinatesInitialized_;

    std::vector<FaceVarScalarDataInfo> faceVarScalarDataInfo_;
    std::vector<FaceVarVectorDataInfo> faceVarVectorDataInfo_;

271
272
273
    std::vector<FaceFieldScalarDataInfo> faceFieldScalarDataInfo_;
    std::vector<FaceFieldVectorDataInfo> faceFieldVectorDataInfo_;

274
275
276
277
278
};

} // end namespace Dumux

#endif