GCC Code Coverage Report

Directory: gridformat/
File: gridformat/traits/dune.hpp
Date: 2024-11-10 16:24:00
Exec Total Coverage
Lines: 355 379 93.7%
Functions: 397 398 99.7%
Branches: 440 798 55.1%
Line Branch Exec Source
1 // SPDX-FileCopyrightText: 2022-2023 Dennis Gläser <dennis.glaeser@iws.uni-stuttgart.de>
2 // SPDX-License-Identifier: MIT
3 /*!
4 * \file
5 * \ingroup PredefinedTraits
6 * \brief Traits specializations for <a href="https://gitlab.dune-project.org/core/dune-grid">dune grid views</a>
7 */
8 #ifndef GRIDFORMAT_TRAITS_DUNE_HPP_
9 #define GRIDFORMAT_TRAITS_DUNE_HPP_
10
11 #include <ranges>
12 #include <cassert>
13 #include <algorithm>
14 #include <tuple>
15 // dune seems to not explicitly include this but uses std::int64_t.
16 // With gcc-13 this leads to an error, maybe before gcc-13 the header
17 // was included by some other standard library header...
18 #include <cstdint>
19
20 #ifdef GRIDFORMAT_IGNORE_DUNE_WARNINGS
21 #pragma GCC diagnostic push
22 #pragma GCC diagnostic ignored "-Wold-style-cast"
23 #pragma GCC diagnostic ignored "-Wuseless-cast"
24 #pragma GCC diagnostic ignored "-Warray-bounds"
25 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
26 #pragma GCC diagnostic ignored "-Wnull-dereference"
27 #pragma GCC diagnostic ignored "-Wzero-as-null-pointer-constant"
28 #endif // GRIDFORMAT_IGNORE_DUNE_WARNINGS
29 #include <dune/common/fmatrix.hh>
30 #include <dune/common/fvector.hh>
31 #include <dune/geometry/type.hh>
32 #include <dune/grid/common/gridview.hh>
33 #include <dune/grid/common/gridenums.hh>
34 #include <dune/grid/common/gridfactory.hh>
35 #include <dune/grid/yaspgrid.hh>
36 #ifdef GRIDFORMAT_IGNORE_DUNE_WARNINGS
37 #pragma GCC diagnostic pop
38 #endif // GRIDFORMAT_IGNORE_DUNE_WARNINGS
39
40 #include <gridformat/common/ranges.hpp>
41 #include <gridformat/common/exceptions.hpp>
42 #include <gridformat/common/type_traits.hpp>
43
44 #include <gridformat/grid/cell_type.hpp>
45 #include <gridformat/grid/traits.hpp>
46
47 namespace GridFormat::Traits {
48
49 template<typename T, int R, int C>
50 struct StaticSize<Dune::FieldMatrix<T, R, C>> {
51 static constexpr std::size_t value = R;
52 };
53
54 #ifndef DOXYGEN
55 namespace DuneDetail {
56
57 1144 inline int map_corner_index(const Dune::GeometryType& gt, int i) {
58
2/2
✓ Branch 1 taken 568 times.
✓ Branch 2 taken 576 times.
 1144 if (gt.isQuadrilateral()) {
59
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 568 times.
 568 assert(i < 4);
60 static constexpr int map[4] = {0, 1, 3, 2};
61 568 return map[i];
62 }
63
1/2
✓ Branch 1 taken 576 times.
✗ Branch 2 not taken.
 576 if (gt.isHexahedron()) {
64
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 576 times.
 576 assert(i < 8);
65 static constexpr int map[8] = {0, 1, 3, 2, 4, 5, 7, 6};
66 576 return map[i];
67 }
68 ✗ return i;
69 }
70
71 257 inline constexpr GridFormat::CellType cell_type(const Dune::GeometryType& gt) {
72
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 257 times.
 257 if (gt.isVertex()) return GridFormat::CellType::vertex;
73
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 257 times.
 257 if (gt.isLine()) return GridFormat::CellType::segment;
74
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 257 times.
 257 if (gt.isTriangle()) return GridFormat::CellType::triangle;
75
2/2
✓ Branch 1 taken 233 times.
✓ Branch 2 taken 24 times.
 257 if (gt.isQuadrilateral()) return GridFormat::CellType::quadrilateral;
76
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
 24 if (gt.isTetrahedron()) return GridFormat::CellType::tetrahedron;
77
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 if (gt.isHexahedron()) return GridFormat::CellType::hexahedron;
78 ✗ throw NotImplemented("Unknown Dune::GeometryType");
79 }
80
81 template<typename GridView>
82 using Element = typename GridView::template Codim<0>::Entity;
83
84 template<typename GridView>
85 using Vertex = typename GridView::template Codim<GridView::dimension>::Entity;
86
87 } // namespace DuneDetail
88 #endif // DOXYGEN
89
90
91 template<typename Traits>
92 struct Points<Dune::GridView<Traits>> {
93 652 static decltype(auto) get(const Dune::GridView<Traits>& grid_view) {
94 using GV = Dune::GridView<Traits>;
95 static constexpr int vertex_codim = GV::dimension;
96 return std::ranges::subrange(
97
1/2
✓ Branch 1 taken 428 times.
✗ Branch 2 not taken.
 654 grid_view.template begin<vertex_codim, Dune::InteriorBorder_Partition>(),
98
1/2
✓ Branch 1 taken 428 times.
✗ Branch 2 not taken.
 1304 grid_view.template end<vertex_codim, Dune::InteriorBorder_Partition>()
99
1/2
✓ Branch 1 taken 428 times.
✗ Branch 2 not taken.
 1304 );
100 }
101 };
102
103 template<typename Traits>
104 struct Cells<Dune::GridView<Traits>> {
105 2352 static decltype(auto) get(const Dune::GridView<Traits>& grid_view) {
106 return std::ranges::subrange(
107
1/2
✓ Branch 1 taken 1713 times.
✗ Branch 2 not taken.
 2534 grid_view.template begin<0, Dune::Interior_Partition>(),
108
1/2
✓ Branch 1 taken 1713 times.
✗ Branch 2 not taken.
 4704 grid_view.template end<0, Dune::Interior_Partition>()
109
1/2
✓ Branch 1 taken 1713 times.
✗ Branch 2 not taken.
 4704 );
110 }
111 };
112
113 template<typename Traits>
114 struct NumberOfPoints<Dune::GridView<Traits>> {
115 786 static auto get(const Dune::GridView<Traits>& grid_view) {
116 static constexpr int point_codim = Dune::GridView<Traits>::dimension;
117
2/2
✓ Branch 2 taken 111 times.
✓ Branch 3 taken 372 times.
 786 if (grid_view.comm().size() == 1)
118 222 return static_cast<std::size_t>(grid_view.size(point_codim));
119 return static_cast<std::size_t>(
120
2/4
✓ Branch 1 taken 372 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 372 times.
✗ Branch 5 not taken.
 564 Ranges::size(Points<Dune::GridView<Traits>>::get(grid_view))
121 564 );
122 }
123 };
124
125 template<typename Traits>
126 struct NumberOfCells<Dune::GridView<Traits>> {
127 2992 static auto get(const Dune::GridView<Traits>& grid_view) {
128
2/2
✓ Branch 2 taken 881 times.
✓ Branch 3 taken 1038 times.
 2992 if (grid_view.comm().size() == 1)
129 1762 return static_cast<std::size_t>(grid_view.size(0));
130 return static_cast<std::size_t>(
131
2/4
✓ Branch 1 taken 1038 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1038 times.
✗ Branch 5 not taken.
 1230 Ranges::size(Cells<Dune::GridView<Traits>>::get(grid_view))
132 1230 );
133 }
134 };
135
136 template<typename Traits>
137 struct NumberOfCellPoints<Dune::GridView<Traits>, DuneDetail::Element<Dune::GridView<Traits>>> {
138 396 static auto get(const Dune::GridView<Traits>&,
139 const DuneDetail::Element<Dune::GridView<Traits>>& cell) {
140 396 return cell.subEntities(Dune::GridView<Traits>::dimension);
141 }
142 };
143
144 template<typename Traits>
145 struct CellPoints<Dune::GridView<Traits>, DuneDetail::Element<Dune::GridView<Traits>>> {
146 1436 static decltype(auto) get(const Dune::GridView<Traits>&,
147 const DuneDetail::Element<Dune::GridView<Traits>>& element) {
148 static constexpr int dim = Dune::GridView<Traits>::dimension;
149
4/8
✓ Branch 1 taken 818 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 818 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 818 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 818 times.
✗ Branch 11 not taken.
 3616 return std::views::iota(unsigned{0}, element.subEntities(dim)) | std::views::transform([&] (int i) {
150
4/8
✓ Branch 1 taken 976 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 976 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 168 times.
✗ Branch 9 not taken.
✓ Branch 12 taken 168 times.
✗ Branch 13 not taken.
 2288 return element.template subEntity<dim>(DuneDetail::map_corner_index(element.type(), i));
151 2872 });
152 }
153 };
154
155 template<typename Traits>
156 struct CellType<Dune::GridView<Traits>, DuneDetail::Element<Dune::GridView<Traits>>> {
157 414 static decltype(auto) get(const Dune::GridView<Traits>&,
158 const DuneDetail::Element<Dune::GridView<Traits>>& element) {
159
2/4
✓ Branch 1 taken 257 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 257 times.
✗ Branch 5 not taken.
 414 return DuneDetail::cell_type(element.type());
160 }
161 };
162
163 template<typename Traits>
164 struct PointCoordinates<Dune::GridView<Traits>, DuneDetail::Vertex<Dune::GridView<Traits>>> {
165 504 static decltype(auto) get(const Dune::GridView<Traits>&,
166 const DuneDetail::Vertex<Dune::GridView<Traits>>& vertex) {
167
2/4
✓ Branch 1 taken 324 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 324 times.
✗ Branch 5 not taken.
 504 return vertex.geometry().center();
168 }
169 };
170
171 template<typename Traits>
172 struct PointId<Dune::GridView<Traits>, DuneDetail::Vertex<Dune::GridView<Traits>>> {
173 1216 static decltype(auto) get(const Dune::GridView<Traits>& grid_view,
174 const DuneDetail::Vertex<Dune::GridView<Traits>>& vertex) {
175 1216 return grid_view.indexSet().index(vertex);
176 }
177 };
178
179
180 #ifndef DOXYGEN
181 namespace DuneDetail {
182
183 template<typename T>
184 struct IsDuneYaspGrid : public std::false_type {};
185
186 template<int dim, typename Coords>
187 struct IsDuneYaspGrid<Dune::YaspGrid<dim, Coords>> : public std::true_type {
188 static constexpr bool uses_tp_coords = std::same_as<
189 Coords,
190 Dune::TensorProductCoordinates<typename Coords::ctype, dim>
191 >;
192 };
193
194 template<typename GridView>
195 inline constexpr bool is_yasp_grid_view = IsDuneYaspGrid<typename GridView::Grid>::value;
196
197 template<typename GridView> requires(is_yasp_grid_view<GridView>)
198 inline constexpr bool uses_tensor_product_coords = IsDuneYaspGrid<typename GridView::Grid>::uses_tp_coords;
199
200 template<typename ctype, int dim>
201 190 auto spacing_in(int direction, const Dune::EquidistantCoordinates<ctype, dim>& coords) {
202 190 return coords.meshsize(direction, 0);
203 }
204 template<typename ctype, int dim>
205 10 auto spacing_in(int direction, const Dune::EquidistantOffsetCoordinates<ctype, dim>& coords) {
206 10 return coords.meshsize(direction, 0);
207 }
208
209 } // namespace DuneDetail
210 #endif // DOXYGEN
211
212
213 // Register YaspGrid as structured grid
214 template<typename Traits> requires(DuneDetail::is_yasp_grid_view<Dune::GridView<Traits>>)
215 struct Extents<Dune::GridView<Traits>> {
216 208 static auto get(const Dune::GridView<Traits>& grid_view) {
217
4/8
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 120 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 120 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 120 times.
✗ Branch 11 not taken.
 208 const auto level = std::ranges::begin(Cells<Dune::GridView<Traits>>::get(grid_view))->level();
218
2/4
✓ Branch 1 taken 76 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 76 times.
✗ Branch 5 not taken.
 208 const auto& grid_level = *grid_view.grid().begin(level);
219 208 const auto& interior_grid = grid_level.interior[0];
220 208 const auto& gc = *interior_grid.dataBegin();
221
222 static constexpr int dim = Traits::Grid::dimension;
223 std::array<std::size_t, dim> result;
224
2/2
✓ Branch 0 taken 284 times.
✓ Branch 1 taken 120 times.
 712 for (int i = 0; i < Traits::Grid::dimension; ++i)
225
1/2
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
 504 result[i] = gc.max(i) - gc.min(i) + 1;
226 208 return result;
227 }
228 };
229
230 template<typename Traits, typename Entity> requires(DuneDetail::is_yasp_grid_view<Dune::GridView<Traits>>)
231 struct Location<Dune::GridView<Traits>, Entity> {
232 12384 static auto get(const Dune::GridView<Traits>& grid_view, const Entity& entity) {
233
3/6
✓ Branch 1 taken 6192 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6192 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6192 times.
✗ Branch 8 not taken.
 12384 auto const& grid_level = *grid_view.grid().begin(entity.level());
234 12384 auto const& interior = grid_level.interior[0];
235 12384 auto const& extent_bounds = *interior.dataBegin();
236
237
1/2
✓ Branch 3 taken 6192 times.
✗ Branch 4 not taken.
 12384 auto result = entity.impl().transformingsubiterator().coord();
238
2/2
✓ Branch 0 taken 17659 times.
✓ Branch 1 taken 6192 times.
 47702 for (int i = 0; i < Dune::GridView<Traits>::dimension; ++i)
239 35318 result[i] -= extent_bounds.min(i);
240 12384 return result;
241 }
242 };
243
244 template<typename Traits> requires(DuneDetail::is_yasp_grid_view<Dune::GridView<Traits>>)
245 struct Origin<Dune::GridView<Traits>> {
246 40 static auto get(const Dune::GridView<Traits>& grid_view) {
247
4/8
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 20 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 20 times.
✗ Branch 11 not taken.
 40 const auto level = std::ranges::begin(Cells<Dune::GridView<Traits>>::get(grid_view))->level();
248
2/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
 40 auto const& grid_level = *grid_view.grid().begin(level);
249 40 auto const& interior = grid_level.interior[0];
250 40 auto const& extent_bounds = *interior.dataBegin();
251
252 std::array<typename Traits::Grid::ctype, Traits::Grid::dimension> result;
253
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 20 times.
 140 for (int i = 0; i < Traits::Grid::dimension; ++i)
254
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 100 result[i] = grid_level.coords.coordinate(i, extent_bounds.min(i));
255 40 return result;
256 }
257 };
258
259 template<typename Traits>
260 requires(DuneDetail::is_yasp_grid_view<Dune::GridView<Traits>> and
261 !DuneDetail::uses_tensor_product_coords<Dune::GridView<Traits>>) // spacing not uniquely defined for tpcoords
262 struct Spacing<Dune::GridView<Traits>> {
263 44 static auto get(const Dune::GridView<Traits>& grid_view) {
264
4/8
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
 44 const auto level = std::ranges::begin(Cells<Dune::GridView<Traits>>::get(grid_view))->level();
265
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 44 auto const& grid_level = *grid_view.grid().begin(level);
266
267 std::array<typename Traits::Grid::ctype, Traits::Grid::dimension> result;
268
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 4 times.
 154 for (int i = 0; i < Traits::Grid::dimension; ++i)
269
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 110 result[i] = DuneDetail::spacing_in(i, grid_level.coords);
270 44 return result;
271 }
272 };
273
274 template<typename Traits> requires(DuneDetail::is_yasp_grid_view<Dune::GridView<Traits>>)
275 struct Ordinates<Dune::GridView<Traits>> {
276 37 static auto get(const Dune::GridView<Traits>& grid_view, int direction) {
277
4/8
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 32 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 32 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 32 times.
✗ Branch 11 not taken.
 37 const auto level = std::ranges::begin(Cells<Dune::GridView<Traits>>::get(grid_view))->level();
278 37 auto const& grid_level = *grid_view.grid().begin(level);
279 37 auto const& interior = grid_level.interior[0];
280 37 auto const& extent_bounds = *interior.dataBegin();
281
282 37 const auto num_point_ordinates = extent_bounds.max(direction) - extent_bounds.min(direction) + 2;
283
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 37 std::vector<typename Traits::Grid::ctype> ordinates(num_point_ordinates);
284
2/2
✓ Branch 0 taken 144 times.
✓ Branch 1 taken 32 times.
 196 for (int i = 0; i < num_point_ordinates; ++i)
285 159 ordinates[i] = grid_level.coords.coordinate(direction, extent_bounds.min(direction) + i);
286 37 return ordinates;
287 }
288 };
289
290 } // namespace GridFormat::Traits
291
292
293 // Higher-order output is only available with dune-localfunctions
294 #if GRIDFORMAT_HAVE_DUNE_LOCALFUNCTIONS
295
296 #include <algorithm>
297 #include <iterator>
298 #include <utility>
299 #include <ranges>
300 #include <cassert>
301 #include <type_traits>
302 #include <unordered_map>
303 #include <map>
304
305 #ifdef GRIDFORMAT_IGNORE_DUNE_WARNINGS
306 #pragma GCC diagnostic push
307 #pragma GCC diagnostic ignored "-Wold-style-cast"
308 #pragma GCC diagnostic ignored "-Wtype-limits"
309 #endif // GRIDFORMAT_IGNORE_DUNE_WARNINGS
310 #include <dune/geometry/referenceelements.hh>
311 #include <dune/grid/common/mcmgmapper.hh>
312 #include <dune/localfunctions/lagrange/equidistantpoints.hh>
313 #ifdef GRIDFORMAT_IGNORE_DUNE_WARNINGS
314 #pragma GCC diagnostic pop
315 #endif // GRIDFORMAT_IGNORE_DUNE_WARNINGS
316
317 #include <gridformat/common/reserved_vector.hpp>
318 #include <gridformat/common/type_traits.hpp>
319 #include <gridformat/common/precision.hpp>
320 #include <gridformat/common/concepts.hpp>
321 #include <gridformat/common/field.hpp>
322 #include <gridformat/grid.hpp>
323
324
325 namespace GridFormat {
326
327 // create this alias s.t. we can explicitly
328 // address the original Dune namespace
329 namespace DUNE = Dune;
330
331 namespace Dune {
332
333 #ifndef DOXYGEN
334 namespace LagrangeDetail {
335
336 2508 inline int dune_to_gfmt_sub_entity(const DUNE::GeometryType& gt, int i, int codim) {
337
2/2
✓ Branch 1 taken 92 times.
✓ Branch 2 taken 2416 times.
 2508 if (gt.isTriangle()) {
338
2/2
✓ Branch 0 taken 44 times.
✓ Branch 1 taken 48 times.
 92 if (codim == 1) {
339
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 44 times.
 44 assert(i < 3);
340 static constexpr int map[3] = {0, 2, 1};
341 44 return map[i];
342 }
343 }
344
2/2
✓ Branch 1 taken 588 times.
✓ Branch 2 taken 1876 times.
 2464 if (gt.isQuadrilateral()) {
345
2/2
✓ Branch 0 taken 252 times.
✓ Branch 1 taken 336 times.
 588 if (codim == 2) {
346
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 252 times.
 252 assert(i < 4);
347 static constexpr int map[4] = {0, 1, 3, 2};
348 252 return map[i];
349 }
350
2/2
✓ Branch 0 taken 300 times.
✓ Branch 1 taken 36 times.
 336 if (codim == 1) {
351
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 300 times.
 300 assert(i < 4);
352 static constexpr int map[4] = {3, 1, 0, 2};
353 300 return map[i];
354 }
355 }
356
2/2
✓ Branch 1 taken 252 times.
✓ Branch 2 taken 1660 times.
 1912 if (gt.isTetrahedron()) {
357
2/2
✓ Branch 0 taken 160 times.
✓ Branch 1 taken 92 times.
 252 if (codim == 2) {
358
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 160 times.
 160 assert(i < 6);
359 static constexpr int map[6] = {0, 2, 1, 3, 4, 5};
360 160 return map[i];
361 }
362
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 68 times.
 92 if (codim == 1) {
363
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
 24 assert(i < 4);
364 static constexpr int map[4] = {3, 0, 2, 1};
365 24 return map[i];
366 }
367 }
368
2/2
✓ Branch 1 taken 1576 times.
✓ Branch 2 taken 152 times.
 1728 if (gt.isHexahedron()) {
369
2/2
✓ Branch 0 taken 236 times.
✓ Branch 1 taken 1340 times.
 1576 if (codim == 3) {
370
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 236 times.
 236 assert(i < 8);
371 static constexpr int map[8] = {0, 1, 3, 2, 4, 5, 7, 6};
372 236 return map[i];
373 }
374
2/2
✓ Branch 0 taken 808 times.
✓ Branch 1 taken 532 times.
 1340 if (codim == 2) {
375
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 808 times.
 808 assert(i < 12);
376 static constexpr int map[12] = {8, 9, 11, 10, 3, 1, 0, 2, 7, 5, 4, 6};
377 808 return map[i];
378 }
379 }
380 684 return i;
381 }
382
383 1818 inline constexpr GridFormat::CellType cell_type(const DUNE::GeometryType& gt) {
384
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1818 times.
 1818 if (gt.isLine()) return GridFormat::CellType::lagrange_segment;
385
2/2
✓ Branch 1 taken 438 times.
✓ Branch 2 taken 1380 times.
 1818 if (gt.isTriangle()) return GridFormat::CellType::lagrange_triangle;
386
2/2
✓ Branch 1 taken 318 times.
✓ Branch 2 taken 1062 times.
 1380 if (gt.isQuadrilateral()) return GridFormat::CellType::lagrange_quadrilateral;
387
2/2
✓ Branch 1 taken 1026 times.
✓ Branch 2 taken 36 times.
 1062 if (gt.isTetrahedron()) return GridFormat::CellType::lagrange_tetrahedron;
388
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 if (gt.isHexahedron()) return GridFormat::CellType::lagrange_hexahedron;
389 ✗ throw NotImplemented("Unsupported Dune::GeometryType");
390 }
391
392 // Exposes the lagrange points of a geometry in gridformat ordering
393 template<typename GridView>
394 class LocalPoints {
395 // reserve space for third-order hexahedra
396 static constexpr std::size_t reserved_size = 64;
397
398 public:
399 using PointSet = DUNE::EquidistantPointSet<typename GridView::ctype, GridView::dimension>;
400 using Point = typename PointSet::LagrangePoint;
401
402 60 explicit LocalPoints(unsigned int order)
403 60 : _points(order)
404 60 {}
405
406 60 void build(const DUNE::GeometryType& geo_type) {
407
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 36 times.
 60 if (geo_type.dim() != GridView::dimension)
408 ✗ throw ValueError("Dimension of given geometry does not match the grid");
409 60 _points.build(geo_type);
410 60 _setup_sorted_indices(geo_type);
411 60 }
412
413 6672 std::size_t size() const {
414 6672 return _points.size();
415 }
416
417 const Point& at(std::size_t i) const {
418 return _points[_sorted_indices.at(i)];
419 }
420
421 6372 std::ranges::range auto get() const {
422
2/4
✓ Branch 1 taken 3336 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3336 times.
✗ Branch 5 not taken.
 44272 return _sorted_indices | std::views::transform([&] (unsigned int i) {
423 34428 return _points[i];
424 12744 });
425 }
426
427 private:
428 60 void _setup_sorted_indices(const DUNE::GeometryType& geo_type) {
429
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 60 std::ranges::copy(
430
1/2
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
 60 std::views::iota(unsigned{0}, static_cast<unsigned int>(_points.size())),
431 60 std::back_inserter(_sorted_indices)
432 );
433 1574 std::ranges::sort(_sorted_indices, [&] (unsigned int i1, unsigned int i2) {
434 1810 const auto& key1 = _points[i1].localKey();
435 1810 const auto& key2 = _points[i2].localKey();
436 using LagrangeDetail::dune_to_gfmt_sub_entity;
437
8/8
✓ Branch 2 taken 322 times.
✓ Branch 3 taken 170 times.
✓ Branch 6 taken 46 times.
✓ Branch 7 taken 26 times.
✓ Branch 10 taken 788 times.
✓ Branch 11 taken 310 times.
✓ Branch 14 taken 98 times.
✓ Branch 15 taken 50 times.
 1810 if (key1.codim() == key2.codim()) {
438 1254 const auto mapped1 = dune_to_gfmt_sub_entity(geo_type, key1.subEntity(), key1.codim());
439 1254 const auto mapped2 = dune_to_gfmt_sub_entity(geo_type, key2.subEntity(), key2.codim());
440
8/8
✓ Branch 0 taken 42 times.
✓ Branch 1 taken 280 times.
✓ Branch 4 taken 6 times.
✓ Branch 5 taken 40 times.
✓ Branch 8 taken 96 times.
✓ Branch 9 taken 692 times.
✓ Branch 12 taken 14 times.
✓ Branch 13 taken 84 times.
 1254 return mapped1 == mapped2 ? key1.index() < key2.index() : mapped1 < mapped2;
441 }
442 556 return _points[i1].localKey().codim() > _points[i2].localKey().codim();
443 });
444 60 }
445
446 PointSet _points;
447 ReservedVector<unsigned int, reserved_size> _sorted_indices;
448 };
449
450 class PointMapper {
451 public:
452 template<typename GridView>
453 60 explicit PointMapper(const GridView& grid_view) {
454
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 60 _codim_to_global_indices.resize(GridView::dimension + 1);
455
2/2
✓ Branch 0 taken 120 times.
✓ Branch 1 taken 36 times.
 264 for (int codim = 0; codim < GridView::dimension + 1; ++codim)
456
2/4
✓ Branch 2 taken 120 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 120 times.
✗ Branch 6 not taken.
 204 _codim_to_global_indices[codim].resize(grid_view.size(codim));
457 60 }
458
459 struct Key {
460 unsigned int codim;
461 std::size_t global_index;
462 unsigned int sub_index;
463 };
464
465 34428 bool contains(const Key& key) const {
466 34428 const auto& entity_dofs = _codim_to_global_indices[key.codim][key.global_index];
467
2/2
✓ Branch 1 taken 11562 times.
✓ Branch 2 taken 22866 times.
 34428 if (entity_dofs.size() <= key.sub_index)
468 11562 return false;
469 22866 return entity_dofs[key.sub_index] != -1;
470 }
471
472 11562 void insert(const Key& key, std::size_t index) {
473 11562 auto& entity_dofs = _codim_to_global_indices[key.codim][key.global_index];
474
1/2
✓ Branch 1 taken 11562 times.
✗ Branch 2 not taken.
 11562 if (entity_dofs.size() < key.sub_index + 1)
475
1/2
✓ Branch 1 taken 11562 times.
✗ Branch 2 not taken.
 11562 entity_dofs.resize(key.sub_index + 1, -1);
476 11562 entity_dofs[key.sub_index] = static_cast<std::int64_t>(index);
477 11562 }
478
479 22866 std::size_t get(const Key& key) const {
480
1/2
✗ Branch 3 not taken.
✓ Branch 4 taken 22866 times.
 22866 assert(_codim_to_global_indices[key.codim][key.global_index].size() > key.sub_index);
481 22866 auto index = _codim_to_global_indices[key.codim][key.global_index][key.sub_index];
482
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 22866 times.
 22866 assert(index != -1);
483 22866 return static_cast<std::size_t>(index);
484 }
485
486 void clear() {
487 _codim_to_global_indices.clear();
488 }
489
490 private:
491 std::vector< // codim
492 std::vector< // entity index
493 ReservedVector<std::int64_t, 20> // sub-entity index to global indices
494 >
495 > _codim_to_global_indices;
496 };
497
498 } // namespace LagrangeDetail
499 #endif // DOXYGEN
500
501
502 /*!
503 * \ingroup PredefinedTraits
504 * \brief Exposes a `Dune::GridView` as a grid composed of lagrange cells with the given order.
505 * Can be used to conveniently write `Dune::Functions` into grid files.
506 * \note This is only available if dune-localfunctions is found on the system.
507 */
508 template<typename GV>
509 class LagrangePolynomialGrid {
510 using Element = typename GV::Grid::template Codim<0>::Entity;
511 using ElementGeometry = typename Element::Geometry;
512 using GlobalCoordinate = typename ElementGeometry::GlobalCoordinate;
513 using LocalCoordinate = typename ElementGeometry::LocalCoordinate;
514
515 using LocalPoints = LagrangeDetail::LocalPoints<GV>;
516 using Mapper = DUNE::MultipleCodimMultipleGeomTypeMapper<GV>;
517 static constexpr int dim = GV::dimension;
518
519 template<typename Coord>
520 struct P {
521 std::size_t index;
522 Coord coordinates;
523 };
524
525 using _GlobalPoint = P<GlobalCoordinate>;
526 using _LocalPoint = P<LocalCoordinate>;
527
528 public:
529 using GridView = GV;
530 using Position = GlobalCoordinate;
531 using Cell = Element;
532 using Point = _GlobalPoint;
533 using LocalPoint = _LocalPoint;
534
535 36 explicit LagrangePolynomialGrid(const GridView& grid_view, unsigned int order = 1)
536 36 : _grid_view{grid_view}
537 36 , _order{order} {
538
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
 36 if (order < 1)
539 ✗ throw InvalidState("Order must be >= 1");
540
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 36 update(grid_view);
541 36 }
542
543 60 void update(const GridView& grid_view) {
544 60 clear();
545 60 _grid_view = grid_view;
546 60 _make_codim_mappers();
547 60 _update_local_points();
548 60 _update_mesh();
549 60 }
550
551 84 void clear() {
552 84 _codim_to_mapper.clear();
553 84 _local_points.clear();
554 84 _points.clear();
555 84 _cells.clear();
556 84 }
557
558 2238 std::size_t number_of_cells() const {
559
2/2
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 1446 times.
 2238 return _cells.empty() ? 0 : Traits::NumberOfCells<GridView>::get(_grid_view);
560 }
561
562 2238 std::size_t number_of_points() const {
563 2238 return _points.size();
564 }
565
566 6672 std::size_t number_of_points(const Element& element) const {
567
2/4
✓ Branch 1 taken 3636 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3636 times.
✗ Branch 5 not taken.
 6672 return _local_points.at(element.type()).size();
568 }
569
570 432 std::ranges::range auto cells() const {
571 432 return Traits::Cells<GridView>::get(_grid_view);
572 }
573
574 108 std::ranges::range auto points() const {
575
1/2
✓ Branch 2 taken 72 times.
✗ Branch 3 not taken.
 108 return std::views::iota(std::size_t{0}, _points.size())
576
2/4
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
 15189 | std::views::transform([&] (std::size_t idx) {
577 19821 return Point{idx, _points[idx]};
578 216 });
579 }
580
581 16380 std::ranges::range auto points(const Element& e) const {
582 16380 const auto& corners = _cells[_codim_to_mapper[0].index(e)];
583
1/2
✓ Branch 2 taken 8790 times.
✗ Branch 3 not taken.
 16380 return std::views::iota(std::size_t{0}, corners.size())
584
2/4
✓ Branch 1 taken 8790 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8790 times.
✗ Branch 5 not taken.
 127960 | std::views::transform([&] (std::size_t i) {
585 return Point{
586 90420 .index = corners[i],
587 90420 .coordinates = _points[corners[i]]
588 180840 };
589 32760 });
590 }
591
592 const GridView& grid_view() const {
593 return _grid_view;
594 }
595
596 private:
597 60 void _update_local_points() {
598
4/6
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
✓ Branch 11 taken 36 times.
✓ Branch 12 taken 36 times.
 120 for (const auto& gt : _grid_view.indexSet().types(0)) {
599
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 60 _local_points.emplace(gt, _order);
600
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
 60 _local_points.at(gt).build(gt);
601 }
602 60 }
603
604 60 void _make_codim_mappers() {
605 60 _codim_to_mapper.reserve(dim + 1);
606
2/4
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 36 times.
✗ Branch 6 not taken.
 60 _codim_to_mapper.emplace_back(Mapper{_grid_view, DUNE::mcmgLayout(DUNE::Codim<0>{})});
607 if constexpr (int(GridView::dimension) >= 1)
608
2/4
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 36 times.
✗ Branch 6 not taken.
 60 _codim_to_mapper.emplace_back(Mapper{_grid_view, DUNE::mcmgLayout(DUNE::Codim<1>{})});
609 if constexpr (int(GridView::dimension) >= 2)
610
2/4
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 36 times.
✗ Branch 6 not taken.
 60 _codim_to_mapper.emplace_back(Mapper{_grid_view, DUNE::mcmgLayout(DUNE::Codim<2>{})});
611 if constexpr (int(GridView::dimension) == 3)
612
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 24 _codim_to_mapper.emplace_back(Mapper{_grid_view, DUNE::mcmgLayout(DUNE::Codim<3>{})});
613 60 }
614
615 60 void _update_mesh() {
616
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 60 LagrangeDetail::PointMapper point_mapper{_grid_view};
617 60 std::size_t dof_index = 0;
618
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
 60 _cells.resize(_grid_view.size(0));
619
10/16
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 36 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2928 times.
✓ Branch 11 taken 408 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 3360 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 408 times.
✓ Branch 17 taken 2964 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 2928 times.
✓ Branch 20 taken 12 times.
 12864 for (const auto& element : Traits::Cells<GridView>::get(_grid_view)) {
620
3/6
✓ Branch 1 taken 3336 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3336 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3336 times.
✗ Branch 8 not taken.
 6372 const auto& local_points = _local_points.at(element.type()).get();
621
1/2
✓ Branch 1 taken 3336 times.
✗ Branch 2 not taken.
 6372 const auto& element_geometry = element.geometry();
622
1/2
✓ Branch 2 taken 3336 times.
✗ Branch 3 not taken.
 6372 const auto element_index = _codim_to_mapper[0].index(element);
623
2/4
✓ Branch 2 taken 3336 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3336 times.
✗ Branch 6 not taken.
 6372 _cells[element_index].reserve(local_points.size());
624
5/8
✓ Branch 1 taken 3336 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3336 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 34428 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 34428 times.
✓ Branch 12 taken 3336 times.
 78700 for (const auto& local_point : local_points) {
625 65956 const auto& localKey = local_point.localKey();
626 197868 const typename LagrangeDetail::PointMapper::Key key{
627 65956 .codim = localKey.codim(),
628
1/2
✓ Branch 4 taken 34428 times.
✗ Branch 5 not taken.
 131912 .global_index = _codim_to_mapper[localKey.codim()].subIndex(
629 65956 element, localKey.subEntity(), localKey.codim()
630 ),
631 65956 .sub_index = localKey.index()
632 };
633
3/4
✓ Branch 1 taken 34428 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 11562 times.
✓ Branch 4 taken 22866 times.
 65956 if (!point_mapper.contains(key)) {
634
1/2
✓ Branch 1 taken 11562 times.
✗ Branch 2 not taken.
 21472 point_mapper.insert(key, dof_index);
635
1/2
✓ Branch 2 taken 11562 times.
✗ Branch 3 not taken.
 21472 _cells[element_index].push_back(dof_index);
636
2/4
✓ Branch 2 taken 11562 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 11562 times.
✗ Branch 6 not taken.
 21472 _points.push_back(element_geometry.global(local_point.point()));
637 21472 dof_index++;
638 } else {
639
2/4
✓ Branch 2 taken 22866 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 22866 times.
✗ Branch 6 not taken.
 44484 _cells[element_index].push_back(point_mapper.get(key));
640 }
641 }
642 }
643 60 }
644
645 GridView _grid_view;
646 unsigned int _order;
647 std::vector<Mapper> _codim_to_mapper;
648 std::map<DUNE::GeometryType, LocalPoints> _local_points;
649 std::vector<GlobalCoordinate> _points;
650 std::vector<std::vector<std::size_t>> _cells;
651 };
652
653 namespace Traits {
654
655 template<typename GV>
656 struct GridView {
657 using type = GV;
658 static const auto& get(const GV& gv) {
659 return gv;
660 }
661 };
662
663 template<typename GV>
664 struct GridView<LagrangePolynomialGrid<GV>> {
665 using type = GV;
666 static const auto& get(const LagrangePolynomialGrid<GV>& grid) {
667 return grid.grid_view();
668 }
669 };
670
671 } // namespace Traits
672 } // namespace Dune
673
674
675 namespace Traits {
676
677 template<typename GridView>
678 struct Points<Dune::LagrangePolynomialGrid<GridView>> {
679 108 static auto get(const Dune::LagrangePolynomialGrid<GridView>& mesh) {
680 108 return mesh.points();
681 }
682 };
683
684 template<typename GridView>
685 struct Cells<Dune::LagrangePolynomialGrid<GridView>> {
686 432 static auto get(const Dune::LagrangePolynomialGrid<GridView>& mesh) {
687 432 return mesh.cells();
688 }
689 };
690
691 template<typename GridView>
692 struct NumberOfPoints<Dune::LagrangePolynomialGrid<GridView>> {
693 2106 static auto get(const Dune::LagrangePolynomialGrid<GridView>& mesh) {
694 2106 return mesh.number_of_points();
695 }
696 };
697
698 template<typename GridView>
699 struct NumberOfCells<Dune::LagrangePolynomialGrid<GridView>> {
700 2166 static auto get(const Dune::LagrangePolynomialGrid<GridView>& mesh) {
701 2166 return mesh.number_of_cells();
702 }
703 };
704
705 template<typename GridView>
706 struct NumberOfCellPoints<Dune::LagrangePolynomialGrid<GridView>,
707 typename Dune::LagrangePolynomialGrid<GridView>::Cell> {
708 6672 static auto get(const Dune::LagrangePolynomialGrid<GridView>& mesh,
709 const typename Dune::LagrangePolynomialGrid<GridView>::Cell& cell) {
710 6672 return mesh.number_of_points(cell);
711 }
712 };
713
714 template<typename GridView>
715 struct CellPoints<Dune::LagrangePolynomialGrid<GridView>,
716 typename Dune::LagrangePolynomialGrid<GridView>::Cell> {
717 3336 static auto get(const Dune::LagrangePolynomialGrid<GridView>& mesh,
718 const typename Dune::LagrangePolynomialGrid<GridView>::Cell& cell) {
719 3336 return mesh.points(cell);
720 }
721 };
722
723 template<typename GridView>
724 struct CellType<Dune::LagrangePolynomialGrid<GridView>,
725 typename Dune::LagrangePolynomialGrid<GridView>::Cell> {
726 3336 static auto get(const Dune::LagrangePolynomialGrid<GridView>&,
727 const typename Dune::LagrangePolynomialGrid<GridView>::Cell& cell) {
728
2/4
✓ Branch 1 taken 1818 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1818 times.
✗ Branch 5 not taken.
 3336 return Dune::LagrangeDetail::cell_type(cell.type());
729 }
730 };
731
732 template<typename GridView>
733 struct PointCoordinates<Dune::LagrangePolynomialGrid<GridView>,
734 typename Dune::LagrangePolynomialGrid<GridView>::Point> {
735 11562 static const auto& get(const Dune::LagrangePolynomialGrid<GridView>&,
736 const typename Dune::LagrangePolynomialGrid<GridView>::Point& point) {
737 11562 return point.coordinates;
738 }
739 };
740
741 template<typename GridView>
742 struct PointId<Dune::LagrangePolynomialGrid<GridView>,
743 typename Dune::LagrangePolynomialGrid<GridView>::Point> {
744 57552 static auto get(const Dune::LagrangePolynomialGrid<GridView>&,
745 const typename Dune::LagrangePolynomialGrid<GridView>::Point& point) {
746 57552 return point.index;
747 }
748 };
749
750 } // namespace Traits
751
752
753 namespace Dune {
754 namespace Concepts {
755
756 template<typename T, typename GridView>
757 concept Function = requires(const T& f, const typename GridView::template Codim<0>::Entity& element) {
758 { localFunction(f) };
759 { localFunction(f).bind(element) };
760 { localFunction(f)(element.geometry().center()) };
761 };
762
763 } // namespace Concepts
764
765 #ifndef DOXYGEN
766 namespace FunctionDetail {
767
768 template<typename Function, typename GridView>
769 using RangeType = std::invoke_result_t<
770 typename std::remove_cvref_t<Function>::LocalFunction,
771 typename GridView::template Codim<0>::Entity::Geometry::LocalCoordinate
772 >;
773
774 template<typename Function, typename GridView>
775 using RangeScalar = FieldScalar<RangeType<Function, GridView>>;
776
777 } // FunctionDetail
778 #endif // DOXYGEN
779
780 /*!
781 * \ingroup PredefinedTraits
782 * \brief Implements the field interface for a Dune::Function defined on a (wrapped) Dune::GridView.
783 * \note Takes ownership of the given function if constructed with an rvalue reference, otherwise it
784 * stores a reference.
785 */
786 template<typename _Function, typename Grid, GridFormat::Concepts::Scalar T>
787 requires(Concepts::Function<_Function, typename Dune::Traits::GridView<Grid>::type>)
788 class FunctionField : public GridFormat::Field {
789 using Function = std::remove_cvref_t<_Function>;
790 using GridView = typename Dune::Traits::GridView<Grid>::type;
791 using Element = typename GridView::template Codim<0>::Entity;
792 using ElementGeometry = typename Element::Geometry;
793
794 static constexpr bool is_higher_order = std::is_same_v<Grid, LagrangePolynomialGrid<GridView>>;
795
796 public:
797 template<typename F>
798 requires(std::is_same_v<std::remove_cvref_t<F>, Function>)
799 336 explicit FunctionField(F&& function,
800 const Grid& grid,
801 const Precision<T>& = {},
802 bool cellwise_constant = false)
803 672 : _function{std::forward<F>(function)}
804 336 , _grid{grid}
805 336 , _cellwise_constant{cellwise_constant} {
806 if constexpr (requires { function.basis(); }) {
807 static_assert(std::is_same_v<typename Function::Basis::GridView, GridView>);
808 if (&function.basis().gridView().grid() != &Dune::Traits::GridView<Grid>::get(_grid).grid())
809 throw ValueError("Function and mesh do not use the same underlying grid");
810 }
811 336 }
812
813 private:
814 3888 MDLayout _layout() const override {
815 return MDLayout{{
816
3/4
✓ Branch 0 taken 918 times.
✓ Branch 1 taken 1026 times.
✓ Branch 3 taken 918 times.
✗ Branch 4 not taken.
 3888 _cellwise_constant ? GridFormat::Traits::NumberOfCells<Grid>::get(_grid)
817
1/2
✓ Branch 1 taken 1026 times.
✗ Branch 2 not taken.
 2052 : GridFormat::Traits::NumberOfPoints<Grid>::get(_grid)
818
3/6
✓ Branch 1 taken 1944 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1944 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1944 times.
✗ Branch 8 not taken.
 3888 }}.template with_sub_layout_from<FunctionDetail::RangeType<Function, GridView>>();
819 }
820
821 1332 DynamicPrecision _precision() const override {
822 1332 return {Precision<T>{}};
823 }
824
825 336 Serialization _serialized() const override {
826
1/2
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
 336 const auto layout = _layout();
827
1/2
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
 336 const auto num_entries = layout.number_of_entries();
828
4/6
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 96 times.
✓ Branch 4 taken 72 times.
✓ Branch 6 taken 96 times.
✗ Branch 7 not taken.
 336 const auto num_entries_per_value = layout.dimension() == 1 ? 1 : layout.number_of_entries(1);
829
830
1/2
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
 336 Serialization result(num_entries*sizeof(T));
831
1/2
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
 336 auto out_data = result.template as_span_of<T>();
832
833 using GridFormat::Traits::Cells;
834
2/2
✓ Branch 0 taken 84 times.
✓ Branch 1 taken 84 times.
 336 if (_cellwise_constant) {
835 168 std::size_t count = 0;
836
1/2
✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
 168 auto local_function = localFunction(_function);
837
12/17
✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 84 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 84 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 5856 times.
✓ Branch 11 taken 1152 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 3636 times.
✓ Branch 14 taken 3396 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 3348 times.
✓ Branch 17 taken 2964 times.
✓ Branch 18 taken 2928 times.
✓ Branch 19 taken 2940 times.
✓ Branch 20 taken 12 times.
 26808 for (const auto& element : Cells<Grid>::get(_grid)) {
838
1/2
✓ Branch 1 taken 6972 times.
✗ Branch 2 not taken.
 13944 local_function.bind(element);
839
1/2
✓ Branch 1 taken 6972 times.
✗ Branch 2 not taken.
 13944 const auto& elem_geo = element.geometry();
840
2/4
✓ Branch 1 taken 6972 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6972 times.
✗ Branch 5 not taken.
 13944 const auto& local_pos = elem_geo.local(elem_geo.center());
841 13944 std::size_t offset = (count++)*num_entries_per_value;
842
3/6
✓ Branch 1 taken 6972 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3636 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 708 times.
✗ Branch 8 not taken.
 13944 _copy_values(local_function(local_pos), out_data, offset);
843 }
844 48 } else {
845
1/2
✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
 168 _fill_point_values(out_data, num_entries_per_value);
846 }
847
848 336 return result;
849 336 }
850
851 void _fill_point_values(std::span<T> out_data, std::size_t num_entries_per_value) const
852 requires(!is_higher_order) {
853 using GridFormat::Traits::Cells;
854 using GridFormat::Traits::CellPoints;
855 using GridFormat::Traits::PointCoordinates;
856 using GridFormat::Traits::PointId;
857 using GridFormat::Traits::NumberOfPoints;
858
859 auto local_function = localFunction(_function);
860 auto point_id_to_running_idx = make_point_id_map(_grid);
861 std::vector<bool> handled(NumberOfPoints<Grid>::get(_grid), false);
862
863 std::ranges::for_each(Cells<Grid>::get(_grid), [&] <typename C> (const C& element) {
864 const auto& element_geometry = element.geometry();
865 local_function.bind(element);
866 std::ranges::for_each(CellPoints<Grid, Element>::get(_grid, element), [&] <typename P> (const P& point) {
867 const auto point_id = PointId<Grid, P>::get(_grid, point);
868 const auto running_idx = point_id_to_running_idx.at(point_id);
869 if (!handled[running_idx]) {
870 const auto local_pos = element_geometry.local(PointCoordinates<Grid, P>::get(_grid, point));
871 std::size_t offset = running_idx*num_entries_per_value;
872 _copy_values(local_function(local_pos), out_data, offset);
873 }
874 handled[running_idx] = true;
875 });
876 });
877 }
878
879 168 void _fill_point_values(std::span<T> out_data, std::size_t num_entries_per_value) const
880 requires(is_higher_order) {
881 using GridFormat::Traits::Cells;
882 using GridFormat::Traits::CellPoints;
883
1/2
✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
 168 auto local_function = localFunction(_function);
884
1/2
✓ Branch 2 taken 84 times.
✗ Branch 3 not taken.
 168 std::vector<bool> handled(_grid.number_of_points(), false);
885
2/4
✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 84 times.
✗ Branch 5 not taken.
 8256 std::ranges::for_each(Cells<Grid>::get(_grid), [&] <typename C> (const C& element) {
886
16/32
✓ Branch 1 taken 318 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 336 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 738 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1026 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 18 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 18 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 18 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 36 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 36 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 36 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 438 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 438 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 438 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 1026 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 1026 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 1026 times.
✗ Branch 47 not taken.
 6972 const auto& element_geometry = element.geometry();
887
16/32
✓ Branch 1 taken 318 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 336 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 738 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1026 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 18 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 18 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 18 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 36 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 36 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 36 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 438 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 438 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 438 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 1026 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 1026 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 1026 times.
✗ Branch 47 not taken.
 6972 local_function.bind(element);
888
32/64
✓ Branch 1 taken 318 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 318 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 336 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 336 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 738 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 738 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1026 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1026 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 18 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 18 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 18 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 18 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 18 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 18 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 36 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 36 times.
✗ Branch 47 not taken.
✓ Branch 49 taken 36 times.
✗ Branch 50 not taken.
✓ Branch 52 taken 36 times.
✗ Branch 53 not taken.
✓ Branch 55 taken 36 times.
✗ Branch 56 not taken.
✓ Branch 58 taken 36 times.
✗ Branch 59 not taken.
✓ Branch 61 taken 438 times.
✗ Branch 62 not taken.
✓ Branch 64 taken 438 times.
✗ Branch 65 not taken.
✓ Branch 67 taken 438 times.
✗ Branch 68 not taken.
✓ Branch 70 taken 438 times.
✗ Branch 71 not taken.
✓ Branch 73 taken 438 times.
✗ Branch 74 not taken.
✓ Branch 76 taken 438 times.
✗ Branch 77 not taken.
✓ Branch 79 taken 1026 times.
✗ Branch 80 not taken.
✓ Branch 82 taken 1026 times.
✗ Branch 83 not taken.
✓ Branch 85 taken 1026 times.
✗ Branch 86 not taken.
✓ Branch 88 taken 1026 times.
✗ Branch 89 not taken.
✓ Branch 91 taken 1026 times.
✗ Branch 92 not taken.
✓ Branch 94 taken 1026 times.
✗ Branch 95 not taken.
 150484 std::ranges::for_each(_grid.points(element), [&] <typename P> (const P& point) {
889
48/64
✓ Branch 1 taken 3074 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1769 times.
✓ Branch 5 taken 1305 times.
✓ Branch 8 taken 4088 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 2353 times.
✓ Branch 12 taken 1735 times.
✓ Branch 15 taken 5674 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 2773 times.
✓ Branch 19 taken 2901 times.
✓ Branch 22 taken 11628 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 3016 times.
✓ Branch 26 taken 8612 times.
✓ Branch 29 taken 174 times.
✗ Branch 30 not taken.
✓ Branch 32 taken 117 times.
✓ Branch 33 taken 57 times.
✓ Branch 36 taken 174 times.
✗ Branch 37 not taken.
✓ Branch 39 taken 117 times.
✓ Branch 40 taken 57 times.
✓ Branch 43 taken 174 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 117 times.
✓ Branch 47 taken 57 times.
✓ Branch 50 taken 1188 times.
✗ Branch 51 not taken.
✓ Branch 53 taken 701 times.
✓ Branch 54 taken 487 times.
✓ Branch 57 taken 1188 times.
✗ Branch 58 not taken.
✓ Branch 60 taken 701 times.
✓ Branch 61 taken 487 times.
✓ Branch 64 taken 1188 times.
✗ Branch 65 not taken.
✓ Branch 67 taken 701 times.
✓ Branch 68 taken 487 times.
✓ Branch 71 taken 2774 times.
✗ Branch 72 not taken.
✓ Branch 74 taken 1121 times.
✓ Branch 75 taken 1653 times.
✓ Branch 78 taken 2774 times.
✗ Branch 79 not taken.
✓ Branch 81 taken 1121 times.
✓ Branch 82 taken 1653 times.
✓ Branch 85 taken 2774 times.
✗ Branch 86 not taken.
✓ Branch 88 taken 1121 times.
✓ Branch 89 taken 1653 times.
✓ Branch 92 taken 11628 times.
✗ Branch 93 not taken.
✓ Branch 95 taken 3016 times.
✓ Branch 96 taken 8612 times.
✓ Branch 99 taken 11628 times.
✗ Branch 100 not taken.
✓ Branch 102 taken 3016 times.
✓ Branch 103 taken 8612 times.
✓ Branch 106 taken 11628 times.
✗ Branch 107 not taken.
✓ Branch 109 taken 3016 times.
✓ Branch 110 taken 8612 times.
 71756 if (!handled[point.index]) {
890
16/32
✓ Branch 1 taken 1769 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2353 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2773 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3016 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 117 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 117 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 117 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 701 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 701 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 701 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 1121 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 1121 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 1121 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 3016 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 3016 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 3016 times.
✗ Branch 47 not taken.
 24776 const auto local_pos = element_geometry.local(point.coordinates);
891 24776 std::size_t offset = point.index*num_entries_per_value;
892
27/52
✓ Branch 1 taken 1769 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 2353 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1652 times.
✓ Branch 9 taken 1121 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 1652 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 3016 times.
✓ Branch 14 taken 1652 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 117 times.
✗ Branch 18 not taken.
✓ Branch 21 taken 117 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 117 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 117 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 117 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 701 times.
✗ Branch 34 not taken.
✓ Branch 37 taken 701 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 701 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 701 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 701 times.
✗ Branch 47 not taken.
✓ Branch 49 taken 1121 times.
✗ Branch 50 not taken.
✓ Branch 53 taken 1121 times.
✗ Branch 54 not taken.
✓ Branch 56 taken 1121 times.
✗ Branch 57 not taken.
✓ Branch 59 taken 1121 times.
✗ Branch 60 not taken.
✓ Branch 62 taken 1121 times.
✗ Branch 63 not taken.
✓ Branch 65 taken 3016 times.
✗ Branch 66 not taken.
✓ Branch 69 taken 3016 times.
✗ Branch 70 not taken.
✓ Branch 72 taken 3016 times.
✗ Branch 73 not taken.
✓ Branch 75 taken 3016 times.
✗ Branch 76 not taken.
✓ Branch 78 taken 3016 times.
✗ Branch 79 not taken.
 24776 _copy_values(local_function(local_pos), out_data, offset);
893 }
894
16/32
✓ Branch 1 taken 3074 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 4088 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 5674 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 11628 times.
✗ Branch 14 not taken.
✓ Branch 17 taken 174 times.
✗ Branch 18 not taken.
✓ Branch 21 taken 174 times.
✗ Branch 22 not taken.
✓ Branch 25 taken 174 times.
✗ Branch 26 not taken.
✓ Branch 29 taken 1188 times.
✗ Branch 30 not taken.
✓ Branch 33 taken 1188 times.
✗ Branch 34 not taken.
✓ Branch 37 taken 1188 times.
✗ Branch 38 not taken.
✓ Branch 41 taken 2774 times.
✗ Branch 42 not taken.
✓ Branch 45 taken 2774 times.
✗ Branch 46 not taken.
✓ Branch 49 taken 2774 times.
✗ Branch 50 not taken.
✓ Branch 53 taken 11628 times.
✗ Branch 54 not taken.
✓ Branch 57 taken 11628 times.
✗ Branch 58 not taken.
✓ Branch 61 taken 11628 times.
✗ Branch 62 not taken.
 71756 handled[point.index] = true;
895 });
896 5856 });
897 168 }
898
899 template<std::ranges::range R>
900 76958 void _copy_values(R&& range, std::span<T> out, std::size_t& offset) const {
901
1/2
✓ Branch 1 taken 38479 times.
✗ Branch 2 not taken.
 278664 std::ranges::for_each(range, [&] (const auto& entry) {
902 100853 _copy_values(entry, out, offset);
903 });
904 76958 }
905
906 template<GridFormat::Concepts::Scalar S>
907 188244 void _copy_values(const S value, std::span<T> out, std::size_t& offset) const {
908 188244 out[offset++] = static_cast<T>(value);
909 188244 }
910
911 _Function _function;
912 const Grid& _grid;
913 bool _cellwise_constant;
914 };
915
916
917 template<typename F, typename G, typename T = FunctionDetail::RangeScalar<F, typename Traits::GridView<G>::type>>
918 requires(std::is_lvalue_reference_v<F>)
919 FunctionField(F&&, const G&, const Precision<T>& = {}, bool = false) -> FunctionField<
920 std::add_lvalue_reference_t<std::add_const_t<std::remove_reference_t<F>>>, G, T
921 >;
922
923 template<typename F, typename G, typename T = FunctionDetail::RangeScalar<F, typename Traits::GridView<G>::type>>
924 requires(!std::is_lvalue_reference_v<F>)
925 FunctionField(F&&, const G&, const Precision<T>& = {}, bool = false) -> FunctionField<std::remove_cvref_t<F>, G, T>;
926
927
928 #ifndef DOXYGEN
929 template<typename T> struct IsLagrangeGrid : public std::false_type {};
930 template<typename GV> struct IsLagrangeGrid<LagrangePolynomialGrid<GV>> : public std::true_type {};
931
932 namespace FunctionDetail {
933 template<typename F, typename W, typename T>
934 336 void set_function(F&& f, W& w, const std::string& name, const Precision<T>& prec, bool is_cellwise) {
935
2/2
✓ Branch 0 taken 84 times.
✓ Branch 1 taken 84 times.
 336 if (is_cellwise)
936
1/2
✓ Branch 4 taken 84 times.
✗ Branch 5 not taken.
 168 w.set_cell_field(name, FunctionField{std::forward<F>(f), w.grid(), prec, true});
937 else
938
1/2
✓ Branch 4 taken 84 times.
✗ Branch 5 not taken.
 168 w.set_point_field(name, FunctionField{std::forward<F>(f), w.grid(), prec, false});
939 336 }
940
941 template<typename F, typename W>
942 288 void set_function(F&& f, W& w, const std::string& name, bool is_cellwise) {
943 using GV = Dune::Traits::GridView<typename W::Grid>::type;
944 using T = RangeScalar<F, GV>;
945
1/2
✓ Branch 2 taken 144 times.
✗ Branch 3 not taken.
 288 set_function(std::forward<F>(f), w, name, Precision<T>{}, is_cellwise);
946 288 }
947 } // namespace FunctionDetail
948 #endif //DOXYGEN
949
950 /*!
951 * \ingroup PredefinedTraits
952 * \brief Insert the given Dune::Function to the writer as point field.
953 * \note This requires the Writer to have been constructed with a LagrangePolynomialGrid.
954 */
955 template<typename Function, typename Writer>
956 144 void set_point_function(Function&& f, Writer& writer, const std::string& name) {
957 144 FunctionDetail::set_function(std::forward<Function>(f), writer, name, false);
958 144 }
959
960 /*!
961 * \ingroup PredefinedTraits
962 * \brief Insert the given Dune::Function to the writer as point field with the given precision.
963 * \note This requires the Writer to have been constructed with a LagrangePolynomialGrid.
964 */
965 template<typename Function, typename Writer, GridFormat::Concepts::Scalar T>
966 24 void set_point_function(Function&& f, Writer& writer, const std::string& name, const Precision<T>& prec) {
967 24 FunctionDetail::set_function(std::forward<Function>(f), writer, name, prec, false);
968 24 }
969
970 /*!
971 * \ingroup PredefinedTraits
972 * \brief Insert the given Dune::Function to the writer as cell field.
973 * \note This requires the Writer to have been constructed with a LagrangePolynomialGrid.
974 */
975 template<typename Writer, typename Function>
976 144 void set_cell_function(Function&& f, Writer& writer, const std::string& name) {
977 144 FunctionDetail::set_function(std::forward<Function>(f), writer, name, true);
978 144 }
979
980 /*!
981 * \ingroup PredefinedTraits
982 * \brief Insert the given Dune::Function to the writer as cell field with the given precision.
983 * \note This requires the Writer to have been constructed with a LagrangePolynomialGrid.
984 */
985 template<typename Writer, typename Function, GridFormat::Concepts::Scalar T>
986 24 void set_cell_function(Function&& f, Writer& writer, const std::string& name, const Precision<T>& prec) {
987 24 FunctionDetail::set_function(std::forward<Function>(f), writer, name, prec, true);
988 24 }
989
990 } // namespace Dune
991 } // namespace GridFormat
992
993 #else // GRIDFORMAT_HAVE_DUNE_LOCALFUNCTIONS
994
995 #ifndef DOXYGEN
996 namespace GridFormat::Dune {
997
998 template<typename... T>
999 class LagrangePolynomialGrid {
1000 template<typename... Ts> struct False { static constexpr bool value = false; };
1001 public:
1002 template<typename... Args>
1003 LagrangePolynomialGrid(Args&&... args) {
1004 static_assert(False<Args...>::value, "Dune-localfunctions required for higher-order output");
1005 }
1006 };
1007
1008 } // namespace GridFormat::Dune
1009 #endif // DOXYGEN
1010
1011 #endif // GRIDFORMAT_HAVE_DUNE_LOCALFUNCTIONS
1012
1013 namespace GridFormat::Dune {
1014
1015 12 inline constexpr ::Dune::GeometryType to_dune_geometry_type(const CellType& ct) {
1016 namespace DGT = ::Dune::GeometryTypes;
1017
1/15
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 12 times.
✗ 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 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 12 switch (ct) {
1018 ✗ case CellType::vertex: return DGT::vertex;
1019 ✗ case CellType::segment: return DGT::line;
1020 ✗ case CellType::triangle: return DGT::triangle;
1021 12 case CellType::pixel: return DGT::quadrilateral;
1022 ✗ case CellType::quadrilateral: return DGT::quadrilateral;
1023 ✗ case CellType::tetrahedron: return DGT::tetrahedron;
1024 ✗ case CellType::hexahedron: return DGT::hexahedron;
1025 ✗ case CellType::voxel: return DGT::hexahedron;
1026 ✗ case CellType::polygon: throw NotImplemented("No conversion from polygon to Dune::GeometryType");
1027 ✗ case CellType::lagrange_segment: throw NotImplemented("Cannot map higher-order cells to Dune::GeometryType");
1028 ✗ case CellType::lagrange_triangle: throw NotImplemented("Cannot map higher-order cells to Dune::GeometryType");
1029 ✗ case CellType::lagrange_quadrilateral: throw NotImplemented("Cannot map higher-order cells to Dune::GeometryType");
1030 ✗ case CellType::lagrange_tetrahedron: throw NotImplemented("Cannot map higher-order cells to Dune::GeometryType");
1031 ✗ case CellType::lagrange_hexahedron: throw NotImplemented("Cannot map higher-order cells to Dune::GeometryType");
1032 ✗ default: throw NotImplemented("Unknown cell type.");
1033 }
1034 }
1035
1036 template<std::integral TargetType, std::integral T>
1037 12 auto to_dune(const CellType& ct, const std::vector<T>& corners) {
1038
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 auto gt = to_dune_geometry_type(ct);
1039
1/2
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
 12 std::vector<TargetType> reordered(corners.size());
1040
1041 // voxels/pixels map to hexes/quads, but reordering has to be skipped
1042
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 12 if (ct != CellType::pixel && ct != CellType::voxel) {
1043 ✗ std::ranges::copy(
1044 ✗ std::views::iota(std::size_t{0}, corners.size())
1045 ✗ | std::views::transform([&] (std::size_t i) {
1046 ✗ return corners[GridFormat::Traits::DuneDetail::map_corner_index(gt, i)];
1047 }),
1048 reordered.begin()
1049 );
1050 ✗ } else {
1051
1/2
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
 12 std::ranges::copy(corners, reordered.begin());
1052 }
1053
1054 24 return std::make_tuple(std::move(gt), std::move(reordered));
1055 12 }
1056
1057 /*!
1058 * \ingroup PredefinedTraits
1059 * \brief Adapter around a Dune::GridFactory to be compatible with GridFormat::Concepts::GridFactory.
1060 * Can be used to export a grid from a reader directly into a Dune::GridFactory. For instance:
1061 * \code{.cpp}
1062 * GridFormat::Reader reader; reader.open(filename);
1063 * Dune::GridFactory<DuneGrid> factory;
1064 * {
1065 * GridFormat::Dune::GridFactoryAdapter adapter{factory};
1066 * reader.export_grid(adapter);
1067 * }
1068 * // ... use dune grid factory
1069 * \endcode
1070 */
1071 template<typename Grid>
1072 class GridFactoryAdapter {
1073 public:
1074 static constexpr int space_dim = Grid::dimensionworld;
1075 using ctype = typename Grid::ctype;
1076 using DuneFactory = ::Dune::GridFactory<Grid>;
1077
1078 1 GridFactoryAdapter(DuneFactory& factory)
1079 1 : _factory{factory}
1080 1 {}
1081
1082 template<std::size_t _space_dim>
1083 20 void insert_point(const std::array<ctype, _space_dim>& point) {
1084 20 ::Dune::FieldVector<ctype, space_dim> p;
1085
2/4
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
 20 std::copy_n(point.begin(), space_dim, p.begin());
1086
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 _factory.insertVertex(p);
1087 20 }
1088
1089 12 void insert_cell(const CellType& ct, const std::vector<std::size_t>& corners) {
1090
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 const auto [dune_gt, dune_corners] = to_dune<unsigned int>(ct, corners);
1091
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 _factory.insertElement(dune_gt, dune_corners);
1092 12 }
1093
1094 private:
1095 DuneFactory& _factory;
1096 };
1097
1098 } // namespace GridFormat::Dune
1099
1100 #endif // GRIDFORMAT_TRAITS_DUNE_HPP_
1101