Compadre 1.7.2
Loading...
Searching...
No Matches
Compadre_Targets.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Compadre: COMpatible PArticle Discretization and REmap Toolkit
4//
5// Copyright 2018 NTESS and the Compadre contributors.
6// SPDX-License-Identifier: BSD-2-Clause
7// *****************************************************************************
8// @HEADER
9#ifndef _COMPADRE_TARGETS_HPP_
10#define _COMPADRE_TARGETS_HPP_
11
12#include "Compadre_GMLS.hpp"
13#include "Compadre_Manifold_Functions.hpp"
14
15namespace Compadre {
16
17/*! \brief Evaluates a polynomial basis with a target functional applied to each member of the basis
18 \param data [in] - GMLSBasisData struct
19 \param teamMember [in] - Kokkos::TeamPolicy member type (created by parallel_for)
20 \param delta [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large is the _basis_multipler*the dimension of the polynomial basis.
21 \param thread_workspace [in/out] - scratch space that is allocated so that each team has its own copy. Must be at least as large is the _poly_order*_global_dimensions.
22 \param P_target_row [out] - 1D Kokkos View where the evaluation of the polynomial basis is stored
23*/
24template <typename TargetData>
25KOKKOS_INLINE_FUNCTION
26void computeTargetFunctionals(const TargetData& data, const member_type& teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, device_unmanaged_matrix_right_type P_target_row) {
27
28 // check if VectorOfScalarClonesTaylorPolynomial is used with a scalar sampling functional other than PointSample
29 if (data._dimensions > 1) {
30 compadre_kernel_assert_debug(
31 (data._reconstruction_space!=ReconstructionSpace::VectorOfScalarClonesTaylorPolynomial
32 || data._data_sampling_multiplier!=0
33 || (data._reconstruction_space==ReconstructionSpace::VectorOfScalarClonesTaylorPolynomial
34 && data._polynomial_sampling_functional==PointSample))
35 && "data._reconstruction_space(VectorOfScalar clones incompatible with scalar output sampling functional which is not a PointSample");
36 }
37
38 // determine if additional evaluation sites are requested by user and handled by target operations
39 bool additional_evaluation_sites_need_handled =
40 (data._additional_pc._source_coordinates.extent(0) > 0) ? true : false; // additional evaluation sites are specified
41
42 const int target_index = data._initial_index_for_batch + teamMember.league_rank();
43
44 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, P_target_row.extent(0)), [&] (const int j) {
45 Kokkos::parallel_for(Kokkos::ThreadVectorRange(teamMember, P_target_row.extent(1)),
46 [&] (const int& k) {
47 P_target_row(j,k) = 0;
48 });
49 });
50 for (int j = 0; j < delta.extent_int(0); ++j) {
51 delta(j) = 0;
52 }
53 for (int j = 0; j < thread_workspace.extent_int(0); ++j) {
54 thread_workspace(j) = 0;
55 }
56 teamMember.team_barrier();
57
58 // not const b.c. of gcc 7.2 issue
59 int target_NP = GMLS::getNP(data._poly_order, data._dimensions, data._reconstruction_space);
60 const int num_evaluation_sites = data._additional_pc._nla.getNumberOfNeighborsDevice(target_index) + 1;
61
62 for (size_t i=0; i<data._operations.size(); ++i) {
63
64 bool additional_evaluation_sites_handled = false; // target operations that can handle these sites should flip this flag
65
66 bool operation_handled = true;
67
68 // USER defined targets should be added to this file
69 // if the USER defined targets don't catch this operation, then operation_handled will be false
70 #include "Compadre_USER_StandardTargetFunctionals.hpp"
71
72 // if the user didn't handle the operation, we pass it along to the toolkit's targets
73 if (!operation_handled) {
74
75 if (data._reconstruction_space == ReconstructionSpace::ScalarTaylorPolynomial) {
76
77 /*
78 * Beginning of ScalarTaylorPolynomial basis
79 */
80
81 if (data._operations(i) == TargetOperation::ScalarPointEvaluation || (data._operations(i) == TargetOperation::VectorPointEvaluation && data._dimensions == 1) /* vector is a scalar in 1D */) {
82 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
83 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, data._dimensions, data._poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
84 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
85 for (int k=0; k<target_NP; ++k) {
86 P_target_row(offset, k) = delta(k);
87 }
88 });
89 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
90 } else if (data._operations(i) == TargetOperation::LaplacianOfScalarPointEvaluation) {
91 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
92 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
93 switch (data._dimensions) {
94 case 3:
95 P_target_row(offset, 4) = std::pow(data._epsilons(target_index), -2);
96 P_target_row(offset, 6) = std::pow(data._epsilons(target_index), -2);
97 P_target_row(offset, 9) = std::pow(data._epsilons(target_index), -2);
98 break;
99 case 2:
100 P_target_row(offset, 3) = std::pow(data._epsilons(target_index), -2);
101 P_target_row(offset, 5) = std::pow(data._epsilons(target_index), -2);
102 break;
103 default:
104 P_target_row(offset, 2) = std::pow(data._epsilons(target_index), -2);
105 }
106 });
107 } else if (data._operations(i) == TargetOperation::GradientOfScalarPointEvaluation) {
108 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
109 for (int d=0; d<data._dimensions; ++d) {
110 int offset = data._d_ss.getTargetOffsetIndex(i, 0, d, j);
111 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
112 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, d /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
113 }
114 });
115 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
116 } else if (data._operations(i) == TargetOperation::PartialXOfScalarPointEvaluation) {
117 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
118 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
119 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
120 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 0 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
121 });
122 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
123 } else if (data._operations(i) == TargetOperation::PartialYOfScalarPointEvaluation) {
124 compadre_kernel_assert_release(data._dimensions>1 && "PartialYOfScalarPointEvaluation requested for dim < 2");
125 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
126 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
127 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
128 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 1 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
129 });
130 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
131 } else if (data._operations(i) == TargetOperation::PartialZOfScalarPointEvaluation) {
132 compadre_kernel_assert_release(data._dimensions>2 && "PartialZOfScalarPointEvaluation requested for dim < 3");
133 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
134 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
135 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
136 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 2 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
137 });
138 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
139 }
140 // staggered gradient w/ edge integrals known analytically, using a basis
141 // of potentials
142 else if (data._operations(i) == TargetOperation::DivergenceOfVectorPointEvaluation
143 && data._polynomial_sampling_functional == StaggeredEdgeAnalyticGradientIntegralSample) {
144 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
145 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
146 switch (data._dimensions) {
147 case 3:
148 P_target_row(offset, 4) = std::pow(data._epsilons(target_index), -2);
149 P_target_row(offset, 6) = std::pow(data._epsilons(target_index), -2);
150 P_target_row(offset, 9) = std::pow(data._epsilons(target_index), -2);
151 break;
152 case 2:
153 P_target_row(offset, 3) = std::pow(data._epsilons(target_index), -2);
154 P_target_row(offset, 5) = std::pow(data._epsilons(target_index), -2);
155 break;
156 default:
157 P_target_row(offset, 2) = std::pow(data._epsilons(target_index), -2);
158 }
159 });
160 } else if (data._operations(i) == TargetOperation::CellAverageEvaluation ||
161 data._operations(i) == TargetOperation::CellIntegralEvaluation) {
162 compadre_kernel_assert_debug(data._local_dimensions==2 &&
163 "CellAverageEvaluation and CellIntegralEvaluation only support 2d or 3d with 2d manifold");
164 const double factorial[15] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
165 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
166
167 // Calculate basis matrix for NON MANIFOLD problems
168 double cutoff_p = data._epsilons(target_index);
169 int alphax, alphay;
170 double alphaf;
171
172 double triangle_coords[3*3]; //data._global_dimensions*3
173 device_unmanaged_matrix_right_type triangle_coords_matrix(triangle_coords, data._global_dimensions, 3);
174
175 for (int j=0; j<data._global_dimensions; ++j) {
176 // midpoint
177 triangle_coords_matrix(j, 0) = data._pc.getTargetCoordinate(target_index, j);
178 }
179
180 size_t num_vertices = (data._target_extra_data(target_index, data._target_extra_data.extent(1)-1)
181 != data._target_extra_data(target_index, data._target_extra_data.extent(1)-1))
182 ? (data._target_extra_data.extent(1) / data._global_dimensions) - 1 :
183 (data._target_extra_data.extent(1) / data._global_dimensions);
184
185 XYZ relative_coord;
186 // loop over each two vertices
187 double entire_cell_area = 0.0;
188 for (size_t v=0; v<num_vertices; ++v) {
189 int v1 = v;
190 int v2 = (v+1) % num_vertices;
191
192 for (int j=0; j<data._global_dimensions; ++j) {
193 triangle_coords_matrix(j,1) = data._target_extra_data(target_index, v1*data._global_dimensions+j)
194 - triangle_coords_matrix(j,0);
195 triangle_coords_matrix(j,2) = data._target_extra_data(target_index, v2*data._global_dimensions+j)
196 - triangle_coords_matrix(j,0);
197 }
198 // triangle_coords now has:
199 // (midpoint_x, midpoint_y, midpoint_z,
200 // v1_x-midpoint_x, v1_y-midpoint_y, v1_z-midpoint_z,
201 // v2_x-midpoint_x, v2_y-midpoint_y, v2_z-midpoint_z);
202 auto T=triangle_coords_matrix;
203 for (int quadrature = 0; quadrature<data._qm.getNumberOfQuadraturePoints(); ++quadrature) {
204 double transformed_qp[2] = {0,0};
205 for (int j=0; j<data._global_dimensions; ++j) {
206 for (int k=1; k<3; ++k) {
207 transformed_qp[j] += T(j,k)*data._qm.getSite(quadrature, k-1);
208 }
209 transformed_qp[j] += T(j,0);
210 }
211 // half the norm of the cross-product is the area of the triangle
212 // so scaling is area / reference area (0.5) = the norm of the cross-product
213 double G_determinant = getAreaFromVectors(teamMember,
214 Kokkos::subview(T, Kokkos::ALL(), 1), Kokkos::subview(T, Kokkos::ALL(), 2));
215
216 for (int j=0; j<data._global_dimensions; ++j) {
217 relative_coord[j] = transformed_qp[j] - data._pc.getTargetCoordinate(target_index, j); // shift quadrature point by target site
218 }
219
220 int k = 0;
221 for (int n = 0; n <= data._poly_order; n++){
222 for (alphay = 0; alphay <= n; alphay++){
223 alphax = n - alphay;
224 alphaf = factorial[alphax]*factorial[alphay];
225 double val_to_sum = G_determinant * ( data._qm.getWeight(quadrature)
226 * std::pow(relative_coord.x/cutoff_p,alphax)
227 * std::pow(relative_coord.y/cutoff_p,alphay)/alphaf);
228 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
229 if (quadrature==0 && v==0) P_target_row(offset, k) = val_to_sum;
230 else P_target_row(offset, k) += val_to_sum;
231 });
232 k++;
233 }
234 }
235 entire_cell_area += G_determinant * data._qm.getWeight(quadrature);
236 }
237 }
238 if (data._operations(i) == TargetOperation::CellAverageEvaluation) {
239 int k = 0;
240 for (int n = 0; n <= data._poly_order; n++){
241 for (alphay = 0; alphay <= n; alphay++){
242 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
243 P_target_row(offset, k) /= entire_cell_area;
244 });
245 k++;
246 }
247 }
248 }
249 } else {
250 compadre_kernel_assert_release((false) && "Functionality not yet available.");
251 }
252
253 /*
254 * End of ScalarTaylorPolynomial basis
255 */
256
257 } else if (data._reconstruction_space == ReconstructionSpace::VectorTaylorPolynomial) {
258
259 /*
260 * Beginning of VectorTaylorPolynomial basis
261 */
262
263 if (data._operations(i) == TargetOperation::ScalarPointEvaluation || (data._operations(i) == TargetOperation::VectorPointEvaluation && data._dimensions == 1) /* vector is a scalar in 1D */) {
264 // copied from ScalarTaylorPolynomial
265 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
266 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, data._dimensions, data._poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
267 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
268 for (int k=0; k<target_NP; ++k) {
269 P_target_row(offset, k) = delta(k);
270 }
271 });
272 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
273 } else if (data._operations(i) == TargetOperation::VectorPointEvaluation) {
274 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int e) {
275 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, data._dimensions, data._poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, e);
276 for (int m=0; m<data._sampling_multiplier; ++m) {
277 int output_components = data._basis_multiplier;
278 for (int c=0; c<output_components; ++c) {
279 int offset = data._d_ss.getTargetOffsetIndex(i, m /*in*/, c /*out*/, e/*additional*/);
280 // for the case where data._sampling_multiplier is > 1,
281 // this approach relies on c*target_NP being equivalent to P_target_row(offset, j) where offset is
282 // data._d_ss.getTargetOffsetIndex(i, m /*in*/, c /*out*/, e/*additional*/)*data._basis_multiplier*target_NP;
283 for (int j=0; j<target_NP; ++j) {
284 P_target_row(offset, c*target_NP + j) = delta(j);
285 }
286 }
287 }
288 });
289 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
290 } else if ( (data._operations(i) == TargetOperation::GradientOfScalarPointEvaluation) && (data._polynomial_sampling_functional == StaggeredEdgeIntegralSample) ) {
291 // when using staggered edge integral sample with vector basis, the gradient of scalar point evaluation
292 // is just the vector point evaluation
293 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int e) {
294 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, data._dimensions, data._poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, e);
295 for (int m=0; m<data._sampling_multiplier; ++m) {
296 int output_components = data._basis_multiplier;
297 for (int c=0; c<output_components; ++c) {
298 int offset = data._d_ss.getTargetOffsetIndex(i, m /*in*/, c /*out*/, e/*additional*/);
299 // for the case where data._sampling_multiplier is > 1,
300 // this approach relies on c*target_NP being equivalent to P_target_row(offset, j) where offset is
301 // data._d_ss.getTargetOffsetIndex(i, m /*in*/, c /*out*/, e/*additional*/)*data._basis_multiplier*target_NP;
302 for (int j=0; j<target_NP; ++j) {
303 P_target_row(offset, c*target_NP + j) = delta(j);
304 }
305 }
306 }
307 });
308 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
309 } else if (data._operations(i) == TargetOperation::GradientOfScalarPointEvaluation) {
310 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int e) {
311 int output_components = data._basis_multiplier;
312 for (int m2=0; m2<output_components; ++m2) { // output components 2
313 int offset = data._d_ss.getTargetOffsetIndex(i, 0 /*in*/, m2 /*out*/, e);
314 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, m2 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, e);
315 for (int j=0; j<target_NP; ++j) {
316 P_target_row(offset, j) = delta(j);
317 }
318 }
319 });
320 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
321 } else if (data._operations(i) == TargetOperation::GradientOfVectorPointEvaluation) {
322 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int e) {
323 for (int m0=0; m0<data._sampling_multiplier; ++m0) { // input components
324 int output_components = data._basis_multiplier;
325 for (int m1=0; m1<output_components; ++m1) { // output components 1
326 for (int m2=0; m2<output_components; ++m2) { // output components 2
327 int offset = data._d_ss.getTargetOffsetIndex(i, m0 /*in*/, m1*output_components+m2 /*out*/, e);
328 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, m2 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, e);
329 for (int j=0; j<target_NP; ++j) {
330 P_target_row(offset, m1*target_NP + j) = delta(j);
331 }
332 }
333 }
334 }
335 });
336 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
337 } else if (data._operations(i) == TargetOperation::DivergenceOfVectorPointEvaluation) {
338 if (data._polynomial_sampling_functional == StaggeredEdgeIntegralSample) {
339 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
340 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);;
341 switch (data._dimensions) {
342 case 3:
343 P_target_row(offset, 1) = std::pow(data._epsilons(target_index), -1);
344 P_target_row(offset, target_NP + 2) = std::pow(data._epsilons(target_index), -1);
345 P_target_row(offset, 2*target_NP + 3) = std::pow(data._epsilons(target_index), -1);
346 break;
347 case 2:
348 P_target_row(offset, 1) = std::pow(data._epsilons(target_index), -1);
349 P_target_row(offset, target_NP + 2) = std::pow(data._epsilons(target_index), -1);
350 break;
351 default:
352 P_target_row(offset, 1) = std::pow(data._epsilons(target_index), -1);
353 }
354 });
355 } else {
356 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int e) {
357 for (int m=0; m<data._sampling_multiplier; ++m) {
358 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, m /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, e);
359 int offset = data._d_ss.getTargetOffsetIndex(i, m /*in*/, 0 /*out*/, e/*additional*/);
360 for (int j=0; j<target_NP; ++j) {
361 P_target_row(offset, m*target_NP + j) = delta(j);
362 }
363 }
364 });
365 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
366 }
367 } else if (data._operations(i) == TargetOperation::CurlOfVectorPointEvaluation) {
368 if (data._dimensions==3) {
369 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
370 // output component 0
371 // u_{2,y} - u_{1,z}
372 {
373 int offset = data._d_ss.getTargetOffsetIndex(i, 0 /*in*/, 0 /*out*/, 0/*additional*/);
374 // role of input 0 on component 0 of curl
375 // (no contribution)
376
377 offset = data._d_ss.getTargetOffsetIndex(i, 1 /*in*/, 0 /*out*/, 0/*additional*/);
378 // role of input 1 on component 0 of curl
379 // -u_{1,z}
380 P_target_row(offset, target_NP + 3) = -std::pow(data._epsilons(target_index), -1);
381
382 offset = data._d_ss.getTargetOffsetIndex(i, 2 /*in*/, 0 /*out*/, 0/*additional*/);
383 // role of input 2 on component 0 of curl
384 // u_{2,y}
385 P_target_row(offset, 2*target_NP + 2) = std::pow(data._epsilons(target_index), -1);
386 }
387
388 // output component 1
389 // -u_{2,x} + u_{0,z}
390 {
391 int offset = data._d_ss.getTargetOffsetIndex(i, 0 /*in*/, 1 /*out*/, 0/*additional*/);
392 // role of input 0 on component 1 of curl
393 // u_{0,z}
394 P_target_row(offset, 3) = std::pow(data._epsilons(target_index), -1);
395
396 // offset = data._d_ss.getTargetOffsetIndex(i, 1 /*in*/, 1 /*out*/, 0/*additional*/);
397 // role of input 1 on component 1 of curl
398 // (no contribution)
399
400 offset = data._d_ss.getTargetOffsetIndex(i, 2 /*in*/, 1 /*out*/, 0/*additional*/);
401 // role of input 2 on component 1 of curl
402 // -u_{2,x}
403 P_target_row(offset, 2*target_NP + 1) = -std::pow(data._epsilons(target_index), -1);
404 }
405
406 // output component 2
407 // u_{1,x} - u_{0,y}
408 {
409 int offset = data._d_ss.getTargetOffsetIndex(i, 0 /*in*/, 2 /*out*/, 0/*additional*/);
410 // role of input 0 on component 1 of curl
411 // -u_{0,y}
412 P_target_row(offset, 2) = -std::pow(data._epsilons(target_index), -1);
413
414 offset = data._d_ss.getTargetOffsetIndex(i, 1 /*in*/, 2 /*out*/, 0/*additional*/);
415 // role of input 1 on component 1 of curl
416 // u_{1,x}
417 P_target_row(offset, target_NP + 1) = std::pow(data._epsilons(target_index), -1);
418
419 // offset = data._d_ss.getTargetOffsetIndex(i, 2 /*in*/, 2 /*out*/, 0/*additional*/);
420 // role of input 2 on component 1 of curl
421 // (no contribution)
422 }
423 });
424 } else if (data._dimensions==2) {
425 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
426 // output component 0
427 // u_{1,y}
428 {
429 int offset = data._d_ss.getTargetOffsetIndex(i, 0 /*in*/, 0 /*out*/, 0/*additional*/);
430 // role of input 0 on component 0 of curl
431 // (no contribution)
432
433 offset = data._d_ss.getTargetOffsetIndex(i, 1 /*in*/, 0 /*out*/, 0/*additional*/);
434 // role of input 1 on component 0 of curl
435 // -u_{1,z}
436 P_target_row(offset, target_NP + 2) = std::pow(data._epsilons(target_index), -1);
437 }
438
439 // output component 1
440 // -u_{0,x}
441 {
442 int offset = data._d_ss.getTargetOffsetIndex(i, 0 /*in*/, 1 /*out*/, 0/*additional*/);
443 // role of input 0 on component 1 of curl
444 // u_{0,z}
445 P_target_row(offset, 1) = -std::pow(data._epsilons(target_index), -1);
446
447 //offset = data._d_ss.getTargetOffsetIndex(i, 1 /*in*/, 1 /*out*/, 0/*additional*/);
448 // role of input 1 on component 1 of curl
449 // (no contribution)
450 }
451 });
452 }
453 } else {
454 compadre_kernel_assert_release((false) && "Functionality not yet available.");
455 }
456
457 /*
458 * End of VectorTaylorPolynomial basis
459 */
460
461 } else if (data._reconstruction_space == ReconstructionSpace::VectorOfScalarClonesTaylorPolynomial) {
462
463 /*
464 * Beginning of VectorOfScalarClonesTaylorPolynomial basis
465 */
466
467 if (data._operations(i) == TargetOperation::ScalarPointEvaluation || (data._operations(i) == TargetOperation::VectorPointEvaluation && data._dimensions == 1) /* vector is a scalar in 1D */) {
468 // copied from ScalarTaylorPolynomial
469 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
470 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, data._dimensions, data._poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
471 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
472 for (int k=0; k<target_NP; ++k) {
473 P_target_row(offset, k) = delta(k);
474 }
475 });
476 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
477 } else if (data._operations(i) == TargetOperation::VectorPointEvaluation) {
478 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int e) {
479 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, data._dimensions, data._poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, e);
480 for (int m=0; m<data._sampling_multiplier; ++m) {
481 for (int c=0; c<data._data_sampling_multiplier; ++c) {
482 int offset = data._d_ss.getTargetOffsetIndex(i, c /*in*/, c /*out*/, e/*additional*/);
483 for (int j=0; j<target_NP; ++j) {
484 P_target_row(offset, j) = delta(j);
485 }
486 }
487 }
488 });
489 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
490 } else if (data._operations(i) == TargetOperation::LaplacianOfScalarPointEvaluation) {
491 // copied from ScalarTaylorPolynomial
492 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
493 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
494 switch (data._dimensions) {
495 case 3:
496 P_target_row(offset, 4) = std::pow(data._epsilons(target_index), -2);
497 P_target_row(offset, 6) = std::pow(data._epsilons(target_index), -2);
498 P_target_row(offset, 9) = std::pow(data._epsilons(target_index), -2);
499 break;
500 case 2:
501 P_target_row(offset, 3) = std::pow(data._epsilons(target_index), -2);
502 P_target_row(offset, 5) = std::pow(data._epsilons(target_index), -2);
503 break;
504 default:
505 P_target_row(offset, 2) = std::pow(data._epsilons(target_index), -2);
506 }
507 });
508 } else if (data._operations(i) == TargetOperation::GradientOfScalarPointEvaluation) {
509 // copied from ScalarTaylorPolynomial
510 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
511 for (int d=0; d<data._dimensions; ++d) {
512 int offset = data._d_ss.getTargetOffsetIndex(i, 0, d, j);
513 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
514 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, d /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
515 }
516 });
517 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
518 } else if (data._operations(i) == TargetOperation::GradientOfVectorPointEvaluation) {
519 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
520 for (int m0=0; m0<data._dimensions; ++m0) { // input components
521 for (int m1=0; m1<data._dimensions; ++m1) { // output components 1
522 for (int m2=0; m2<data._dimensions; ++m2) { // output componets 2
523 int offset = data._d_ss.getTargetOffsetIndex(i, m0 /*in*/, m1*data._dimensions+m2 /*out*/, j);
524 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
525 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, m2 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
526 }
527 }
528 }
529 });
530 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
531 } else if (data._operations(i) == TargetOperation::PartialXOfScalarPointEvaluation) {
532 // copied from ScalarTaylorPolynomial
533 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
534 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
535 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
536 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 0 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
537 });
538 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
539 } else if (data._operations(i) == TargetOperation::PartialYOfScalarPointEvaluation) {
540 // copied from ScalarTaylorPolynomial
541 if (data._dimensions>1) {
542 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
543 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
544 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
545 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 1 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
546 });
547 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
548 }
549 } else if (data._operations(i) == TargetOperation::PartialZOfScalarPointEvaluation) {
550 // copied from ScalarTaylorPolynomial
551 if (data._dimensions>2) {
552 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
553 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
554 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
555 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 2 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
556 });
557 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
558 }
559 } else if (data._operations(i) == TargetOperation::DivergenceOfVectorPointEvaluation) {
560 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
561 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
562 for (int j=0; j<target_NP; ++j) {
563 P_target_row(offset, j) = 0;
564 }
565
566 P_target_row(offset, 1) = std::pow(data._epsilons(target_index), -1);
567
568 if (data._dimensions>1) {
569 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, 0);
570 for (int j=0; j<target_NP; ++j) {
571 P_target_row(offset, j) = 0;
572 }
573 P_target_row(offset, 2) = std::pow(data._epsilons(target_index), -1);
574 }
575
576 if (data._dimensions>2) {
577 offset = data._d_ss.getTargetOffsetIndex(i, 2, 0, 0);
578 for (int j=0; j<target_NP; ++j) {
579 P_target_row(offset, j) = 0;
580 }
581 P_target_row(offset, 3) = std::pow(data._epsilons(target_index), -1);
582 }
583 });
584 } else if (data._operations(i) == TargetOperation::CurlOfVectorPointEvaluation) {
585 // comments based on taking curl of vector [u_{0},u_{1},u_{2}]^T
586 // with as e.g., u_{1,z} being the partial derivative with respect to z of
587 // u_{1}
588 if (data._dimensions==3) {
589 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
590 // output component 0
591 // u_{2,y} - u_{1,z}
592 {
593 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
594 for (int j=0; j<target_NP; ++j) {
595 P_target_row(offset, j) = 0;
596 }
597 // role of input 0 on component 0 of curl
598 // (no contribution)
599
600 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, 0);
601 for (int j=0; j<target_NP; ++j) {
602 P_target_row(offset, j) = 0;
603 }
604 // role of input 1 on component 0 of curl
605 // -u_{1,z}
606 P_target_row(offset, 3) = -std::pow(data._epsilons(target_index), -1);
607
608 offset = data._d_ss.getTargetOffsetIndex(i, 2, 0, 0);
609 for (int j=0; j<target_NP; ++j) {
610 P_target_row(offset, j) = 0;
611 }
612 // role of input 2 on component 0 of curl
613 // u_{2,y}
614 P_target_row(offset, 2) = std::pow(data._epsilons(target_index), -1);
615 }
616
617 // output component 1
618 // -u_{2,x} + u_{0,z}
619 {
620 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, 0);
621 for (int j=0; j<target_NP; ++j) {
622 P_target_row(offset, j) = 0;
623 }
624 // role of input 0 on component 1 of curl
625 // u_{0,z}
626 P_target_row(offset, 3) = std::pow(data._epsilons(target_index), -1);
627
628 offset = data._d_ss.getTargetOffsetIndex(i, 1, 1, 0);
629 for (int j=0; j<target_NP; ++j) {
630 P_target_row(offset, j) = 0;
631 }
632 // role of input 1 on component 1 of curl
633 // (no contribution)
634
635 offset = data._d_ss.getTargetOffsetIndex(i, 2, 1, 0);
636 for (int j=0; j<target_NP; ++j) {
637 P_target_row(offset, j) = 0;
638 }
639 // role of input 2 on component 1 of curl
640 // -u_{2,x}
641 P_target_row(offset, 1) = -std::pow(data._epsilons(target_index), -1);
642 }
643
644 // output component 2
645 // u_{1,x} - u_{0,y}
646 {
647 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 2, 0);
648 for (int j=0; j<target_NP; ++j) {
649 P_target_row(offset, j) = 0;
650 }
651 // role of input 0 on component 1 of curl
652 // -u_{0,y}
653 P_target_row(offset, 2) = -std::pow(data._epsilons(target_index), -1);
654
655 offset = data._d_ss.getTargetOffsetIndex(i, 1, 2, 0);
656 for (int j=0; j<target_NP; ++j) {
657 P_target_row(offset, j) = 0;
658 }
659 // role of input 1 on component 1 of curl
660 // u_{1,x}
661 P_target_row(offset, 1) = std::pow(data._epsilons(target_index), -1);
662
663 offset = data._d_ss.getTargetOffsetIndex(i, 2, 2, 0);
664 for (int j=0; j<target_NP; ++j) {
665 P_target_row(offset, j) = 0;
666 }
667 // role of input 2 on component 1 of curl
668 // (no contribution)
669 }
670 });
671 } else if (data._dimensions==2) {
672 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
673 // output component 0
674 // u_{1,y}
675 {
676 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
677 for (int j=0; j<target_NP; ++j) {
678 P_target_row(offset, j) = 0;
679 }
680 // role of input 0 on component 0 of curl
681 // (no contribution)
682
683 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, 0);
684 for (int j=0; j<target_NP; ++j) {
685 P_target_row(offset, j) = 0;
686 }
687 // role of input 1 on component 0 of curl
688 // -u_{1,z}
689 P_target_row(offset, 2) = std::pow(data._epsilons(target_index), -1);
690 }
691
692 // output component 1
693 // -u_{0,x}
694 {
695 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, 0);
696 for (int j=0; j<target_NP; ++j) {
697 P_target_row(offset, j) = 0;
698 }
699 // role of input 0 on component 1 of curl
700 // u_{0,z}
701 P_target_row(offset, 1) = -std::pow(data._epsilons(target_index), -1);
702
703 offset = data._d_ss.getTargetOffsetIndex(i, 1, 1, 0);
704 for (int j=0; j<target_NP; ++j) {
705 P_target_row(offset, j) = 0;
706 }
707 // role of input 1 on component 1 of curl
708 // (no contribution)
709 }
710 });
711 }
712 } else {
713 compadre_kernel_assert_release((false) && "Functionality not yet available.");
714 }
715
716 /*
717 * End of VectorOfScalarClonesTaylorPolynomial basis
718 */
719
720 } else if (data._reconstruction_space == ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial) {
721
722 /*
723 * Beginning of DivergenceFreeVectorTaylorPolynomial basis
724 */
725
726 if (data._operations(i) == TargetOperation::VectorPointEvaluation) {
727 // copied from VectorTaylorPolynomial
728 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int e) {
729 for (int m0=0; m0<data._sampling_multiplier; ++m0) {
730 for (int m1=0; m1<data._sampling_multiplier; ++m1) {
731
732 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(m1+1) /* target is neighbor, but also which component */, 1 /*alpha*/, data._dimensions, data._poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
733
734 int offset = data._d_ss.getTargetOffsetIndex(i, m0 /*in*/, m1 /*out*/, e /*no additional*/);
735 for (int j=0; j<target_NP; ++j) {
736 P_target_row(offset, j) = delta(j);
737 }
738 }
739 }
740 });
741 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
742 } else if (data._operations(i) == TargetOperation::CurlOfVectorPointEvaluation) {
743 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int e) {
744 for (int m0=0; m0<data._sampling_multiplier; ++m0) { // input components
745 for (int m1=0; m1<data._sampling_multiplier; ++m1) { // output components
746 int offset = data._d_ss.getTargetOffsetIndex(i, m0 /*in*/, m1 /*out*/, e /*no additional*/);
747 if (data._dimensions==3) {
748 switch (m1) {
749 case 0:
750 // output component 0
751 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
752 // u2y
753 for (int j=0; j<target_NP; ++j) {
754 P_target_row(offset, j) = delta(j);
755 }
756 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 2 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
757 // -u1z
758 for (int j=0; j<target_NP; ++j) {
759 P_target_row(offset, j) -= delta(j);
760 }
761 // u2y - u1z
762 break;
763 case 1:
764 // output component 1
765 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
766 // -u2x
767 for (int j=0; j<target_NP; ++j) {
768 P_target_row(offset, j) = -delta(j);
769 }
770 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 2 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
771 // u0z
772 for (int j=0; j<target_NP; ++j) {
773 P_target_row(offset, j) += delta(j);
774 }
775 // -u2x + u0z
776 break;
777 default:
778 // output component 2
779 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
780 // u1x
781 for (int j=0; j<target_NP; ++j) {
782 P_target_row(offset, j) = delta(j);
783 }
784 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
785 // -u0y
786 for (int j=0; j<target_NP; ++j) {
787 P_target_row(offset, j) -= delta(j);
788 }
789 // u1x - u0y
790 break;
791 }
792 } else {
793 if (m1==0) {
794 // curl results in 1D output
795 P_target_row(offset, 2) = -std::pow(data._epsilons(target_index), -1);
796 P_target_row(offset, 3) = std::pow(data._epsilons(target_index), -1);
797
798 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
799 // u1x
800 for (int j=0; j<target_NP; ++j) {
801 P_target_row(offset, j) = delta(j);
802 }
803 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
804 // -u0y
805 for (int j=0; j<target_NP; ++j) {
806 P_target_row(offset, j) -= delta(j);
807 }
808 // u1x - u0y
809 }
810 }
811 }
812 }
813 });
814 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
815 } else if (data._operations(i) == TargetOperation::CurlCurlOfVectorPointEvaluation) {
816 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int e) {
817 for (int m0=0; m0<data._sampling_multiplier; ++m0) { // input components
818 for (int m1=0; m1<data._sampling_multiplier; ++m1) { // output components
819 int offset = data._d_ss.getTargetOffsetIndex(i, m0 /*in*/, m1 /*out*/, e /*no additional*/);
820 if (data._dimensions == 3) {
821 switch (m1) {
822 // manually compute the output components
823 case 0:
824 // output component 0
825 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 1 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
826 // u1xy
827 for (int j=0; j<target_NP; ++j) {
828 P_target_row(offset, j) = delta(j);
829 }
830 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 1 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
831 // -u0yy
832 for (int j=0; j<target_NP; ++j) {
833 P_target_row(offset, j) -= delta(j);
834 }
835 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 2 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
836 // u2xz
837 for (int j=0; j<target_NP; ++j) {
838 P_target_row(offset, j) += delta(j);
839 }
840 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 2 /*partial_direction_1*/, 2 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
841 // -u0zz
842 for (int j=0; j<target_NP; ++j) {
843 P_target_row(offset, j) -= delta(j);
844 }
845 // u1xy - u0yy + u2xz - u0zz
846 break;
847 case 1:
848 // output component 1
849 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 0 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
850 // -u1xx
851 for (int j=0; j<target_NP; ++j) {
852 P_target_row(offset, j) = -delta(j);
853 }
854 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 0 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
855 // u0yx
856 for (int j=0; j<target_NP; ++j) {
857 P_target_row(offset, j) += delta(j);
858 }
859 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 2 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
860 // u2yz
861 for (int j=0; j<target_NP; ++j) {
862 P_target_row(offset, j) += delta(j);
863 }
864 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 2 /*partial_direction_1*/, 2 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
865 // -u1zz
866 for (int j=0; j<target_NP; ++j) {
867 P_target_row(offset, j) -= delta(j);
868 }
869 // -u1xx + u0yx + u2yz - u1zz
870 break;
871 default:
872 // output component 2
873 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 0 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
874 // -u2xx
875 for (int j=0; j<target_NP; ++j) {
876 P_target_row(offset, j) = -delta(j);
877 }
878 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 2 /*partial_direction_1*/, 0 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
879 // u0zx
880 for (int j=0; j<target_NP; ++j) {
881 P_target_row(offset, j) += delta(j);
882 }
883 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 1 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
884 // -u2yy
885 for (int j=0; j<target_NP; ++j) {
886 P_target_row(offset, j) -= delta(j);
887 }
888 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 2 /*partial_direction_1*/, 1 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
889 // u1zy
890 for (int j=0; j<target_NP; ++j) {
891 P_target_row(offset, j) += delta(j);
892 }
893 // -u2xx + u0zx - u2yy + u1zy
894 break;
895 }
896 }
897 if (data._dimensions == 2) {
898 // uses curl curl u = grad ( div (u) ) - laplace (u)
899 switch (m1) {
900 case 0:
901 // output component 0
902 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 0 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
903 // u0xx
904 for (int j=0; j<target_NP; ++j) {
905 P_target_row(offset, j) = delta(j);
906 }
907 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 0 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
908 // u1yx
909 for (int j=0; j<target_NP; ++j) {
910 P_target_row(offset, j) += delta(j);
911 }
912 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 0 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
913 // -u0xx
914 for (int j=0; j<target_NP; ++j) {
915 P_target_row(offset, j) -= delta(j);
916 }
917 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 1 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
918 // -u0yy
919 for (int j=0; j<target_NP; ++j) {
920 P_target_row(offset, j) -= delta(j);
921 }
922 // u0xx + u1yx - u0xx - u0yy
923 break;
924 case 1:
925 // output component 1
926 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 1 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
927 // u0xy
928 for (int j=0; j<target_NP; ++j) {
929 P_target_row(offset, j) = delta(j);
930 }
931 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 1 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
932 // u1yy
933 for (int j=0; j<target_NP; ++j) {
934 P_target_row(offset, j) += delta(j);
935 }
936 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 0 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
937 // -u1xx
938 for (int j=0; j<target_NP; ++j) {
939 P_target_row(offset, j) -= delta(j);
940 }
941 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 1 /*partial_direction_2*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
942 // -u1yy
943 for (int j=0; j<target_NP; ++j) {
944 P_target_row(offset, j) -= delta(j);
945 }
946 // u0xy + u1yy - u1xx - u1yy
947 break;
948 }
949 }
950 }
951 }
952 });
953 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
954 } else if (data._operations(i) == TargetOperation::GradientOfVectorPointEvaluation) {
955 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int e) {
956 for (int m0=0; m0<data._sampling_multiplier; ++m0) { // input components
957 for (int m1=0; m1<data._sampling_multiplier; ++m1) { // output components 1
958 for (int m2=0; m2<data._sampling_multiplier; ++m2) { // output components 2
959 int offset = data._d_ss.getTargetOffsetIndex(i, m0 /*in*/, m1*data._sampling_multiplier+m2 /*out*/, e);
960 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(m1+1) /* target is neighbor */, 1 /*alpha*/, m2 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
961 for (int j=0; j<target_NP; ++j) {
962 P_target_row(offset, j) = delta(j);
963 }
964 }
965 }
966 }
967 });
968 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
969 } else {
970 compadre_kernel_assert_release((false) && "Functionality not yet available.");
971 }
972
973 /*
974 * End of DivergenceFreeVectorTaylorPolynomial basis
975 */
976
977 } else if (data._reconstruction_space == ReconstructionSpace::BernsteinPolynomial) {
978
979 /*
980 * Beginning of BernsteinPolynomial basis
981 */
982 // TODO: Copied from ScalarTaylorPolynomial section. Could be defined once in ScalarTaylorPolynomial if refactored.
983 if (data._operations(i) == TargetOperation::ScalarPointEvaluation || (data._operations(i) == TargetOperation::VectorPointEvaluation && data._dimensions == 1) /* vector is a scalar in 1D */) {
984 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
985 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, data._dimensions, data._poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::BernsteinPolynomial, PointSample, j);
986 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
987 for (int k=0; k<target_NP; ++k) {
988 P_target_row(offset, k) = delta(k);
989 }
990 });
991 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
992 } else if (data._operations(i) == TargetOperation::GradientOfScalarPointEvaluation) {
993 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
994 for (int d=0; d<data._dimensions; ++d) {
995 int offset = data._d_ss.getTargetOffsetIndex(i, 0, d, j);
996 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
997 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, d /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::BernsteinPolynomial, PointSample, j);
998 }
999 });
1000 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
1001 } else if (data._operations(i) == TargetOperation::PartialXOfScalarPointEvaluation) {
1002 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
1003 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
1004 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
1005 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 0 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::BernsteinPolynomial, PointSample, j);
1006 });
1007 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
1008 } else if (data._operations(i) == TargetOperation::PartialYOfScalarPointEvaluation) {
1009 compadre_kernel_assert_release(data._dimensions>1 && "PartialYOfScalarPointEvaluation requested for dim < 2");
1010 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
1011 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
1012 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
1013 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 1 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::BernsteinPolynomial, PointSample, j);
1014 });
1015 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
1016 } else if (data._operations(i) == TargetOperation::PartialZOfScalarPointEvaluation) {
1017 compadre_kernel_assert_release(data._dimensions>2 && "PartialZOfScalarPointEvaluation requested for dim < 3");
1018 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
1019 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
1020 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
1021 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 2 /*partial_direction*/, data._dimensions, data._poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::BernsteinPolynomial, PointSample, j);
1022 });
1023 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
1024 } else {
1025 compadre_kernel_assert_release((false) && "Functionality not yet available.");
1026 }
1027
1028 /*
1029 * End of BernsteinPolynomial basis
1030 */
1031
1032 } else {
1033 compadre_kernel_assert_release((false) && "Functionality not yet available.");
1034 }
1035
1036 (void)additional_evaluation_sites_handled; // handle unused variable warning/error
1037 (void)additional_evaluation_sites_need_handled; // handle unused variable warning/error
1038 compadre_kernel_assert_release(((additional_evaluation_sites_need_handled && additional_evaluation_sites_handled) || (!additional_evaluation_sites_need_handled)) && "Auxiliary evaluation coordinates are specified by user, but are calling a target operation that can not handle evaluating additional sites.");
1039 } // !operation_handled
1040 }
1041 teamMember.team_barrier();
1042}
1043
1044/*! \brief Evaluates a polynomial basis for the curvature with a gradient target functional applied
1045
1046 data._operations is used by this function which is set through a modifier function
1047
1048 \param data [in] - GMLSBasisData struct
1049 \param teamMember [in] - Kokkos::TeamPolicy member type (created by parallel_for)
1050 \param delta [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large is the _basis_multipler*the dimension of the polynomial basis.
1051 \param thread_workspace [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large as the _curvature_poly_order*the spatial dimension of the polynomial basis.
1052 \param P_target_row [out] - 1D Kokkos View where the evaluation of the polynomial basis is stored
1053 \param V [in] - orthonormal basis matrix size _dimensions * _dimensions whose first _dimensions-1 columns are an approximation of the tangent plane
1054*/
1055template <typename TargetData>
1056KOKKOS_INLINE_FUNCTION
1057void computeCurvatureFunctionals(const TargetData& data, const member_type& teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, device_unmanaged_matrix_right_type P_target_row, const device_unmanaged_matrix_right_type* V, const local_index_type local_neighbor_index = -1) {
1058
1059 compadre_kernel_assert_release((thread_workspace.extent_int(0)>=(data._curvature_poly_order+1)*data._local_dimensions) && "Workspace thread_workspace not large enough.");
1060
1061 const int target_index = data._initial_index_for_batch + teamMember.league_rank();
1062
1063 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, P_target_row.extent(0)), [&] (const int j) {
1064 Kokkos::parallel_for(Kokkos::ThreadVectorRange(teamMember, P_target_row.extent(1)),
1065 [&] (const int& k) {
1066 P_target_row(j,k) = 0;
1067 });
1068 });
1069 for (int j = 0; j < delta.extent_int(0); ++j) {
1070 delta(j) = 0;
1071 }
1072 for (int j = 0; j < thread_workspace.extent_int(0); ++j) {
1073 thread_workspace(j) = 0;
1074 }
1075 teamMember.team_barrier();
1076
1077 // not const b.c. of gcc 7.2 issue
1078 int manifold_NP = GMLS::getNP(data._curvature_poly_order, data._dimensions-1, ReconstructionSpace::ScalarTaylorPolynomial);
1079 for (size_t i=0; i<data._curvature_support_operations.size(); ++i) {
1080 if (data._curvature_support_operations(i) == TargetOperation::ScalarPointEvaluation) {
1081 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1082 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1083 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, local_neighbor_index, 0 /*alpha*/, data._dimensions-1, data._curvature_poly_order, false /*bool on only specific order*/, V, ReconstructionSpace::ScalarTaylorPolynomial, PointSample);
1084 for (int j=0; j<manifold_NP; ++j) {
1085 P_target_row(offset, j) = delta(j);
1086 }
1087 });
1088 } else if (data._curvature_support_operations(i) == TargetOperation::GradientOfScalarPointEvaluation) {
1089 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1090 //int offset = i*manifold_NP;
1091 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1092 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, local_neighbor_index, 0 /*alpha*/, 0 /*partial_direction*/, data._dimensions-1, data._curvature_poly_order, false /*specific order only*/, V, ReconstructionSpace::ScalarTaylorPolynomial, PointSample);
1093 for (int j=0; j<manifold_NP; ++j) {
1094 P_target_row(offset, j) = delta(j);
1095 }
1096 if (data._dimensions>2) { // _dimensions-1 > 1
1097 //offset = (i+1)*manifold_NP;
1098 offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, 0);
1099 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, local_neighbor_index, 0 /*alpha*/, 1 /*partial_direction*/, data._dimensions-1, data._curvature_poly_order, false /*specific order only*/, V, ReconstructionSpace::ScalarTaylorPolynomial, PointSample);
1100 for (int j=0; j<manifold_NP; ++j) {
1101 P_target_row(offset, j) = delta(j);
1102 }
1103 }
1104 });
1105 } else {
1106 compadre_kernel_assert_release((false) && "Functionality not yet available.");
1107 }
1108 }
1109 teamMember.team_barrier();
1110}
1111
1112/*! \brief Evaluates a polynomial basis with a target functional applied, using information from the manifold curvature
1113
1114 data._operations is used by this function which is set through a modifier function
1115
1116 \param data [in] - GMLSBasisData struct
1117 \param teamMember [in] - Kokkos::TeamPolicy member type (created by parallel_for)
1118 \param delta [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large is the _basis_multipler*the dimension of the polynomial basis.
1119 \param thread_workspace [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large as the _curvature_poly_order*the spatial dimension of the polynomial basis.
1120 \param P_target_row [out] - 1D Kokkos View where the evaluation of the polynomial basis is stored
1121 \param V [in] - orthonormal basis matrix size _dimensions * _dimensions whose first _dimensions-1 columns are an approximation of the tangent plane
1122 \param curvature_coefficients [in] - polynomial coefficients for curvature
1123*/
1124template <typename TargetData>
1125KOKKOS_INLINE_FUNCTION
1126void computeTargetFunctionalsOnManifold(const TargetData& data, const member_type& teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, device_unmanaged_matrix_right_type P_target_row, device_unmanaged_matrix_right_type V, device_unmanaged_vector_type curvature_coefficients) {
1127
1128 compadre_kernel_assert_release((thread_workspace.extent_int(0)>=(data._poly_order+1)*data._local_dimensions) && "Workspace thread_workspace not large enough.");
1129
1130 // only designed for 2D manifold embedded in 3D space
1131 const int target_index = data._initial_index_for_batch + teamMember.league_rank();
1132 // not const b.c. of gcc 7.2 issue
1133 int target_NP = GMLS::getNP(data._poly_order, data._dimensions-1, data._reconstruction_space);
1134
1135 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, P_target_row.extent(0)), [&] (const int j) {
1136 Kokkos::parallel_for(Kokkos::ThreadVectorRange(teamMember, P_target_row.extent(1)),
1137 [&] (const int& k) {
1138 P_target_row(j,k) = 0;
1139 });
1140 });
1141 for (int j = 0; j < delta.extent_int(0); ++j) {
1142 delta(j) = 0;
1143 }
1144 for (int j = 0; j < thread_workspace.extent_int(0); ++j) {
1145 thread_workspace(j) = 0;
1146 }
1147 teamMember.team_barrier();
1148
1149 // determine if additional evaluation sites are requested by user and handled by target operations
1150 bool additional_evaluation_sites_need_handled =
1151 (data._additional_pc._source_coordinates.extent(0) > 0) ? true : false; // additional evaluation sites are specified
1152
1153 const int num_evaluation_sites = data._additional_pc._nla.getNumberOfNeighborsDevice(target_index) + 1;
1154
1155 for (size_t i=0; i<data._operations.size(); ++i) {
1156
1157 bool additional_evaluation_sites_handled = false; // target operations that can handle these sites should flip this flag
1158
1159 bool operation_handled = true;
1160
1161 // USER defined targets on the manifold should be added to this file
1162 // if the USER defined targets don't catch this operation, then operation_handled will be false
1163 #include "Compadre_USER_ManifoldTargetFunctionals.hpp"
1164
1165 // if the user didn't handle the operation, we pass it along to the toolkit's targets
1166 if (!operation_handled) {
1167 if (data._dimensions>2) {
1168 if (data._operations(i) == TargetOperation::ScalarPointEvaluation) {
1169 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
1170 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, data._dimensions-1, data._poly_order, false /*bool on only specific order*/, &V, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
1171 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
1172 for (int k=0; k<target_NP; ++k) {
1173 P_target_row(offset, k) = delta(k);
1174 }
1175 });
1176 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
1177 } else if (data._operations(i) == TargetOperation::VectorPointEvaluation) {
1178 // vector basis
1179 if (data._reconstruction_space_rank == 1) {
1180 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int k) {
1181 // output component 0
1182 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, data._dimensions-1, data._poly_order, false /*bool on only specific order*/, &V, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, k);
1183 for (int m=0; m<data._sampling_multiplier; ++m) {
1184 int output_components = data._basis_multiplier;
1185 for (int c=0; c<output_components; ++c) {
1186 int offset = data._d_ss.getTargetOffsetIndex(i, m /*in*/, c /*out*/, k/*additional*/);
1187 // for the case where data._sampling_multiplier is > 1,
1188 // this approach relies on c*target_NP being equivalent to P_target_row(offset, j) where offset is
1189 // data._d_ss.getTargetOffsetIndex(i, m /*in*/, c /*out*/, e/*additional*/)*data._basis_multiplier*target_NP;
1190 for (int j=0; j<target_NP; ++j) {
1191 P_target_row(offset, c*target_NP + j) = delta(j);
1192 }
1193 }
1194 }
1195 });
1196 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
1197 // scalar basis times number of components in the vector
1198 } else if (data._reconstruction_space_rank == 0) {
1199 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int k) {
1200 // output component 0
1201 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, data._dimensions-1, data._poly_order, false /*bool on only specific order*/, &V, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, k);
1202 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, k);
1203 for (int j=0; j<target_NP; ++j) {
1204 P_target_row(offset, j) = delta(j);
1205 }
1206 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, k);
1207 for (int j=0; j<target_NP; ++j) {
1208 P_target_row(offset, j) = 0;
1209 }
1210
1211 // output component 1
1212 offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, k);
1213 for (int j=0; j<target_NP; ++j) {
1214 P_target_row(offset, j) = 0;
1215 }
1216 offset = data._d_ss.getTargetOffsetIndex(i, 1, 1, k);
1217 for (int j=0; j<target_NP; ++j) {
1218 P_target_row(offset, j) = delta(j);
1219 }
1220 });
1221 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
1222 } else {
1223 compadre_kernel_assert_release((false) && "Functionality not yet available.");
1224 }
1225
1226 } else if (data._operations(i) == TargetOperation::LaplacianOfScalarPointEvaluation) {
1227 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1228
1229 double h = data._epsilons(target_index);
1230 double a1=0, a2=0, a3=0, a4=0, a5=0;
1231 if (data._curvature_poly_order > 0) {
1232 a1 = curvature_coefficients(1);
1233 a2 = curvature_coefficients(2);
1234 }
1235 if (data._curvature_poly_order > 1) {
1236 a3 = curvature_coefficients(3);
1237 a4 = curvature_coefficients(4);
1238 a5 = curvature_coefficients(5);
1239 }
1240 double den = (h*h + a1*a1 + a2*a2);
1241
1242 // Gaussian Curvature sanity check
1243 //double K_curvature = ( - a4*a4 + a3*a5) / den / den;
1244 //std::cout << "Gaussian curvature is: " << K_curvature << std::endl;
1245
1246
1247 const int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1248 for (int j=0; j<target_NP; ++j) {
1249 P_target_row(offset, j) = 0;
1250 }
1251 // scaled
1252 if (data._poly_order > 0 && data._curvature_poly_order > 1) {
1253 P_target_row(offset, 1) = (-a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5))/den/den/(h*h);
1254 P_target_row(offset, 2) = (-a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5))/den/den/(h*h);
1255 }
1256 if (data._poly_order > 1 && data._curvature_poly_order > 0) {
1257 P_target_row(offset, 3) = (h*h+a2*a2)/den/(h*h);
1258 P_target_row(offset, 4) = -2*a1*a2/den/(h*h);
1259 P_target_row(offset, 5) = (h*h+a1*a1)/den/(h*h);
1260 }
1261
1262 });
1263 } else if (data._operations(i) == TargetOperation::ChainedStaggeredLaplacianOfScalarPointEvaluation) {
1264 if (data._reconstruction_space == ReconstructionSpace::VectorTaylorPolynomial) {
1265 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1266
1267 double h = data._epsilons(target_index);
1268 double a1=0, a2=0, a3=0, a4=0, a5=0;
1269 if (data._curvature_poly_order > 0) {
1270 a1 = curvature_coefficients(1);
1271 a2 = curvature_coefficients(2);
1272 }
1273 if (data._curvature_poly_order > 1) {
1274 a3 = curvature_coefficients(3);
1275 a4 = curvature_coefficients(4);
1276 a5 = curvature_coefficients(5);
1277 }
1278 double den = (h*h + a1*a1 + a2*a2);
1279
1280 double c0a = -a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5)/den/den/h;
1281 double c1a = (h*h+a2*a2)/den/h;
1282 double c2a = -a1*a2/den/h;
1283
1284 double c0b = -a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5)/den/den/h;
1285 double c1b = -a1*a2/den/h;
1286 double c2b = (h*h+a1*a1)/den/h;
1287
1288 // 1st input component
1289 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1290 for (int j=0; j<target_NP; ++j) {
1291 P_target_row(offset, j) = 0;
1292 P_target_row(offset, target_NP + j) = 0;
1293 }
1294 P_target_row(offset, 0) = c0a;
1295 P_target_row(offset, 1) = c1a;
1296 P_target_row(offset, 2) = c2a;
1297 P_target_row(offset, target_NP + 0) = c0b;
1298 P_target_row(offset, target_NP + 1) = c1b;
1299 P_target_row(offset, target_NP + 2) = c2b;
1300 });
1301 } else if (data._reconstruction_space == ReconstructionSpace::ScalarTaylorPolynomial) {
1302 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1303
1304 double h = data._epsilons(target_index);
1305 double a1=0, a2=0, a3=0, a4=0, a5=0;
1306 if (data._curvature_poly_order > 0) {
1307 a1 = curvature_coefficients(1);
1308 a2 = curvature_coefficients(2);
1309 }
1310 if (data._curvature_poly_order > 1) {
1311 a3 = curvature_coefficients(3);
1312 a4 = curvature_coefficients(4);
1313 a5 = curvature_coefficients(5);
1314 }
1315 double den = (h*h + a1*a1 + a2*a2);
1316
1317 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1318 for (int j=0; j<target_NP; ++j) {
1319 P_target_row(offset, j) = 0;
1320 }
1321
1322 // verified
1323 if (data._poly_order > 0 && data._curvature_poly_order > 1) {
1324 P_target_row(offset, 1) = (-a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5))/den/den/(h*h);
1325 P_target_row(offset, 2) = (-a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5))/den/den/(h*h);
1326 }
1327 if (data._poly_order > 1 && data._curvature_poly_order > 0) {
1328 P_target_row(offset, 3) = (h*h+a2*a2)/den/(h*h);
1329 P_target_row(offset, 4) = -2*a1*a2/den/(h*h);
1330 P_target_row(offset, 5) = (h*h+a1*a1)/den/(h*h);
1331 }
1332
1333 });
1334 } else {
1335 compadre_kernel_assert_release((false) && "Functionality not yet available.");
1336 }
1337 } else if (data._operations(i) == TargetOperation::VectorLaplacianPointEvaluation) {
1338 // vector basis
1339 if (data._reconstruction_space_rank == 1) {
1340 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1341
1342 double h = data._epsilons(target_index);
1343 double a1=0, a2=0, a3=0, a4=0, a5=0;
1344 if (data._curvature_poly_order > 0) {
1345 a1 = curvature_coefficients(1);
1346 a2 = curvature_coefficients(2);
1347 }
1348 if (data._curvature_poly_order > 1) {
1349 a3 = curvature_coefficients(3);
1350 a4 = curvature_coefficients(4);
1351 a5 = curvature_coefficients(5);
1352 }
1353 double den = (h*h + a1*a1 + a2*a2);
1354
1355 for (int j=0; j<target_NP; ++j) {
1356 delta(j) = 0;
1357 }
1358
1359 // 1/Sqrt[Det[G[r, s]]])*Div[Sqrt[Det[G[r, s]]]*Inv[G]*P
1360 if (data._poly_order > 0 && data._curvature_poly_order > 1) {
1361 delta(1) = (-a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4+(h*h+a1*a1)*a5))/den/den/(h*h);
1362 delta(2) = (-a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4+(h*h+a1*a1)*a5))/den/den/(h*h);
1363 }
1364 if (data._poly_order > 1 && data._curvature_poly_order > 0) {
1365 delta(3) = (h*h+a2*a2)/den/(h*h);
1366 delta(4) = -2*a1*a2/den/(h*h);
1367 delta(5) = (h*h+a1*a1)/den/(h*h);
1368 }
1369
1370 // output component 0
1371 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1372 for (int j=0; j<target_NP; ++j) {
1373 P_target_row(offset, j) = delta(j);
1374 P_target_row(offset, target_NP + j) = 0;
1375 }
1376 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, 0);
1377 for (int j=0; j<target_NP; ++j) {
1378 P_target_row(offset, j) = 0;
1379 P_target_row(offset, target_NP + j) = 0;
1380 }
1381
1382 // output component 1
1383 offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, 0);
1384 for (int j=0; j<target_NP; ++j) {
1385 P_target_row(offset, j) = 0;
1386 P_target_row(offset, target_NP + j) = 0;
1387 }
1388 offset = data._d_ss.getTargetOffsetIndex(i, 1, 1, 0);
1389 for (int j=0; j<target_NP; ++j) {
1390 P_target_row(offset, j) = 0;
1391 P_target_row(offset, target_NP + j) = delta(j);
1392 }
1393
1394 });
1395 // scalar basis times number of components in the vector
1396 } else if (data._reconstruction_space_rank == 0) {
1397 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1398
1399 double h = data._epsilons(target_index);
1400 double a1=0, a2=0, a3=0, a4=0, a5=0;
1401 if (data._curvature_poly_order > 0) {
1402 a1 = curvature_coefficients(1);
1403 a2 = curvature_coefficients(2);
1404 }
1405 if (data._curvature_poly_order > 1) {
1406 a3 = curvature_coefficients(3);
1407 a4 = curvature_coefficients(4);
1408 a5 = curvature_coefficients(5);
1409 }
1410 double den = (h*h + a1*a1 + a2*a2);
1411
1412 for (int j=0; j<target_NP; ++j) {
1413 delta(j) = 0;
1414 }
1415
1416 // 1/Sqrt[Det[G[r, s]]])*Div[Sqrt[Det[G[r, s]]]*Inv[G]*P
1417 if (data._poly_order > 0 && data._curvature_poly_order > 1) {
1418 delta(1) = (-a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4+(h*h+a1*a1)*a5))/den/den/(h*h);
1419 delta(2) = (-a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4+(h*h+a1*a1)*a5))/den/den/(h*h);
1420 }
1421 if (data._poly_order > 1 && data._curvature_poly_order > 0) {
1422 delta(3) = (h*h+a2*a2)/den/(h*h);
1423 delta(4) = -2*a1*a2/den/(h*h);
1424 delta(5) = (h*h+a1*a1)/den/(h*h);
1425 }
1426
1427 // output component 0
1428 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1429 for (int j=0; j<target_NP; ++j) {
1430 P_target_row(offset, j) = delta(j);
1431 }
1432 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, 0);
1433 for (int j=0; j<target_NP; ++j) {
1434 P_target_row(offset, j) = 0;
1435 }
1436
1437 // output component 1
1438 offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, 0);
1439 for (int j=0; j<target_NP; ++j) {
1440 P_target_row(offset, j) = 0;
1441 }
1442 offset = data._d_ss.getTargetOffsetIndex(i, 1, 1, 0);
1443 for (int j=0; j<target_NP; ++j) {
1444 P_target_row(offset, j) = delta(j);
1445 }
1446 });
1447 } else {
1448 compadre_kernel_assert_release((false) && "Functionality not yet available.");
1449 }
1450 } else if (data._operations(i) == TargetOperation::GradientOfScalarPointEvaluation) {
1451 if (data._reconstruction_space_rank == 0
1452 && (data._polynomial_sampling_functional == PointSample
1453 || data._polynomial_sampling_functional == ManifoldVectorPointSample)) {
1454 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1455
1456 double h = data._epsilons(target_index);
1457 double a1 = curvature_coefficients(1);
1458 double a2 = curvature_coefficients(2);
1459
1460 double q1 = (h*h + a2*a2)/(h*h*h + h*a1*a1 + h*a2*a2);
1461 double q2 = -(a1*a2)/(h*h*h + h*a1*a1 + h*a2*a2);
1462 double q3 = (h*h + a1*a1)/(h*h*h + h*a1*a1 + h*a2*a2);
1463
1464 double t1a = q1*1 + q2*0;
1465 double t2a = q1*0 + q2*1;
1466
1467 double t1b = q2*1 + q3*0;
1468 double t2b = q2*0 + q3*1;
1469
1470 // scaled
1471 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1472 for (int j=0; j<target_NP; ++j) {
1473 P_target_row(offset, j) = 0;
1474 }
1475 if (data._poly_order > 0 && data._curvature_poly_order > 0) {
1476 P_target_row(offset, 1) = t1a + t2a;
1477 P_target_row(offset, 2) = 0;
1478 }
1479
1480 offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, 0);
1481 for (int j=0; j<target_NP; ++j) {
1482 P_target_row(offset, j) = 0;
1483 }
1484 if (data._poly_order > 0 && data._curvature_poly_order > 0) {
1485 P_target_row(offset, 1) = 0;
1486 P_target_row(offset, 2) = t1b + t2b;
1487 }
1488
1489 });
1490 // staggered gradient w/ edge integrals performed by numerical integration
1491 // with a vector basis
1492 } else if (data._reconstruction_space_rank == 1
1493 && data._polynomial_sampling_functional
1494 == StaggeredEdgeIntegralSample) {
1495 compadre_kernel_assert_release((false) && "Functionality not yet available.");
1496
1497 // staggered gradient w/ edge integrals known analytically, using a basis
1498 // of potentials
1499 } else if (data._reconstruction_space_rank == 0
1500 && data._polynomial_sampling_functional
1501 == StaggeredEdgeAnalyticGradientIntegralSample) {
1502 compadre_kernel_assert_release((false) && "Functionality not yet available.");
1503
1504 } else {
1505 compadre_kernel_assert_release((false) && "Functionality not yet available.");
1506 }
1507 } else if (data._operations(i) == TargetOperation::DivergenceOfVectorPointEvaluation) {
1508 // vector basis
1509 if (data._reconstruction_space_rank == 1
1510 && data._polynomial_sampling_functional == ManifoldVectorPointSample) {
1511 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1512
1513 double h = data._epsilons(target_index);
1514 double a1=0, a2=0, a3=0, a4=0, a5=0;
1515 if (data._curvature_poly_order > 0) {
1516 a1 = curvature_coefficients(1);
1517 a2 = curvature_coefficients(2);
1518 }
1519 if (data._curvature_poly_order > 1) {
1520 a3 = curvature_coefficients(3);
1521 a4 = curvature_coefficients(4);
1522 a5 = curvature_coefficients(5);
1523 }
1524 double den = (h*h + a1*a1 + a2*a2);
1525
1526 // 1/Sqrt[Det[G[r, s]]])*Div[Sqrt[Det[G[r, s]]]*P
1527 // i.e. P recovers G^{-1}*grad of scalar
1528 double c0a = (a1*a3+a2*a4)/(h*den);
1529 double c1a = 1./h;
1530 double c2a = 0;
1531
1532 double c0b = (a1*a4+a2*a5)/(h*den);
1533 double c1b = 0;
1534 double c2b = 1./h;
1535
1536 // 1st input component
1537 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1538 for (int j=0; j<target_NP; ++j) {
1539 P_target_row(offset, j) = 0;
1540 P_target_row(offset, target_NP + j) = 0;
1541 }
1542 P_target_row(offset, 0) = c0a;
1543 P_target_row(offset, 1) = c1a;
1544 P_target_row(offset, 2) = c2a;
1545
1546 // 2nd input component
1547 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, 0);
1548 for (int j=0; j<target_NP; ++j) {
1549 P_target_row(offset, j) = 0;
1550 P_target_row(offset, target_NP + j) = 0;
1551 }
1552 P_target_row(offset, target_NP + 0) = c0b;
1553 P_target_row(offset, target_NP + 1) = c1b;
1554 P_target_row(offset, target_NP + 2) = c2b;
1555 });
1556 // scalar basis times number of components in the vector
1557 } else if (data._reconstruction_space_rank == 0
1558 && data._polynomial_sampling_functional == ManifoldVectorPointSample) {
1559 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1560
1561 double h = data._epsilons(target_index);
1562 double a1=0, a2=0, a3=0, a4=0, a5=0;
1563 if (data._curvature_poly_order > 0) {
1564 a1 = curvature_coefficients(1);
1565 a2 = curvature_coefficients(2);
1566 }
1567 if (data._curvature_poly_order > 1) {
1568 a3 = curvature_coefficients(3);
1569 a4 = curvature_coefficients(4);
1570 a5 = curvature_coefficients(5);
1571 }
1572 double den = (h*h + a1*a1 + a2*a2);
1573
1574 // 1/Sqrt[Det[G[r, s]]])*Div[Sqrt[Det[G[r, s]]]*P
1575 // i.e. P recovers G^{-1}*grad of scalar
1576 double c0a = (a1*a3+a2*a4)/(h*den);
1577 double c1a = 1./h;
1578 double c2a = 0;
1579
1580 double c0b = (a1*a4+a2*a5)/(h*den);
1581 double c1b = 0;
1582 double c2b = 1./h;
1583
1584 // 1st input component
1585 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1586 for (int j=0; j<target_NP; ++j) {
1587 P_target_row(offset, j) = 0;
1588 }
1589 P_target_row(offset, 0) = c0a;
1590 P_target_row(offset, 1) = c1a;
1591 P_target_row(offset, 2) = c2a;
1592
1593 // 2nd input component
1594 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, 0);
1595 for (int j=0; j<target_NP; ++j) {
1596 P_target_row(offset, j) = 0;
1597 }
1598 P_target_row(offset, 0) = c0b;
1599 P_target_row(offset, 1) = c1b;
1600 P_target_row(offset, 2) = c2b;
1601 });
1602 // staggered divergence acting on vector polynomial space
1603 } else if (data._reconstruction_space_rank == 1
1604 && data._polynomial_sampling_functional
1605 == StaggeredEdgeIntegralSample) {
1606
1607 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1608
1609 double h = data._epsilons(target_index);
1610 double a1=0, a2=0, a3=0, a4=0, a5=0;
1611 if (data._curvature_poly_order > 0) {
1612 a1 = curvature_coefficients(1);
1613 a2 = curvature_coefficients(2);
1614 }
1615 if (data._curvature_poly_order > 1) {
1616 a3 = curvature_coefficients(3);
1617 a4 = curvature_coefficients(4);
1618 a5 = curvature_coefficients(5);
1619 }
1620 double den = (h*h + a1*a1 + a2*a2);
1621
1622 // 1/Sqrt[Det[G[r, s]]])*Div[Sqrt[Det[G[r, s]]]*Inv[G].P
1623 // i.e. P recovers grad of scalar
1624 double c0a = -a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5)/den/den/h;
1625 double c1a = (h*h+a2*a2)/den/h;
1626 double c2a = -a1*a2/den/h;
1627
1628 double c0b = -a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5)/den/den/h;
1629 double c1b = -a1*a2/den/h;
1630 double c2b = (h*h+a1*a1)/den/h;
1631
1632 // 1st input component
1633 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1634 for (int j=0; j<target_NP; ++j) {
1635 P_target_row(offset, j) = 0;
1636 P_target_row(offset, target_NP + j) = 0;
1637 }
1638 P_target_row(offset, 0) = c0a;
1639 P_target_row(offset, 1) = c1a;
1640 P_target_row(offset, 2) = c2a;
1641 P_target_row(offset, target_NP + 0) = c0b;
1642 P_target_row(offset, target_NP + 1) = c1b;
1643 P_target_row(offset, target_NP + 2) = c2b;
1644
1645 });
1646 } else {
1647 compadre_kernel_assert_release((false) && "Functionality not yet available.");
1648 }
1649 } else if (data._operations(i) == TargetOperation::GaussianCurvaturePointEvaluation) {
1650 double h = data._epsilons(target_index);
1651
1652 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int k) {
1653 XYZ relative_coord;
1654 if (k > 0) {
1655 for (int d=0; d<data._dimensions-1; ++d) {
1656 // k indexing is for evaluation site, which includes target site
1657 // the k-1 converts to the local index for ADDITIONAL evaluation sites
1658 relative_coord[d] = data._additional_pc.getNeighborCoordinate(target_index, k-1, d, &V);
1659 relative_coord[d] -= data._pc.getTargetCoordinate(target_index, d, &V);
1660 }
1661 } else {
1662 for (int j=0; j<3; ++j) relative_coord[j] = 0;
1663 }
1664
1665 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, k);
1666 P_target_row(offset, 0) = GaussianCurvature(curvature_coefficients, h, relative_coord.x, relative_coord.y);
1667 });
1668 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
1669 } else if (data._operations(i) == TargetOperation::CurlOfVectorPointEvaluation) {
1670 double h = data._epsilons(target_index);
1671
1672 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int k) {
1673 int alphax, alphay;
1674 XYZ relative_coord;
1675 if (k > 0) {
1676 for (int d=0; d<data._dimensions-1; ++d) {
1677 // k indexing is for evaluation site, which includes target site
1678 // the k-1 converts to the local index for ADDITIONAL evaluation sites
1679 relative_coord[d] = data._additional_pc.getNeighborCoordinate(target_index, k-1, d, &V);
1680 relative_coord[d] -= data._pc.getTargetCoordinate(target_index, d, &V);
1681 }
1682 } else {
1683 for (int j=0; j<3; ++j) relative_coord[j] = 0;
1684 }
1685
1686 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, k);
1687 int index = 0;
1688 for (int n = 0; n <= data._poly_order; n++){
1689 for (alphay = 0; alphay <= n; alphay++){
1690 alphax = n - alphay;
1691 P_target_row(offset, index) = SurfaceCurlOfScalar(curvature_coefficients, h, relative_coord.x, relative_coord.y, alphax, alphay, 0);
1692 index++;
1693 }
1694 }
1695
1696 offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, k);
1697 index = 0;
1698 for (int n = 0; n <= data._poly_order; n++){
1699 for (alphay = 0; alphay <= n; alphay++){
1700 alphax = n - alphay;
1701 P_target_row(offset, index) = SurfaceCurlOfScalar(curvature_coefficients, h, relative_coord.x, relative_coord.y, alphax, alphay, 1);
1702 index++;
1703 }
1704 }
1705 });
1706 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
1707 } else if (data._operations(i) == TargetOperation::CellAverageEvaluation ||
1708 data._operations(i) == TargetOperation::CellIntegralEvaluation) {
1709 compadre_kernel_assert_debug(data._local_dimensions==2 &&
1710 "CellAverageEvaluation and CellIntegralEvaluation only support 2d or 3d with 2d manifold");
1711 const double factorial[15] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
1712 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1713
1714 double cutoff_p = data._epsilons(target_index);
1715 int alphax, alphay;
1716 double alphaf;
1717
1718 // global dimension cannot be determined in a constexpr way, so we use a largest case scenario
1719 // of dimensions 3 for _global_dimension
1720 double G_data[3*3]; //data._global_dimensions*3
1721 double triangle_coords[3*3]; //data._global_dimensions*3
1722 for (int j=0; j<data._global_dimensions*3; ++j) G_data[j] = 0;
1723 for (int j=0; j<data._global_dimensions*3; ++j) triangle_coords[j] = 0;
1724 // 3 is for # vertices in sub-triangle
1725 device_unmanaged_matrix_right_type G(G_data, data._global_dimensions, 3);
1726 device_unmanaged_matrix_right_type triangle_coords_matrix(triangle_coords, data._global_dimensions, 3);
1727
1728 double radius = 0.0;
1729 for (int j=0; j<data._global_dimensions; ++j) {
1730 // midpoint
1731 triangle_coords_matrix(j, 0) = data._pc.getTargetCoordinate(target_index, j);
1732 radius += triangle_coords_matrix(j, 0)*triangle_coords_matrix(j, 0);
1733 }
1734 radius = std::sqrt(radius);
1735
1736 // NaN in entry (data._global_dimensions) is a convention for indicating fewer vertices
1737 // for this cell and NaN is checked by entry!=entry
1738 int num_vertices = 0;
1739 for (int j=0; j<data._target_extra_data.extent_int(1); ++j) {
1740 auto val = data._target_extra_data(target_index, j);
1741 if (val != val) {
1742 break;
1743 } else {
1744 num_vertices++;
1745 }
1746 }
1747 num_vertices = num_vertices / data._global_dimensions;
1748 auto T = triangle_coords_matrix;
1749
1750 // loop over each two vertices
1751 XYZ relative_coord;
1752 double entire_cell_area = 0.0;
1753 for (int v=0; v<num_vertices; ++v) {
1754 int v1 = v;
1755 int v2 = (v+1) % num_vertices;
1756
1757 for (int j=0; j<data._global_dimensions; ++j) {
1758 triangle_coords_matrix(j,1) = data._target_extra_data(target_index,
1759 v1*data._global_dimensions+j)
1760 - triangle_coords_matrix(j,0);
1761 triangle_coords_matrix(j,2) = data._target_extra_data(target_index,
1762 v2*data._global_dimensions+j)
1763 - triangle_coords_matrix(j,0);
1764 }
1765
1766 // triangle_coords now has:
1767 // (midpoint_x, midpoint_y, midpoint_z,
1768 // v1_x-midpoint_x, v1_y-midpoint_y, v1_z-midpoint_z,
1769 // v2_x-midpoint_x, v2_y-midpoint_y, v2_z-midpoint_z);
1770 for (int quadrature = 0; quadrature<data._qm.getNumberOfQuadraturePoints(); ++quadrature) {
1771 double unscaled_transformed_qp[3] = {0,0,0};
1772 double scaled_transformed_qp[3] = {0,0,0};
1773 for (int j=0; j<data._global_dimensions; ++j) {
1774 for (int k=1; k<3; ++k) { // 3 is for # vertices in subtriangle
1775 // uses vertex-midpoint as one direction
1776 // and other vertex-midpoint as other direction
1777 unscaled_transformed_qp[j] += T(j,k)*data._qm.getSite(quadrature, k-1);
1778 }
1779 // adds back on shift by midpoint
1780 unscaled_transformed_qp[j] += T(j,0);
1781 }
1782
1783 // project onto the sphere
1784 double transformed_qp_norm = 0;
1785 for (int j=0; j<data._global_dimensions; ++j) {
1786 transformed_qp_norm += unscaled_transformed_qp[j]*unscaled_transformed_qp[j];
1787 }
1788 transformed_qp_norm = std::sqrt(transformed_qp_norm);
1789
1790 // project back onto sphere
1791 for (int j=0; j<data._global_dimensions; ++j) {
1792 scaled_transformed_qp[j] = unscaled_transformed_qp[j] * radius / transformed_qp_norm;
1793 }
1794
1795 // u_qp = midpoint + r_qp[1]*(v_1-midpoint) + r_qp[2]*(v_2-midpoint)
1796 // s_qp = u_qp * radius / norm(u_qp)
1797 //
1798 // so G(:,i) is \partial{s_qp}/ \partial{r_qp[i]}
1799 // where r_qp is reference quadrature point (R^2 in 2D manifold in R^3)
1800 //
1801 // G(:,i) = radius * ( \partial{u_qp}/\partial{r_qp[i]} * (\sum_m u_qp[k]^2)^{-1/2}
1802 // + u_qp * \partial{(\sum_m u_qp[k]^2)^{-1/2}}/\partial{r_qp[i]} )
1803 //
1804 // = radius * ( T(:,i)/norm(u_qp) + u_qp*(-1/2)*(\sum_m u_qp[k]^2)^{-3/2}
1805 // *2*(\sum_k u_qp[k]*\partial{u_qp[k]}/\partial{r_qp[i]}) )
1806 //
1807 // = radius * ( T(:,i)/norm(u_qp) + u_qp*(-1/2)*(\sum_m u_qp[k]^2)^{-3/2}
1808 // *2*(\sum_k u_qp[k]*T(k,i)) )
1809 //
1810 // NOTE: we do not multiply G by radius before determining area from vectors,
1811 // so we multiply by radius**2 after calculation
1812 double qp_norm_sq = transformed_qp_norm*transformed_qp_norm;
1813 for (int j=0; j<data._global_dimensions; ++j) {
1814 G(j,1) = T(j,1)/transformed_qp_norm;
1815 G(j,2) = T(j,2)/transformed_qp_norm;
1816 for (int k=0; k<data._global_dimensions; ++k) {
1817 G(j,1) += unscaled_transformed_qp[j]*(-0.5)*std::pow(qp_norm_sq,-1.5)
1818 *2*(unscaled_transformed_qp[k]*T(k,1));
1819 G(j,2) += unscaled_transformed_qp[j]*(-0.5)*std::pow(qp_norm_sq,-1.5)
1820 *2*(unscaled_transformed_qp[k]*T(k,2));
1821 }
1822 }
1823 double G_determinant = getAreaFromVectors(teamMember,
1824 Kokkos::subview(G, Kokkos::ALL(), 1), Kokkos::subview(G, Kokkos::ALL(), 2));
1825 G_determinant *= radius*radius;
1826 XYZ qp = XYZ(scaled_transformed_qp[0], scaled_transformed_qp[1], scaled_transformed_qp[2]);
1827 for (int j=0; j<data._local_dimensions; ++j) {
1828 relative_coord[j] = data._pc.convertGlobalToLocalCoordinate(qp,j,V)
1829 - data._pc.getTargetCoordinate(target_index,j,&V);
1830 // shift quadrature point by target site
1831 }
1832
1833 int k = 0;
1834 for (int n = 0; n <= data._poly_order; n++){
1835 for (alphay = 0; alphay <= n; alphay++){
1836 alphax = n - alphay;
1837 alphaf = factorial[alphax]*factorial[alphay];
1838 double val_to_sum = G_determinant * (data._qm.getWeight(quadrature)
1839 * std::pow(relative_coord.x/cutoff_p,alphax)
1840 * std::pow(relative_coord.y/cutoff_p,alphay) / alphaf);
1841 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1842 if (quadrature==0 && v==0) P_target_row(offset, k) = val_to_sum;
1843 else P_target_row(offset, k) += val_to_sum;
1844 });
1845 k++;
1846 }
1847 }
1848 entire_cell_area += G_determinant * data._qm.getWeight(quadrature);
1849 }
1850 }
1851 if (data._operations(i) == TargetOperation::CellAverageEvaluation) {
1852 int k = 0;
1853 for (int n = 0; n <= data._poly_order; n++){
1854 for (alphay = 0; alphay <= n; alphay++){
1855 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1856 P_target_row(offset, k) /= entire_cell_area;
1857 });
1858 k++;
1859 }
1860 }
1861 }
1862 } else if (data._operations(i) == TargetOperation::FaceNormalAverageEvaluation ||
1863 data._operations(i) == TargetOperation::FaceNormalIntegralEvaluation ||
1864 data._operations(i) == TargetOperation::EdgeTangentAverageEvaluation ||
1865 data._operations(i) == TargetOperation::EdgeTangentIntegralEvaluation) {
1866 compadre_kernel_assert_debug(data._local_dimensions==2 &&
1867 "FaceNormalIntegralSample, EdgeTangentIntegralSample, FaceNormalAverageSample, \
1868 and EdgeTangentAverageSample only support 2d or 3d with 2d manifold");
1869 compadre_kernel_assert_debug(data._qm.getDimensionOfQuadraturePoints()==1 \
1870 && "Only 1d quadrature may be specified for edge integrals");
1871 compadre_kernel_assert_debug(data._qm.getNumberOfQuadraturePoints()>=1 \
1872 && "Quadrature points not generated");
1873 compadre_kernel_assert_debug(data._target_extra_data.extent(0)>0 && "Extra data used but not set.");
1874 const double factorial[15] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
1875
1876 double cutoff_p = data._epsilons(target_index);
1877 int alphax, alphay;
1878 double alphaf;
1879
1880 /*
1881 * requires quadrature points defined on an edge, not a target/source edge (spoke)
1882 *
1883 * data._source_extra_data will contain the endpoints (2 for 2D, 3 for 3D) and then the unit normals
1884 * (e0_x, e0_y, e1_x, e1_y, n_x, n_y, t_x, t_y)
1885 */
1886
1887 int quadrature_point_loop = data._qm.getNumberOfQuadraturePoints();
1888
1889 const int TWO = 2; // used because of # of vertices on an edge
1890 double G_data[3*TWO]; // max(2,3)*TWO
1891 double edge_coords[3*TWO];
1892 for (int j=0; j<data._global_dimensions*TWO; ++j) G_data[j] = 0;
1893 for (int j=0; j<data._global_dimensions*TWO; ++j) edge_coords[j] = 0;
1894 // 2 is for # vertices on an edge
1895 device_unmanaged_matrix_right_type G(G_data, data._global_dimensions, TWO);
1896 device_unmanaged_matrix_right_type edge_coords_matrix(edge_coords, data._global_dimensions, TWO);
1897
1898 // neighbor coordinate is assumed to be midpoint
1899 // could be calculated, but is correct for sphere
1900 // and also for non-manifold problems
1901 // uses given midpoint, rather than computing the midpoint from vertices
1902 double radius = 0.0;
1903 // this midpoint should lie on the sphere, or this will be the wrong radius
1904 for (int j=0; j<data._global_dimensions; ++j) {
1905 edge_coords_matrix(j, 0) = data._target_extra_data(target_index, j);
1906 edge_coords_matrix(j, 1) = data._target_extra_data(target_index, data._global_dimensions + j) - edge_coords_matrix(j, 0);
1907 radius += edge_coords_matrix(j, 0)*edge_coords_matrix(j, 0);
1908 }
1909 radius = std::sqrt(radius);
1910
1911 // edge_coords now has:
1912 // (v0_x, v0_y, v1_x-v0_x, v1_y-v0_y)
1913 auto E = edge_coords_matrix;
1914
1915 // get arc length of edge on manifold
1916 double theta = 0.0;
1917 if (data._problem_type == ProblemType::MANIFOLD) {
1918 XYZ a = {E(0,0), E(1,0), E(2,0)};
1919 XYZ b = {E(0,1)+E(0,0), E(1,1)+E(1,0), E(2,1)+E(2,0)};
1920 double a_dot_b = a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
1921 double norm_a_cross_b = getAreaFromVectors(teamMember, a, b);
1922 theta = std::atan(norm_a_cross_b / a_dot_b);
1923 }
1924
1925 for (int c=0; c<data._local_dimensions; ++c) {
1926 int input_component = (data._sampling_multiplier==1 && data._reconstruction_space_rank==1) ? 0 : c;
1927 int offset = data._d_ss.getTargetOffsetIndex(i, input_component /*in*/, 0 /*out*/, 0/*additional*/);
1928 int column_offset = (data._reconstruction_space_rank==1) ? c*target_NP : 0;
1929 for (int j=0; j<target_NP; ++j) {
1930 P_target_row(offset, column_offset + j) = 0;
1931 }
1932 }
1933
1934 // loop
1935 double entire_edge_length = 0.0;
1936 XYZ relative_coord;
1937 for (int quadrature = 0; quadrature<quadrature_point_loop; ++quadrature) {
1938
1939 double G_determinant = 1.0;
1940 double scaled_transformed_qp[3] = {0,0,0};
1941 double unscaled_transformed_qp[3] = {0,0,0};
1942 for (int j=0; j<data._global_dimensions; ++j) {
1943 unscaled_transformed_qp[j] += E(j,1)*data._qm.getSite(quadrature, 0);
1944 // adds back on shift by endpoint
1945 unscaled_transformed_qp[j] += E(j,0);
1946 }
1947
1948 // project onto the sphere
1949 if (data._problem_type == ProblemType::MANIFOLD) {
1950 // unscaled_transformed_qp now lives on cell edge, but if on manifold,
1951 // not directly on the sphere, just near by
1952
1953 // normalize to project back onto sphere
1954 double transformed_qp_norm = 0;
1955 for (int j=0; j<data._global_dimensions; ++j) {
1956 transformed_qp_norm += unscaled_transformed_qp[j]*unscaled_transformed_qp[j];
1957 }
1958 transformed_qp_norm = std::sqrt(transformed_qp_norm);
1959 // transformed_qp made radius in length
1960 for (int j=0; j<data._global_dimensions; ++j) {
1961 scaled_transformed_qp[j] = unscaled_transformed_qp[j] * radius / transformed_qp_norm;
1962 }
1963
1964 G_determinant = radius * theta;
1965 XYZ qp = XYZ(scaled_transformed_qp[0], scaled_transformed_qp[1], scaled_transformed_qp[2]);
1966 for (int j=0; j<data._local_dimensions; ++j) {
1967 relative_coord[j] = data._pc.convertGlobalToLocalCoordinate(qp,j,V)
1968 - data._pc.getTargetCoordinate(target_index,j,&V);
1969 // shift quadrature point by target site
1970 }
1971 relative_coord[2] = 0;
1972 } else { // not on a manifold, but still integrated
1973 XYZ endpoints_difference = {E(0,1), E(1,1), 0};
1974 G_determinant = data._pc.EuclideanVectorLength(endpoints_difference, 2);
1975 for (int j=0; j<data._local_dimensions; ++j) {
1976 relative_coord[j] = unscaled_transformed_qp[j]
1977 - data._pc.getTargetCoordinate(target_index,j,&V);
1978 // shift quadrature point by target site
1979 }
1980 relative_coord[2] = 0;
1981 }
1982
1983 // get normal or tangent direction (ambient)
1984 XYZ direction;
1985 if (data._operations(i) == TargetOperation::FaceNormalIntegralEvaluation
1986 || data._operations(i) == FaceNormalAverageEvaluation) {
1987 for (int j=0; j<data._global_dimensions; ++j) {
1988 // normal direction
1989 direction[j] = data._target_extra_data(target_index, 2*data._global_dimensions + j);
1990 }
1991 } else {
1992 if (data._problem_type == ProblemType::MANIFOLD) {
1993 // generate tangent from outward normal direction of the sphere and edge normal
1994 XYZ k = {scaled_transformed_qp[0], scaled_transformed_qp[1], scaled_transformed_qp[2]};
1995 //XYZ k = {data._pc.getTargetCoordinate(target_index, 0), data._pc.getTargetCoordinate(target_index, 1),data._pc.getTargetCoordinate(target_index, 2)};
1996 double k_norm = std::sqrt(k[0]*k[0]+k[1]*k[1]+k[2]*k[2]);
1997 k[0] = k[0]/k_norm; k[1] = k[1]/k_norm; k[2] = k[2]/k_norm;
1998 XYZ n = {data._target_extra_data(target_index, 2*data._global_dimensions + 0),
1999 data._target_extra_data(target_index, 2*data._global_dimensions + 1),
2000 data._target_extra_data(target_index, 2*data._global_dimensions + 2)};
2001
2002 double norm_k_cross_n = getAreaFromVectors(teamMember, k, n);
2003 direction[0] = (k[1]*n[2] - k[2]*n[1]) / norm_k_cross_n;
2004 direction[1] = (k[2]*n[0] - k[0]*n[2]) / norm_k_cross_n;
2005 direction[2] = (k[0]*n[1] - k[1]*n[0]) / norm_k_cross_n;
2006 } else {
2007 for (int j=0; j<data._global_dimensions; ++j) {
2008 // tangent direction
2009 direction[j] = data._target_extra_data(target_index, 3*data._global_dimensions + j);
2010 }
2011 }
2012 }
2013
2014 // convert direction to local chart (for manifolds)
2015 XYZ local_direction;
2016 if (data._problem_type == ProblemType::MANIFOLD) {
2017 for (int j=0; j<data._local_dimensions; ++j) {
2018 // Project ambient normal direction onto local chart basis as a local direction.
2019 // Using V alone to provide vectors only gives tangent vectors at
2020 // the target site. This could result in accuracy < 3rd order.
2021
2022 local_direction[j] = data._pc.convertGlobalToLocalCoordinate(direction,j,V);
2023 }
2024 }
2025
2026 // if sampling multiplier is 1 && vector, then vector is [(u_x, u_y)]
2027 // if sampling multiplier is 2 && vector, then vector is [(u_x, 0), (0, u_y)]
2028 // if sampling multiplier is 1 && scalar, then vector is [u_x][u_y]
2029 // if sampling multiplier is 2 && scalar, then vector is [(u_x, 0), (0, u_y)]
2030 for (int c=0; c<data._local_dimensions; ++c) {
2031 int input_component = (data._sampling_multiplier==1 && data._reconstruction_space_rank==1) ? 0 : c;
2032 int offset = data._d_ss.getTargetOffsetIndex(i, input_component, 0 /*out*/, 0/*additional*/);
2033 int column_offset = (data._reconstruction_space_rank==1) ? c*target_NP : 0;
2034 int k = 0;
2035 for (int n = 0; n <= data._poly_order; n++){
2036 for (alphay = 0; alphay <= n; alphay++){
2037 alphax = n - alphay;
2038 alphaf = factorial[alphax]*factorial[alphay];
2039
2040 // local evaluation of vector [0,p] or [p,0]
2041 double v0, v1;
2042 v0 = (c==0) ? std::pow(relative_coord.x/cutoff_p,alphax)
2043 *std::pow(relative_coord.y/cutoff_p,alphay)/alphaf : 0;
2044 v1 = (c==1) ? std::pow(relative_coord.x/cutoff_p,alphax)
2045 *std::pow(relative_coord.y/cutoff_p,alphay)/alphaf : 0;
2046
2047 // either n*v or t*v
2048 double dot_product = 0.0;
2049 if (data._problem_type == ProblemType::MANIFOLD) {
2050 // alternate option for projection
2051 dot_product = local_direction[0]*v0 + local_direction[1]*v1;
2052 } else {
2053 dot_product = direction[0]*v0 + direction[1]*v1;
2054 }
2055
2056 // multiply by quadrature weight
2057 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
2058 if (quadrature==0 && c==0) P_target_row(offset, column_offset + k) =
2059 dot_product * data._qm.getWeight(quadrature) * G_determinant;
2060 else P_target_row(offset, column_offset + k) +=
2061 dot_product * data._qm.getWeight(quadrature) * G_determinant;
2062 });
2063 k++;
2064 }
2065 }
2066 }
2067 entire_edge_length += G_determinant * data._qm.getWeight(quadrature);
2068 } // end of quadrature loop
2069 if (data._operations(i) == TargetOperation::FaceNormalAverageEvaluation
2070 || data._operations(i) == TargetOperation::EdgeTangentAverageEvaluation) {
2071 for (int c=0; c<data._local_dimensions; ++c) {
2072 int k = 0;
2073 //int offset = data._d_ss.getTargetOffsetIndex(i, c, 0, 0);
2074 //int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
2075 //int input_component = (data._sampling_multiplier==1) ? 0 : c;
2076 int input_component = (data._sampling_multiplier==1 && data._reconstruction_space_rank==1) ? 0 : c;
2077 //int input_component = c;
2078 int offset = data._d_ss.getTargetOffsetIndex(i, input_component, 0 /*out*/, 0/*additional*/);
2079 int column_offset = (data._reconstruction_space_rank == 1) ? c*target_NP : 0;
2080 for (int n = 0; n <= data._poly_order; n++){
2081 for (alphay = 0; alphay <= n; alphay++){
2082 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
2083 P_target_row(offset, column_offset + k) /= entire_edge_length;
2084 });
2085 k++;
2086 }
2087 }
2088 }
2089 }
2090 } else {
2091 compadre_kernel_assert_release((false) && "Functionality not yet available.");
2092 }
2093 } else if (data._dimensions==2) { // 1D manifold in 2D problem
2094 if (data._operations(i) == TargetOperation::ScalarPointEvaluation) {
2095 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (const int j) {
2096 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, data._dimensions-1, data._poly_order, false /*bool on only specific order*/, &V, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
2097 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
2098 for (int k=0; k<target_NP; ++k) {
2099 P_target_row(offset, k) = delta(k);
2100 }
2101 });
2102 additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
2103 } else {
2104 compadre_kernel_assert_release((false) && "Functionality not yet available.");
2105 }
2106 } else {
2107 compadre_kernel_assert_release((false) && "Functionality not yet available.");
2108 }
2109 (void)additional_evaluation_sites_handled; // handle unused variable warning/error
2110 (void)additional_evaluation_sites_need_handled; // handle unused variable warning/error
2111 compadre_kernel_assert_release(((additional_evaluation_sites_need_handled && additional_evaluation_sites_handled) || (!additional_evaluation_sites_need_handled)) && "Auxiliary evaluation coordinates are specified by user, but are calling a target operation that can not handle evaluating additional sites.");
2112 } // !operation_handled
2113 }
2114 teamMember.team_barrier();
2115}
2116
2117} // Compadre
2118#endif
Compadre_GMLS.hpp
Compadre_Manifold_Functions.hpp
compadre_kernel_assert_debug
#define compadre_kernel_assert_debug(condition)
Definition Compadre_Typedefs.hpp:195
compadre_kernel_assert_release
#define compadre_kernel_assert_release(condition)
compadre_kernel_assert_release is similar to compadre_assert_release, but is a call on the device,...
Definition Compadre_Typedefs.hpp:175
Compadre_USER_ManifoldTargetFunctionals.hpp
if
if(some_conditions_for_a_user_defined_operation)
Definition Compadre_USER_ManifoldTargetFunctionals.hpp:18
Compadre_USER_StandardTargetFunctionals.hpp
Compadre::GMLS::getNP
static KOKKOS_INLINE_FUNCTION int getNP(const int m, const int dimension=3, const ReconstructionSpace r_space=ReconstructionSpace::ScalarTaylorPolynomial)
Returns size of the basis for a given polynomial order and dimension General to dimension 1....
Definition Compadre_GMLS.hpp:512
Compadre::SurfaceCurlOfScalar
KOKKOS_INLINE_FUNCTION double SurfaceCurlOfScalar(const device_unmanaged_vector_type a_, const double h, const double u1, const double u2, int x_pow, int y_pow, const int component)
Surface curl at any point in the local chart.
Definition Compadre_Manifold_Functions.hpp:104
Compadre::computeTargetFunctionalsOnManifold
KOKKOS_INLINE_FUNCTION void computeTargetFunctionalsOnManifold(const TargetData &data, const member_type &teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, device_unmanaged_matrix_right_type P_target_row, device_unmanaged_matrix_right_type V, device_unmanaged_vector_type curvature_coefficients)
Evaluates a polynomial basis with a target functional applied, using information from the manifold cu...
Definition Compadre_Targets.hpp:1126
Compadre::StaggeredEdgeIntegralSample
constexpr SamplingFunctional StaggeredEdgeIntegralSample
Samples consist of the result of integrals of a vector dotted with the tangent along edges between ne...
Definition Compadre_Operators.hpp:193
Compadre::PointSample
constexpr SamplingFunctional PointSample
Available sampling functionals.
Definition Compadre_Operators.hpp:180
Compadre::member_type
team_policy::member_type member_type
Definition Compadre_Typedefs.hpp:61
Compadre::scratch_vector_type
Kokkos::View< double *, device_scratch, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_vector_type
Definition Compadre_Typedefs.hpp:76
Compadre::DivergenceFreeVectorTaylorPolynomial
@ DivergenceFreeVectorTaylorPolynomial
Divergence-free vector polynomial basis.
Definition Compadre_Operators.hpp:114
Compadre::VectorTaylorPolynomial
@ VectorTaylorPolynomial
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...
Definition Compadre_Operators.hpp:109
Compadre::ScalarTaylorPolynomial
@ ScalarTaylorPolynomial
Scalar polynomial basis centered at the target site and scaled by sum of basis powers e....
Definition Compadre_Operators.hpp:107
Compadre::VectorOfScalarClonesTaylorPolynomial
@ VectorOfScalarClonesTaylorPolynomial
Scalar basis reused as many times as there are components in the vector resulting in a much cheaper p...
Definition Compadre_Operators.hpp:112
Compadre::StaggeredEdgeAnalyticGradientIntegralSample
constexpr SamplingFunctional StaggeredEdgeAnalyticGradientIntegralSample
Analytical integral of a gradient source vector is just a difference of the scalar source at neighbor...
Definition Compadre_Operators.hpp:190
Compadre::computeTargetFunctionals
KOKKOS_INLINE_FUNCTION void computeTargetFunctionals(const TargetData &data, const member_type &teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, device_unmanaged_matrix_right_type P_target_row)
Evaluates a polynomial basis with a target functional applied to each member of the basis.
Definition Compadre_Targets.hpp:26
Compadre::GaussianCurvature
KOKKOS_INLINE_FUNCTION double GaussianCurvature(const device_unmanaged_vector_type a_, const double h, const double u1, const double u2)
Gaussian curvature K at any point in the local chart.
Definition Compadre_Manifold_Functions.hpp:58
Compadre::VectorPointSample
constexpr SamplingFunctional VectorPointSample
Point evaluations of the entire vector source function.
Definition Compadre_Operators.hpp:183
Compadre::computeCurvatureFunctionals
KOKKOS_INLINE_FUNCTION void computeCurvatureFunctionals(const TargetData &data, const member_type &teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, device_unmanaged_matrix_right_type P_target_row, const device_unmanaged_matrix_right_type *V, const local_index_type local_neighbor_index=-1)
Evaluates a polynomial basis for the curvature with a gradient target functional applied.
Definition Compadre_Targets.hpp:1057
Compadre::LaplacianOfScalarPointEvaluation
@ LaplacianOfScalarPointEvaluation
Point evaluation of the laplacian of a scalar (could be on a manifold or not)
Definition Compadre_Operators.hpp:26
Compadre::GradientOfScalarPointEvaluation
@ GradientOfScalarPointEvaluation
Point evaluation of the gradient of a scalar.
Definition Compadre_Operators.hpp:30
Compadre::CurlOfVectorPointEvaluation
@ CurlOfVectorPointEvaluation
Point evaluation of the curl of a vector (results in a vector)
Definition Compadre_Operators.hpp:36
Compadre::PartialYOfScalarPointEvaluation
@ PartialYOfScalarPointEvaluation
Point evaluation of the partial with respect to y of a scalar.
Definition Compadre_Operators.hpp:42
Compadre::CurlCurlOfVectorPointEvaluation
@ CurlCurlOfVectorPointEvaluation
Point evaluation of the curl of a curl of a vector (results in a vector)
Definition Compadre_Operators.hpp:38
Compadre::GradientOfVectorPointEvaluation
@ GradientOfVectorPointEvaluation
Point evaluation of the gradient of a vector (results in a matrix, NOT CURRENTLY IMPLEMENTED)
Definition Compadre_Operators.hpp:32
Compadre::PartialZOfScalarPointEvaluation
@ PartialZOfScalarPointEvaluation
Point evaluation of the partial with respect to z of a scalar.
Definition Compadre_Operators.hpp:44
Compadre::DivergenceOfVectorPointEvaluation
@ DivergenceOfVectorPointEvaluation
Point evaluation of the divergence of a vector (results in a scalar)
Definition Compadre_Operators.hpp:34
Compadre::CellIntegralEvaluation
@ CellIntegralEvaluation
Integral values over cell using quadrature Supported on 2D faces in 3D problems (manifold) and 2D cel...
Definition Compadre_Operators.hpp:55
Compadre::CellAverageEvaluation
@ CellAverageEvaluation
Average values of a cell using quadrature Supported on 2D faces in 3D problems (manifold) and 2D cell...
Definition Compadre_Operators.hpp:52
Compadre::VectorPointEvaluation
@ VectorPointEvaluation
Point evaluation of a vector (reconstructs entire vector at once, requiring a ReconstructionSpace hav...
Definition Compadre_Operators.hpp:24
Compadre::ScalarPointEvaluation
@ ScalarPointEvaluation
Point evaluation of a scalar.
Definition Compadre_Operators.hpp:21
Compadre::PartialXOfScalarPointEvaluation
@ PartialXOfScalarPointEvaluation
Point evaluation of the partial with respect to x of a scalar.
Definition Compadre_Operators.hpp:40
Compadre::device_unmanaged_matrix_right_type
Kokkos::View< double **, layout_right, device_execution_space, Kokkos::MemoryTraits< Kokkos::Unmanaged > > device_unmanaged_matrix_right_type
Definition Compadre_Typedefs.hpp:92
Compadre::device_unmanaged_vector_type
Kokkos::View< double *, device_execution_space, Kokkos::MemoryTraits< Kokkos::Unmanaged > > device_unmanaged_vector_type
Definition Compadre_Typedefs.hpp:96
Compadre::getAreaFromVectors
KOKKOS_INLINE_FUNCTION double getAreaFromVectors(const member_type &teamMember, view_type_1 v1, view_type_2 v2)
Definition Compadre_Utilities.hpp:32
Generated by doxygen 1.9.8