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