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 for (index_t component = 0;
159 component < sourceValuesPtr.extent(1); ++component) {
160 targetValuesPtr(predicate_index, component) +=
161 weights[j] * sourceValuesPtr(nodeId, component);
162 }
163 }
164 targetStatusPtr(predicate_index) = statusOutside;
165 });
166 Kokkos::Profiling::popRegion(); // Compute projection and interpolation
167 // coefficient
168};
169
181template <typename ExecSpace>
182void ComputeProjectionOn2DSkin(Transfer<ExecSpace, 2> &transfer,
183 bool extrapol = false) {
184 Kokkos::Profiling::pushRegion("FiniteElements::ComputeProjectionOn2DSkin");
185 using MemorySpace = typename ExecSpace::memory_space;
186 using Point = ArborX::Point<2, coordinates_t>;
187
188 ExecSpace execSpace{};
189
190 auto sourcePointsPtr = transfer.sourcePoints;
191 auto sourceValuesPtr = transfer.sourceValues;
192 auto connValPtr = transfer.connValues;
193 auto connOffPtr = transfer.connOffsets;
194 auto CellTypesPtr = transfer.cellTypes;
195
196 auto targetPointsPtr = transfer.targetPoints;
197 auto targetValuesPtr = transfer.targetValues;
198 auto targetStatusPtr = transfer.targetStatus;
199 auto nbtargetPoints = targetPointsPtr.extent_int(0);
200
201 auto skinFacesPtr = transfer.skinFaces;
202 auto skinParentsPtr = transfer.skinParents;
203
204 auto statusOutside =
205 (extrapol) ? TransferStatus::EXTRAP : TransferStatus::CLAMP;
206
207 auto targetElement = Kokkos::View<int_t *, MemorySpace>(
208 Kokkos::view_alloc(execSpace, Kokkos::WithoutInitializing,
209 "Target to elem"),
210 nbtargetPoints);
211
212 Kokkos::Profiling::pushRegion(
213 "Compute projection and interpolation coefficient");
214 Kokkos::parallel_for(
215 "Retrieve closest bar",
216 Kokkos::RangePolicy<ExecSpace>(execSpace, 0,
217 targetPointsPtr.extent_int(0)),
218 KOKKOS_LAMBDA(const int_t &i) {
219 if (targetStatusPtr(i) == TransferStatus::OUTSIDE) {
220 Kokkos::Array<fp_t, MaxNodesPerElt> warr;
221 Kokkos::Array<coordinates_t, MaxNodesPerElt * 2> Xarr;
222 Kokkos::Array<coordinates_t, 2> tpa;
223 Kokkos::View<coordinates_t *, ExecSpace> tp(tpa.data(), 2);
224 const Point target = {targetPointsPtr(i, 0), targetPointsPtr(i, 1)};
225 fp_t mindistsqr = fp_consts::max();
226 int_t closestId = -1;
227 Point closest;
228 for (int_t s = 0; s < skinFacesPtr.extent_int(0); s++) {
229 auto p1 = sourcePointsPtr(skinFacesPtr(s, 0));
230 auto p2 = sourcePointsPtr(skinFacesPtr(s, 1));
231 auto cur = closest_point_to_edge(target, p1, p2);
232 coordinates_t dx = cur[0] - target[0];
233 coordinates_t dy = cur[1] - target[1];
234 coordinates_t curdistsqr = dx * dx + dy * dy;
235 if (curdistsqr < mindistsqr) {
236 closest = cur;
237 mindistsqr = curdistsqr;
238 closestId = s;
239 }
240 }
241 if (!extrapol) {
242 targetPointsPtr(i, 0) = closest[0];
243 targetPointsPtr(i, 1) = closest[1];
244 }
245 for (int_t j = 0; j < 2; j++) {
246 tp(j) = targetPointsPtr(i, j);
247 }
248 targetElement(i) = skinParentsPtr(closestId);
249 const int_t curElem = targetElement(i);
250 const auto cellType = CellTypesPtr(curElem);
251 auto localConn = Kokkos::subview(
252 connValPtr,
253 Kokkos::make_pair(connOffPtr(curElem), connOffPtr(curElem + 1)));
254 const int_t nbConnNodes =
255 connOffPtr(curElem + 1) - connOffPtr(curElem);
256 Kokkos::View<fp_t *, ExecSpace> weights(warr.data(), nbConnNodes);
257 Kokkos::View<coordinates_t **, ExecSpace> Xcoor(Xarr.data(),
258 nbConnNodes, 2);
259 for (int_t j = 0; j < nbConnNodes; ++j) {
260 const auto nodeId = localConn(j);
261 for (int_t k = 0; k < 2; ++k) {
262 Xcoor(j, k) = sourcePointsPtr(nodeId)[k];
263 }
264 }
265 auto isInside = FiniteElements::ApplyNewtonOnElement<ExecSpace, 2>(
266 cellType, Xcoor, tp, weights, true);
267 for (int_t j = 0; j < nbConnNodes; ++j) {
268 const auto nodeId = localConn(j);
269 for (index_t component = 0;
270 component < sourceValuesPtr.extent(1); ++component) {
271 targetValuesPtr(i, component) +=
272 weights[j] * sourceValuesPtr(nodeId, component);
273 }
274 }
275 targetStatusPtr(i) = statusOutside;
276 }
277 });
278 Kokkos::Profiling::popRegion(); // Compute projection and interpolation
279 // coefficient
280};
281
293template <typename ExecSpace>
294void ComputeProjectionOn1DSkin(Transfer<ExecSpace, 1> &transfer,
295 bool extrapol = false) {
296 Kokkos::Profiling::pushRegion("FiniteElement::ComputeProjectionOn1DSkin");
297 using MemorySpace = typename ExecSpace::memory_space;
298 using Point = ArborX::Point<1, coordinates_t>;
299
300 ExecSpace execSpace{};
301
302 auto sourcePointsPtr = transfer.sourcePoints;
303 auto sourceValuesPtr = transfer.sourceValues;
304 auto connValPtr = transfer.connValues;
305 auto connOffPtr = transfer.connOffsets;
306 auto CellTypesPtr = transfer.cellTypes;
307
308 auto targetPointsPtr = transfer.targetPoints;
309 auto targetValuesPtr = transfer.targetValues;
310 auto targetStatusPtr = transfer.targetStatus;
311 auto nbtargetPoints = targetPointsPtr.extent_int(0);
312
313 auto skinFacesPtr = transfer.skinFaces;
314 auto skinParentsPtr = transfer.skinParents;
315
316 auto statusOutside =
317 (extrapol) ? TransferStatus::EXTRAP : TransferStatus::CLAMP;
318
319 auto targetElement = Kokkos::View<int_t *, MemorySpace>(
320 Kokkos::view_alloc(execSpace, Kokkos::WithoutInitializing,
321 "Target to elem"),
322 nbtargetPoints);
323
324 Kokkos::Profiling::pushRegion(
325 "Compute projection and interpolation coefficient");
326 Kokkos::parallel_for(
327 "Retrieve closest point",
328 Kokkos::RangePolicy<ExecSpace>(execSpace, 0,
329 targetPointsPtr.extent_int(0)),
330 KOKKOS_LAMBDA(const int &i) {
331 if (targetStatusPtr(i) == TransferStatus::OUTSIDE) {
332 Kokkos::Array<fp_t, MaxNodesPerElt> warr;
333 Kokkos::Array<coordinates_t, MaxNodesPerElt> Xarr;
334 Kokkos::Array<coordinates_t, 1> tpa;
335 Kokkos::View<coordinates_t *, ExecSpace> tp(tpa.data(), 1);
336 const Point target = {targetPointsPtr(i, 0)};
337 fp_t mindist = fp_consts::max();
338 int_t closestId = -1;
339 Point closest;
340 for (int_t s = 0; s < skinFacesPtr.extent_int(0); s++) {
341 auto p = sourcePointsPtr(skinFacesPtr(s, 0));
342 const auto curdist = Kokkos::abs(p[0] - target[0]);
343 if (curdist < mindist) {
344 closest = p;
345 mindist = curdist;
346 closestId = s;
347 }
348 }
349 if (!extrapol) {
350 targetPointsPtr(i, 0) = closest[0];
351 }
352 for (int_t j = 0; j < 1; j++) {
353 tp(j) = targetPointsPtr(i, j);
354 }
355 targetElement(i) = skinParentsPtr(closestId);
356 const int_t curElem = targetElement(i);
357 const auto cellType = CellTypesPtr(curElem);
358 auto localConn = Kokkos::subview(
359 connValPtr,
360 Kokkos::make_pair(connOffPtr(curElem), connOffPtr(curElem + 1)));
361 int nbConnNodes = connOffPtr(curElem + 1) - connOffPtr(curElem);
362 Kokkos::View<fp_t *, ExecSpace> weights(warr.data(), nbConnNodes);
363 Kokkos::View<coordinates_t **, ExecSpace> Xcoor(Xarr.data(),
364 nbConnNodes, 1);
365 for (int_t j = 0; j < nbConnNodes; ++j) {
366 const auto nodeId = localConn(j);
367 for (int_t k = 0; k < 1; ++k) {
368 Xcoor(j, k) = sourcePointsPtr(nodeId)[k];
369 }
370 }
371 auto isInside = ApplyNewtonOnElement<ExecSpace, 1>(cellType, Xcoor,
372 tp, weights, true);
373 for (int_t j = 0; j < nbConnNodes; ++j) {
374 const auto nodeId = localConn(j);
375 for (index_t component = 0;
376 component < sourceValuesPtr.extent(1); ++component) {
377 targetValuesPtr(i, component) +=
378 weights[j] * sourceValuesPtr(nodeId, component);
379 }
380 }
381 targetStatusPtr(i) = statusOutside;
382 }
383 });
384 Kokkos::Profiling::popRegion(); // Compute projection and interpolation
385 // coefficient
386};
387
388} // namespace FiniteElements
389} // 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:182
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:294
PACMAN::FiniteElements::ApplyNewtonOnElement
KOKKOS_FUNCTION bool ApplyNewtonOnElement(const cell_t type, const Kokkos::View< coordinates_t **, ExecSpace > Xcoor, const Kokkos::View< coordinates_t[Dim], ExecSpace > targetPoint, Kokkos::View< fp_t *, ExecSpace > weights, const bool forceEvaluation=false)
Dispatch Newton-based FE inversion to the correct element type.
Definition NewtonKokkos.hpp:190
PACMAN::fp_consts::max
KOKKOS_INLINE_FUNCTION fp_t max(void)
Definition types.hpp:89
PACMAN::fp_t
double fp_t
Definition types.hpp:15
PACMAN::coordinates_t
double coordinates_t
Definition types.hpp:16
PACMAN::index_t
uint32_t index_t
Definition types.hpp:20
PACMAN::int_t
int32_t int_t
Definition types.hpp:18
PACMAN::FiniteElements::PointCloudNearest
Predicate wrapper for nearest-neighbor queries over point clouds.
Definition ArborXCallbacks.hpp:169
PACMAN::FiniteElements::PointTriangleProjectionExtrapol
ArborX callback variant for extrapolation mode.
Definition ArborXCallbacks.hpp:136
PACMAN::FiniteElements::PointTriangleProjection
Project outside target points onto nearest 3D skin triangle.
Definition ArborXCallbacks.hpp:91