version 3.8
intersectionwriter.hh
Go to the documentation of this file.
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//
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>
24
25namespace Dumux::Detail {
26
33template<typename GV>
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{
40public:
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
53 GridIntersectionIterator(const GV& gv, bool end = false)
54 : gridView_(gv), eIt_(end ? gridView_.template end<0>() : gridView_.template begin<0>())
55 {
56 if (eIt_ != gridView_.template end<0>())
57 iIt_ = IntersectionIterator(gridView_.ibegin(*eIt_));
58 }
59
61 {
62 if constexpr (std::is_lvalue_reference_v<decltype(*std::declval<IntersectionIterator>())>)
63 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 intersection_ = *iIt_;
70 return intersection_;
71 }
72 }
73
74 bool equals(const DerivedType& other) const
75 {
76 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 bool mePassedTheEnd = eIt_ == gridView_.template end<0>();
82 bool otherPassedTheEnd = other.eIt_ == other.gridView_.template end<0>();
83
84 // both passed-the-end => consider them equal
85 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 void increment()
97 {
98 ++iIt_;
99 if (iIt_ == gridView_.iend(*eIt_))
100 {
101 iIt_ = IntersectionIterator();
102 ++eIt_;
103 if (eIt_ != gridView_.template end<0>())
104 iIt_ = IntersectionIterator(gridView_.ibegin(*eIt_));
105 }
106 }
107
108private:
109 const GV gridView_;
110 ElementIterator eIt_;
112 mutable Intersection intersection_;
113};
114
119template<class GridView>
121{
122public:
123 static constexpr auto dimCell = GridView::dimension-1;
124 using Cell = typename GridView::Intersection;
126 using Corner = Dune::VTK::Corner<Cell>;
127 using CornerIterator = Dune::VTK::CornerIterator<CellIterator>;
128 using Point = Corner;
130 using ConnectivityWriter = Dune::VTK::NonConformingConnectivityWriter<Cell>;
131 using Communication = typename GridView::CollectiveCommunication;
132
133 explicit NonConformingIntersectionIteratorFactory(const GridView& gv)
134 : gridView_(gv) {}
135
137 { return CellIterator(gridView_); }
138
140 { return CellIterator(gridView_, true); }
141
143 { return CornerIterator(beginCells(), endCells()); }
144
146 { return CornerIterator(endCells()); }
147
149 { return beginCorners(); }
150
152 { return endCorners(); }
153
155 { return ConnectivityWriter(); }
156
157 const Communication& comm() const
158 { return gridView_.comm(); }
159
160private:
161 GridView gridView_;
162};
163
168template<class GridView, class Mapper, class F>
170{
171 using Intersection = typename GridView::Intersection;
172public:
173 using Traits = Dune::VTK::SkeletonFunctionTraits<GridView, typename GridView::ctype>;
174
175 SkeletonFunction(const GridView& gv, const Mapper& mapper, const F& field)
176 : gv_(gv)
177 , mapper_(mapper)
178 , field_(field)
179 , 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 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 components_ = field[0].size();
201 }
202 }
203
205 unsigned dimRange() const { return components_; }
206
207 void evaluate(const typename Traits::Cell& intersection,
208 const typename Traits::Domain& xl,
209 typename Traits::Range& result) const
210 {
211 assert(intersection.conforming());
212 result.resize(components_);
213 const auto idx = mapper_(intersection, gv_);
214
215 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 return field_(intersection, idx);
223 }
224 else
225 {
226 if constexpr (Dune::IsIndexable<std::decay_t<decltype(std::declval<F>()[0])>>{})
227 return field_[idx][i];
228 else
229 return field_[idx];
230 }
231 };
232
233 for (int i = 0; i < components_; ++i)
234 result[i] = accessEntry(i);
235 }
236
237private:
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
250namespace Dumux {
255template<class GridView>
258, public Dune::VTK::BasicWriter<Detail::NonConformingIntersectionIteratorFactory<GridView>>
259{
261 using Base = Dune::VTK::BasicWriter<Factory>;
262
263public:
264 ConformingIntersectionWriter(const GridView& gridView, const std::string& paramGroup = "")
265 : Factory(gridView), Base(static_cast<const Factory&>(*this)), gridView_(gridView)
266 {
267 static bool addProcessRank = getParamFromGroup<bool>(paramGroup, "Vtk.AddProcessRank");
268 if (addProcessRank)
269 {
270 auto getRank = [rank = gridView_.comm().rank()](const auto& is, const auto idx)
271 { return rank; };
272
273 auto mapper = getStandardMapper();
274 using SF = Detail::SkeletonFunction<GridView, decltype(mapper), decltype(getRank)>;
275 auto fun = std::make_shared<SF>(gridView_, mapper, getRank);
276 addCellData(fun, "processRank");
277 }
278 }
279
280 using Base::addCellData;
281
282 static auto getStandardMapper()
283 {
284 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 return std::make_shared<Detail::SkeletonFunction<GridView, decltype(mapper), decltype(f)>>(gridView_, mapper, f);
291 }
292
293 template<class Func>
294 void addCellData(const std::shared_ptr<Func>& p, const std::string& name)
295 {
296 addCellData(std::shared_ptr<typename Base::FunctionWriter>
297 (new Dune::VTK::SkeletonFunctionWriter<Func>(p, name)));
298 }
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 void addField(const F& field, const std::string& name)
308 {
309 auto mapper = getStandardMapper();
310 addCellData(makeSkeletonFunction(field, mapper), name);
311 }
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 return Base::write(name, outputType);
331 }
332
333private:
334 GridView gridView_;
335};
336
337} // end namespace Dumux
338
339#endif
Conforming intersection writer.
Definition: intersectionwriter.hh:259
void addCellData(const std::shared_ptr< Func > &p, const std::string &name)
Definition: intersectionwriter.hh:294
void addPointData(const std::shared_ptr< Func > &p, const std::string &name)
Definition: intersectionwriter.hh:316
ConformingIntersectionWriter(const GridView &gridView, const std::string &paramGroup="")
Definition: intersectionwriter.hh:264
static auto getStandardMapper()
Definition: intersectionwriter.hh:282
void addCellData(Func *p, const std::string &name)
Definition: intersectionwriter.hh:301
auto makeSkeletonFunction(const F &f, const Mapper &mapper=getStandardMapper()) const
Definition: intersectionwriter.hh:288
void addField(const F &field, const std::string &name)
Definition: intersectionwriter.hh:307
void addPointData(Func *p, const std::string &name)
Definition: intersectionwriter.hh:323
std::string write(const std::string &name, Dune::VTK::OutputType outputType=Dune::VTK::OutputType::ascii)
Definition: intersectionwriter.hh:328
Iterate over the GridViews boundary intersections This will visit all intersections for which boundar...
Definition: intersectionwriter.hh:39
typename GV::Intersection Intersection
Definition: intersectionwriter.hh:42
typename std::iterator_traits< ElementIterator >::difference_type DifferenceType
Definition: intersectionwriter.hh:46
const Intersection & dereference() const
Definition: intersectionwriter.hh:60
typename GV::IntersectionIterator IntersectionIterator
Definition: intersectionwriter.hh:45
GridIntersectionIterator(const GV &gv, bool end=false)
Construct a GridIntersectionIterator If end == true, construct an end iterator for the given gridview...
Definition: intersectionwriter.hh:53
void increment()
Definition: intersectionwriter.hh:96
typename GV::template Codim< 0 >::Iterator ElementIterator
Definition: intersectionwriter.hh:44
const Intersection Value
Definition: intersectionwriter.hh:43
bool equals(const DerivedType &other) const
Definition: intersectionwriter.hh:74
Non conforming intersection iterator factory.
Definition: intersectionwriter.hh:121
PointIterator beginPoints() const
Definition: intersectionwriter.hh:148
CornerIterator PointIterator
Definition: intersectionwriter.hh:129
Corner Point
Definition: intersectionwriter.hh:128
PointIterator endPoints() const
Definition: intersectionwriter.hh:151
static constexpr auto dimCell
Definition: intersectionwriter.hh:123
GridIntersectionIterator< GridView > CellIterator
Definition: intersectionwriter.hh:125
Dune::VTK::Corner< Cell > Corner
Definition: intersectionwriter.hh:126
Dune::VTK::NonConformingConnectivityWriter< Cell > ConnectivityWriter
Definition: intersectionwriter.hh:130
typename GridView::CollectiveCommunication Communication
Definition: intersectionwriter.hh:131
NonConformingIntersectionIteratorFactory(const GridView &gv)
Definition: intersectionwriter.hh:133
CornerIterator endCorners() const
Definition: intersectionwriter.hh:145
CornerIterator beginCorners() const
Definition: intersectionwriter.hh:142
typename GridView::Intersection Cell
Definition: intersectionwriter.hh:124
CellIterator beginCells() const
Definition: intersectionwriter.hh:136
ConnectivityWriter makeConnectivity() const
Definition: intersectionwriter.hh:154
CellIterator endCells() const
Definition: intersectionwriter.hh:139
Dune::VTK::CornerIterator< CellIterator > CornerIterator
Definition: intersectionwriter.hh:127
const Communication & comm() const
Definition: intersectionwriter.hh:157
Skeleton function for intersection writer.
Definition: intersectionwriter.hh:170
Dune::VTK::SkeletonFunctionTraits< GridView, typename GridView::ctype > Traits
Definition: intersectionwriter.hh:173
SkeletonFunction(const GridView &gv, const Mapper &mapper, const F &field)
Definition: intersectionwriter.hh:175
unsigned dimRange() const
return number of components
Definition: intersectionwriter.hh:205
void evaluate(const typename Traits::Cell &intersection, const typename Traits::Domain &xl, typename Traits::Range &result) const
Definition: intersectionwriter.hh:207
Distance implementation details.
Definition: cvfelocalresidual.hh:25
double accessEntry(const Field &f, int mycomp, int i)
Definition: function.hh:34
Definition: adapt.hh:17
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.