9#ifndef GRIDFORMAT_TRAITS_DEAL_II_HPP_
10#define GRIDFORMAT_TRAITS_DEAL_II_HPP_
18#include <deal.II/grid/tria.h>
19#include <deal.II/distributed/tria.h>
20#include <deal.II/distributed/fully_distributed_tria.h>
22#include <gridformat/common/exceptions.hpp>
23#include <gridformat/common/filtered_range.hpp>
25#include <gridformat/grid/cell_type.hpp>
26#include <gridformat/grid/traits.hpp>
34template<
int dim,
int space_dim>
using Tria = dealii::Triangulation<dim, space_dim>;
35template<
int dim,
int space_dim>
using DTria = dealii::parallel::distributed::Triangulation<dim, space_dim>;
36template<
int dim,
int space_dim>
using FDTria = dealii::parallel::fullydistributed::Triangulation<dim, space_dim>;
38template<
typename T>
struct IsTriangulation :
public std::false_type {};
39template<
typename T>
struct IsParallelTriangulation :
public std::false_type {};
41template<
int d,
int sd>
struct IsTriangulation<Tria<d, sd>> :
public std::true_type {};
42template<
int d,
int sd>
struct IsTriangulation<DTria<d, sd>> :
public std::true_type {};
43template<
int d,
int sd>
struct IsTriangulation<FDTria<d, sd>> :
public std::true_type {};
44template<
int d,
int sd>
struct IsParallelTriangulation<DTria<d, sd>> :
public std::true_type {};
45template<
int d,
int sd>
struct IsParallelTriangulation<FDTria<d, sd>> :
public std::true_type {};
47template<
typename T>
inline constexpr bool is_triangulation = IsTriangulation<T>::value;
48template<
typename T>
inline constexpr bool is_parallel_triangulation = IsParallelTriangulation<T>::value;
58template<
typename Iterator>
59class ForwardIteratorWrapper
60:
public ForwardIteratorFacade<ForwardIteratorWrapper<Iterator>,
61 typename Iterator::value_type,
62 const typename Iterator::value_type&,
63 typename Iterator::pointer,
64 typename Iterator::difference_type> {
65 using Reference =
const typename Iterator::value_type&;
68 ForwardIteratorWrapper() =
default;
69 ForwardIteratorWrapper(Iterator it) : _it{it} {}
77 friend bool operator==(
const ForwardIteratorWrapper& self,
const ForwardIteratorWrapper<I>& other) {
78 return self._it == other._get_it();
82 friend bool operator!=(
const ForwardIteratorWrapper& self,
const ForwardIteratorWrapper<I>& other) {
83 return self._it != other._get_it();
86 const Iterator& _get_it()
const {
91 friend IteratorAccess;
93 Reference _dereference()
const {
return *_it; }
94 bool _is_equal(
const ForwardIteratorWrapper& other) {
return _it == other._it; }
95 void _increment() { ++_it; }
101ForwardIteratorWrapper(It&& it) -> ForwardIteratorWrapper<std::remove_cvref_t<It>>;
103template<
typename T>
requires(is_triangulation<T>)
104using Cell =
typename T::active_cell_iterator::AccessorType;
106template<
typename T>
requires(is_triangulation<T>)
107using Point =
typename T::active_vertex_iterator::AccessorType;
109const std::vector<int>& cell_corners_in_gridformat_order(
unsigned int cell_dimension,
110 unsigned int number_of_cell_corners) {
111 if (cell_dimension == 0) {
112 assert(number_of_cell_corners == 1);
113 static const std::vector<int> corners{0};
115 }
else if (cell_dimension == 1) {
116 assert(number_of_cell_corners == 2);
117 static const std::vector<int> corners{0, 1};
119 }
else if (cell_dimension == 2) {
120 if (number_of_cell_corners == 3) {
121 static const std::vector<int> corners{0, 1, 2};
123 }
else if (number_of_cell_corners == 4) {
124 static const std::vector<int> corners{0, 1, 3, 2};
127 }
else if (cell_dimension == 3) {
128 if (number_of_cell_corners == 4) {
129 static const std::vector<int> corners{0, 1, 2, 3};
131 }
else if (number_of_cell_corners == 8) {
132 static const std::vector<int> corners{0, 1, 3, 2, 4, 5, 7, 6};
137 throw NotImplemented(
138 "deal.ii cell corner indices for cell of dimension " + std::to_string(cell_dimension) +
139 " and " + std::to_string(number_of_cell_corners) +
" corners"
142 static const std::vector<int> v{};
152template<
typename T>
requires(DealII::is_triangulation<T>)
154 static std::ranges::range
decltype(
auto) get(
const T& grid) {
155 return std::ranges::subrange(
156 DealII::ForwardIteratorWrapper{grid.begin_active_vertex()},
157 DealII::ForwardIteratorWrapper{grid.end_vertex()}
162template<
typename T>
requires(DealII::is_triangulation<T>)
164 static std::ranges::range
auto get(
const T& grid) {
165 if constexpr (DealII::is_parallel_triangulation<T>)
166 return Ranges::filter_by(
167 [] (
const auto& cell) {
return cell.is_locally_owned(); },
168 std::ranges::subrange(
169 DealII::ForwardIteratorWrapper{grid.begin_active()},
170 DealII::ForwardIteratorWrapper{grid.end()}
174 return std::ranges::subrange(
175 DealII::ForwardIteratorWrapper{grid.begin_active()},
176 DealII::ForwardIteratorWrapper{grid.end()}
181template<
typename T>
requires(DealII::is_triangulation<T>)
182struct CellType<T, DealII::Cell<T>> {
183 static GridFormat::CellType get(
const T&,
const DealII::Cell<T>& cell) {
184 using CT = GridFormat::CellType;
185 static constexpr std::array cubes{CT::vertex, CT::segment, CT::quadrilateral, CT::hexahedron};
186 static constexpr std::array simplices{CT::vertex, CT::segment, CT::triangle, CT::tetrahedron};
187 const auto& ref_cell = cell.reference_cell();
188 if (ref_cell.is_hyper_cube())
189 return cubes[ref_cell.get_dimension()];
190 else if (ref_cell.is_simplex())
191 return simplices[ref_cell.get_dimension()];
192 throw NotImplemented(
"CellType only implemented for hypercubes & simplices");
196template<
typename T>
requires(DealII::is_triangulation<T>)
197struct CellPoints<T, DealII::Cell<T>> {
198 static std::ranges::range
decltype(
auto) get(
const T&,
const DealII::Cell<T>& cell) {
199 return DealII::cell_corners_in_gridformat_order(cell.reference_cell().get_dimension(), cell.n_vertices())
200 | std::views::transform([&] (
int i) {
return *cell.vertex_iterator(i); });
204template<
typename T>
requires(DealII::is_triangulation<T>)
205struct PointId<T, DealII::Point<T>> {
206 static auto get(
const T&,
const DealII::Point<T>& point) {
207 return point.index();
211template<
typename T>
requires(DealII::is_triangulation<T>)
212struct PointCoordinates<T, DealII::Point<T>> {
213 static std::ranges::range
decltype(
auto) get(
const T&,
const DealII::Point<T>& point) {
214 static_assert(T::dimension >= 1 && T::dimension <= 3);
215 const auto& center = point.center();
216 if constexpr (T::dimension == 1)
217 return std::array{center[0]};
218 else if constexpr (T::dimension == 2)
219 return std::array{center[0], center[1]};
220 else if constexpr (T::dimension == 3)
221 return std::array{center[0], center[1], center[2]};
225template<
typename T>
requires(DealII::is_triangulation<T>)
226struct NumberOfPoints<T> {
227 static std::integral
auto get(
const T& grid) {
228 return grid.n_used_vertices();
232template<
typename T>
requires(DealII::is_triangulation<T>)
233struct NumberOfCells<T> {
234 static std::integral
auto get(
const T& grid) {
235 if constexpr (DealII::is_parallel_triangulation<T>)
236 return grid.n_locally_owned_active_cells();
238 return grid.n_active_cells();
242template<
typename T>
requires(DealII::is_triangulation<T>)
243struct NumberOfCellPoints<T, DealII::Cell<T>> {
244 static std::integral
auto get(
const T&,
const DealII::Cell<T>& cell) {
245 return cell.n_vertices();