GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/io/staggeredvtkoutputmodule.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 73 77 94.8%
Functions: 225 357 63.0%
Branches: 59 130 45.4%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup InputOutput
10 * \brief A VTK output module to simplify writing dumux simulation data to VTK format. Specialization for staggered grids with dofs on faces.
11 */
12 #ifndef DUMUX_STAGGERED_VTK_OUTPUT_MODULE_HH
13 #define DUMUX_STAGGERED_VTK_OUTPUT_MODULE_HH
14
15 #include <dune/common/fvector.hh>
16
17 #include <dumux/io/vtkoutputmodule.hh>
18 #include <dumux/io/pointcloudvtkwriter.hh>
19 #include <dumux/io/vtksequencewriter.hh>
20 #include <dumux/discretization/staggered/freeflow/velocityoutput.hh>
21
22 namespace Dumux {
23
24 template<class Scalar, class GlobalPosition>
25 class PointCloudVtkWriter;
26
27 /*!
28 * \ingroup InputOutput
29 * \brief A VTK output module to simplify writing dumux simulation data to VTK format
30 * Specialization for staggered grids with dofs on faces.
31 *
32 * \tparam GridVariables The grid variables
33 * \tparam SolutionVector The solution vector
34 */
35 template<class GridVariables, class SolutionVector>
36 class StaggeredVtkOutputModule
37 : public VtkOutputModule<GridVariables, SolutionVector>
38 {
39 using ParentType = VtkOutputModule<GridVariables, SolutionVector>;
40 using GridGeometry = typename GridVariables::GridGeometry;
41 using GridView = typename GridGeometry::GridView;
42 using Scalar = typename GridVariables::Scalar;
43 using FaceVariables = typename GridVariables::GridFaceVariables::FaceVariables;
44 using FVElementGeometry = typename GridGeometry::LocalView;
45 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
46
47 using Element = typename GridView::template Codim<0>::Entity;
48 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
49
50 struct FaceVarScalarDataInfo { std::function<Scalar(const FaceVariables&)> get; std::string name; };
51 struct FaceVarVectorDataInfo { std::function<GlobalPosition(const SubControlVolumeFace& scvf, const FaceVariables&)> get; std::string name; };
52
53 ✗ struct FaceFieldScalarDataInfo
54 {
55 FaceFieldScalarDataInfo(const std::vector<Scalar>& f, const std::string& n) : data(f), name(n) {}
56 const std::vector<Scalar>& data;
57 const std::string name;
58 };
59
60
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 9 struct FaceFieldVectorDataInfo
61 {
62
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
 18 FaceFieldVectorDataInfo(const std::vector<GlobalPosition>& f, const std::string& n) : data(f), name(n) {}
63 const std::vector<GlobalPosition>& data;
64 const std::string name;
65 };
66
67 public:
68
69 template<class Sol>
70 99 StaggeredVtkOutputModule(const GridVariables& gridVariables,
71 const Sol& sol,
72 const std::string& name,
73 const std::string& paramGroup = "",
74 Dune::VTK::DataMode dm = Dune::VTK::conforming,
75 bool verbose = true)
76 : ParentType(gridVariables, sol, name, paramGroup, dm, verbose)
77 198 , faceWriter_(std::make_shared<PointCloudVtkWriter<Scalar, GlobalPosition>>(coordinates_))
78
3/6
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 75 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 75 times.
✗ Branch 9 not taken.
 198 , sequenceWriter_(faceWriter_, name + "-face", "","",
79
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 99 gridVariables.curGridVolVars().problem().gridGeometry().gridView().comm().rank(),
80
4/8
✓ Branch 2 taken 75 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 75 times.
✓ Branch 7 taken 75 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 75 times.
✗ Branch 11 not taken.
 198 gridVariables.curGridVolVars().problem().gridGeometry().gridView().comm().size() )
81
82 {
83 static_assert(std::is_same<Sol, SolutionVector>::value, "Make sure that sol has the same type as SolutionVector."
84 "Use StaggeredVtkOutputModule<GridVariables, decltype(sol)> when calling the constructor.");
85
86 // enable velocity output per default
87
3/6
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 75 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 75 times.
✗ Branch 8 not taken.
 198 this->addVelocityOutput(std::make_shared<StaggeredFreeFlowVelocityOutput<GridVariables, SolutionVector>>(gridVariables, sol));
88
1/2
✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
 99 writeFaceVars_ = getParamFromGroup<bool>(paramGroup, "Vtk.WriteFaceData", false);
89 99 coordinatesInitialized_ = false;
90 99 }
91
92 //////////////////////////////////////////////////////////////////////////////////////////////
93 //! Methods to conveniently add face variables
94 //! Do not call these methods after initialization
95 //////////////////////////////////////////////////////////////////////////////////////////////
96
97 //! Add a scalar valued field
98 //! \param v The field to be added
99 //! \param name The name of the vtk field
100 void addFaceField(const std::vector<Scalar>& v, const std::string& name)
101 {
102 if (v.size() == this->gridGeometry().gridView().size(1))
103 faceFieldScalarDataInfo_.emplace_back(v, name);
104 else
105 DUNE_THROW(Dune::RangeError, "Size mismatch of added field!");
106 }
107
108 //! Add a vector valued field
109 //! \param v The field to be added
110 //! \param name The name of the vtk field
111 9 void addFaceField(const std::vector<GlobalPosition>& v, const std::string& name)
112 {
113
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 if (v.size() == this->gridGeometry().gridView().size(1))
114 9 faceFieldVectorDataInfo_.emplace_back(v, name);
115 else
116 ✗ DUNE_THROW(Dune::RangeError, "Size mismatch of added field!");
117 9 }
118
119 //! Add a scalar-valued faceVarible
120 //! \param f A function taking a FaceVariables object and returning the desired scalar
121 //! \param name The name of the vtk field
122 50 void addFaceVariable(std::function<Scalar(const FaceVariables&)>&& f, const std::string& name)
123 {
124 50 faceVarScalarDataInfo_.push_back(FaceVarScalarDataInfo{f, name});
125 150 }
126
127 //! Add a vector-valued faceVarible
128 //! \param f A function taking a SubControlVolumeFace and FaceVariables object and returning the desired vector
129 //! \param name The name of the vtk field
130 50 void addFaceVariable(std::function<GlobalPosition(const SubControlVolumeFace& scvf, const FaceVariables&)>&& f, const std::string& name)
131 {
132 50 faceVarVectorDataInfo_.push_back(FaceVarVectorDataInfo{f, name});
133 150 }
134
135 //! Write the values to vtp files
136 //! \param time The current time
137 //! \param type The output type
138 1379 void write(double time, Dune::VTK::OutputType type = Dune::VTK::ascii)
139 {
140 1379 ParentType::write(time, type);
141
2/2
✓ Branch 0 taken 87 times.
✓ Branch 1 taken 1268 times.
 1379 if(writeFaceVars_)
142 111 getFaceDataAndWrite_(time);
143 1379 }
144
145
146 private:
147
148 //! Update the coordinates (the face centers)
149 111 void updateCoordinates_()
150 {
151 111 coordinates_.resize(this->gridGeometry().numFaceDofs());
152
2/2
✓ Branch 2 taken 25167 times.
✓ Branch 3 taken 24975 times.
 103101 for(auto&& facet : facets(this->gridGeometry().gridView()))
153 {
154 51495 const int dofIdxGlobal = this->gridGeometry().gridView().indexSet().index(facet);
155 102990 coordinates_[dofIdxGlobal] = facet.geometry().center();
156 }
157 111 coordinatesInitialized_ = true;
158 111 }
159
160 //! Gathers all face-related data and invokes the face vtk-writer using these data.
161 //! \param time The current time
162 111 void getFaceDataAndWrite_(const Scalar time)
163 {
164 111 const auto numPoints = this->gridGeometry().numFaceDofs();
165
166 // make sure not to iterate over the same dofs twice
167
1/2
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
 111 std::vector<bool> dofVisited(numPoints, false);
168
169 // get fields for all primary coordinates and variables
170
1/2
✓ Branch 0 taken 87 times.
✗ Branch 1 not taken.
 111 if(!coordinatesInitialized_)
171
1/2
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
 111 updateCoordinates_();
172
173 // prepare some containers to store the relevant data
174 111 std::vector<std::vector<Scalar>> faceVarScalarData;
175 111 std::vector<std::vector<GlobalPosition>> faceVarVectorData;
176
177
1/2
✓ Branch 0 taken 87 times.
✗ Branch 1 not taken.
 111 if(!faceVarScalarDataInfo_.empty())
178
2/4
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 87 times.
✗ Branch 5 not taken.
 222 faceVarScalarData.resize(faceVarScalarDataInfo_.size(), std::vector<Scalar>(numPoints));
179
180
1/2
✓ Branch 0 taken 87 times.
✗ Branch 1 not taken.
 111 if(!faceVarVectorDataInfo_.empty())
181
2/4
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 87 times.
✗ Branch 5 not taken.
 222 faceVarVectorData.resize(faceVarVectorDataInfo_.size(), std::vector<GlobalPosition>(numPoints));
182
183
1/2
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
 111 auto fvGeometry = localView(this->gridGeometry());
184
1/2
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
 111 auto elemFaceVars = localView(this->gridVariables().curGridFaceVars());
185
3/8
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23787 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 23700 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 73011 for (const auto& element : elements(this->gridGeometry().gridView(), Dune::Partitions::interior))
186 {
187
1/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 23700 times.
 24300 if (!faceVarScalarDataInfo_.empty() || !faceVarVectorDataInfo_.empty())
188 {
189 24300 fvGeometry.bind(element);
190 24300 elemFaceVars.bindElement(element, fvGeometry, this->sol());
191
192
4/4
✓ Branch 0 taken 44745 times.
✓ Branch 1 taken 50055 times.
✓ Branch 2 taken 94800 times.
✓ Branch 3 taken 23700 times.
 121500 for (auto&& scvf : scvfs(fvGeometry))
193 {
194 97200 const auto dofIdxGlobal = scvf.dofIndex();
195
2/2
✓ Branch 0 taken 44745 times.
✓ Branch 1 taken 50055 times.
 97200 if(dofVisited[dofIdxGlobal])
196 45705 continue;
197
198 51495 dofVisited[dofIdxGlobal] = true;
199
200 51495 const auto& faceVars = elemFaceVars[scvf];
201
202 // get the scalar-valued data
203
2/2
✓ Branch 0 taken 50055 times.
✓ Branch 1 taken 50055 times.
 102990 for (std::size_t i = 0; i < faceVarScalarDataInfo_.size(); ++i)
204
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 50055 times.
 102990 faceVarScalarData[i][dofIdxGlobal] = faceVarScalarDataInfo_[i].get(faceVars);
205
206 // get the vector-valued data
207
2/2
✓ Branch 0 taken 50055 times.
✓ Branch 1 taken 50055 times.
 102990 for (std::size_t i = 0; i < faceVarVectorDataInfo_.size(); ++i)
208
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 50055 times.
 102990 faceVarVectorData[i][dofIdxGlobal] = faceVarVectorDataInfo_[i].get(scvf, faceVars);
209 }
210 }
211 }
212
213 // transfer the data to the point writer
214
1/2
✓ Branch 0 taken 87 times.
✗ Branch 1 not taken.
 111 if(!faceVarScalarDataInfo_.empty())
215
2/2
✓ Branch 0 taken 87 times.
✓ Branch 1 taken 87 times.
 222 for (std::size_t i = 0; i < faceVarScalarDataInfo_.size(); ++i)
216
1/2
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
 111 faceWriter_->addPointData(faceVarScalarData[i], faceVarScalarDataInfo_[i].name);
217
218
1/2
✓ Branch 0 taken 87 times.
✗ Branch 1 not taken.
 111 if(!faceVarVectorDataInfo_.empty())
219
2/2
✓ Branch 0 taken 87 times.
✓ Branch 1 taken 87 times.
 222 for (std::size_t i = 0; i < faceVarVectorDataInfo_.size(); ++i)
220
1/2
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
 111 faceWriter_->addPointData(faceVarVectorData[i], faceVarVectorDataInfo_[i].name);
221
222 // account for the custom fields
223
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 87 times.
 111 for(const auto& field: faceFieldScalarDataInfo_)
224 ✗ faceWriter_->addPointData(field.data, field.name);
225
226
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 87 times.
 111 for(const auto& field: faceFieldVectorDataInfo_)
227 ✗ faceWriter_->addPointData(field.data, field.name);
228
229 // write for the current time step
230
1/2
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
 111 sequenceWriter_.write(time);
231
232 // clear coordinates to save some memory
233
1/2
✓ Branch 0 taken 87 times.
✗ Branch 1 not taken.
 111 coordinates_.clear();
234
1/2
✓ Branch 0 taken 87 times.
✗ Branch 1 not taken.
 111 coordinates_.shrink_to_fit();
235 111 coordinatesInitialized_ = false;
236 111 }
237
238
239 std::shared_ptr<PointCloudVtkWriter<Scalar, GlobalPosition>> faceWriter_;
240
241 VTKSequenceWriter<PointCloudVtkWriter<Scalar, GlobalPosition>> sequenceWriter_;
242
243 bool writeFaceVars_;
244
245 std::vector<GlobalPosition> coordinates_;
246 bool coordinatesInitialized_;
247
248 std::vector<FaceVarScalarDataInfo> faceVarScalarDataInfo_;
249 std::vector<FaceVarVectorDataInfo> faceVarVectorDataInfo_;
250
251 std::vector<FaceFieldScalarDataInfo> faceFieldScalarDataInfo_;
252 std::vector<FaceFieldVectorDataInfo> faceFieldVectorDataInfo_;
253
254 };
255
256 } // end namespace Dumux
257
258 #endif
259