1 #ifndef _COMPADRE_TARGETS_HPP_
2 #define _COMPADRE_TARGETS_HPP_
16 template <
typename TargetData>
17 KOKKOS_INLINE_FUNCTION
21 if (data._dimensions > 1) {
24 || data._data_sampling_multiplier!=0
26 && data._polynomial_sampling_functional==
PointSample))
27 &&
"data._reconstruction_space(VectorOfScalar clones incompatible with scalar output sampling functional which is not a PointSample");
31 bool additional_evaluation_sites_need_handled =
32 (data._additional_pc._source_coordinates.extent(0) > 0) ?
true :
false;
34 const int target_index = data._initial_index_for_batch + teamMember.league_rank();
36 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, P_target_row.extent(0)), [&] (
const int j) {
37 Kokkos::parallel_for(Kokkos::ThreadVectorRange(teamMember, P_target_row.extent(1)),
39 P_target_row(j,k) = 0;
42 for (
int j = 0; j < delta.extent_int(0); ++j) {
45 for (
int j = 0; j < thread_workspace.extent_int(0); ++j) {
46 thread_workspace(j) = 0;
48 teamMember.team_barrier();
51 int target_NP =
GMLS::getNP(data._poly_order, data._dimensions, data._reconstruction_space);
52 const int num_evaluation_sites = data._additional_pc._nla.getNumberOfNeighborsDevice(target_index) + 1;
54 for (
size_t i=0; i<data._operations.size(); ++i) {
56 bool additional_evaluation_sites_handled =
false;
58 bool operation_handled =
true;
65 if (!operation_handled) {
74 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
75 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, j);
76 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
77 for (
int k=0; k<target_NP; ++k) {
78 P_target_row(offset, k) = delta(k);
81 additional_evaluation_sites_handled =
true;
83 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
84 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
85 switch (data._dimensions) {
87 P_target_row(offset, 4) = std::pow(data._epsilons(target_index), -2);
88 P_target_row(offset, 6) = std::pow(data._epsilons(target_index), -2);
89 P_target_row(offset, 9) = std::pow(data._epsilons(target_index), -2);
92 P_target_row(offset, 3) = std::pow(data._epsilons(target_index), -2);
93 P_target_row(offset, 5) = std::pow(data._epsilons(target_index), -2);
96 P_target_row(offset, 2) = std::pow(data._epsilons(target_index), -2);
100 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
101 for (
int d=0;
d<data._dimensions; ++
d) {
102 int offset = data._d_ss.getTargetOffsetIndex(i, 0,
d, j);
103 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
104 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 , 1 ,
d , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, j);
107 additional_evaluation_sites_handled =
true;
109 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
110 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
111 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
112 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 , 1 , 0 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, j);
114 additional_evaluation_sites_handled =
true;
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 , 1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, j);
122 additional_evaluation_sites_handled =
true;
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 , 1 , 2 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, j);
130 additional_evaluation_sites_handled =
true;
136 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
137 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
138 switch (data._dimensions) {
140 P_target_row(offset, 4) = std::pow(data._epsilons(target_index), -2);
141 P_target_row(offset, 6) = std::pow(data._epsilons(target_index), -2);
142 P_target_row(offset, 9) = std::pow(data._epsilons(target_index), -2);
145 P_target_row(offset, 3) = std::pow(data._epsilons(target_index), -2);
146 P_target_row(offset, 5) = std::pow(data._epsilons(target_index), -2);
149 P_target_row(offset, 2) = std::pow(data._epsilons(target_index), -2);
155 "CellAverageEvaluation and CellIntegralEvaluation only support 2d or 3d with 2d manifold");
156 const double factorial[15] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
157 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
160 double cutoff_p = data._epsilons(target_index);
164 double triangle_coords[3*3];
167 for (
int j=0; j<data._global_dimensions; ++j) {
169 triangle_coords_matrix(j, 0) = data._pc.getTargetCoordinate(target_index, j);
172 size_t num_vertices = (data._target_extra_data(target_index, data._target_extra_data.extent(1)-1)
173 != data._target_extra_data(target_index, data._target_extra_data.extent(1)-1))
174 ? (data._target_extra_data.extent(1) / data._global_dimensions) - 1 :
175 (data._target_extra_data.extent(1) / data._global_dimensions);
179 double entire_cell_area = 0.0;
180 for (
size_t v=0; v<num_vertices; ++v) {
182 int v2 = (v+1) % num_vertices;
184 for (
int j=0; j<data._global_dimensions; ++j) {
185 triangle_coords_matrix(j,1) = data._target_extra_data(target_index, v1*data._global_dimensions+j)
186 - triangle_coords_matrix(j,0);
187 triangle_coords_matrix(j,2) = data._target_extra_data(target_index, v2*data._global_dimensions+j)
188 - triangle_coords_matrix(j,0);
194 auto T=triangle_coords_matrix;
195 for (
int quadrature = 0; quadrature<data._qm.getNumberOfQuadraturePoints(); ++quadrature) {
196 double transformed_qp[2] = {0,0};
197 for (
int j=0; j<data._global_dimensions; ++j) {
198 for (
int k=1; k<3; ++k) {
199 transformed_qp[j] += T(j,k)*data._qm.getSite(quadrature, k-1);
201 transformed_qp[j] += T(j,0);
206 Kokkos::subview(T, Kokkos::ALL(), 1), Kokkos::subview(T, Kokkos::ALL(), 2));
208 for (
int j=0; j<data._global_dimensions; ++j) {
209 relative_coord[j] = transformed_qp[j] - data._pc.getTargetCoordinate(target_index, j);
213 for (
int n = 0; n <= data._poly_order; n++){
214 for (alphay = 0; alphay <= n; alphay++){
216 alphaf = factorial[alphax]*factorial[alphay];
217 double val_to_sum = G_determinant * ( data._qm.getWeight(quadrature)
218 * std::pow(relative_coord.
x/cutoff_p,alphax)
219 * std::pow(relative_coord.
y/cutoff_p,alphay)/alphaf);
220 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
221 if (quadrature==0 && v==0) P_target_row(offset, k) = val_to_sum;
222 else P_target_row(offset, k) += val_to_sum;
227 entire_cell_area += G_determinant * data._qm.getWeight(quadrature);
232 for (
int n = 0; n <= data._poly_order; n++){
233 for (alphay = 0; alphay <= n; alphay++){
234 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
235 P_target_row(offset, k) /= entire_cell_area;
257 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
258 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, j);
259 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
260 for (
int k=0; k<target_NP; ++k) {
261 P_target_row(offset, k) = delta(k);
264 additional_evaluation_sites_handled =
true;
266 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int e) {
267 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, e);
268 for (
int m=0; m<data._sampling_multiplier; ++m) {
269 int output_components = data._basis_multiplier;
270 for (
int c=0; c<output_components; ++c) {
271 int offset = data._d_ss.getTargetOffsetIndex(i, m , c , e);
275 for (
int j=0; j<target_NP; ++j) {
276 P_target_row(offset, c*target_NP + j) = delta(j);
281 additional_evaluation_sites_handled =
true;
285 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int e) {
286 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, e);
287 for (
int m=0; m<data._sampling_multiplier; ++m) {
288 int output_components = data._basis_multiplier;
289 for (
int c=0; c<output_components; ++c) {
290 int offset = data._d_ss.getTargetOffsetIndex(i, m , c , e);
294 for (
int j=0; j<target_NP; ++j) {
295 P_target_row(offset, c*target_NP + j) = delta(j);
300 additional_evaluation_sites_handled =
true;
302 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int e) {
303 int output_components = data._basis_multiplier;
304 for (
int m2=0; m2<output_components; ++m2) {
305 int offset = data._d_ss.getTargetOffsetIndex(i, 0 , m2 , e);
306 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 , 1 , m2 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, e);
307 for (
int j=0; j<target_NP; ++j) {
308 P_target_row(offset, j) = delta(j);
312 additional_evaluation_sites_handled =
true;
314 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int e) {
315 for (
int m0=0; m0<data._sampling_multiplier; ++m0) {
316 int output_components = data._basis_multiplier;
317 for (
int m1=0; m1<output_components; ++m1) {
318 for (
int m2=0; m2<output_components; ++m2) {
319 int offset = data._d_ss.getTargetOffsetIndex(i, m0 , m1*output_components+m2 , e);
320 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 , 1 , m2 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, e);
321 for (
int j=0; j<target_NP; ++j) {
322 P_target_row(offset, m1*target_NP + j) = delta(j);
328 additional_evaluation_sites_handled =
true;
331 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
332 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);;
333 switch (data._dimensions) {
335 P_target_row(offset, 1) = std::pow(data._epsilons(target_index), -1);
336 P_target_row(offset, target_NP + 2) = std::pow(data._epsilons(target_index), -1);
337 P_target_row(offset, 2*target_NP + 3) = std::pow(data._epsilons(target_index), -1);
340 P_target_row(offset, 1) = std::pow(data._epsilons(target_index), -1);
341 P_target_row(offset, target_NP + 2) = std::pow(data._epsilons(target_index), -1);
344 P_target_row(offset, 1) = std::pow(data._epsilons(target_index), -1);
348 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int e) {
349 for (
int m=0; m<data._sampling_multiplier; ++m) {
350 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 , 1 , m , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, e);
351 int offset = data._d_ss.getTargetOffsetIndex(i, m , 0 , e);
352 for (
int j=0; j<target_NP; ++j) {
353 P_target_row(offset, m*target_NP + j) = delta(j);
357 additional_evaluation_sites_handled =
true;
360 if (data._dimensions==3) {
361 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
365 int offset = data._d_ss.getTargetOffsetIndex(i, 0 , 0 , 0);
369 offset = data._d_ss.getTargetOffsetIndex(i, 1 , 0 , 0);
372 P_target_row(offset, target_NP + 3) = -std::pow(data._epsilons(target_index), -1);
374 offset = data._d_ss.getTargetOffsetIndex(i, 2 , 0 , 0);
377 P_target_row(offset, 2*target_NP + 2) = std::pow(data._epsilons(target_index), -1);
383 int offset = data._d_ss.getTargetOffsetIndex(i, 0 , 1 , 0);
386 P_target_row(offset, 3) = std::pow(data._epsilons(target_index), -1);
392 offset = data._d_ss.getTargetOffsetIndex(i, 2 , 1 , 0);
395 P_target_row(offset, 2*target_NP + 1) = -std::pow(data._epsilons(target_index), -1);
401 int offset = data._d_ss.getTargetOffsetIndex(i, 0 , 2 , 0);
404 P_target_row(offset, 2) = -std::pow(data._epsilons(target_index), -1);
406 offset = data._d_ss.getTargetOffsetIndex(i, 1 , 2 , 0);
409 P_target_row(offset, target_NP + 1) = std::pow(data._epsilons(target_index), -1);
416 }
else if (data._dimensions==2) {
417 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
421 int offset = data._d_ss.getTargetOffsetIndex(i, 0 , 0 , 0);
425 offset = data._d_ss.getTargetOffsetIndex(i, 1 , 0 , 0);
428 P_target_row(offset, target_NP + 2) = std::pow(data._epsilons(target_index), -1);
434 int offset = data._d_ss.getTargetOffsetIndex(i, 0 , 1 , 0);
437 P_target_row(offset, 1) = -std::pow(data._epsilons(target_index), -1);
461 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
462 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, j);
463 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
464 for (
int k=0; k<target_NP; ++k) {
465 P_target_row(offset, k) = delta(k);
468 additional_evaluation_sites_handled =
true;
470 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int e) {
471 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, e);
472 for (
int m=0; m<data._sampling_multiplier; ++m) {
473 for (
int c=0; c<data._data_sampling_multiplier; ++c) {
474 int offset = data._d_ss.getTargetOffsetIndex(i, c , c , e);
475 for (
int j=0; j<target_NP; ++j) {
476 P_target_row(offset, j) = delta(j);
481 additional_evaluation_sites_handled =
true;
484 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
485 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
486 switch (data._dimensions) {
488 P_target_row(offset, 4) = std::pow(data._epsilons(target_index), -2);
489 P_target_row(offset, 6) = std::pow(data._epsilons(target_index), -2);
490 P_target_row(offset, 9) = std::pow(data._epsilons(target_index), -2);
493 P_target_row(offset, 3) = std::pow(data._epsilons(target_index), -2);
494 P_target_row(offset, 5) = std::pow(data._epsilons(target_index), -2);
497 P_target_row(offset, 2) = std::pow(data._epsilons(target_index), -2);
502 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
503 for (
int d=0;
d<data._dimensions; ++
d) {
504 int offset = data._d_ss.getTargetOffsetIndex(i, 0,
d, j);
505 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
506 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 , 1 ,
d , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, j);
509 additional_evaluation_sites_handled =
true;
511 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
512 for (
int m0=0; m0<data._dimensions; ++m0) {
513 for (
int m1=0; m1<data._dimensions; ++m1) {
514 for (
int m2=0; m2<data._dimensions; ++m2) {
515 int offset = data._d_ss.getTargetOffsetIndex(i, m0 , m1*data._dimensions+m2 , j);
516 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
517 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 , 1 , m2 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, j);
522 additional_evaluation_sites_handled =
true;
525 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
526 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
527 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
528 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 , 1 , 0 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, j);
530 additional_evaluation_sites_handled =
true;
533 if (data._dimensions>1) {
534 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
535 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
536 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
537 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 , 1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, j);
539 additional_evaluation_sites_handled =
true;
543 if (data._dimensions>2) {
544 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
545 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
546 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
547 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 , 1 , 2 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, j);
549 additional_evaluation_sites_handled =
true;
552 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
553 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
554 for (
int j=0; j<target_NP; ++j) {
555 P_target_row(offset, j) = 0;
558 P_target_row(offset, 1) = std::pow(data._epsilons(target_index), -1);
560 if (data._dimensions>1) {
561 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, 0);
562 for (
int j=0; j<target_NP; ++j) {
563 P_target_row(offset, j) = 0;
565 P_target_row(offset, 2) = std::pow(data._epsilons(target_index), -1);
568 if (data._dimensions>2) {
569 offset = data._d_ss.getTargetOffsetIndex(i, 2, 0, 0);
570 for (
int j=0; j<target_NP; ++j) {
571 P_target_row(offset, j) = 0;
573 P_target_row(offset, 3) = std::pow(data._epsilons(target_index), -1);
580 if (data._dimensions==3) {
581 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
585 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
586 for (
int j=0; j<target_NP; ++j) {
587 P_target_row(offset, j) = 0;
592 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, 0);
593 for (
int j=0; j<target_NP; ++j) {
594 P_target_row(offset, j) = 0;
598 P_target_row(offset, 3) = -std::pow(data._epsilons(target_index), -1);
600 offset = data._d_ss.getTargetOffsetIndex(i, 2, 0, 0);
601 for (
int j=0; j<target_NP; ++j) {
602 P_target_row(offset, j) = 0;
606 P_target_row(offset, 2) = std::pow(data._epsilons(target_index), -1);
612 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, 0);
613 for (
int j=0; j<target_NP; ++j) {
614 P_target_row(offset, j) = 0;
618 P_target_row(offset, 3) = std::pow(data._epsilons(target_index), -1);
620 offset = data._d_ss.getTargetOffsetIndex(i, 1, 1, 0);
621 for (
int j=0; j<target_NP; ++j) {
622 P_target_row(offset, j) = 0;
627 offset = data._d_ss.getTargetOffsetIndex(i, 2, 1, 0);
628 for (
int j=0; j<target_NP; ++j) {
629 P_target_row(offset, j) = 0;
633 P_target_row(offset, 1) = -std::pow(data._epsilons(target_index), -1);
639 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 2, 0);
640 for (
int j=0; j<target_NP; ++j) {
641 P_target_row(offset, j) = 0;
645 P_target_row(offset, 2) = -std::pow(data._epsilons(target_index), -1);
647 offset = data._d_ss.getTargetOffsetIndex(i, 1, 2, 0);
648 for (
int j=0; j<target_NP; ++j) {
649 P_target_row(offset, j) = 0;
653 P_target_row(offset, 1) = std::pow(data._epsilons(target_index), -1);
655 offset = data._d_ss.getTargetOffsetIndex(i, 2, 2, 0);
656 for (
int j=0; j<target_NP; ++j) {
657 P_target_row(offset, j) = 0;
663 }
else if (data._dimensions==2) {
664 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
668 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
669 for (
int j=0; j<target_NP; ++j) {
670 P_target_row(offset, j) = 0;
675 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, 0);
676 for (
int j=0; j<target_NP; ++j) {
677 P_target_row(offset, j) = 0;
681 P_target_row(offset, 2) = std::pow(data._epsilons(target_index), -1);
687 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, 0);
688 for (
int j=0; j<target_NP; ++j) {
689 P_target_row(offset, j) = 0;
693 P_target_row(offset, 1) = -std::pow(data._epsilons(target_index), -1);
695 offset = data._d_ss.getTargetOffsetIndex(i, 1, 1, 0);
696 for (
int j=0; j<target_NP; ++j) {
697 P_target_row(offset, j) = 0;
720 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int e) {
721 for (
int m0=0; m0<data._sampling_multiplier; ++m0) {
722 for (
int m1=0; m1<data._sampling_multiplier; ++m1) {
724 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(m1+1) , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
726 int offset = data._d_ss.getTargetOffsetIndex(i, m0 , m1 , e );
727 for (
int j=0; j<target_NP; ++j) {
728 P_target_row(offset, j) = delta(j);
733 additional_evaluation_sites_handled =
true;
735 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int e) {
736 for (
int m0=0; m0<data._sampling_multiplier; ++m0) {
737 for (
int m1=0; m1<data._sampling_multiplier; ++m1) {
738 int offset = data._d_ss.getTargetOffsetIndex(i, m0 , m1 , e );
739 if (data._dimensions==3) {
743 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) , 1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
745 for (
int j=0; j<target_NP; ++j) {
746 P_target_row(offset, j) = delta(j);
748 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) , 1 , 2 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
750 for (
int j=0; j<target_NP; ++j) {
751 P_target_row(offset, j) -= delta(j);
757 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) , 1 , 0 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
759 for (
int j=0; j<target_NP; ++j) {
760 P_target_row(offset, j) = -delta(j);
762 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) , 1 , 2 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
764 for (
int j=0; j<target_NP; ++j) {
765 P_target_row(offset, j) += delta(j);
771 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) , 1 , 0 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
773 for (
int j=0; j<target_NP; ++j) {
774 P_target_row(offset, j) = delta(j);
776 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) , 1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
778 for (
int j=0; j<target_NP; ++j) {
779 P_target_row(offset, j) -= delta(j);
787 P_target_row(offset, 2) = -std::pow(data._epsilons(target_index), -1);
788 P_target_row(offset, 3) = std::pow(data._epsilons(target_index), -1);
790 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) , 1 , 0 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
792 for (
int j=0; j<target_NP; ++j) {
793 P_target_row(offset, j) = delta(j);
795 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) , 1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
797 for (
int j=0; j<target_NP; ++j) {
798 P_target_row(offset, j) -= delta(j);
806 additional_evaluation_sites_handled =
true;
808 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int e) {
809 for (
int m0=0; m0<data._sampling_multiplier; ++m0) {
810 for (
int m1=0; m1<data._sampling_multiplier; ++m1) {
811 int offset = data._d_ss.getTargetOffsetIndex(i, m0 , m1 , e );
812 if (data._dimensions == 3) {
817 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) , 1 , 0 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
819 for (
int j=0; j<target_NP; ++j) {
820 P_target_row(offset, j) = delta(j);
822 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) , 1 , 1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
824 for (
int j=0; j<target_NP; ++j) {
825 P_target_row(offset, j) -= delta(j);
827 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) , 1 , 0 , 2 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
829 for (
int j=0; j<target_NP; ++j) {
830 P_target_row(offset, j) += delta(j);
832 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) , 1 , 2 , 2 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
834 for (
int j=0; j<target_NP; ++j) {
835 P_target_row(offset, j) -= delta(j);
841 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) , 1 , 0 , 0 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
843 for (
int j=0; j<target_NP; ++j) {
844 P_target_row(offset, j) = -delta(j);
846 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) , 1 , 1 , 0 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
848 for (
int j=0; j<target_NP; ++j) {
849 P_target_row(offset, j) += delta(j);
851 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) , 1 , 1 , 2 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
853 for (
int j=0; j<target_NP; ++j) {
854 P_target_row(offset, j) += delta(j);
856 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) , 1 , 2 , 2 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
858 for (
int j=0; j<target_NP; ++j) {
859 P_target_row(offset, j) -= delta(j);
865 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) , 1 , 0 , 0 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
867 for (
int j=0; j<target_NP; ++j) {
868 P_target_row(offset, j) = -delta(j);
870 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) , 1 , 2 , 0 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
872 for (
int j=0; j<target_NP; ++j) {
873 P_target_row(offset, j) += delta(j);
875 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) , 1 , 1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
877 for (
int j=0; j<target_NP; ++j) {
878 P_target_row(offset, j) -= delta(j);
880 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) , 1 , 2 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
882 for (
int j=0; j<target_NP; ++j) {
883 P_target_row(offset, j) += delta(j);
889 if (data._dimensions == 2) {
894 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) , 1 , 0 , 0 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
896 for (
int j=0; j<target_NP; ++j) {
897 P_target_row(offset, j) = delta(j);
899 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) , 1 , 1 , 0 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
901 for (
int j=0; j<target_NP; ++j) {
902 P_target_row(offset, j) += delta(j);
904 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) , 1 , 0 , 0 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
906 for (
int j=0; j<target_NP; ++j) {
907 P_target_row(offset, j) -= delta(j);
909 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) , 1 , 1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
911 for (
int j=0; j<target_NP; ++j) {
912 P_target_row(offset, j) -= delta(j);
918 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) , 1 , 0 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
920 for (
int j=0; j<target_NP; ++j) {
921 P_target_row(offset, j) = delta(j);
923 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) , 1 , 1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
925 for (
int j=0; j<target_NP; ++j) {
926 P_target_row(offset, j) += delta(j);
928 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) , 1 , 0 , 0 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
930 for (
int j=0; j<target_NP; ++j) {
931 P_target_row(offset, j) -= delta(j);
933 calcHessianPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) , 1 , 1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
935 for (
int j=0; j<target_NP; ++j) {
936 P_target_row(offset, j) -= delta(j);
945 additional_evaluation_sites_handled =
true;
947 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int e) {
948 for (
int m0=0; m0<data._sampling_multiplier; ++m0) {
949 for (
int m1=0; m1<data._sampling_multiplier; ++m1) {
950 for (
int m2=0; m2<data._sampling_multiplier; ++m2) {
951 int offset = data._d_ss.getTargetOffsetIndex(i, m0 , m1*data._sampling_multiplier+m2 , e);
952 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -(m1+1) , 1 , m2 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial,
VectorPointSample, e);
953 for (
int j=0; j<target_NP; ++j) {
954 P_target_row(offset, j) = delta(j);
960 additional_evaluation_sites_handled =
true;
976 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
977 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::BernsteinPolynomial,
PointSample, j);
978 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
979 for (
int k=0; k<target_NP; ++k) {
980 P_target_row(offset, k) = delta(k);
983 additional_evaluation_sites_handled =
true;
985 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
986 for (
int d=0;
d<data._dimensions; ++
d) {
987 int offset = data._d_ss.getTargetOffsetIndex(i, 0,
d, j);
988 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
989 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 , 1 ,
d , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::BernsteinPolynomial,
PointSample, j);
992 additional_evaluation_sites_handled =
true;
994 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
995 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
996 auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
997 calcGradientPij<TargetData>(data, teamMember, row.data(), thread_workspace.data(), target_index, -1 , 1 , 0 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::BernsteinPolynomial,
PointSample, j);
999 additional_evaluation_sites_handled =
true;
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 , 1 , 1 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::BernsteinPolynomial,
PointSample, j);
1007 additional_evaluation_sites_handled =
true;
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 , 1 , 2 , data._dimensions, data._poly_order,
false , NULL ,
ReconstructionSpace::BernsteinPolynomial,
PointSample, j);
1015 additional_evaluation_sites_handled =
true;
1028 compadre_kernel_assert_release(((additional_evaluation_sites_need_handled && additional_evaluation_sites_handled) || (!additional_evaluation_sites_need_handled)) &&
"Auxiliary evaluation coordinates are specified by user, but are calling a target operation that can not handle evaluating additional sites.");
1031 teamMember.team_barrier();
1045 template <
typename TargetData>
1046 KOKKOS_INLINE_FUNCTION
1049 compadre_kernel_assert_release((thread_workspace.extent_int(0)>=(data._curvature_poly_order+1)*data._local_dimensions) &&
"Workspace thread_workspace not large enough.");
1051 const int target_index = data._initial_index_for_batch + teamMember.league_rank();
1053 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, P_target_row.extent(0)), [&] (
const int j) {
1054 Kokkos::parallel_for(Kokkos::ThreadVectorRange(teamMember, P_target_row.extent(1)),
1055 [&] (
const int& k) {
1056 P_target_row(j,k) = 0;
1059 for (
int j = 0; j < delta.extent_int(0); ++j) {
1062 for (
int j = 0; j < thread_workspace.extent_int(0); ++j) {
1063 thread_workspace(j) = 0;
1065 teamMember.team_barrier();
1069 for (
size_t i=0; i<data._curvature_support_operations.size(); ++i) {
1071 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1072 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1073 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, local_neighbor_index, 0 , data._dimensions-1, data._curvature_poly_order,
false , V,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample);
1074 for (
int j=0; j<manifold_NP; ++j) {
1075 P_target_row(offset, j) = delta(j);
1079 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1081 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1082 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, local_neighbor_index, 0 , 0 , data._dimensions-1, data._curvature_poly_order,
false , V,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample);
1083 for (
int j=0; j<manifold_NP; ++j) {
1084 P_target_row(offset, j) = delta(j);
1086 if (data._dimensions>2) {
1088 offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, 0);
1089 calcGradientPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, local_neighbor_index, 0 , 1 , data._dimensions-1, data._curvature_poly_order,
false , V,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample);
1090 for (
int j=0; j<manifold_NP; ++j) {
1091 P_target_row(offset, j) = delta(j);
1099 teamMember.team_barrier();
1114 template <
typename TargetData>
1115 KOKKOS_INLINE_FUNCTION
1118 compadre_kernel_assert_release((thread_workspace.extent_int(0)>=(data._poly_order+1)*data._local_dimensions) &&
"Workspace thread_workspace not large enough.");
1121 const int target_index = data._initial_index_for_batch + teamMember.league_rank();
1123 int target_NP =
GMLS::getNP(data._poly_order, data._dimensions-1, data._reconstruction_space);
1125 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, P_target_row.extent(0)), [&] (
const int j) {
1126 Kokkos::parallel_for(Kokkos::ThreadVectorRange(teamMember, P_target_row.extent(1)),
1127 [&] (
const int& k) {
1128 P_target_row(j,k) = 0;
1131 for (
int j = 0; j < delta.extent_int(0); ++j) {
1134 for (
int j = 0; j < thread_workspace.extent_int(0); ++j) {
1135 thread_workspace(j) = 0;
1137 teamMember.team_barrier();
1140 bool additional_evaluation_sites_need_handled =
1141 (data._additional_pc._source_coordinates.extent(0) > 0) ?
true :
false;
1143 const int num_evaluation_sites = data._additional_pc._nla.getNumberOfNeighborsDevice(target_index) + 1;
1145 for (
size_t i=0; i<data._operations.size(); ++i) {
1147 bool additional_evaluation_sites_handled =
false;
1149 bool operation_handled =
true;
1156 if (!operation_handled) {
1157 if (data._dimensions>2) {
1159 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
1160 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 , 1 , data._dimensions-1, data._poly_order,
false , &V,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, j);
1161 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
1162 for (
int k=0; k<target_NP; ++k) {
1163 P_target_row(offset, k) = delta(k);
1166 additional_evaluation_sites_handled =
true;
1169 if (data._reconstruction_space_rank == 1) {
1170 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int k) {
1172 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 , 1 , data._dimensions-1, data._poly_order,
false , &V,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, k);
1173 for (
int m=0; m<data._sampling_multiplier; ++m) {
1174 int output_components = data._basis_multiplier;
1175 for (
int c=0; c<output_components; ++c) {
1176 int offset = data._d_ss.getTargetOffsetIndex(i, m , c , k);
1180 for (
int j=0; j<target_NP; ++j) {
1181 P_target_row(offset, c*target_NP + j) = delta(j);
1186 additional_evaluation_sites_handled =
true;
1188 }
else if (data._reconstruction_space_rank == 0) {
1189 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int k) {
1191 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 , 1 , data._dimensions-1, data._poly_order,
false , &V,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, k);
1192 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, k);
1193 for (
int j=0; j<target_NP; ++j) {
1194 P_target_row(offset, j) = delta(j);
1196 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, k);
1197 for (
int j=0; j<target_NP; ++j) {
1198 P_target_row(offset, j) = 0;
1202 offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, k);
1203 for (
int j=0; j<target_NP; ++j) {
1204 P_target_row(offset, j) = 0;
1206 offset = data._d_ss.getTargetOffsetIndex(i, 1, 1, k);
1207 for (
int j=0; j<target_NP; ++j) {
1208 P_target_row(offset, j) = delta(j);
1211 additional_evaluation_sites_handled =
true;
1217 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1219 double h = data._epsilons(target_index);
1220 double a1=0, a2=0, a3=0, a4=0, a5=0;
1221 if (data._curvature_poly_order > 0) {
1222 a1 = curvature_coefficients(1);
1223 a2 = curvature_coefficients(2);
1225 if (data._curvature_poly_order > 1) {
1226 a3 = curvature_coefficients(3);
1227 a4 = curvature_coefficients(4);
1228 a5 = curvature_coefficients(5);
1230 double den = (h*h + a1*a1 + a2*a2);
1237 const int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1238 for (
int j=0; j<target_NP; ++j) {
1239 P_target_row(offset, j) = 0;
1242 if (data._poly_order > 0 && data._curvature_poly_order > 1) {
1243 P_target_row(offset, 1) = (-a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5))/den/den/(h*h);
1244 P_target_row(offset, 2) = (-a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5))/den/den/(h*h);
1246 if (data._poly_order > 1 && data._curvature_poly_order > 0) {
1247 P_target_row(offset, 3) = (h*h+a2*a2)/den/(h*h);
1248 P_target_row(offset, 4) = -2*a1*a2/den/(h*h);
1249 P_target_row(offset, 5) = (h*h+a1*a1)/den/(h*h);
1255 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1257 double h = data._epsilons(target_index);
1258 double a1=0, a2=0, a3=0, a4=0, a5=0;
1259 if (data._curvature_poly_order > 0) {
1260 a1 = curvature_coefficients(1);
1261 a2 = curvature_coefficients(2);
1263 if (data._curvature_poly_order > 1) {
1264 a3 = curvature_coefficients(3);
1265 a4 = curvature_coefficients(4);
1266 a5 = curvature_coefficients(5);
1268 double den = (h*h + a1*a1 + a2*a2);
1270 double c0a = -a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5)/den/den/h;
1271 double c1a = (h*h+a2*a2)/den/h;
1272 double c2a = -a1*a2/den/h;
1274 double c0b = -a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5)/den/den/h;
1275 double c1b = -a1*a2/den/h;
1276 double c2b = (h*h+a1*a1)/den/h;
1279 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1280 for (
int j=0; j<target_NP; ++j) {
1281 P_target_row(offset, j) = 0;
1282 P_target_row(offset, target_NP + j) = 0;
1284 P_target_row(offset, 0) = c0a;
1285 P_target_row(offset, 1) = c1a;
1286 P_target_row(offset, 2) = c2a;
1287 P_target_row(offset, target_NP + 0) = c0b;
1288 P_target_row(offset, target_NP + 1) = c1b;
1289 P_target_row(offset, target_NP + 2) = c2b;
1292 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1294 double h = data._epsilons(target_index);
1295 double a1=0, a2=0, a3=0, a4=0, a5=0;
1296 if (data._curvature_poly_order > 0) {
1297 a1 = curvature_coefficients(1);
1298 a2 = curvature_coefficients(2);
1300 if (data._curvature_poly_order > 1) {
1301 a3 = curvature_coefficients(3);
1302 a4 = curvature_coefficients(4);
1303 a5 = curvature_coefficients(5);
1305 double den = (h*h + a1*a1 + a2*a2);
1307 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1308 for (
int j=0; j<target_NP; ++j) {
1309 P_target_row(offset, j) = 0;
1313 if (data._poly_order > 0 && data._curvature_poly_order > 1) {
1314 P_target_row(offset, 1) = (-a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5))/den/den/(h*h);
1315 P_target_row(offset, 2) = (-a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5))/den/den/(h*h);
1317 if (data._poly_order > 1 && data._curvature_poly_order > 0) {
1318 P_target_row(offset, 3) = (h*h+a2*a2)/den/(h*h);
1319 P_target_row(offset, 4) = -2*a1*a2/den/(h*h);
1320 P_target_row(offset, 5) = (h*h+a1*a1)/den/(h*h);
1329 if (data._reconstruction_space_rank == 1) {
1330 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1332 double h = data._epsilons(target_index);
1333 double a1=0, a2=0, a3=0, a4=0, a5=0;
1334 if (data._curvature_poly_order > 0) {
1335 a1 = curvature_coefficients(1);
1336 a2 = curvature_coefficients(2);
1338 if (data._curvature_poly_order > 1) {
1339 a3 = curvature_coefficients(3);
1340 a4 = curvature_coefficients(4);
1341 a5 = curvature_coefficients(5);
1343 double den = (h*h + a1*a1 + a2*a2);
1345 for (
int j=0; j<target_NP; ++j) {
1350 if (data._poly_order > 0 && data._curvature_poly_order > 1) {
1351 delta(1) = (-a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4+(h*h+a1*a1)*a5))/den/den/(h*h);
1352 delta(2) = (-a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4+(h*h+a1*a1)*a5))/den/den/(h*h);
1354 if (data._poly_order > 1 && data._curvature_poly_order > 0) {
1355 delta(3) = (h*h+a2*a2)/den/(h*h);
1356 delta(4) = -2*a1*a2/den/(h*h);
1357 delta(5) = (h*h+a1*a1)/den/(h*h);
1361 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1362 for (
int j=0; j<target_NP; ++j) {
1363 P_target_row(offset, j) = delta(j);
1364 P_target_row(offset, target_NP + j) = 0;
1366 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, 0);
1367 for (
int j=0; j<target_NP; ++j) {
1368 P_target_row(offset, j) = 0;
1369 P_target_row(offset, target_NP + j) = 0;
1373 offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, 0);
1374 for (
int j=0; j<target_NP; ++j) {
1375 P_target_row(offset, j) = 0;
1376 P_target_row(offset, target_NP + j) = 0;
1378 offset = data._d_ss.getTargetOffsetIndex(i, 1, 1, 0);
1379 for (
int j=0; j<target_NP; ++j) {
1380 P_target_row(offset, j) = 0;
1381 P_target_row(offset, target_NP + j) = delta(j);
1386 }
else if (data._reconstruction_space_rank == 0) {
1387 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1389 double h = data._epsilons(target_index);
1390 double a1=0, a2=0, a3=0, a4=0, a5=0;
1391 if (data._curvature_poly_order > 0) {
1392 a1 = curvature_coefficients(1);
1393 a2 = curvature_coefficients(2);
1395 if (data._curvature_poly_order > 1) {
1396 a3 = curvature_coefficients(3);
1397 a4 = curvature_coefficients(4);
1398 a5 = curvature_coefficients(5);
1400 double den = (h*h + a1*a1 + a2*a2);
1402 for (
int j=0; j<target_NP; ++j) {
1407 if (data._poly_order > 0 && data._curvature_poly_order > 1) {
1408 delta(1) = (-a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4+(h*h+a1*a1)*a5))/den/den/(h*h);
1409 delta(2) = (-a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4+(h*h+a1*a1)*a5))/den/den/(h*h);
1411 if (data._poly_order > 1 && data._curvature_poly_order > 0) {
1412 delta(3) = (h*h+a2*a2)/den/(h*h);
1413 delta(4) = -2*a1*a2/den/(h*h);
1414 delta(5) = (h*h+a1*a1)/den/(h*h);
1418 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1419 for (
int j=0; j<target_NP; ++j) {
1420 P_target_row(offset, j) = delta(j);
1422 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, 0);
1423 for (
int j=0; j<target_NP; ++j) {
1424 P_target_row(offset, j) = 0;
1428 offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, 0);
1429 for (
int j=0; j<target_NP; ++j) {
1430 P_target_row(offset, j) = 0;
1432 offset = data._d_ss.getTargetOffsetIndex(i, 1, 1, 0);
1433 for (
int j=0; j<target_NP; ++j) {
1434 P_target_row(offset, j) = delta(j);
1441 if (data._reconstruction_space_rank == 0
1442 && (data._polynomial_sampling_functional ==
PointSample
1444 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1446 double h = data._epsilons(target_index);
1447 double a1 = curvature_coefficients(1);
1448 double a2 = curvature_coefficients(2);
1450 double q1 = (h*h + a2*a2)/(h*h*h + h*a1*a1 + h*a2*a2);
1451 double q2 = -(a1*a2)/(h*h*h + h*a1*a1 + h*a2*a2);
1452 double q3 = (h*h + a1*a1)/(h*h*h + h*a1*a1 + h*a2*a2);
1454 double t1a = q1*1 + q2*0;
1455 double t2a = q1*0 + q2*1;
1457 double t1b = q2*1 + q3*0;
1458 double t2b = q2*0 + q3*1;
1461 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1462 for (
int j=0; j<target_NP; ++j) {
1463 P_target_row(offset, j) = 0;
1465 if (data._poly_order > 0 && data._curvature_poly_order > 0) {
1466 P_target_row(offset, 1) = t1a + t2a;
1467 P_target_row(offset, 2) = 0;
1470 offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, 0);
1471 for (
int j=0; j<target_NP; ++j) {
1472 P_target_row(offset, j) = 0;
1474 if (data._poly_order > 0 && data._curvature_poly_order > 0) {
1475 P_target_row(offset, 1) = 0;
1476 P_target_row(offset, 2) = t1b + t2b;
1482 }
else if (data._reconstruction_space_rank == 1
1483 && data._polynomial_sampling_functional
1489 }
else if (data._reconstruction_space_rank == 0
1490 && data._polynomial_sampling_functional
1499 if (data._reconstruction_space_rank == 1
1501 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1503 double h = data._epsilons(target_index);
1504 double a1=0, a2=0, a3=0, a4=0, a5=0;
1505 if (data._curvature_poly_order > 0) {
1506 a1 = curvature_coefficients(1);
1507 a2 = curvature_coefficients(2);
1509 if (data._curvature_poly_order > 1) {
1510 a3 = curvature_coefficients(3);
1511 a4 = curvature_coefficients(4);
1512 a5 = curvature_coefficients(5);
1514 double den = (h*h + a1*a1 + a2*a2);
1518 double c0a = (a1*a3+a2*a4)/(h*den);
1522 double c0b = (a1*a4+a2*a5)/(h*den);
1527 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1528 for (
int j=0; j<target_NP; ++j) {
1529 P_target_row(offset, j) = 0;
1530 P_target_row(offset, target_NP + j) = 0;
1532 P_target_row(offset, 0) = c0a;
1533 P_target_row(offset, 1) = c1a;
1534 P_target_row(offset, 2) = c2a;
1537 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, 0);
1538 for (
int j=0; j<target_NP; ++j) {
1539 P_target_row(offset, j) = 0;
1540 P_target_row(offset, target_NP + j) = 0;
1542 P_target_row(offset, target_NP + 0) = c0b;
1543 P_target_row(offset, target_NP + 1) = c1b;
1544 P_target_row(offset, target_NP + 2) = c2b;
1547 }
else if (data._reconstruction_space_rank == 0
1549 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1551 double h = data._epsilons(target_index);
1552 double a1=0, a2=0, a3=0, a4=0, a5=0;
1553 if (data._curvature_poly_order > 0) {
1554 a1 = curvature_coefficients(1);
1555 a2 = curvature_coefficients(2);
1557 if (data._curvature_poly_order > 1) {
1558 a3 = curvature_coefficients(3);
1559 a4 = curvature_coefficients(4);
1560 a5 = curvature_coefficients(5);
1562 double den = (h*h + a1*a1 + a2*a2);
1566 double c0a = (a1*a3+a2*a4)/(h*den);
1570 double c0b = (a1*a4+a2*a5)/(h*den);
1575 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1576 for (
int j=0; j<target_NP; ++j) {
1577 P_target_row(offset, j) = 0;
1579 P_target_row(offset, 0) = c0a;
1580 P_target_row(offset, 1) = c1a;
1581 P_target_row(offset, 2) = c2a;
1584 offset = data._d_ss.getTargetOffsetIndex(i, 1, 0, 0);
1585 for (
int j=0; j<target_NP; ++j) {
1586 P_target_row(offset, j) = 0;
1588 P_target_row(offset, 0) = c0b;
1589 P_target_row(offset, 1) = c1b;
1590 P_target_row(offset, 2) = c2b;
1593 }
else if (data._reconstruction_space_rank == 1
1594 && data._polynomial_sampling_functional
1597 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1599 double h = data._epsilons(target_index);
1600 double a1=0, a2=0, a3=0, a4=0, a5=0;
1601 if (data._curvature_poly_order > 0) {
1602 a1 = curvature_coefficients(1);
1603 a2 = curvature_coefficients(2);
1605 if (data._curvature_poly_order > 1) {
1606 a3 = curvature_coefficients(3);
1607 a4 = curvature_coefficients(4);
1608 a5 = curvature_coefficients(5);
1610 double den = (h*h + a1*a1 + a2*a2);
1614 double c0a = -a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5)/den/den/h;
1615 double c1a = (h*h+a2*a2)/den/h;
1616 double c2a = -a1*a2/den/h;
1618 double c0b = -a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5)/den/den/h;
1619 double c1b = -a1*a2/den/h;
1620 double c2b = (h*h+a1*a1)/den/h;
1623 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1624 for (
int j=0; j<target_NP; ++j) {
1625 P_target_row(offset, j) = 0;
1626 P_target_row(offset, target_NP + j) = 0;
1628 P_target_row(offset, 0) = c0a;
1629 P_target_row(offset, 1) = c1a;
1630 P_target_row(offset, 2) = c2a;
1631 P_target_row(offset, target_NP + 0) = c0b;
1632 P_target_row(offset, target_NP + 1) = c1b;
1633 P_target_row(offset, target_NP + 2) = c2b;
1640 double h = data._epsilons(target_index);
1642 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int k) {
1645 for (
int d=0;
d<data._dimensions-1; ++
d) {
1648 relative_coord[
d] = data._additional_pc.getNeighborCoordinate(target_index, k-1,
d, &V);
1649 relative_coord[
d] -= data._pc.getTargetCoordinate(target_index,
d, &V);
1652 for (
int j=0; j<3; ++j) relative_coord[j] = 0;
1655 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, k);
1656 P_target_row(offset, 0) =
GaussianCurvature(curvature_coefficients, h, relative_coord.
x, relative_coord.
y);
1658 additional_evaluation_sites_handled =
true;
1660 double h = data._epsilons(target_index);
1662 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int k) {
1666 for (
int d=0;
d<data._dimensions-1; ++
d) {
1669 relative_coord[
d] = data._additional_pc.getNeighborCoordinate(target_index, k-1,
d, &V);
1670 relative_coord[
d] -= data._pc.getTargetCoordinate(target_index,
d, &V);
1673 for (
int j=0; j<3; ++j) relative_coord[j] = 0;
1676 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, k);
1678 for (
int n = 0; n <= data._poly_order; n++){
1679 for (alphay = 0; alphay <= n; alphay++){
1680 alphax = n - alphay;
1681 P_target_row(offset, index) =
SurfaceCurlOfScalar(curvature_coefficients, h, relative_coord.
x, relative_coord.
y, alphax, alphay, 0);
1686 offset = data._d_ss.getTargetOffsetIndex(i, 0, 1, k);
1688 for (
int n = 0; n <= data._poly_order; n++){
1689 for (alphay = 0; alphay <= n; alphay++){
1690 alphax = n - alphay;
1691 P_target_row(offset, index) =
SurfaceCurlOfScalar(curvature_coefficients, h, relative_coord.
x, relative_coord.
y, alphax, alphay, 1);
1696 additional_evaluation_sites_handled =
true;
1700 "CellAverageEvaluation and CellIntegralEvaluation only support 2d or 3d with 2d manifold");
1701 const double factorial[15] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
1702 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, 0);
1704 double cutoff_p = data._epsilons(target_index);
1711 double triangle_coords[3*3];
1712 for (
int j=0; j<data._global_dimensions*3; ++j) G_data[j] = 0;
1713 for (
int j=0; j<data._global_dimensions*3; ++j) triangle_coords[j] = 0;
1718 double radius = 0.0;
1719 for (
int j=0; j<data._global_dimensions; ++j) {
1721 triangle_coords_matrix(j, 0) = data._pc.getTargetCoordinate(target_index, j);
1722 radius += triangle_coords_matrix(j, 0)*triangle_coords_matrix(j, 0);
1724 radius = std::sqrt(radius);
1728 int num_vertices = 0;
1729 for (
int j=0; j<data._target_extra_data.extent_int(1); ++j) {
1730 auto val = data._target_extra_data(target_index, j);
1737 num_vertices = num_vertices / data._global_dimensions;
1738 auto T = triangle_coords_matrix;
1742 double entire_cell_area = 0.0;
1743 for (
int v=0; v<num_vertices; ++v) {
1745 int v2 = (v+1) % num_vertices;
1747 for (
int j=0; j<data._global_dimensions; ++j) {
1748 triangle_coords_matrix(j,1) = data._target_extra_data(target_index,
1749 v1*data._global_dimensions+j)
1750 - triangle_coords_matrix(j,0);
1751 triangle_coords_matrix(j,2) = data._target_extra_data(target_index,
1752 v2*data._global_dimensions+j)
1753 - triangle_coords_matrix(j,0);
1760 for (
int quadrature = 0; quadrature<data._qm.getNumberOfQuadraturePoints(); ++quadrature) {
1761 double unscaled_transformed_qp[3] = {0,0,0};
1762 double scaled_transformed_qp[3] = {0,0,0};
1763 for (
int j=0; j<data._global_dimensions; ++j) {
1764 for (
int k=1; k<3; ++k) {
1767 unscaled_transformed_qp[j] += T(j,k)*data._qm.getSite(quadrature, k-1);
1770 unscaled_transformed_qp[j] += T(j,0);
1774 double transformed_qp_norm = 0;
1775 for (
int j=0; j<data._global_dimensions; ++j) {
1776 transformed_qp_norm += unscaled_transformed_qp[j]*unscaled_transformed_qp[j];
1778 transformed_qp_norm = std::sqrt(transformed_qp_norm);
1781 for (
int j=0; j<data._global_dimensions; ++j) {
1782 scaled_transformed_qp[j] = unscaled_transformed_qp[j] * radius / transformed_qp_norm;
1802 double qp_norm_sq = transformed_qp_norm*transformed_qp_norm;
1803 for (
int j=0; j<data._global_dimensions; ++j) {
1804 G(j,1) = T(j,1)/transformed_qp_norm;
1805 G(j,2) = T(j,2)/transformed_qp_norm;
1806 for (
int k=0; k<data._global_dimensions; ++k) {
1807 G(j,1) += unscaled_transformed_qp[j]*(-0.5)*std::pow(qp_norm_sq,-1.5)
1808 *2*(unscaled_transformed_qp[k]*T(k,1));
1809 G(j,2) += unscaled_transformed_qp[j]*(-0.5)*std::pow(qp_norm_sq,-1.5)
1810 *2*(unscaled_transformed_qp[k]*T(k,2));
1814 Kokkos::subview(G, Kokkos::ALL(), 1), Kokkos::subview(G, Kokkos::ALL(), 2));
1815 G_determinant *= radius*radius;
1816 XYZ qp =
XYZ(scaled_transformed_qp[0], scaled_transformed_qp[1], scaled_transformed_qp[2]);
1817 for (
int j=0; j<data._local_dimensions; ++j) {
1818 relative_coord[j] = data._pc.convertGlobalToLocalCoordinate(qp,j,V)
1819 - data._pc.getTargetCoordinate(target_index,j,&V);
1824 for (
int n = 0; n <= data._poly_order; n++){
1825 for (alphay = 0; alphay <= n; alphay++){
1826 alphax = n - alphay;
1827 alphaf = factorial[alphax]*factorial[alphay];
1828 double val_to_sum = G_determinant * (data._qm.getWeight(quadrature)
1829 * std::pow(relative_coord.
x/cutoff_p,alphax)
1830 * std::pow(relative_coord.
y/cutoff_p,alphay) / alphaf);
1831 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1832 if (quadrature==0 && v==0) P_target_row(offset, k) = val_to_sum;
1833 else P_target_row(offset, k) += val_to_sum;
1838 entire_cell_area += G_determinant * data._qm.getWeight(quadrature);
1843 for (
int n = 0; n <= data._poly_order; n++){
1844 for (alphay = 0; alphay <= n; alphay++){
1845 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1846 P_target_row(offset, k) /= entire_cell_area;
1857 "FaceNormalIntegralSample, EdgeTangentIntegralSample, FaceNormalAverageSample, \
1858 and EdgeTangentAverageSample only support 2d or 3d with 2d manifold");
1860 &&
"Only 1d quadrature may be specified for edge integrals");
1862 &&
"Quadrature points not generated");
1864 const double factorial[15] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
1866 double cutoff_p = data._epsilons(target_index);
1877 int quadrature_point_loop = data._qm.getNumberOfQuadraturePoints();
1880 double G_data[3*TWO];
1881 double edge_coords[3*TWO];
1882 for (
int j=0; j<data._global_dimensions*TWO; ++j) G_data[j] = 0;
1883 for (
int j=0; j<data._global_dimensions*TWO; ++j) edge_coords[j] = 0;
1892 double radius = 0.0;
1894 for (
int j=0; j<data._global_dimensions; ++j) {
1895 edge_coords_matrix(j, 0) = data._target_extra_data(target_index, j);
1896 edge_coords_matrix(j, 1) = data._target_extra_data(target_index, data._global_dimensions + j) - edge_coords_matrix(j, 0);
1897 radius += edge_coords_matrix(j, 0)*edge_coords_matrix(j, 0);
1899 radius = std::sqrt(radius);
1903 auto E = edge_coords_matrix;
1908 XYZ a = {E(0,0), E(1,0), E(2,0)};
1909 XYZ b = {E(0,1)+E(0,0), E(1,1)+E(1,0), E(2,1)+E(2,0)};
1910 double a_dot_b = a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
1912 theta = std::atan(norm_a_cross_b / a_dot_b);
1915 for (
int c=0; c<data._local_dimensions; ++c) {
1916 int input_component = (data._sampling_multiplier==1 && data._reconstruction_space_rank==1) ? 0 : c;
1917 int offset = data._d_ss.getTargetOffsetIndex(i, input_component , 0 , 0);
1918 int column_offset = (data._reconstruction_space_rank==1) ? c*target_NP : 0;
1919 for (
int j=0; j<target_NP; ++j) {
1920 P_target_row(offset, column_offset + j) = 0;
1925 double entire_edge_length = 0.0;
1927 for (
int quadrature = 0; quadrature<quadrature_point_loop; ++quadrature) {
1929 double G_determinant = 1.0;
1930 double scaled_transformed_qp[3] = {0,0,0};
1931 double unscaled_transformed_qp[3] = {0,0,0};
1932 for (
int j=0; j<data._global_dimensions; ++j) {
1933 unscaled_transformed_qp[j] += E(j,1)*data._qm.getSite(quadrature, 0);
1935 unscaled_transformed_qp[j] += E(j,0);
1944 double transformed_qp_norm = 0;
1945 for (
int j=0; j<data._global_dimensions; ++j) {
1946 transformed_qp_norm += unscaled_transformed_qp[j]*unscaled_transformed_qp[j];
1948 transformed_qp_norm = std::sqrt(transformed_qp_norm);
1950 for (
int j=0; j<data._global_dimensions; ++j) {
1951 scaled_transformed_qp[j] = unscaled_transformed_qp[j] * radius / transformed_qp_norm;
1954 G_determinant = radius * theta;
1955 XYZ qp =
XYZ(scaled_transformed_qp[0], scaled_transformed_qp[1], scaled_transformed_qp[2]);
1956 for (
int j=0; j<data._local_dimensions; ++j) {
1957 relative_coord[j] = data._pc.convertGlobalToLocalCoordinate(qp,j,V)
1958 - data._pc.getTargetCoordinate(target_index,j,&V);
1961 relative_coord[2] = 0;
1963 XYZ endpoints_difference = {E(0,1), E(1,1), 0};
1964 G_determinant = data._pc.EuclideanVectorLength(endpoints_difference, 2);
1965 for (
int j=0; j<data._local_dimensions; ++j) {
1966 relative_coord[j] = unscaled_transformed_qp[j]
1967 - data._pc.getTargetCoordinate(target_index,j,&V);
1970 relative_coord[2] = 0;
1977 for (
int j=0; j<data._global_dimensions; ++j) {
1979 direction[j] = data._target_extra_data(target_index, 2*data._global_dimensions + j);
1984 XYZ k = {scaled_transformed_qp[0], scaled_transformed_qp[1], scaled_transformed_qp[2]};
1986 double k_norm = std::sqrt(k[0]*k[0]+k[1]*k[1]+k[2]*k[2]);
1987 k[0] = k[0]/k_norm; k[1] = k[1]/k_norm; k[2] = k[2]/k_norm;
1988 XYZ n = {data._target_extra_data(target_index, 2*data._global_dimensions + 0),
1989 data._target_extra_data(target_index, 2*data._global_dimensions + 1),
1990 data._target_extra_data(target_index, 2*data._global_dimensions + 2)};
1993 direction[0] = (k[1]*n[2] - k[2]*n[1]) / norm_k_cross_n;
1994 direction[1] = (k[2]*n[0] - k[0]*n[2]) / norm_k_cross_n;
1995 direction[2] = (k[0]*n[1] - k[1]*n[0]) / norm_k_cross_n;
1997 for (
int j=0; j<data._global_dimensions; ++j) {
1999 direction[j] = data._target_extra_data(target_index, 3*data._global_dimensions + j);
2005 XYZ local_direction;
2007 for (
int j=0; j<data._local_dimensions; ++j) {
2012 local_direction[j] = data._pc.convertGlobalToLocalCoordinate(direction,j,V);
2020 for (
int c=0; c<data._local_dimensions; ++c) {
2021 int input_component = (data._sampling_multiplier==1 && data._reconstruction_space_rank==1) ? 0 : c;
2022 int offset = data._d_ss.getTargetOffsetIndex(i, input_component, 0 , 0);
2023 int column_offset = (data._reconstruction_space_rank==1) ? c*target_NP : 0;
2025 for (
int n = 0; n <= data._poly_order; n++){
2026 for (alphay = 0; alphay <= n; alphay++){
2027 alphax = n - alphay;
2028 alphaf = factorial[alphax]*factorial[alphay];
2032 v0 = (c==0) ? std::pow(relative_coord.
x/cutoff_p,alphax)
2033 *std::pow(relative_coord.
y/cutoff_p,alphay)/alphaf : 0;
2034 v1 = (c==1) ? std::pow(relative_coord.
x/cutoff_p,alphax)
2035 *std::pow(relative_coord.
y/cutoff_p,alphay)/alphaf : 0;
2038 double dot_product = 0.0;
2041 dot_product = local_direction[0]*v0 + local_direction[1]*v1;
2043 dot_product = direction[0]*v0 + direction[1]*v1;
2047 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
2048 if (quadrature==0 && c==0) P_target_row(offset, column_offset + k) =
2049 dot_product * data._qm.getWeight(quadrature) * G_determinant;
2050 else P_target_row(offset, column_offset + k) +=
2051 dot_product * data._qm.getWeight(quadrature) * G_determinant;
2057 entire_edge_length += G_determinant * data._qm.getWeight(quadrature);
2061 for (
int c=0; c<data._local_dimensions; ++c) {
2066 int input_component = (data._sampling_multiplier==1 && data._reconstruction_space_rank==1) ? 0 : c;
2068 int offset = data._d_ss.getTargetOffsetIndex(i, input_component, 0 , 0);
2069 int column_offset = (data._reconstruction_space_rank == 1) ? c*target_NP : 0;
2070 for (
int n = 0; n <= data._poly_order; n++){
2071 for (alphay = 0; alphay <= n; alphay++){
2072 Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
2073 P_target_row(offset, column_offset + k) /= entire_edge_length;
2083 }
else if (data._dimensions==2) {
2085 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [&] (
const int j) {
2086 calcPij<TargetData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, -1 , 1 , data._dimensions-1, data._poly_order,
false , &V,
ReconstructionSpace::ScalarTaylorPolynomial,
PointSample, j);
2087 int offset = data._d_ss.getTargetOffsetIndex(i, 0, 0, j);
2088 for (
int k=0; k<target_NP; ++k) {
2089 P_target_row(offset, k) = delta(k);
2092 additional_evaluation_sites_handled =
true;
2099 compadre_kernel_assert_release(((additional_evaluation_sites_need_handled && additional_evaluation_sites_handled) || (!additional_evaluation_sites_need_handled)) &&
"Auxiliary evaluation coordinates are specified by user, but are calling a target operation that can not handle evaluating additional sites.");
2102 teamMember.team_barrier();
Average value of vector dotted with tangent directions Supported on 1D edges in 3D problems (2D-manif...
Point evaluation of Gaussian curvature.
Point evaluation of a scalar.
Scalar basis reused as many times as there are components in the vector resulting in a much cheaper p...
#define compadre_kernel_assert_release(condition)
compadre_kernel_assert_release is similar to compadre_assert_release, but is a call on the device...
constexpr SamplingFunctional ManifoldVectorPointSample
Point evaluations of the entire vector source function (but on a manifold, so it includes a transform...
Point evaluation of the gradient of a scalar.
team_policy::member_type member_type
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.
Point evaluation of the laplacian of a scalar (could be on a manifold or not)
Point evaluation of the partial with respect to y of a scalar.
Point evaluation of the partial with respect to z of a scalar.
Integral value of vector dotted with tangent directions Supported on 1D edges in 3D problems (2D-mani...
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...
KOKKOS_INLINE_FUNCTION double getAreaFromVectors(const member_type &teamMember, view_type_1 v1, view_type_2 v2)
Point evaluation of the curl of a vector (results in a vector)
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...
Point evaluation of the curl of a curl of a vector (results in a vector)
constexpr SamplingFunctional StaggeredEdgeAnalyticGradientIntegralSample
Analytical integral of a gradient source vector is just a difference of the scalar source at neighbor...
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...
Average value of vector dotted with normal direction Supported on 1D edges in 3D problems (2D-manifol...
Solve GMLS problem on a manifold (will use QR or SVD to solve the resultant GMLS problem dependent on...
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...
Point evaluation of the partial with respect to x of a scalar.
constexpr SamplingFunctional StaggeredEdgeIntegralSample
Samples consist of the result of integrals of a vector dotted with the tangent along edges between ne...
KOKKOS_INLINE_FUNCTION void computeCurvatureFunctionals(const TargetData &data, const member_type &teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P_target_row, const scratch_matrix_right_type *V, const local_index_type local_neighbor_index=-1)
Evaluates a polynomial basis for the curvature with a gradient target functional applied.
Kokkos::View< double **, layout_right, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_matrix_right_type
Kokkos::View< double *, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_vector_type
Scalar polynomial basis centered at the target site and scaled by sum of basis powers e...
import subprocess import os import re import math import sys import argparse d d d
Point evaluation of the divergence of a vector (results in a scalar)
Average values of a cell using quadrature Supported on 2D faces in 3D problems (manifold) and 2D cell...
Point evaluation of the laplacian of each component of a vector.
Divergence-free vector polynomial basis.
constexpr SamplingFunctional PointSample
Available sampling functionals.
Integral value of vector dotted with normal direction Supported on 1D edges in 3D problems (2D-manifo...
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.
Integral values over cell using quadrature Supported on 2D faces in 3D problems (manifold) and 2D cel...
constexpr SamplingFunctional VectorPointSample
Point evaluations of the entire vector source function.
Point evaluation of the gradient of a vector (results in a matrix, NOT CURRENTLY IMPLEMENTED) ...
Point evaluation of the chained staggered Laplacian acting on VectorTaylorPolynomial basis + Staggere...
Point evaluation of a vector (reconstructs entire vector at once, requiring a ReconstructionSpace hav...
#define compadre_kernel_assert_debug(condition)
Bernstein polynomial basis.