PACMAN 0.1.0
Portable Algorithms for Coupling, Mapping, and Adaptive iNterpolation
Loading...
Searching...
No Matches
SkinProjection.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 <ArborX_Point.hpp>
9#include <ArborX_Triangle.hpp>
10#include <Kokkos_Core.hpp>
11#include <limits>
12
13#include "common/types.hpp"
14#include "finite_elements/utils/ArborXCallbacks.hpp"
15
16namespace PACMAN {
17namespace FiniteElements {
18
26KOKKOS_INLINE_FUNCTION ArborX::Point<2, coordinates_t>
27closest_point_to_edge(const ArborX::Point<2, coordinates_t> &p,
28 const ArborX::Point<2, coordinates_t> &a,
29 const ArborX::Point<2, coordinates_t> &b) {
30 const auto ab = b - a;
31 const auto ap = p - a;
32 coordinates_t t =
33 (ap[0] * ab[0] + ap[1] * ab[1]) / (ab[0] * ab[0] + ab[1] * ab[1]);
34 auto tclamp = Kokkos::max(0.0, Kokkos::min(1.0, t));
35 ArborX::Point<2, coordinates_t> ret = {a[0] + ab[0] * tclamp,
36 a[1] + ab[1] * tclamp};
37 return ret;
38};
39
51template <typename ExecSpace>
52void ComputeProjectionOn3DSkin(Transfer<ExecSpace, 3> &transfer,
53 bool extrapol = false) {
54 Kokkos::Profiling::pushRegion("FiniteElements::ComputeProjectionOn3DSkin");
55 using MemorySpace = typename ExecSpace::memory_space;
56 using Triangle = ArborX::Triangle<3, coordinates_t>;
57 using Point = ArborX::Point<3, coordinates_t>;
58
59 ExecSpace execSpace{};
60
61 auto sourcePointsPtr = transfer.sourcePoints;
62 auto sourceValuesPtr = transfer.sourceValues;
63 auto connValPtr = transfer.connValues;
64 auto connOffPtr = transfer.connOffsets;
65 auto CellTypesPtr = transfer.cellTypes;
66
67 auto targetPointsPtr = transfer.targetPoints;
68 auto targetValuesPtr = transfer.targetValues;
69 auto targetStatusPtr = transfer.targetStatus;
70 auto nbtargetPoints = targetPointsPtr.extent_int(0);
71
72 auto skinFacesPtr = transfer.skinFaces;
73 auto skinParentsPtr = transfer.skinParents;
74
75 auto nbTri = transfer.skinFaces.extent_int(0);
76 Kokkos::View<Triangle *, MemorySpace> skinFacesView(
77 Kokkos::view_alloc(execSpace, Kokkos::WithoutInitializing,
78 "Skin ArborX Triangles"),
79 nbTri);
80 Kokkos::parallel_for(
81 "Fill ArborX triangles",
82 Kokkos::RangePolicy<ExecSpace>(execSpace, 0, nbTri),
83 KOKKOS_LAMBDA(const int_t &i) {
84 for (int_t d = 0; d < 3; d++) {
85 skinFacesView(i).a[d] = sourcePointsPtr(skinFacesPtr(i, 0))[d];
86 skinFacesView(i).b[d] = sourcePointsPtr(skinFacesPtr(i, 1))[d];
87 skinFacesView(i).c[d] = sourcePointsPtr(skinFacesPtr(i, 2))[d];
88 }
89 });
90 ArborX::BoundingVolumeHierarchy skinBVH(
91 execSpace, ArborX::Experimental::attach_indices(skinFacesView));
92 Kokkos::Profiling::popRegion(); // Build triangle BoundingVolumeHierarchy
93
94 Kokkos::Profiling::pushRegion("Compute query of nearest skin");
95
96 Kokkos::View<int_t *, MemorySpace> values("values", 0);
97 Kokkos::View<offset_t *, MemorySpace> offsets("offsets", 0);
98
99 auto targetElement = Kokkos::View<int_t *, MemorySpace>(
100 Kokkos::view_alloc(execSpace, Kokkos::WithoutInitializing,
101 "Target to elem"),
102 nbtargetPoints);
103
104 Kokkos::Profiling::pushRegion("Point projection on triangle");
105 PointCloudNearest<MemorySpace, 3> pcn{targetPointsPtr};
106 if (extrapol) {
107 skinBVH.query(execSpace, pcn,
108 PointTriangleProjectionExtrapol<ExecSpace, 3>{
109 targetElement, targetStatusPtr, skinParentsPtr},
110 values, offsets);
111 } else {
112 skinBVH.query(
113 execSpace, pcn,
114 PointTriangleProjection<ExecSpace, 3>{targetPointsPtr, targetElement,
115 targetStatusPtr, skinParentsPtr},
116 values, offsets);
117 }
118 Kokkos::Profiling::popRegion(); // Point projection on triangle
119
120 auto statusOutside =
121 (extrapol) ? TransferStatus::EXTRAP : TransferStatus::CLAMP;
122
123 Kokkos::Profiling::pushRegion(
124 "Compute projection and interpolation coefficient");
125 Kokkos::parallel_for(
126 "Compute projection values ",
127 Kokkos::RangePolicy<ExecSpace>(execSpace, 0, values.extent_int(0)),
128 KOKKOS_LAMBDA(const int_t &i) {
129 auto predicate_index = values(i);
130 Kokkos::Array<fp_t, MaxNodesPerElt> warr;
131 Kokkos::Array<coordinates_t, MaxNodesPerElt * 3> Xarr;
132 Kokkos::Array<coordinates_t, 3> tpa;
133 Kokkos::View<coordinates_t *, ExecSpace> tp(tpa.data(), 3);
134 for (int_t j = 0; j < 3; j++) {
135 tp(j) = targetPointsPtr(predicate_index, j);
136 }
137 // predicate = target point data
138 // primitive intersected bounding box finite element data
139 const int_t curElem = targetElement(predicate_index);
140 const auto cellType = CellTypesPtr(curElem);
141 auto localConn = Kokkos::subview(
142 connValPtr,
143 Kokkos::make_pair(connOffPtr(curElem), connOffPtr(curElem + 1)));
144 int nbConnNodes = connOffPtr(curElem + 1) - connOffPtr(curElem);
145 Kokkos::View<fp_t *, ExecSpace> weights(warr.data(), nbConnNodes);
146 Kokkos::View<coordinates_t **, ExecSpace> Xcoor(Xarr.data(),
147 nbConnNodes, 3);
148 for (int_t j = 0; j < nbConnNodes; ++j) {
149 auto nodeId = localConn(j);
150 for (int_t k = 0; k < 3; ++k) {
151 Xcoor(j, k) = sourcePointsPtr(nodeId)[k];
152 }
153 }
154 auto isInside = ApplyNewtonOnElement<ExecSpace, 3>(cellType, Xcoor, tp,
155 weights, true);
156 for (int_t j = 0; j < nbConnNodes; ++j) {
157 auto nodeId = localConn(j);
158 targetValuesPtr(predicate_index) +=
159 weights[j] * sourceValuesPtr(nodeId);
160 }
161 targetStatusPtr(predicate_index) = statusOutside;
162 });
163 Kokkos::Profiling::popRegion(); // Compute projection and interpolation
164 // coefficient
165};
166
178template <typename ExecSpace>
179void ComputeProjectionOn2DSkin(Transfer<ExecSpace, 2> &transfer,
180 bool extrapol = false) {
181 Kokkos::Profiling::pushRegion("FiniteElements::ComputeProjectionOn2DSkin");
182 using MemorySpace = typename ExecSpace::memory_space;
183 using Point = ArborX::Point<2, coordinates_t>;
184
185 ExecSpace execSpace{};
186
187 auto sourcePointsPtr = transfer.sourcePoints;
188 auto sourceValuesPtr = transfer.sourceValues;
189 auto connValPtr = transfer.connValues;
190 auto connOffPtr = transfer.connOffsets;
191 auto CellTypesPtr = transfer.cellTypes;
192
193 auto targetPointsPtr = transfer.targetPoints;
194 auto targetValuesPtr = transfer.targetValues;
195 auto targetStatusPtr = transfer.targetStatus;
196 auto nbtargetPoints = targetPointsPtr.extent_int(0);
197
198 auto skinFacesPtr = transfer.skinFaces;
199 auto skinParentsPtr = transfer.skinParents;
200
201 auto statusOutside =
202 (extrapol) ? TransferStatus::EXTRAP : TransferStatus::CLAMP;
203
204 auto targetElement = Kokkos::View<int_t *, MemorySpace>(
205 Kokkos::view_alloc(execSpace, Kokkos::WithoutInitializing,
206 "Target to elem"),
207 nbtargetPoints);
208
209 Kokkos::Profiling::pushRegion(
210 "Compute projection and interpolation coefficient");
211 Kokkos::parallel_for(
212 "Retrieve closest bar",
213 Kokkos::RangePolicy<ExecSpace>(execSpace, 0,
214 targetPointsPtr.extent_int(0)),
215 KOKKOS_LAMBDA(const int_t &i) {
216 if (targetStatusPtr(i) == TransferStatus::OUTSIDE) {
217 Kokkos::Array<fp_t, MaxNodesPerElt> warr;
218 Kokkos::Array<coordinates_t, MaxNodesPerElt * 2> Xarr;
219 Kokkos::Array<coordinates_t, 2> tpa;
220 Kokkos::View<coordinates_t *, ExecSpace> tp(tpa.data(), 2);
221 const Point target = {targetPointsPtr(i, 0), targetPointsPtr(i, 1)};
222 fp_t mindistsqr = fp_consts::max();
223 int_t closestId = -1;
224 Point closest;
225 for (int_t s = 0; s < skinFacesPtr.extent_int(0); s++) {
226 auto p1 = sourcePointsPtr(skinFacesPtr(s, 0));
227 auto p2 = sourcePointsPtr(skinFacesPtr(s, 1));
228 auto cur = closest_point_to_edge(target, p1, p2);
229 coordinates_t dx = cur[0] - target[0];
230 coordinates_t dy = cur[1] - target[1];
231 coordinates_t curdistsqr = dx * dx + dy * dy;
232 if (curdistsqr < mindistsqr) {
233 closest = cur;
234 mindistsqr = curdistsqr;
235 closestId = s;
236 }
237 }
238 if (!extrapol) {
239 targetPointsPtr(i, 0) = closest[0];
240 targetPointsPtr(i, 1) = closest[1];
241 }
242 for (int_t j = 0; j < 2; j++) {
243 tp(j) = targetPointsPtr(i, j);
244 }
245 targetElement(i) = skinParentsPtr(closestId);
246 const int_t curElem = targetElement(i);
247 const auto cellType = CellTypesPtr(curElem);
248 auto localConn = Kokkos::subview(
249 connValPtr,
250 Kokkos::make_pair(connOffPtr(curElem), connOffPtr(curElem + 1)));
251 const int_t nbConnNodes =
252 connOffPtr(curElem + 1) - connOffPtr(curElem);
253 Kokkos::View<fp_t *, ExecSpace> weights(warr.data(), nbConnNodes);
254 Kokkos::View<coordinates_t **, ExecSpace> Xcoor(Xarr.data(),
255 nbConnNodes, 2);
256 for (int_t j = 0; j < nbConnNodes; ++j) {
257 const auto nodeId = localConn(j);
258 for (int_t k = 0; k < 2; ++k) {
259 Xcoor(j, k) = sourcePointsPtr(nodeId)[k];
260 }
261 }
262 auto isInside = FiniteElements::ApplyNewtonOnElement<ExecSpace, 2>(
263 cellType, Xcoor, tp, weights, true);
264 for (int_t j = 0; j < nbConnNodes; ++j) {
265 const auto nodeId = localConn(j);
266 targetValuesPtr(i) += weights[j] * sourceValuesPtr(nodeId);
267 }
268 targetStatusPtr(i) = statusOutside;
269 }
270 });
271 Kokkos::Profiling::popRegion(); // Compute projection and interpolation
272 // coefficient
273};
274
286template <typename ExecSpace>
287void ComputeProjectionOn1DSkin(Transfer<ExecSpace, 1> &transfer,
288 bool extrapol = false) {
289 Kokkos::Profiling::pushRegion("FiniteElement::ComputeProjectionOn1DSkin");
290 using MemorySpace = typename ExecSpace::memory_space;
291 using Point = ArborX::Point<1, coordinates_t>;
292
293 ExecSpace execSpace{};
294
295 auto sourcePointsPtr = transfer.sourcePoints;
296 auto sourceValuesPtr = transfer.sourceValues;
297 auto connValPtr = transfer.connValues;
298 auto connOffPtr = transfer.connOffsets;
299 auto CellTypesPtr = transfer.cellTypes;
300
301 auto targetPointsPtr = transfer.targetPoints;
302 auto targetValuesPtr = transfer.targetValues;
303 auto targetStatusPtr = transfer.targetStatus;
304 auto nbtargetPoints = targetPointsPtr.extent_int(0);
305
306 auto skinFacesPtr = transfer.skinFaces;
307 auto skinParentsPtr = transfer.skinParents;
308
309 auto statusOutside =
310 (extrapol) ? TransferStatus::EXTRAP : TransferStatus::CLAMP;
311
312 auto targetElement = Kokkos::View<int_t *, MemorySpace>(
313 Kokkos::view_alloc(execSpace, Kokkos::WithoutInitializing,
314 "Target to elem"),
315 nbtargetPoints);
316
317 Kokkos::Profiling::pushRegion(
318 "Compute projection and interpolation coefficient");
319 Kokkos::parallel_for(
320 "Retrieve closest point",
321 Kokkos::RangePolicy<ExecSpace>(execSpace, 0,
322 targetPointsPtr.extent_int(0)),
323 KOKKOS_LAMBDA(const int &i) {
324 if (targetStatusPtr(i) == TransferStatus::OUTSIDE) {
325 Kokkos::Array<fp_t, MaxNodesPerElt> warr;
326 Kokkos::Array<coordinates_t, MaxNodesPerElt> Xarr;
327 Kokkos::Array<coordinates_t, 1> tpa;
328 Kokkos::View<coordinates_t *, ExecSpace> tp(tpa.data(), 1);
329 const Point target = {targetPointsPtr(i, 0)};
330 fp_t mindist = fp_consts::max();
331 int_t closestId = -1;
332 Point closest;
333 for (int_t s = 0; s < skinFacesPtr.extent_int(0); s++) {
334 auto p = sourcePointsPtr(skinFacesPtr(s, 0));
335 const auto curdist = Kokkos::abs(p[0] - target[0]);
336 if (curdist < mindist) {
337 closest = p;
338 mindist = curdist;
339 closestId = s;
340 }
341 }
342 if (!extrapol) {
343 targetPointsPtr(i, 0) = closest[0];
344 }
345 for (int_t j = 0; j < 1; j++) {
346 tp(j) = targetPointsPtr(i, j);
347 }
348 targetElement(i) = skinParentsPtr(closestId);
349 const int_t curElem = targetElement(i);
350 const auto cellType = CellTypesPtr(curElem);
351 auto localConn = Kokkos::subview(
352 connValPtr,
353 Kokkos::make_pair(connOffPtr(curElem), connOffPtr(curElem + 1)));
354 int nbConnNodes = connOffPtr(curElem + 1) - connOffPtr(curElem);
355 Kokkos::View<fp_t *, ExecSpace> weights(warr.data(), nbConnNodes);
356 Kokkos::View<coordinates_t **, ExecSpace> Xcoor(Xarr.data(),
357 nbConnNodes, 1);
358 for (int_t j = 0; j < nbConnNodes; ++j) {
359 const auto nodeId = localConn(j);
360 for (int_t k = 0; k < 1; ++k) {
361 Xcoor(j, k) = sourcePointsPtr(nodeId)[k];
362 }
363 }
364 auto isInside = ApplyNewtonOnElement<ExecSpace, 1>(cellType, Xcoor,
365 tp, weights, true);
366 for (int_t j = 0; j < nbConnNodes; ++j) {
367 const auto nodeId = localConn(j);
368 targetValuesPtr(i) += weights[j] * sourceValuesPtr(nodeId);
369 }
370 targetStatusPtr(i) = statusOutside;
371 }
372 });
373 Kokkos::Profiling::popRegion(); // Compute projection and interpolation
374 // coefficient
375};
376
377} // namespace FiniteElements
378} // namespace PACMAN
ArborXCallbacks.hpp
ArborX predicate/callback helpers for PACMAN finite-elements kernels.
PACMAN::FiniteElements::closest_point_to_edge
KOKKOS_INLINE_FUNCTION ArborX::Point< 2, coordinates_t > closest_point_to_edge(const ArborX::Point< 2, coordinates_t > &p, const ArborX::Point< 2, coordinates_t > &a, const ArborX::Point< 2, coordinates_t > &b)
Compute the closest point on a 2D edge segment.
Definition SkinProjection.hpp:27
PACMAN::FiniteElements::ComputeProjectionOn3DSkin
void ComputeProjectionOn3DSkin(Transfer< ExecSpace, 3 > &transfer, bool extrapol=false)
Project outside target points on a 3D skin and evaluate FE values.
Definition SkinProjection.hpp:52
PACMAN::FiniteElements::ComputeProjectionOn2DSkin
void ComputeProjectionOn2DSkin(Transfer< ExecSpace, 2 > &transfer, bool extrapol=false)
Project outside target points on a 2D skin and evaluate FE values.
Definition SkinProjection.hpp:179
PACMAN::FiniteElements::ComputeProjectionOn1DSkin
void ComputeProjectionOn1DSkin(Transfer< ExecSpace, 1 > &transfer, bool extrapol=false)
Project outside target points on a 1D skin and evaluate FE values.
Definition SkinProjection.hpp:287
PACMAN::fp_consts::max
KOKKOS_INLINE_FUNCTION fp_t max(void)
Definition types.hpp:89
PACMAN::is_kokkos_view
Definition concepts.hpp:17