PACMAN 0.1.0
Portable Algorithms for Coupling, Mapping, and Adaptive iNterpolation
Loading...
Searching...
No Matches
ComputeLinearSkin.hpp
Go to the documentation of this file.
1//
2// This file is subject to the terms and conditions defined in
3// file 'LICENSE', which is part of this source code package.
4//
5
6#pragma once
7
8#include <Kokkos_Core.hpp>
9#include <stdexcept>
10
11#include "common/concepts.hpp"
12#include "common/transfer.hxx"
13#include "common/types.hpp"
14
15namespace PACMAN {
16namespace FiniteElements {
17
23template <KokkosViewRank<1> T, int_t Dim> struct permute {
26 KOKKOS_INLINE_FUNCTION
27 static void permuteCol(T kokkos_sv) {
28 if constexpr (Dim >= 2) {
29 if (kokkos_sv(0) > kokkos_sv(1)) {
30 Kokkos::kokkos_swap(kokkos_sv(0), kokkos_sv(1));
31 }
32 }
33
34 if constexpr (Dim >= 3) {
35 if (kokkos_sv(2) > kokkos_sv(3)) {
36 Kokkos::kokkos_swap(kokkos_sv(2), kokkos_sv(3));
37 }
38 if (kokkos_sv(0) > kokkos_sv(2)) {
39 Kokkos::kokkos_swap(kokkos_sv(0), kokkos_sv(2));
40 }
41 if (kokkos_sv(1) > kokkos_sv(3)) {
42 Kokkos::kokkos_swap(kokkos_sv(1), kokkos_sv(3));
43 }
44 if (kokkos_sv(1) > kokkos_sv(2)) {
45 Kokkos::kokkos_swap(kokkos_sv(1), kokkos_sv(2));
46 }
47 }
48 }
49};
50
67template <typename ExecSpace, int_t Dim>
68void ComputeLinearSkin(Transfer<ExecSpace, Dim> &transfer) {
69 Kokkos::Profiling::pushRegion("Compute Linear Skin");
70
71 using MemorySpace = typename ExecSpace::memory_space;
72 ExecSpace execSpace{};
73
74 // TODO: Rewrite -> Ugly !
75 auto LFPE_val_vec = getLinearFacesPaddedEntries(Dim);
76 auto LFPE_val = Kokkos::View<shortint_t *,
77 Kokkos::DefaultHostExecutionSpace::memory_space,
78 Kokkos::MemoryTraits<Kokkos::Unmanaged>>(
79 LFPE_val_vec.data(), LFPE_val_vec.size());
80 auto linearFacesPadded_val =
81 Kokkos::create_mirror_view_and_copy(execSpace, LFPE_val);
82 auto LFPE_off_vec = getLinearFacesPaddedOffsets(Dim);
83 auto LFPE_off =
84 Kokkos::View<offset_t *, Kokkos::DefaultHostExecutionSpace::memory_space,
85 Kokkos::MemoryTraits<Kokkos::Unmanaged>>(
86 LFPE_off_vec.data(), LFPE_off_vec.size());
87 auto linearFacesPadded_off =
88 Kokkos::create_mirror_view_and_copy(execSpace, LFPE_off);
89
90 auto TF_val_vec = getTriFacesEntries(Dim);
91 auto TF_val = Kokkos::View<shortint_t *,
92 Kokkos::DefaultHostExecutionSpace::memory_space,
93 Kokkos::MemoryTraits<Kokkos::Unmanaged>>(
94 TF_val_vec.data(), TF_val_vec.size());
95 auto TriFaces_val = Kokkos::create_mirror_view_and_copy(execSpace, TF_val);
96 auto TF_off_vec = getTriFacesOffsets(Dim);
97 auto TF_off =
98 Kokkos::View<offset_t *, Kokkos::DefaultHostExecutionSpace::memory_space,
99 Kokkos::MemoryTraits<Kokkos::Unmanaged>>(TF_off_vec.data(),
100 TF_off_vec.size());
101 auto TriFaces_off = Kokkos::create_mirror_view_and_copy(execSpace, TF_off);
102
103 auto TLF_val_vec = getTriLocFacesEntries(Dim);
104 auto TLF_val = Kokkos::View<shortint_t *,
105 Kokkos::DefaultHostExecutionSpace::memory_space,
106 Kokkos::MemoryTraits<Kokkos::Unmanaged>>(
107 TLF_val_vec.data(), TLF_val_vec.size());
108 auto TriLocFaces_val =
109 Kokkos::create_mirror_view_and_copy(execSpace, TLF_val);
110 auto TLF_off_vec = getTriLocFacesOffsets(Dim);
111 auto TLF_off =
112 Kokkos::View<offset_t *, Kokkos::DefaultHostExecutionSpace::memory_space,
113 Kokkos::MemoryTraits<Kokkos::Unmanaged>>(TLF_off_vec.data(),
114 TLF_off_vec.size());
115 auto TriLocFaces_off =
116 Kokkos::create_mirror_view_and_copy(execSpace, TLF_off);
117
118 auto nbNode = transfer.sourcePoints.extent_int(0);
119 auto nbLinFaces = transfer.nbLinearSkinFaces;
120 auto CellTypesPtr = transfer.cellTypes;
121 auto connValPtr = transfer.connValues;
122 auto connOffPtr = transfer.connOffsets;
123
124 auto nbElem = CellTypesPtr.extent(0);
125
126 int_t nodesPerFace = -1;
127 switch (Dim) {
128 case 0:
129 nodesPerFace = 0;
130 break;
131 case 1:
132 nodesPerFace = 1;
133 break;
134 case 2:
135 nodesPerFace = 2;
136 break;
137 case 3:
138 nodesPerFace = 4;
139 break;
140 default:
141 std::cerr << "Invalid Dim " << Dim << std::endl;
142 assert(false);
143 break;
144 }
145
146 auto linFaces = Kokkos::View<int_t **, ExecSpace>(
147 Kokkos::view_alloc(execSpace, Kokkos::WithoutInitializing, "linFaces"),
148 nbLinFaces, nodesPerFace + 1); // last column is CellType
149 auto allFaceId = Kokkos::View<int_t *, ExecSpace>(
150 Kokkos::view_alloc(execSpace, Kokkos::WithoutInitializing, "allFaceId"),
151 nbLinFaces); // global linface id and local element face id
152 auto allFaceLocId = Kokkos::View<int_t *, ExecSpace>(
153 Kokkos::view_alloc(execSpace, Kokkos::WithoutInitializing,
154 "allFaceLocId"),
155 nbLinFaces); // global linface id and local element face id
156 auto allParents = Kokkos::View<int_t *, ExecSpace>(
157 Kokkos::view_alloc(execSpace, Kokkos::WithoutInitializing, "allParents"),
158 nbLinFaces);
159
160 constexpr int_t maxi = std::numeric_limits<int_t>::max();
161
162 Kokkos::Profiling::pushRegion("Fill all faces and sort per row permute");
163 // Avoid branching but very ugly: to change !
164 Kokkos::parallel_scan(
165 "Fill all faces and sort per row : Col permutation",
166 Kokkos::RangePolicy(execSpace, 0, nbElem),
167 KOKKOS_LAMBDA(const index_t &i, index_t &partial_sum, bool is_final) {
168 auto ielType = static_cast<cell_t>(CellTypesPtr(i));
169 auto facesPadded = Kokkos::subview(
170 linearFacesPadded_val,
171 Kokkos::make_pair(linearFacesPadded_off(ielType),
172 linearFacesPadded_off(ielType + 1)));
173 auto localnbFaces = facesPadded.extent_int(0) / nodesPerFace;
174 auto firstIndex = partial_sum;
175 partial_sum += localnbFaces;
176 if (is_final) {
177 auto localConnectivity = Kokkos::subview(
178 connValPtr, Kokkos::make_pair(connOffPtr(i), connOffPtr(i + 1)));
179 for (index_t j = 0; j < localnbFaces; j++) { // loop on element faces
180 auto linFaces_sv =
181 Kokkos::subview(linFaces, firstIndex + j, Kokkos::ALL());
182 allParents(firstIndex + j) = i;
183 allFaceId(firstIndex + j) = firstIndex + j;
184 allFaceLocId(firstIndex + j) = j;
185 for (index_t k = 0; k < nodesPerFace - 1;
186 k++) { // loop on non-padded face connectivity
187 auto index = facesPadded(j * nodesPerFace + k);
188 linFaces_sv(k) = localConnectivity(index);
189 }
190 auto index = facesPadded((j + 1) * nodesPerFace - 1);
191 linFaces_sv(nodesPerFace - 1) =
192 (index == -1) ? maxi : localConnectivity(index);
193 using subviewType = decltype(linFaces_sv);
194 permute<subviewType, Dim>::permuteCol(linFaces_sv);
195 linFaces_sv(nodesPerFace) = ielType;
196 }
197 }
198 });
199 Kokkos::Profiling::popRegion();
200
201 Kokkos::Profiling::pushRegion(
202 "Create permutation vector to sort linFaces by row");
203 auto firstNodeIds = Kokkos::subview(linFaces, Kokkos::ALL(), 0);
204
205 using IdViewType = decltype(firstNodeIds);
206 using CompType = Kokkos::BinOp1D<IdViewType>;
207 using DeviceType = Kokkos::Device<ExecSpace, MemorySpace>;
208
209 Kokkos::BinSort<IdViewType, CompType, DeviceType, int_t> bin_sort(
210 firstNodeIds, CompType(nbLinFaces, 0, nbNode), true /*sort within*/);
211 bin_sort.create_permute_vector();
212 auto permutationVector = bin_sort.get_permute_vector();
213 Kokkos::Profiling::popRegion();
214
215 Kokkos::Profiling::pushRegion("Permute and Compute Nbtris");
216
217 index_t nbTris = 0;
218 Kokkos::parallel_reduce(
219 "Permute and Compute NbTris",
220 Kokkos::RangePolicy(execSpace, 0, nbLinFaces),
221 KOKKOS_LAMBDA(const index_t &i, index_t &sum) {
222 auto index = permutationVector(i);
223
224 // vector who would be at pos "i" in the (first column)-sorted matrix
225 // faces
226 auto vect = Kokkos::subview(linFaces, index, Kokkos::ALL());
227
228 index_t curr_row_occur =
229 1; // Number of occurrences of the row (excluding current)
230 index_t curr = i;
231 // Find the first row starting with vect(0)
232 while (curr > 0 && linFaces(permutationVector(curr), 0) == vect(0)) {
233 --curr;
234 }
235 // In case curr is 0, we do not increase
236 if (linFaces(permutationVector(curr), 0) != vect(0)) {
237 ++curr;
238 }
239 // For all elements starting with vect(0)
240 while (curr <= nbLinFaces - 1 &&
241 linFaces(permutationVector(curr), 0) == vect(0)) {
242 if (curr == i) { // skip current
243 ++curr;
244 continue;
245 }
246 bool diff = false;
247 // Check if identical
248 for (index_t t = 1; t < nodesPerFace + 1; ++t)
249 if (linFaces(permutationVector(curr), t) != vect(t)) {
250 diff = true;
251 break;
252 }
253 if (!diff) {
254 ++curr_row_occur;
255 // if (curr_row_occur == 2)
256 break;
257 }
258 ++curr;
259 }
260 if (curr_row_occur == 2) {
261 allFaceId(index) = -1;
262 } else {
263 auto ielType = (int)CellTypesPtr(allParents(index));
264 auto iTriLocFaces = Kokkos::subview(
265 TriLocFaces_val, Kokkos::make_pair(TriLocFaces_off(ielType),
266 TriLocFaces_off(ielType + 1)));
267 auto iloc = allFaceLocId(index);
268 sum += (iTriLocFaces(iloc + 1) - iTriLocFaces(iloc)) / Dim;
269 }
270 },
271 nbTris);
272 Kokkos::Profiling::popRegion();
273
274 Kokkos::Profiling::pushRegion("Remove");
275
276 // remove the -1 values while saving the relative order of elements (in place)
277 const auto end = Kokkos::Experimental::remove(execSpace, allFaceId, -1);
278 const int dist =
279 Kokkos::Experimental::distance(Kokkos::Experimental::begin(allFaceId),
280 end); // number of not -1
281
282 Kokkos::Profiling::popRegion();
283
284 Kokkos::Profiling::pushRegion("Finalize: compute linear skin");
285
286 auto skinFaces = Kokkos::View<int_t **, ExecSpace>(
287 Kokkos::view_alloc(execSpace, Kokkos::WithoutInitializing, "skinFaces"),
288 nbTris, Dim);
289 auto skinParents = Kokkos::View<int_t *, ExecSpace>(
290 Kokkos::view_alloc(execSpace, Kokkos::WithoutInitializing, "skinParents"),
291 nbTris);
292
293 Kokkos::parallel_scan(
294 "Fill Skin Faces & Parents", Kokkos::RangePolicy(execSpace, 0, dist),
295 KOKKOS_LAMBDA(const index_t &i, index_t &partial_sum, bool is_final) {
296 auto index = allFaceId(i);
297 auto iParent = allParents(index);
298 auto ielType = static_cast<cell_t>(CellTypesPtr(iParent));
299 auto iTriFaces = Kokkos::subview(
300 TriFaces_val, Kokkos::make_pair(TriFaces_off(ielType),
301 TriFaces_off(ielType + 1)));
302 auto iTriLocFaces = Kokkos::subview(
303 TriLocFaces_val, Kokkos::make_pair(TriLocFaces_off(ielType),
304 TriLocFaces_off(ielType + 1)));
305 auto iloc = allFaceLocId(index);
306 auto localnbFaces = (iTriLocFaces(iloc + 1) - iTriLocFaces(iloc)) / Dim;
307 auto firstIndex = partial_sum;
308 partial_sum += localnbFaces;
309 if (is_final) {
310 auto localConnectivity = Kokkos::subview(
311 connValPtr,
312 Kokkos::make_pair(connOffPtr(iParent), connOffPtr(iParent + 1)));
313 for (index_t j = 0; j < localnbFaces; j++) {
314 skinParents(firstIndex + j) = iParent;
315 auto pos = iTriLocFaces(iloc) + j * Dim;
316 for (int_t k = 0; k < Dim; k++) {
317 const auto iconn = iTriFaces(pos + k);
318 skinFaces(firstIndex + j, k) = localConnectivity(iconn);
319 }
320 }
321 }
322 });
323 transfer.skinFaces = skinFaces;
324 transfer.skinParents = skinParents;
325
326 Kokkos::Profiling::popRegion();
327
328 Kokkos::Profiling::popRegion();
329}
330
331} // namespace FiniteElements
332
333} // namespace PACMAN
PACMAN::FiniteElements::getLinearFacesPaddedOffsets
std::vector< offset_t > getLinearFacesPaddedOffsets(const int_t dimension)
Get padded linear-face offset table for a given dimension.
Definition FESkinTools.hpp:475
PACMAN::FiniteElements::getLinearFacesPaddedEntries
std::vector< shortint_t > getLinearFacesPaddedEntries(const int_t dimension)
Get padded linear-face entry table for a given dimension.
Definition FESkinTools.hpp:495
PACMAN::FiniteElements::getTriLocFacesOffsets
std::vector< offset_t > getTriLocFacesOffsets(const int_t dimension)
Get local triangulated-face offset table for a given dimension.
Definition FESkinTools.hpp:1033
PACMAN::FiniteElements::ComputeLinearSkin
void ComputeLinearSkin(Transfer< ExecSpace, Dim > &transfer)
Build tetrahedralized boundary skin from source elements.
Definition ComputeLinearSkin.hpp:68
PACMAN::FiniteElements::getTriFacesEntries
std::vector< shortint_t > getTriFacesEntries(const int_t dimension)
Get triangulated-face entry table for a given dimension.
Definition FESkinTools.hpp:801
PACMAN::FiniteElements::getTriLocFacesEntries
std::vector< shortint_t > getTriLocFacesEntries(const int_t dimension)
Get local triangulated-face entry table for a given dimension.
Definition FESkinTools.hpp:1051
PACMAN::FiniteElements::getTriFacesOffsets
std::vector< offset_t > getTriFacesOffsets(const int_t dimension)
Get triangulated-face offset table for a given dimension.
Definition FESkinTools.hpp:783
PACMAN::cell_t
int32_t cell_t
Definition types.hpp:25
PACMAN::index_t
uint32_t index_t
Definition types.hpp:20
PACMAN::offset_t
int32_t offset_t
Definition types.hpp:21
PACMAN::shortint_t
int32_t shortint_t
Definition types.hpp:19
PACMAN::int_t
int32_t int_t
Definition types.hpp:18
PACMAN::FiniteElements::permute
Small helper to canonicalize face-node ordering for row comparison.
Definition ComputeLinearSkin.hpp:23
PACMAN::FiniteElements::permute::permuteCol
static KOKKOS_INLINE_FUNCTION void permuteCol(T kokkos_sv)
In-place permutation used to keep equivalent faces sortable.
Definition ComputeLinearSkin.hpp:27