GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/io/vtk/intersectionwriter.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 87 90 96.7%
Functions: 103 103 100.0%
Branches: 109 190 57.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 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 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 times.
✗ 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 31 taken 2 times.
✗ Branch 32 not taken.
✓ Branch 44 taken 2 times.
✗ Branch 45 not taken.
✓ Branch 55 taken 2 times.
✗ Branch 56 not taken.
 136 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 926 GridIntersectionIterator(const GV& gv, bool end = false)
54
5/7
✓ Branch 0 taken 530 times.
✓ Branch 1 taken 387 times.
✓ Branch 5 taken 480 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 480 times.
✗ Branch 9 not taken.
✓ Branch 2 taken 9 times.
 1786 : gridView_(gv), eIt_(end ? gridView_.template end<0>() : gridView_.template begin<0>())
55 {
56
5/5
✓ Branch 1 taken 882 times.
✓ Branch 2 taken 44 times.
✓ Branch 3 taken 27 times.
✓ Branch 4 taken 239 times.
✓ Branch 5 taken 280 times.
 1216 if (eIt_ != gridView_.template end<0>())
57
2/4
✓ Branch 1 taken 218 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 200 times.
✗ Branch 5 not taken.
 410 iIt_ = IntersectionIterator(gridView_.ibegin(*eIt_));
58 926 }
59
60 1360692 const Intersection& dereference() const
61 {
62 if constexpr (std::is_lvalue_reference_v<decltype(*std::declval<IntersectionIterator>())>)
63
1/2
✓ Branch 4 taken 329724 times.
✗ Branch 5 not taken.
 2048836 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 2650119 bool equals(const DerivedType& other) const
75 {
76
4/4
✓ Branch 0 taken 2126911 times.
✓ Branch 1 taken 523208 times.
✓ Branch 2 taken 19312 times.
✓ Branch 3 taken 9139 times.
 3311018 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 298743 times.
✗ Branch 2 not taken.
 821951 bool mePassedTheEnd = eIt_ == gridView_.template end<0>();
82
1/2
✓ Branch 1 taken 298743 times.
✗ Branch 2 not taken.
 821951 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 821951 times.
 821951 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 1006920 void increment()
97 {
98
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3600 times.
 1006920 ++iIt_;
99
6/6
✓ Branch 2 taken 913468 times.
✓ Branch 3 taken 4968 times.
✓ Branch 5 taken 205560 times.
✓ Branch 6 taken 448000 times.
✓ Branch 1 taken 86684 times.
✓ Branch 0 taken 2700 times.
 1754364 if (iIt_ == gridView_.iend(*eIt_))
100 {
101 293900 iIt_ = IntersectionIterator();
102 293900 ++eIt_;
103
4/4
✓ Branch 1 taken 86684 times.
✓ Branch 2 taken 207016 times.
✓ Branch 3 taken 2729 times.
✓ Branch 4 taken 27 times.
 296612 if (eIt_ != gridView_.template end<0>())
104 296190 iIt_ = IntersectionIterator(gridView_.ibegin(*eIt_));
105 }
106 1006920 }
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 27 explicit NonConformingIntersectionIteratorFactory(const GridView& gv)
134 28 : gridView_(gv) {}
135
136 383 CellIterator beginCells() const
137
3/5
✓ Branch 1 taken 108 times.
✓ Branch 2 taken 17 times.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✗ Branch 3 not taken.
 223 { return CellIterator(gridView_); }
138
139 543 CellIterator endCells() const
140
1/2
✓ Branch 2 taken 39 times.
✗ Branch 3 not taken.
 223 { return CellIterator(gridView_, true); }
141
142 160 CornerIterator beginCorners() const
143
3/5
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 92 times.
✗ Branch 5 not taken.
✓ Branch 1 taken 88 times.
 172 { return CornerIterator(beginCells(), endCells()); }
144
145
1/2
✓ Branch 3 taken 17 times.
✗ Branch 4 not taken.
 160 CornerIterator endCorners() const
146
2/4
✓ Branch 5 taken 17 times.
✗ Branch 6 not taken.
✓ Branch 1 taken 92 times.
✗ Branch 2 not taken.
 172 { return CornerIterator(endCells()); }
147
148 80 PointIterator beginPoints() const
149
3/5
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 17 times.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✗ Branch 3 not taken.
 80 { return beginCorners(); }
150
151 80 PointIterator endPoints() const
152
1/2
✓ Branch 2 taken 39 times.
✗ Branch 3 not taken.
 80 { return endCorners(); }
153
154 27 ConnectivityWriter makeConnectivity() const
155
1/2
✓ Branch 0 taken 26 times.
✗ Branch 1 not taken.
 27 { return ConnectivityWriter(); }
156
157 214 const Communication& comm() const
158
20/30
✗ Branch 0 not taken.
✓ Branch 1 taken 28 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 24 times.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 32 times.
✓ Branch 10 taken 16 times.
✓ Branch 12 taken 16 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 25 times.
✗ Branch 17 not taken.
✓ Branch 18 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 17 times.
✓ Branch 28 taken 8 times.
✓ Branch 29 taken 9 times.
✗ Branch 6 not taken.
✗ Branch 16 not taken.
✓ Branch 19 taken 12 times.
✗ Branch 11 not taken.
 214 { 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/6
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 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 85 SkeletonFunction(const GridView& gv, const Mapper& mapper, const F& field)
176 30 : gv_(gv)
177 85 , mapper_(mapper)
178 85 , field_(field)
179
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 28 , 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
2/3
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 23 times.
 24 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 24 components_ = field[0].size();
201 }
202 24 }
203
204 //! return number of components
205 85 unsigned dimRange() const { return components_; }
206
207 781968 void evaluate(const typename Traits::Cell& intersection,
208 const typename Traits::Domain& xl,
209 typename Traits::Range& result) const
210 {
211
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 263632 times.
 527264 assert(intersection.conforming());
212 781968 result.resize(components_);
213 781968 const auto idx = mapper_(intersection, gv_);
214
215 1241784 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 389848 return field_(intersection, idx);
223 }
224 else
225 {
226 if constexpr (Dune::IsIndexable<std::decay_t<decltype(std::declval<F>()[0])>>{})
227 410224 return field_[idx][i];
228 else
229 187008 return field_[idx];
230 }
231 };
232
233
3/3
✓ Branch 1 taken 427448 times.
✓ Branch 2 taken 127352 times.
✓ Branch 0 taken 329724 times.
 1769048 for (int i = 0; i < components_; ++i)
234 987080 result[i] = accessEntry(i);
235 781968 }
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 6 times.
✗ Branch 2 not taken.
 26 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 27 ConformingIntersectionWriter(const GridView& gridView, const std::string& paramGroup = "")
265
1/2
✓ Branch 0 taken 26 times.
✗ Branch 1 not taken.
 27 : Factory(gridView), Base(static_cast<const Factory&>(*this)), gridView_(gridView)
266 {
267
3/6
✓ Branch 0 taken 27 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 27 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 27 times.
✗ Branch 7 not taken.
 27 static bool addProcessRank = getParamFromGroup<bool>(paramGroup, "Vtk.AddProcessRank");
268
1/2
✓ Branch 0 taken 27 times.
✗ Branch 1 not taken.
 27 if (addProcessRank)
269 {
270
3/5
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
✓ Branch 3 taken 14 times.
✗ Branch 4 not taken.
✓ Branch 2 taken 12 times.
 27 auto getRank = [rank = gridView_.comm().rank()](const auto& is, const auto idx)
271 102656 { return rank; };
272
273 auto mapper = getStandardMapper();
274 using SF = Detail::SkeletonFunction<GridView, decltype(mapper), decltype(getRank)>;
275
1/2
✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
 27 auto fun = std::make_shared<SF>(gridView_, mapper, getRank);
276
3/8
✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 27 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 54 addCellData(fun, "processRank");
277 27 }
278 27 }
279
280 using Base::addCellData;
281
282
1/2
✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
 27 static auto getStandardMapper()
283 {
284
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 128360 times.
✓ Branch 4 taken 261424 times.
✗ Branch 5 not taken.
 394392 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 58 auto makeSkeletonFunction(const F& f, const Mapper& mapper = getStandardMapper()) const
289 {
290 116 return std::make_shared<Detail::SkeletonFunction<GridView, decltype(mapper), decltype(f)>>(gridView_, mapper, f);
291 }
292
293 template<class Func>
294 168 void addCellData(const std::shared_ptr<Func>& p, const std::string& name)
295 {
296
4/8
✓ Branch 1 taken 85 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 85 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 85 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 85 times.
✗ Branch 9 not taken.
 672 addCellData(std::shared_ptr<typename Base::FunctionWriter>
297
1/4
✓ Branch 3 taken 85 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 168 (new Dune::VTK::SkeletonFunctionWriter<Func>(p, name)));
298 168 }
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 109 void addField(const F& field, const std::string& name)
308 {
309 109 auto mapper = getStandardMapper();
310
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 109 addCellData(makeSkeletonFunction(field, mapper), name);
311 109 }
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 40 std::string write(const std::string& name, Dune::VTK::OutputType outputType = Dune::VTK::OutputType::ascii)
329 {
330
2/4
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
 40 return Base::write(name, outputType);
331 }
332
333 private:
334 GridView gridView_;
335 };
336
337 } // end namespace Dumux
338
339 #endif
340