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