GCC Code Coverage Report

Directory: gridformat/
File: gridformat/traits/dealii.hpp
Date: 2024-11-10 16:24:00
Exec Total Coverage
Lines: 62 73 84.9%
Functions: 96 96 100.0%
Branches: 55 134 41.0%
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
7 * <a href="https://www.dealii.org/current/doxygen/deal.II/group__grid.html">dealii triangulations</a>
8 */
9 #ifndef GRIDFORMAT_TRAITS_DEAL_II_HPP_
10 #define GRIDFORMAT_TRAITS_DEAL_II_HPP_
11
12 #include <cassert>
13 #include <ranges>
14 #include <vector>
15 #include <array>
16 #include <type_traits>
17
18 #include <deal.II/grid/tria.h>
19 #include <deal.II/distributed/tria.h>
20 #include <deal.II/distributed/fully_distributed_tria.h>
21
22 #include <gridformat/common/exceptions.hpp>
23 #include <gridformat/common/filtered_range.hpp>
24
25 #include <gridformat/grid/cell_type.hpp>
26 #include <gridformat/grid/traits.hpp>
27
28
29 namespace GridFormat {
30
31 #ifndef DOXYGEN
32 namespace DealII {
33
34 template<int dim, int space_dim> using Tria = dealii::Triangulation<dim, space_dim>;
35 template<int dim, int space_dim> using DTria = dealii::parallel::distributed::Triangulation<dim, space_dim>;
36 template<int dim, int space_dim> using FDTria = dealii::parallel::fullydistributed::Triangulation<dim, space_dim>;
37
38 template<typename T> struct IsTriangulation : public std::false_type {};
39 template<typename T> struct IsParallelTriangulation : public std::false_type {};
40
41 template<int d, int sd> struct IsTriangulation<Tria<d, sd>> : public std::true_type {};
42 template<int d, int sd> struct IsTriangulation<DTria<d, sd>> : public std::true_type {};
43 template<int d, int sd> struct IsTriangulation<FDTria<d, sd>> : public std::true_type {};
44 template<int d, int sd> struct IsParallelTriangulation<DTria<d, sd>> : public std::true_type {};
45 template<int d, int sd> struct IsParallelTriangulation<FDTria<d, sd>> : public std::true_type {};
46
47 template<typename T> inline constexpr bool is_triangulation = IsTriangulation<T>::value;
48 template<typename T> inline constexpr bool is_parallel_triangulation = IsParallelTriangulation<T>::value;
49
50 // We need to wrap the deal.ii iterators because they do not fulfill the
51 // requirements for the std::ranges concepts. The issue comes from the
52 // std::indirectly_readable concept (https://en.cppreference.com/w/cpp/iterator/indirectly_readable)
53 // which expects that the reference type obtained from Iterator& and const Iterator&
54 // is the same: https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2019/p1878r0.pdf.
55 // However, deal.ii iterators return refs or const refs depending on the constness of the iterator.
56 // For the traits, we only need read-access to cells/points, so we wrap the iterator and expose only
57 // the const interface in order to be compatible with std::ranges.
58 // Note: deal.ii iterators return copies when dereferenced, see e.g. here: https://www.dealii.org/current/doxygen/deal.II/classTriaRawIterator.html#a6e56c256239a9806bd372cfad7fba744
59 // However, Iterator::reference defines a reference, and therefore, we deduce the dereferenced type.
60 // TODO: revisit the above comment, this may no longer be true after reusing typename Iterator::reference
61 // we should revisit our iterator wrappers and hopefully no longer need to wrap dealii iterators
62 template<typename Iterator>
63 class ForwardIteratorWrapper
64 : public ForwardIteratorFacade<ForwardIteratorWrapper<Iterator>,
65 typename Iterator::value_type,
66 decltype(*(std::declval<std::add_const_t<Iterator>>())),
67 typename Iterator::pointer,
68 typename Iterator::difference_type> {
69 using Reference = decltype(*(std::declval<std::add_const_t<Iterator>>()));
70
71 public:
72 ForwardIteratorWrapper() = default;
73 648 ForwardIteratorWrapper(Iterator it) : _it{it} {}
74
75 // we need to explicitly implement the equality operators because with
76 // deal.ii we have compare iterator types that are not convertible into
77 // each other (which is a requirement for the generic implementation provided
78 // in GridFormat::IteratorFacade). For instance, active_cell_iterator has to
79 // be compared against cell_iterator which is used as sentinel.
80 template<typename I>
81 27056 friend bool operator==(const ForwardIteratorWrapper& self, const ForwardIteratorWrapper<I>& other) {
82 27056 return self._it == other._get_it();
83 }
84
85 template<typename I>
86 37040 friend bool operator!=(const ForwardIteratorWrapper& self, const ForwardIteratorWrapper<I>& other) {
87 37040 return self._it != other._get_it();
88 }
89
90 64096 const Iterator& _get_it() const {
91 64096 return _it;
92 }
93
94 private:
95 friend IteratorAccess;
96
97 42972 Reference _dereference() const { return *_it; }
98 bool _is_equal(const ForwardIteratorWrapper& other) { return _it == other._it; }
99 47568 void _increment() { ++_it; }
100
101 Iterator _it;
102 };
103
104 template<typename It>
105 ForwardIteratorWrapper(It&& it) -> ForwardIteratorWrapper<std::remove_cvref_t<It>>;
106
107 template<typename T> requires(is_triangulation<T>)
108 using Cell = typename T::active_cell_iterator::AccessorType;
109
110 template<typename T> requires(is_triangulation<T>)
111 using Point = typename T::active_vertex_iterator::AccessorType;
112
113 2932 const std::vector<int>& cell_corners_in_gridformat_order(unsigned int cell_dimension,
114 unsigned int number_of_cell_corners) {
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2932 times.
 2932 if (cell_dimension == 0) {
116 ✗ assert(number_of_cell_corners == 1);
117 ✗ static const std::vector<int> corners{0};
118 ✗ return corners;
119
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2932 times.
 2932 } else if (cell_dimension == 1) {
120 ✗ assert(number_of_cell_corners == 2);
121 ✗ static const std::vector<int> corners{0, 1};
122 ✗ return corners;
123
2/2
✓ Branch 0 taken 904 times.
✓ Branch 1 taken 2028 times.
 2932 } else if (cell_dimension == 2) {
124
2/2
✓ Branch 0 taken 584 times.
✓ Branch 1 taken 320 times.
 904 if (number_of_cell_corners == 3) { // triangles
125
4/8
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 583 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 586 static const std::vector<int> corners{0, 1, 2};
126 584 return corners;
127
1/2
✓ Branch 0 taken 320 times.
✗ Branch 1 not taken.
 320 } else if (number_of_cell_corners == 4) { // quads
128
4/8
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 317 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 326 static const std::vector<int> corners{0, 1, 3, 2};
129 320 return corners;
130 }
131
1/2
✓ Branch 0 taken 2028 times.
✗ Branch 1 not taken.
 2028 } else if (cell_dimension == 3) {
132
2/2
✓ Branch 0 taken 684 times.
✓ Branch 1 taken 1344 times.
 2028 if (number_of_cell_corners == 4) { // tets
133
4/8
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 683 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 686 static const std::vector<int> corners{0, 1, 2, 3};
134 684 return corners;
135
1/2
✓ Branch 0 taken 1344 times.
✗ Branch 1 not taken.
 1344 } else if (number_of_cell_corners == 8) { // hexes
136
4/8
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1341 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 1350 static const std::vector<int> corners{0, 1, 3, 2, 4, 5, 7, 6};
137 1344 return corners;
138 }
139 }
140
141 throw NotImplemented(
142 ✗ "deal.ii cell corner indices for cell of dimension " + std::to_string(cell_dimension) +
143 ✗ " and " + std::to_string(number_of_cell_corners) + " corners"
144 ✗ );
145 // In order to make compilers happy
146 static const std::vector<int> v{};
147 return v;
148 }
149
150 } // namespace DealII
151 #endif // DOXYGEN
152
153 // Specializations of the traits required for the `UnstructuredGrid` concept for dealii triangulations
154 namespace Traits {
155
156 template<typename T> requires(DealII::is_triangulation<T>)
157 struct Points<T> {
158 128 static std::ranges::range decltype(auto) get(const T& grid) {
159
2/4
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64 times.
✗ Branch 5 not taken.
 256 return std::ranges::subrange(
160
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 128 DealII::ForwardIteratorWrapper{grid.begin_active_vertex()},
161
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 256 DealII::ForwardIteratorWrapper{grid.end_vertex()}
162
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 256 );
163 }
164 };
165
166 template<typename T> requires(DealII::is_triangulation<T>)
167 struct Cells<T> {
168 196 static std::ranges::range auto get(const T& grid) {
169 if constexpr (DealII::is_parallel_triangulation<T>)
170 return Ranges::filter_by(
171 5376 [] (const auto& cell) { return cell.is_locally_owned(); },
172
3/6
✓ Branch 1 taken 62 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 62 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 62 times.
✗ Branch 8 not taken.
 248 std::ranges::subrange(
173
1/2
✓ Branch 1 taken 62 times.
✗ Branch 2 not taken.
 124 DealII::ForwardIteratorWrapper{grid.begin_active()},
174
1/2
✓ Branch 1 taken 62 times.
✗ Branch 2 not taken.
 248 DealII::ForwardIteratorWrapper{grid.end()}
175 )
176
1/2
✓ Branch 1 taken 62 times.
✗ Branch 2 not taken.
 248 );
177 else
178
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
 144 return std::ranges::subrange(
179
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 72 DealII::ForwardIteratorWrapper{grid.begin_active()},
180
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 144 DealII::ForwardIteratorWrapper{grid.end()}
181
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 144 );
182 }
183 };
184
185 template<typename T> requires(DealII::is_triangulation<T>)
186 struct CellType<T, DealII::Cell<T>> {
187 8244 static GridFormat::CellType get(const T&, const DealII::Cell<T>& cell) {
188 using CT = GridFormat::CellType;
189 static constexpr std::array cubes{CT::vertex, CT::segment, CT::quadrilateral, CT::hexahedron};
190 static constexpr std::array simplices{CT::vertex, CT::segment, CT::triangle, CT::tetrahedron};
191
1/2
✓ Branch 1 taken 4122 times.
✗ Branch 2 not taken.
 8244 const auto& ref_cell = cell.reference_cell();
192
2/2
✓ Branch 1 taken 3488 times.
✓ Branch 2 taken 634 times.
 8244 if (ref_cell.is_hyper_cube())
193
1/2
✓ Branch 1 taken 3488 times.
✗ Branch 2 not taken.
 6976 return cubes[ref_cell.get_dimension()];
194
1/2
✓ Branch 1 taken 634 times.
✗ Branch 2 not taken.
 1268 else if (ref_cell.is_simplex())
195
1/2
✓ Branch 1 taken 634 times.
✗ Branch 2 not taken.
 1268 return simplices[ref_cell.get_dimension()];
196 ✗ throw NotImplemented("CellType only implemented for hypercubes & simplices");
197 }
198 };
199
200 template<typename T> requires(DealII::is_triangulation<T>)
201 struct CellPoints<T, DealII::Cell<T>> {
202 5864 static std::ranges::range decltype(auto) get(const T&, const DealII::Cell<T>& cell) {
203 ✗ return DealII::cell_corners_in_gridformat_order(cell.reference_cell().get_dimension(), cell.n_vertices())
204
6/12
✓ Branch 1 taken 2932 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2932 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2932 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2932 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2932 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2932 times.
✗ Branch 17 not taken.
 13868 | std::views::transform([&] (int i) { return *cell.vertex_iterator(i); });
205 }
206 };
207
208 template<typename T> requires(DealII::is_triangulation<T>)
209 struct PointId<T, DealII::Point<T>> {
210 24224 static auto get(const T&, const DealII::Point<T>& point) {
211 24224 return point.index();
212 }
213 };
214
215 template<typename T> requires(DealII::is_triangulation<T>)
216 struct PointCoordinates<T, DealII::Point<T>> {
217 4108 static std::ranges::range decltype(auto) get(const T&, const DealII::Point<T>& point) {
218 static_assert(T::dimension >= 1 && T::dimension <= 3);
219 4108 const auto& center = point.center();
220 if constexpr (T::dimension == 1)
221 return std::array{center[0]};
222 else if constexpr (T::dimension == 2)
223 3536 return std::array{center[0], center[1]};
224 else if constexpr (T::dimension == 3)
225 16120 return std::array{center[0], center[1], center[2]};
226 }
227 };
228
229 template<typename T> requires(DealII::is_triangulation<T>)
230 struct NumberOfPoints<T> {
231 640 static std::integral auto get(const T& grid) {
232 640 return grid.n_used_vertices();
233 }
234 };
235
236 template<typename T> requires(DealII::is_triangulation<T>)
237 struct NumberOfCells<T> {
238 572 static std::integral auto get(const T& grid) {
239 if constexpr (DealII::is_parallel_triangulation<T>)
240 380 return grid.n_locally_owned_active_cells();
241 else
242 192 return grid.n_active_cells();
243 }
244 };
245
246 template<typename T> requires(DealII::is_triangulation<T>)
247 struct NumberOfCellPoints<T, DealII::Cell<T>> {
248 4366 static std::integral auto get(const T&, const DealII::Cell<T>& cell) {
249 4366 return cell.n_vertices();
250 }
251 };
252
253 } // namespace Traits
254 } // namespace GridFormat
255
256 #endif // GRIDFORMAT_TRAITS_DEAL_II_HPP_
257