SimiLie
Loading...
Searching...
No Matches
boundary.hpp
1// SPDX-FileCopyrightText: 2024 Baptiste Legouix
2// SPDX-License-Identifier: AGPL-3.0-or-later
3
4#pragma once
5
6#include <array>
7
8#include <ddc/ddc.hpp>
9
10#include <similie/misc/specialization.hpp>
11
12#include "chain.hpp"
13#include "local_chain.hpp"
14
15namespace sil {
16
17namespace exterior {
18
19namespace detail {
20
21template <class SimplexType>
22struct BoundaryType;
23
24template <std::size_t K, class... Tag>
25struct BoundaryType<Simplex<K, Tag...>>
26{
27 using type = Simplex<K - 1, Tag...>;
28};
29
30template <class SimplexType, class LayoutStridedPolicy, class MemorySpace>
31struct BoundaryType<Chain<SimplexType, LayoutStridedPolicy, MemorySpace>>
32{
33 using type = Chain<typename BoundaryType<SimplexType>::type, LayoutStridedPolicy, MemorySpace>;
34};
35
36} // namespace detail
37
38template <class T>
39using boundary_t = typename detail::BoundaryType<T>::type;
40
41namespace detail {
42
43template <class MemorySpace = Kokkos::HostSpace, class SimplexType>
44KOKKOS_FUNCTION constexpr LocalChain<boundary_t<SimplexType>, Kokkos::LayoutRight, MemorySpace>
45generate_local_half_subchain(
46 typename SimplexType::discrete_element_type elem,
47 typename SimplexType::discrete_vector_type vect,
48 bool negative = false)
49{
50 auto array = ddc::detail::array(vect);
51 LocalChain<boundary_t<SimplexType>, Kokkos::LayoutRight, MemorySpace> chain(elem);
52 auto id_dist = -1;
53 for (std::size_t i = 0; i < SimplexType::dimension(); ++i) {
54 auto array_ = array;
55 auto j = array_.begin() + id_dist + 1;
56 auto id = std::find_if(j, array_.end(), [](int k) { return k != 0; });
57 id_dist = std::distance(array_.begin(), id);
58 *id = 0;
59 typename SimplexType::discrete_vector_type vect_;
60 ddc::detail::array(vect_) = array_;
61 chain += boundary_t<SimplexType>(elem, vect_, (negative + i) % 2);
62 j = id + 1;
63 }
64 return chain;
65}
66
67// TODO Kokkosify
68template <class SimplexType, misc::Specialization<Kokkos::View> AllocationType>
69KOKKOS_FUNCTION constexpr Chain<
70 boundary_t<SimplexType>,
71 typename AllocationType::array_layout,
72 typename AllocationType::memory_space>
73generate_half_subchain(
74 AllocationType allocation,
75 typename SimplexType::discrete_element_type elem,
76 typename SimplexType::discrete_vector_type vect,
77 bool negative = false)
78{
79 auto array = ddc::detail::array(vect);
80 Chain<boundary_t<SimplexType>,
81 typename AllocationType::array_layout,
82 typename AllocationType::memory_space>
83 chain(allocation);
84 auto id_dist = -1;
85 for (std::size_t i = 0; i < SimplexType::dimension(); ++i) {
86 auto array_ = array;
87 auto j = array_.begin() + id_dist + 1;
88 auto id = std::find_if(j, array_.end(), [](int k) { return k != 0; });
89 id_dist = std::distance(array_.begin(), id);
90 *id = 0;
91 typename SimplexType::discrete_vector_type vect_;
92 ddc::detail::array(vect_) = array_;
93 chain += boundary_t<SimplexType>(elem, vect_, (negative + i) % 2);
94 j = id + 1;
95 }
96 return chain;
97}
98
99} // namespace detail
100
101template <misc::Specialization<Kokkos::View> AllocationType, class SimplexType>
102struct Boundary
103{
104 KOKKOS_FUNCTION static Chain<
105 boundary_t<SimplexType>,
106 typename AllocationType::array_layout,
107 typename AllocationType::memory_space>
108 run(AllocationType allocation, SimplexType simplex)
109 {
110 Chain<boundary_t<SimplexType>,
111 typename AllocationType::array_layout,
112 typename AllocationType::memory_space>
113 chain(allocation);
114 detail::generate_half_subchain<SimplexType>(
115 Kokkos::
116 subview(allocation,
117 std::pair<std::size_t, std::size_t>(0, SimplexType::dimension())),
118 simplex.discrete_element(),
119 simplex.discrete_vector(),
120 SimplexType::dimension() % 2);
121 chain += SimplexType::dimension();
122 detail::generate_half_subchain<SimplexType>(
123 Kokkos::
124 subview(allocation,
125 std::pair<std::size_t, std::size_t>(
126 SimplexType::dimension(),
127 2 * SimplexType::dimension())),
128 simplex.discrete_element() + simplex.discrete_vector(),
129 -simplex.discrete_vector());
130 chain += SimplexType::dimension();
131 chain *= (SimplexType::dimension() % 2 ? 1 : -1) * (simplex.negative() ? -1 : 1);
132 return chain;
133 }
134
135 KOKKOS_FUNCTION static Chain<
136 boundary_t<SimplexType>,
137 typename AllocationType::array_layout,
138 typename AllocationType::memory_space>
139 run(AllocationType allocation,
140 Chain<SimplexType,
141 typename AllocationType::array_layout,
142 typename AllocationType::memory_space> chain)
143 {
144 Chain<boundary_t<SimplexType>,
145 typename AllocationType::array_layout,
146 typename AllocationType::memory_space>
147 boundary_chain(allocation);
148 for (auto i = chain.begin(); i < chain.end(); ++i) {
149 std::size_t const distance = Kokkos::Experimental::distance(chain.begin(), i);
150 Boundary::
151 run(Kokkos::
152 subview(allocation,
153 std::pair<std::size_t, std::size_t>(
154 2 * SimplexType::dimension() * distance,
155 2 * SimplexType::dimension() * (distance + 1))),
156 *i);
157 boundary_chain += 2 * SimplexType::dimension();
158 }
159 boundary_chain.optimize();
160 return boundary_chain;
161 }
162};
163
164template <misc::Specialization<Kokkos::View> AllocationType, class SimplexType>
165KOKKOS_FUNCTION Chain<
166 boundary_t<SimplexType>,
167 typename AllocationType::array_layout,
168 typename AllocationType::memory_space>
169boundary(AllocationType allocation, SimplexType simplex)
170{
171 return Boundary<AllocationType, SimplexType>::run(allocation, simplex);
172}
173
174template <class MemorySpace = Kokkos::HostSpace, class SimplexType>
175KOKKOS_FUNCTION LocalChain<boundary_t<SimplexType>, Kokkos::LayoutRight, MemorySpace> boundary(
176 SimplexType simplex)
177{
178 LocalChain<boundary_t<SimplexType>, Kokkos::LayoutRight, MemorySpace> chain(
179 simplex.discrete_element());
180 chain += detail::generate_local_half_subchain<MemorySpace, SimplexType>(
181 simplex.discrete_element(),
182 simplex.discrete_vector(),
183 SimplexType::dimension() % 2);
184 chain += detail::generate_local_half_subchain<MemorySpace, SimplexType>(
185 simplex.discrete_element() + simplex.discrete_vector(),
186 -simplex.discrete_vector());
187 chain *= (SimplexType::dimension() % 2 ? 1 : -1) * (simplex.negative() ? -1 : 1);
188 return chain;
189}
190
191template <misc::Specialization<Kokkos::View> AllocationType, class SimplexType>
192KOKKOS_FUNCTION Chain<
193 boundary_t<SimplexType>,
194 typename AllocationType::array_layout,
195 typename AllocationType::memory_space>
196boundary(
197 AllocationType allocation,
198 Chain<SimplexType,
199 typename AllocationType::array_layout,
200 typename AllocationType::memory_space> chain)
201{
202 return Boundary<AllocationType, SimplexType>::run(allocation, chain);
203}
204
205} // namespace exterior
206
207} // namespace sil
sil::exterior::Chain
Chain class.
Definition chain.hpp:26
sil::exterior::Chain::optimize
KOKKOS_FUNCTION void optimize()
Definition chain.hpp:132
sil::exterior::Simplex
Simplex class.
Definition simplex.hpp:75
sil::exterior::Chain
Chain(Head, Tail...) -> Chain< typename Head::value_type, typename Head::array_layout, typename Head::memory_space >
sil::exterior::boundary_t
typename detail::BoundaryType< T >::type boundary_t
Definition boundary.hpp:39
sil::exterior::boundary
KOKKOS_FUNCTION Chain< boundary_t< SimplexType >, typename AllocationType::array_layout, typename AllocationType::memory_space > boundary(AllocationType allocation, SimplexType simplex)
Definition boundary.hpp:169
sil
The top-level namespace of SimiLie.
Definition csr.hpp:15
sil::exterior::Boundary::run
static KOKKOS_FUNCTION Chain< boundary_t< SimplexType >, typename AllocationType::array_layout, typename AllocationType::memory_space > run(AllocationType allocation, Chain< SimplexType, typename AllocationType::array_layout, typename AllocationType::memory_space > chain)
Definition boundary.hpp:139
sil::exterior::Boundary::run
static KOKKOS_FUNCTION Chain< boundary_t< SimplexType >, typename AllocationType::array_layout, typename AllocationType::memory_space > run(AllocationType allocation, SimplexType simplex)
Definition boundary.hpp:108