GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/io/vtk/intersectionwriter.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 73 78 93.6%
Functions: 85 127 66.9%
Branches: 119 241 49.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-FileCopyrightInfo: 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 Intersection writer
11 */
12 #ifndef DUMUX_IO_VTK_INTERSECTIONWRITER_HH
13 #define DUMUX_IO_VTK_INTERSECTIONWRITER_HH
14
15 #include <memory>
16 #include <string>
17
18 #include <dune/common/typetraits.hh>
19 #include <dune/grid/io/file/vtk.hh>
20 #include <dune/grid/io/file/vtk/basicwriter.hh>
21 #include <dune/grid/io/file/vtk/function.hh>
22 #include <dune/grid/io/file/vtk/skeletonfunction.hh>
23 #include <dumux/common/parameters.hh>
24
25 namespace Dumux::Detail {
26
27 /*!
28 * \ingroup InputOutput
29 * \brief Iterate over the GridViews boundary intersections
30 * This will visit all intersections for which boundary() is true and
31 * neighbor() is false.
32 */
33 template<typename GV>
34
9/22
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 times.
✗ Branch 21 not taken.
 673 class GridIntersectionIterator
35 : public Dune::ForwardIteratorFacade<GridIntersectionIterator<GV>,
36 const typename GV::Intersection,
37 const typename GV::Intersection&,
38 typename std::iterator_traits<typename GV::template Codim<0>::Iterator>::difference_type>
39 {
40 public:
41 using DerivedType = GridIntersectionIterator<GV>;
42 using Intersection = typename GV::Intersection;
43 using Value = const Intersection;
44 using ElementIterator = typename GV::template Codim<0>::Iterator;
45 using IntersectionIterator = typename GV::IntersectionIterator;
46 using DifferenceType = typename std::iterator_traits<ElementIterator>::difference_type;
47
48 /*!
49 * \brief Construct a GridIntersectionIterator
50 * If end == true, construct an end iterator for the given gridview.
51 * Otherwise, construct a begin iterator.
52 */
53 802 GridIntersectionIterator(const GV& gv, bool end = false)
54
5/6
✓ Branch 0 taken 456 times.
✓ Branch 1 taken 337 times.
✓ Branch 2 taken 465 times.
✓ Branch 3 taken 324 times.
✓ Branch 7 taken 480 times.
✗ Branch 8 not taken.
 1582 : gridView_(gv), eIt_(end ? gridView_.template end<0>() : gridView_.template begin<0>())
55 {
56
5/6
✓ Branch 1 taken 802 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
✓ Branch 4 taken 13 times.
✓ Branch 5 taken 200 times.
✓ Branch 6 taken 280 times.
 1000 if (eIt_ != gridView_.template end<0>())
57
2/4
✓ Branch 1 taken 200 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 200 times.
✗ Branch 5 not taken.
 484 iIt_ = IntersectionIterator(gridView_.ibegin(*eIt_));
58 802 }
59
60 8000 const Intersection& dereference() const
61 {
62 if constexpr (std::is_lvalue_reference_v<decltype(*std::declval<IntersectionIterator>())>)
63
2/4
✓ Branch 7 taken 326780 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 326780 times.
✗ Branch 11 not taken.
 3359148 return *iIt_;
64 else
65 {
66 // Some grids only return temporary intersection objects.
67 // If this is the case, story a copy of the intersection here so that
68 // its address can be safely taken later on.
69 16000 intersection_ = *iIt_;
70 8000 return intersection_;
71 }
72 }
73
74 2451141 bool equals(const DerivedType& other) const
75 {
76
4/4
✓ Branch 0 taken 1274733 times.
✓ Branch 1 taken 1176408 times.
✓ Branch 2 taken 526425 times.
✓ Branch 3 taken 1183208 times.
 4160774 if (eIt_ != other.eIt_)
77 return false;
78
79 // this is a bit tricky, since we may not compare iIt_ if we are
80 // passed-the-end
81
1/2
✓ Branch 1 taken 232597 times.
✗ Branch 2 not taken.
 755805 bool mePassedTheEnd = eIt_ == gridView_.template end<0>();
82
1/2
✓ Branch 1 taken 232597 times.
✗ Branch 2 not taken.
 755805 bool otherPassedTheEnd = other.eIt_ == other.gridView_.template end<0>();
83
84 // both passed-the-end => consider them equal
85
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 755805 times.
 755805 if(mePassedTheEnd && otherPassedTheEnd)
86 return true;
87
88 // one passed the end => not equal
89 ✗ if(mePassedTheEnd || otherPassedTheEnd)
90 return false;
91
92 // none passed-the-end => do their iIt_ iterators match?
93 ✗ return iIt_ == other.iIt_;
94 }
95
96 940136 void increment()
97 {
98
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3600 times.
 940136 ++iIt_;
99
6/6
✓ Branch 0 taken 2700 times.
✓ Branch 1 taken 900 times.
✓ Branch 2 taken 865792 times.
✓ Branch 3 taken 70744 times.
✓ Branch 5 taken 205560 times.
✓ Branch 6 taken 448000 times.
 1883872 if (iIt_ == gridView_.iend(*eIt_))
100 {
101 348848 iIt_ = IntersectionIterator();
102 277204 ++eIt_;
103
3/4
✓ Branch 1 taken 71644 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 206251 times.
✓ Branch 4 taken 209 times.
 278228 if (eIt_ != gridView_.template end<0>())
104 350164 iIt_ = IntersectionIterator(gridView_.ibegin(*eIt_));
105 }
106 940136 }
107
108 private:
109 const GV gridView_;
110 ElementIterator eIt_;
111 IntersectionIterator iIt_;
112 mutable Intersection intersection_;
113 };
114
115 /*!
116 * \ingroup InputOutput
117 * \brief Non conforming intersection iterator factory
118 */
119 template<class GridView>
120 ✗ class NonConformingIntersectionIteratorFactory
121 {
122 public:
123 static constexpr auto dimCell = GridView::dimension-1;
124 using Cell = typename GridView::Intersection;
125 using CellIterator = GridIntersectionIterator<GridView>;
126 using Corner = Dune::VTK::Corner<Cell>;
127 using CornerIterator = Dune::VTK::CornerIterator<CellIterator>;
128 using Point = Corner;
129 using PointIterator = CornerIterator;
130 using ConnectivityWriter = Dune::VTK::NonConformingConnectivityWriter<Cell>;
131 using Communication = typename GridView::Communication;
132
133 explicit NonConformingIntersectionIteratorFactory(const GridView& gv)
134 42 : gridView_(gv) {}
135
136 CellIterator beginCells() const
137
3/5
✓ Branch 1 taken 100 times.
✓ Branch 2 taken 13 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
 197 { return CellIterator(gridView_); }
138
139 CellIterator endCells() const
140
1/2
✓ Branch 2 taken 33 times.
✗ Branch 3 not taken.
 197 { return CellIterator(gridView_, true); }
141
142 136 CornerIterator beginCorners() const
143
3/5
✓ Branch 2 taken 80 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 84 times.
✗ Branch 6 not taken.
 140 { return CornerIterator(beginCells(), endCells()); }
144
145 110 CornerIterator endCorners() const
146
3/6
✓ Branch 2 taken 84 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 13 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 13 times.
✗ Branch 11 not taken.
 140 { return CornerIterator(endCells()); }
147
148 PointIterator beginPoints() const
149
3/5
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 13 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
 68 { return beginCorners(); }
150
151 PointIterator endPoints() const
152
1/2
✓ Branch 2 taken 33 times.
✗ Branch 3 not taken.
 68 { return endCorners(); }
153
154 ✗ ConnectivityWriter makeConnectivity() const
155
2/4
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
 42 { return ConnectivityWriter(); }
156
157 const Communication& comm() const
158
21/32
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 24 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 32 times.
✓ Branch 10 taken 16 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 16 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 24 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 18 taken 12 times.
✓ Branch 19 taken 12 times.
✓ Branch 20 taken 8 times.
✓ Branch 21 taken 16 times.
✓ Branch 22 taken 4 times.
✓ Branch 23 taken 8 times.
✓ Branch 24 taken 2 times.
✓ Branch 25 taken 6 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 13 times.
✓ Branch 28 taken 8 times.
✓ Branch 29 taken 6 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 1 times.
 209 { return gridView_.comm(); }
159
160 private:
161 GridView gridView_;
162 };
163
164 /*!
165 * \ingroup InputOutput
166 * \brief Skeleton function for intersection writer
167 */
168 template<class GridView, class Mapper, class F>
169
3/18
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
 3 class SkeletonFunction
170 {
171 using Intersection = typename GridView::Intersection;
172 public:
173 using Traits = Dune::VTK::SkeletonFunctionTraits<GridView, typename GridView::ctype>;
174
175 55 SkeletonFunction(const GridView& gv, const Mapper& mapper, const F& field)
176 : gv_(gv)
177 , mapper_(mapper)
178 , field_(field)
179 41 , components_(1)
180 {
181 if constexpr (std::is_invocable_v<F, Intersection, int>)
182 {
183 if constexpr (Dune::IsIndexable<std::decay_t<decltype(field(std::declval<Intersection>(), 0))>>{})
184 {
185 if constexpr (Dune::IsIndexable<std::decay_t<decltype(field(std::declval<Intersection>(), 0)[0])>>{})
186 DUNE_THROW(Dune::InvalidStateException, "Invalid field type");
187 else
188 {
189 const auto& is = *(intersections(gv, *(elements(gv).begin())).begin());
190 components_ = field(is, mapper_(is, gv_)).size();
191 }
192 }
193 }
194 else if constexpr (Dune::IsIndexable<std::decay_t<decltype(field[0])>>{})
195 {
196
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 18 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
 67 assert(field.size() == gv.size(1));
197 if constexpr (Dune::IsIndexable<std::decay_t<decltype(field[0][0])>>{})
198 DUNE_THROW(Dune::InvalidStateException, "Invalid field type");
199 else
200 33 components_ = field[0].size();
201 }
202 33 }
203
204 //! return number of components
205 71 unsigned dimRange() const { return components_; }
206
207 562936 void evaluate(const typename Traits::Cell& intersection,
208 const typename Traits::Domain& xl,
209 typename Traits::Range& result) const
210 {
211
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 196068 times.
 562936 assert(intersection.conforming());
212
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 562936 result.resize(components_);
213
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 562936 const auto idx = mapper_(intersection, gv_);
214
215 562936 auto accessEntry = [&](auto i)
216 {
217 if constexpr (std::is_invocable_v<F, Intersection, int>)
218 {
219 if constexpr (Dune::IsIndexable<std::decay_t<decltype(field_(intersection, idx))>>{})
220 return field_(intersection, idx)[i];
221 else
222 188800 return field_(intersection, idx);
223 }
224 else
225 {
226 if constexpr (Dune::IsIndexable<std::decay_t<decltype(std::declval<F>()[0])>>{})
227 377200 return field_[idx][i];
228 else
229 371072 return field_[idx];
230 }
231 };
232
233
2/2
✓ Branch 0 taken 375768 times.
✓ Branch 1 taken 281468 times.
 1314472 for (int i = 0; i < components_; ++i)
234
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 2065808 result[i] = accessEntry(i);
235 562936 }
236
237 private:
238 GridView gv_;
239 const Mapper& mapper_;
240
241 // If F is callable we store a copy of it, otherwise we assume that F is a container
242 // stored somewhere else, thus we only keep a reference here
243 std::conditional_t<std::is_invocable_v<F,Intersection, int>, const F, const F&> field_;
244
245 std::size_t components_;
246 };
247
248 } // end namespace Dumux::Detail
249
250 namespace Dumux {
251 /*!
252 * \ingroup InputOutput
253 * \brief Conforming intersection writer
254 */
255 template<class GridView>
256
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 20 class ConformingIntersectionWriter
257 : public Detail::NonConformingIntersectionIteratorFactory<GridView>
258 , public Dune::VTK::BasicWriter<Detail::NonConformingIntersectionIteratorFactory<GridView>>
259 {
260 using Factory = Detail::NonConformingIntersectionIteratorFactory<GridView>;
261 using Base = Dune::VTK::BasicWriter<Factory>;
262
263 public:
264 21 ConformingIntersectionWriter(const GridView& gridView, const std::string& paramGroup = "")
265
3/8
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 61 : Factory(gridView), Base(static_cast<const Factory&>(*this)), gridView_(gridView)
266 {
267
3/6
✓ Branch 0 taken 21 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 21 times.
✗ Branch 7 not taken.
 21 static bool addProcessRank = getParamFromGroup<bool>(paramGroup, "Vtk.AddProcessRank");
268
1/2
✓ Branch 0 taken 21 times.
✗ Branch 1 not taken.
 21 if (addProcessRank)
269 {
270
4/6
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 22 auto getRank = [rank = gridView_.comm().rank()](const auto& is, const auto idx)
271 ✗ { return rank; };
272
273 21 auto mapper = getStandardMapper();
274 using SF = Detail::SkeletonFunction<GridView, decltype(mapper), decltype(getRank)>;
275
1/4
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 42 auto fun = std::make_shared<SF>(gridView_, mapper, getRank);
276
6/14
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 21 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 21 times.
✓ Branch 11 taken 21 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 21 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
 42 addCellData(fun, "processRank");
277 }
278 21 }
279
280 using Base::addCellData;
281
282 static auto getStandardMapper()
283 {
284
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 111112 times.
✓ Branch 4 taken 261424 times.
✗ Branch 5 not taken.
 749872 return [](const auto& is, const GridView& gridView){ return gridView.indexSet().subIndex(is.inside(), is.indexInInside(), 1); };
285 }
286
287 template<class F, class Mapper = decltype(getStandardMapper())>
288 auto makeSkeletonFunction(const F& f, const Mapper& mapper = getStandardMapper()) const
289 {
290
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
 39 return std::make_shared<Detail::SkeletonFunction<GridView, decltype(mapper), decltype(f)>>(gridView_, mapper, f);
291 }
292
293 template<class Func>
294 140 void addCellData(const std::shared_ptr<Func>& p, const std::string& name)
295 {
296
4/14
✓ Branch 2 taken 71 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 71 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 71 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 71 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
 420 addCellData(std::shared_ptr<typename Base::FunctionWriter>
297
2/4
✓ Branch 1 taken 71 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 71 times.
✗ Branch 5 not taken.
 140 (new Dune::VTK::SkeletonFunctionWriter<Func>(p, name)));
298 140 }
299
300 template<class Func>
301 void addCellData(Func* p, const std::string& name)
302 {
303 addCellData(std::shared_ptr<Func>(p), name);
304 }
305
306 template<class F>
307 97 void addField(const F& field, const std::string& name)
308 {
309 97 auto mapper = getStandardMapper();
310
2/6
✓ Branch 2 taken 50 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 97 addCellData(makeSkeletonFunction(field, mapper), name);
311 97 }
312
313 using Base::addPointData;
314
315 template<class Func>
316 void addPointData(const std::shared_ptr<Func>& p, const std::string& name)
317 {
318 addPointData(std::shared_ptr<typename Base::FunctionWriter>
319 (new Dune::VTK::SkeletonFunctionWriter<Func>(p, name)));
320 }
321
322 template<class Func>
323 void addPointData(Func* p, const std::string& name)
324 {
325 addPointData(std::shared_ptr<Func>(p), name);
326 }
327
328 std::string write(const std::string& name, Dune::VTK::OutputType outputType = Dune::VTK::OutputType::ascii)
329 {
330
2/4
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
 34 return Base::write(name, outputType);
331 }
332
333 private:
334 GridView gridView_;
335 };
336
337 } // end namespace Dumux
338
339 #endif
340