Compadre  1.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Compadre_GMLS_Targets.hpp
Go to the documentation of this file.
1 #ifndef _COMPADRE_GMLS_TARGETS_HPP_
2 #define _COMPADRE_GMLS_TARGETS_HPP_
3 
4 #include "Compadre_GMLS.hpp"
5 #include "Compadre_Manifold_Functions.hpp"
6 
7 namespace Compadre {
8 
9 KOKKOS_INLINE_FUNCTION
10 void GMLS::computeTargetFunctionals(const member_type& teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P_target_row) const {
11 
12  // check if VectorOfScalarClonesTaylorPolynomial is used with a scalar sampling functional other than PointSample
13  if (_dimensions > 1) {
14  compadre_kernel_assert_debug(
15  (_reconstruction_space!=ReconstructionSpace::VectorOfScalarClonesTaylorPolynomial
16  || _data_sampling_multiplier!=0
17  || (_reconstruction_space==ReconstructionSpace::VectorOfScalarClonesTaylorPolynomial
18  && _polynomial_sampling_functional==PointSample))
19  && "_reconstruction_space(VectorOfScalar clones incompatible with scalar output sampling functional which is not a PointSample");
20  }
21 
22  // determine if additional evaluation sites are requested by user and handled by target operations
23  bool additional_evaluation_sites_need_handled =
24  (_additional_evaluation_coordinates.extent(0) > 0) ? true : false; // additional evaluation sites are specified
25 
26  const int target_index = _initial_index_for_batch + teamMember.league_rank();
27 
28  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, P_target_row.extent(0)), [&] (const int j) {
29  Kokkos::parallel_for(Kokkos::ThreadVectorRange(teamMember, P_target_row.extent(1)),
30  [=] (const int& k) {
31  P_target_row(j,k) = 0;
32  });
33  });
34  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, delta.extent(0)), [&] (const int j) { delta(j) = 0; });
35  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, thread_workspace.extent(0)), [&] (const int j) { thread_workspace(j) = 0; });
36  teamMember.team_barrier();
37 
38  const int target_NP = this->getNP(_poly_order, _dimensions, _reconstruction_space);
39  const int num_evaluation_sites = getNEvaluationSitesPerTarget(target_index);
40 
41  for (size_t i=0; i<_operations.size(); ++i) {
42 
43  bool additional_evaluation_sites_handled = false; // target operations that can handle these sites should flip this flag
44 
45  bool operation_handled = true;
46 
47  // USER defined targets should be added to this file
48  // if the USER defined targets don't catch this operation, then operation_handled will be false
49  #include "Compadre_USER_StandardTargetFunctionals.hpp"
50 
51  // if the user didn't handle the operation, we pass it along to the toolkit's targets
52  if (!operation_handled) {
53 
54  if (_reconstruction_space == ReconstructionSpace::ScalarTaylorPolynomial) {
55 
56  /*
57  * Beginning of ScalarTaylorPolynomial basis
58  */
59 
60  if (_operations(i) == TargetOperation::ScalarPointEvaluation || (_operations(i) == TargetOperation::VectorPointEvaluation && _dimensions == 1) /* vector is a scalar in 1D */) {
61  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int j) {
62  this->calcPij(teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, _dimensions, _poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
63  int offset = getTargetOffsetIndexDevice(i, 0, 0, j);
64  for (int k=0; k<target_NP; ++k) {
65  P_target_row(offset, k) = delta(k);
66  }
67  });
68  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
69  } else if (_operations(i) == TargetOperation::LaplacianOfScalarPointEvaluation) {
70  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
71  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
72  switch (_dimensions) {
73  case 3:
74  P_target_row(offset, 4) = std::pow(_epsilons(target_index), -2);
75  P_target_row(offset, 6) = std::pow(_epsilons(target_index), -2);
76  P_target_row(offset, 9) = std::pow(_epsilons(target_index), -2);
77  break;
78  case 2:
79  P_target_row(offset, 3) = std::pow(_epsilons(target_index), -2);
80  P_target_row(offset, 5) = std::pow(_epsilons(target_index), -2);
81  break;
82  default:
83  P_target_row(offset, 2) = std::pow(_epsilons(target_index), -2);
84  }
85  });
86  } else if (_operations(i) == TargetOperation::GradientOfScalarPointEvaluation) {
87  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int j) {
88  for (int d=0; d<_dimensions; ++d) {
89  int offset = getTargetOffsetIndexDevice(i, 0, d, j);
90  auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
91  this->calcGradientPij(teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, d /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
92  }
93  });
94  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
95  } else if (_operations(i) == TargetOperation::PartialXOfScalarPointEvaluation) {
96  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int j) {
97  int offset = getTargetOffsetIndexDevice(i, 0, 0, j);
98  auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
99  this->calcGradientPij(teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 0 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
100  });
101  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
102  } else if (_operations(i) == TargetOperation::PartialYOfScalarPointEvaluation) {
103  compadre_kernel_assert_release(_dimensions>1 && "PartialYOfScalarPointEvaluation requested for dim < 2");
104  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int j) {
105  int offset = getTargetOffsetIndexDevice(i, 0, 0, j);
106  auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
107  this->calcGradientPij(teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 1 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
108  });
109  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
110  } else if (_operations(i) == TargetOperation::PartialZOfScalarPointEvaluation) {
111  compadre_kernel_assert_release(_dimensions>2 && "PartialZOfScalarPointEvaluation requested for dim < 3");
112  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int j) {
113  int offset = getTargetOffsetIndexDevice(i, 0, 0, j);
114  auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
115  this->calcGradientPij(teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 2 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
116  });
117  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
118  }
119  // staggered gradient w/ edge integrals known analytically, using a basis
120  // of potentials
121  else if (_operations(i) == TargetOperation::DivergenceOfVectorPointEvaluation
122  && _polynomial_sampling_functional == StaggeredEdgeAnalyticGradientIntegralSample) {
123  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
124  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
125  switch (_dimensions) {
126  case 3:
127  P_target_row(offset, 4) = std::pow(_epsilons(target_index), -2);
128  P_target_row(offset, 6) = std::pow(_epsilons(target_index), -2);
129  P_target_row(offset, 9) = std::pow(_epsilons(target_index), -2);
130  break;
131  case 2:
132  P_target_row(offset, 3) = std::pow(_epsilons(target_index), -2);
133  P_target_row(offset, 5) = std::pow(_epsilons(target_index), -2);
134  break;
135  default:
136  P_target_row(offset, 2) = std::pow(_epsilons(target_index), -2);
137  }
138  });
139  } else if (_operations(i) == TargetOperation::ScalarFaceAverageEvaluation) {
140  //printf("ScalarFaceAverageEvaluation\n");
141  const double factorial[15] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
142  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
143 
144  // Calculate basis matrix for NON MANIFOLD problems
145  double cutoff_p = _epsilons(target_index);
146  int alphax, alphay, alphaz;
147  double alphaf;
148 
149  double triangle_coords[9];
150  scratch_matrix_right_type triangle_coords_matrix(triangle_coords, 3, 3);
151  scratch_vector_type midpoint(delta.data(), 3);
152 
153  getMidpointFromCellVertices(teamMember, midpoint, _target_extra_data, target_index, 3 /*dim*/);
154  for (int j=0; j<3; ++j) {
155  triangle_coords_matrix(j, 0) = midpoint(j);
156  }
157  size_t num_vertices = _target_extra_data.extent(1) / 3;
158 
159  XYZ relative_coord;
160 
161  // loop over each two vertices
162  for (size_t v=0; v<num_vertices; ++v) {
163  int v1 = v;
164  int v2 = (v+1) % num_vertices;
165 
166  for (int j=0; j<3; ++j) {
167  triangle_coords_matrix(j,1) = _target_extra_data(target_index, v1*3+j) - triangle_coords_matrix(j,0);
168  triangle_coords_matrix(j,2) = _target_extra_data(target_index, v2*3+j) - triangle_coords_matrix(j,0);
169  }
170  // triangle_coords now has:
171  // (midpoint_x, midpoint_y, midpoint_z,
172  // v1_x-midpoint_x, v1_y-midpoint_y, v1_z-midpoint_z,
173  // v2_x-midpoint_x, v2_y-midpoint_y, v2_z-midpoint_z);
174  auto T=triangle_coords_matrix;
175  for (int quadrature = 0; quadrature<_qm.getNumberOfQuadraturePoints(); ++quadrature) {
176  double transformed_qp[3] = {0,0,0};
177  for (int j=0; j<3; ++j) {
178  for (int k=1; k<3; ++k) {
179  transformed_qp[j] += T(j,k)*_qm.getSite(quadrature, k-1);
180  }
181  transformed_qp[j] += T(j,0);
182  }
183  // half the norm of the cross-product is the area of the triangle
184  // so scaling is area / reference area (0.5) = the norm of the cross-product
185  double area_scaling = getAreaFromVectors(teamMember,
186  Kokkos::subview(T, Kokkos::ALL(), 1), Kokkos::subview(T, Kokkos::ALL(), 2));
187 
188  for (int j=0; j<3; ++j) {
189  relative_coord[j] = transformed_qp[j] - getTargetCoordinate(target_index, j); // shift quadrature point by target site
190  }
191 
192  int k = 0;
193 
194  const int start_index = 0;
195  for (int n = start_index; n <= _poly_order; n++){
196  for (alphaz = 0; alphaz <= n; alphaz++){
197  int s = n - alphaz;
198  for (alphay = 0; alphay <= s; alphay++){
199  alphax = s - alphay;
200  alphaf = factorial[alphax]*factorial[alphay]*factorial[alphaz];
201  double val_to_sum = (area_scaling * _qm.getWeight(quadrature)
202  * std::pow(relative_coord.x/cutoff_p,alphax)
203  *std::pow(relative_coord.y/cutoff_p,alphay)
204  *std::pow(relative_coord.z/cutoff_p,alphaz)/alphaf) / (0.5 * area_scaling);
205  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
206  if (quadrature==0) P_target_row(offset, k) = val_to_sum;
207  else P_target_row(offset, k) += val_to_sum;
208  });
209  k++;
210  }
211  }
212  }
213  }
214  }
215 
216  } else {
217  compadre_kernel_assert_release((false) && "Functionality not yet available.");
218  }
219 
220  /*
221  * End of ScalarTaylorPolynomial basis
222  */
223 
224  } else if (_reconstruction_space == ReconstructionSpace::VectorTaylorPolynomial) {
225 
226  /*
227  * Beginning of VectorTaylorPolynomial basis
228  */
229 
230  if (_operations(i) == TargetOperation::ScalarPointEvaluation || (_operations(i) == TargetOperation::VectorPointEvaluation && _dimensions == 1) /* vector is a scalar in 1D */) {
231  // copied from ScalarTaylorPolynomial
232  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int j) {
233  this->calcPij(teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, _dimensions, _poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
234  int offset = getTargetOffsetIndexDevice(i, 0, 0, j);
235  for (int k=0; k<target_NP; ++k) {
236  P_target_row(offset, k) = delta(k);
237  }
238  });
239  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
240  } else if (_operations(i) == TargetOperation::VectorPointEvaluation) {
241  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int e) {
242  this->calcPij(teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, _dimensions, _poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, e);
243  for (int m=0; m<_sampling_multiplier; ++m) {
244  int output_components = _basis_multiplier;
245  for (int c=0; c<output_components; ++c) {
246  int offset = getTargetOffsetIndexDevice(i, m /*in*/, c /*out*/, e/*additional*/);
247  // for the case where _sampling_multiplier is > 1,
248  // this approach relies on c*target_NP being equivalent to P_target_row(offset, j) where offset is
249  // getTargetOffsetIndexDevice(i, m /*in*/, c /*out*/, e/*additional*/)*_basis_multiplier*target_NP;
250  for (int j=0; j<target_NP; ++j) {
251  P_target_row(offset, c*target_NP + j) = delta(j);
252  }
253  }
254  }
255  });
256  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
257  } else if ( (_operations(i) == TargetOperation::GradientOfScalarPointEvaluation) && (_polynomial_sampling_functional == StaggeredEdgeIntegralSample) ) {
258  // when using staggered edge integral sample with vector basis, the gradient of scalar point evaluation
259  // is just the vector point evaluation
260  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int e) {
261  this->calcPij(teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, _dimensions, _poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, e);
262  for (int m=0; m<_sampling_multiplier; ++m) {
263  int output_components = _basis_multiplier;
264  for (int c=0; c<output_components; ++c) {
265  int offset = getTargetOffsetIndexDevice(i, m /*in*/, c /*out*/, e/*additional*/);
266  // for the case where _sampling_multiplier is > 1,
267  // this approach relies on c*target_NP being equivalent to P_target_row(offset, j) where offset is
268  // getTargetOffsetIndexDevice(i, m /*in*/, c /*out*/, e/*additional*/)*_basis_multiplier*target_NP;
269  for (int j=0; j<target_NP; ++j) {
270  P_target_row(offset, c*target_NP + j) = delta(j);
271  }
272  }
273  }
274  });
275  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
276  } else if (_operations(i) == TargetOperation::GradientOfScalarPointEvaluation) {
277  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int e) {
278  int output_components = _basis_multiplier;
279  for (int m2=0; m2<output_components; ++m2) { // output components 2
280  int offset = getTargetOffsetIndexDevice(i, 0 /*in*/, m2 /*out*/, e);
281  this->calcGradientPij(teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, m2 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, e);
282  for (int j=0; j<target_NP; ++j) {
283  P_target_row(offset, j) = delta(j);
284  }
285  }
286  });
287  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
288  } else if (_operations(i) == TargetOperation::GradientOfVectorPointEvaluation) {
289  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int e) {
290  for (int m0=0; m0<_sampling_multiplier; ++m0) { // input components
291  int output_components = _basis_multiplier;
292  for (int m1=0; m1<output_components; ++m1) { // output components 1
293  for (int m2=0; m2<output_components; ++m2) { // output components 2
294  int offset = getTargetOffsetIndexDevice(i, m0 /*in*/, m1*output_components+m2 /*out*/, e);
295  this->calcGradientPij(teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, m2 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, e);
296  for (int j=0; j<target_NP; ++j) {
297  P_target_row(offset, m1*target_NP + j) = delta(j);
298  }
299  }
300  }
301  }
302  });
303  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
304  } else if (_operations(i) == TargetOperation::DivergenceOfVectorPointEvaluation) {
305  if (_polynomial_sampling_functional == StaggeredEdgeIntegralSample) {
306  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
307  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);;
308  switch (_dimensions) {
309  case 3:
310  P_target_row(offset, 1) = std::pow(_epsilons(target_index), -1);
311  P_target_row(offset, target_NP + 2) = std::pow(_epsilons(target_index), -1);
312  P_target_row(offset, 2*target_NP + 3) = std::pow(_epsilons(target_index), -1);
313  break;
314  case 2:
315  P_target_row(offset, 1) = std::pow(_epsilons(target_index), -1);
316  P_target_row(offset, target_NP + 2) = std::pow(_epsilons(target_index), -1);
317  break;
318  default:
319  P_target_row(offset, 1) = std::pow(_epsilons(target_index), -1);
320  }
321  });
322  } else {
323  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int e) {
324  for (int m=0; m<_sampling_multiplier; ++m) {
325  this->calcGradientPij(teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, m /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, e);
326  int offset = getTargetOffsetIndexDevice(i, m /*in*/, 0 /*out*/, e/*additional*/);
327  for (int j=0; j<target_NP; ++j) {
328  P_target_row(offset, m*target_NP + j) = delta(j);
329  }
330  }
331  });
332  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
333  }
334  } else if (_operations(i) == TargetOperation::CurlOfVectorPointEvaluation) {
335  if (_dimensions==3) {
336  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
337  // output component 0
338  // u_{2,y} - u_{1,z}
339  {
340  int offset = getTargetOffsetIndexDevice(i, 0 /*in*/, 0 /*out*/, 0/*additional*/);
341  // role of input 0 on component 0 of curl
342  // (no contribution)
343 
344  offset = getTargetOffsetIndexDevice(i, 1 /*in*/, 0 /*out*/, 0/*additional*/);
345  // role of input 1 on component 0 of curl
346  // -u_{1,z}
347  P_target_row(offset, target_NP + 3) = -std::pow(_epsilons(target_index), -1);
348 
349  offset = getTargetOffsetIndexDevice(i, 2 /*in*/, 0 /*out*/, 0/*additional*/);
350  // role of input 2 on component 0 of curl
351  // u_{2,y}
352  P_target_row(offset, 2*target_NP + 2) = std::pow(_epsilons(target_index), -1);
353  }
354 
355  // output component 1
356  // -u_{2,x} + u_{0,z}
357  {
358  int offset = getTargetOffsetIndexDevice(i, 0 /*in*/, 1 /*out*/, 0/*additional*/);
359  // role of input 0 on component 1 of curl
360  // u_{0,z}
361  P_target_row(offset, 3) = std::pow(_epsilons(target_index), -1);
362 
363  // offset = getTargetOffsetIndexDevice(i, 1 /*in*/, 1 /*out*/, 0/*additional*/);
364  // role of input 1 on component 1 of curl
365  // (no contribution)
366 
367  offset = getTargetOffsetIndexDevice(i, 2 /*in*/, 1 /*out*/, 0/*additional*/);
368  // role of input 2 on component 1 of curl
369  // -u_{2,x}
370  P_target_row(offset, 2*target_NP + 1) = -std::pow(_epsilons(target_index), -1);
371  }
372 
373  // output component 2
374  // u_{1,x} - u_{0,y}
375  {
376  int offset = getTargetOffsetIndexDevice(i, 0 /*in*/, 2 /*out*/, 0/*additional*/);
377  // role of input 0 on component 1 of curl
378  // -u_{0,y}
379  P_target_row(offset, 2) = -std::pow(_epsilons(target_index), -1);
380 
381  offset = getTargetOffsetIndexDevice(i, 1 /*in*/, 2 /*out*/, 0/*additional*/);
382  // role of input 1 on component 1 of curl
383  // u_{1,x}
384  P_target_row(offset, target_NP + 1) = std::pow(_epsilons(target_index), -1);
385 
386  // offset = getTargetOffsetIndexDevice(i, 2 /*in*/, 2 /*out*/, 0/*additional*/);
387  // role of input 2 on component 1 of curl
388  // (no contribution)
389  }
390  });
391  } else if (_dimensions==2) {
392  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
393  // output component 0
394  // u_{1,y}
395  {
396  int offset = getTargetOffsetIndexDevice(i, 0 /*in*/, 0 /*out*/, 0/*additional*/);
397  // role of input 0 on component 0 of curl
398  // (no contribution)
399 
400  offset = getTargetOffsetIndexDevice(i, 1 /*in*/, 0 /*out*/, 0/*additional*/);
401  // role of input 1 on component 0 of curl
402  // -u_{1,z}
403  P_target_row(offset, target_NP + 2) = std::pow(_epsilons(target_index), -1);
404  }
405 
406  // output component 1
407  // -u_{0,x}
408  {
409  int offset = getTargetOffsetIndexDevice(i, 0 /*in*/, 1 /*out*/, 0/*additional*/);
410  // role of input 0 on component 1 of curl
411  // u_{0,z}
412  P_target_row(offset, 1) = -std::pow(_epsilons(target_index), -1);
413 
414  //offset = getTargetOffsetIndexDevice(i, 1 /*in*/, 1 /*out*/, 0/*additional*/);
415  // role of input 1 on component 1 of curl
416  // (no contribution)
417  }
418  });
419  }
420  } else {
421  compadre_kernel_assert_release((false) && "Functionality not yet available.");
422  }
423 
424  /*
425  * End of VectorTaylorPolynomial basis
426  */
427 
428  } else if (_reconstruction_space == ReconstructionSpace::VectorOfScalarClonesTaylorPolynomial) {
429 
430  /*
431  * Beginning of VectorOfScalarClonesTaylorPolynomial basis
432  */
433 
434  if (_operations(i) == TargetOperation::ScalarPointEvaluation || (_operations(i) == TargetOperation::VectorPointEvaluation && _dimensions == 1) /* vector is a scalar in 1D */) {
435  // copied from ScalarTaylorPolynomial
436  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int j) {
437  this->calcPij(teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, _dimensions, _poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
438  int offset = getTargetOffsetIndexDevice(i, 0, 0, j);
439  for (int k=0; k<target_NP; ++k) {
440  P_target_row(offset, k) = delta(k);
441  }
442  });
443  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
444  } else if (_operations(i) == TargetOperation::VectorPointEvaluation) {
445  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int e) {
446  this->calcPij(teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, _dimensions, _poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, e);
447  for (int m=0; m<_sampling_multiplier; ++m) {
448  for (int c=0; c<_data_sampling_multiplier; ++c) {
449  int offset = getTargetOffsetIndexDevice(i, c /*in*/, c /*out*/, e/*additional*/);
450  for (int j=0; j<target_NP; ++j) {
451  P_target_row(offset, j) = delta(j);
452  }
453  }
454  }
455  });
456  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
457  } else if (_operations(i) == TargetOperation::LaplacianOfScalarPointEvaluation) {
458  // copied from ScalarTaylorPolynomial
459  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
460  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
461  switch (_dimensions) {
462  case 3:
463  P_target_row(offset, 4) = std::pow(_epsilons(target_index), -2);
464  P_target_row(offset, 6) = std::pow(_epsilons(target_index), -2);
465  P_target_row(offset, 9) = std::pow(_epsilons(target_index), -2);
466  break;
467  case 2:
468  P_target_row(offset, 3) = std::pow(_epsilons(target_index), -2);
469  P_target_row(offset, 5) = std::pow(_epsilons(target_index), -2);
470  break;
471  default:
472  P_target_row(offset, 2) = std::pow(_epsilons(target_index), -2);
473  }
474  });
475  } else if (_operations(i) == TargetOperation::GradientOfScalarPointEvaluation) {
476  // copied from ScalarTaylorPolynomial
477  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int j) {
478  for (int d=0; d<_dimensions; ++d) {
479  int offset = getTargetOffsetIndexDevice(i, 0, d, j);
480  auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
481  this->calcGradientPij(teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, d /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
482  }
483  });
484  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
485  } else if (_operations(i) == TargetOperation::GradientOfVectorPointEvaluation) {
486  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int j) {
487  for (int m0=0; m0<_dimensions; ++m0) { // input components
488  for (int m1=0; m1<_dimensions; ++m1) { // output components 1
489  for (int m2=0; m2<_dimensions; ++m2) { // output componets 2
490  int offset = getTargetOffsetIndexDevice(i, m0 /*in*/, m1*_dimensions+m2 /*out*/, j);
491  auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
492  this->calcGradientPij(teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, m2 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
493  }
494  }
495  }
496  });
497  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
498  } else if (_operations(i) == TargetOperation::PartialXOfScalarPointEvaluation) {
499  // copied from ScalarTaylorPolynomial
500  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int j) {
501  int offset = getTargetOffsetIndexDevice(i, 0, 0, j);
502  auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
503  this->calcGradientPij(teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 0 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
504  });
505  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
506  } else if (_operations(i) == TargetOperation::PartialYOfScalarPointEvaluation) {
507  // copied from ScalarTaylorPolynomial
508  if (_dimensions>1) {
509  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int j) {
510  int offset = getTargetOffsetIndexDevice(i, 0, 0, j);
511  auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
512  this->calcGradientPij(teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 1 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
513  });
514  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
515  }
516  } else if (_operations(i) == TargetOperation::PartialZOfScalarPointEvaluation) {
517  // copied from ScalarTaylorPolynomial
518  if (_dimensions>2) {
519  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int j) {
520  int offset = getTargetOffsetIndexDevice(i, 0, 0, j);
521  auto row = Kokkos::subview(P_target_row, offset, Kokkos::ALL());
522  this->calcGradientPij(teamMember, row.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, 2 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
523  });
524  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
525  }
526  } else if (_operations(i) == TargetOperation::DivergenceOfVectorPointEvaluation) {
527  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
528  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
529  for (int j=0; j<target_NP; ++j) {
530  P_target_row(offset, j) = 0;
531  }
532 
533  P_target_row(offset, 1) = std::pow(_epsilons(target_index), -1);
534 
535  if (_dimensions>1) {
536  offset = getTargetOffsetIndexDevice(i, 1, 0, 0);
537  for (int j=0; j<target_NP; ++j) {
538  P_target_row(offset, j) = 0;
539  }
540  P_target_row(offset, 2) = std::pow(_epsilons(target_index), -1);
541  }
542 
543  if (_dimensions>2) {
544  offset = getTargetOffsetIndexDevice(i, 2, 0, 0);
545  for (int j=0; j<target_NP; ++j) {
546  P_target_row(offset, j) = 0;
547  }
548  P_target_row(offset, 3) = std::pow(_epsilons(target_index), -1);
549  }
550  });
551  } else if (_operations(i) == TargetOperation::CurlOfVectorPointEvaluation) {
552  // comments based on taking curl of vector [u_{0},u_{1},u_{2}]^T
553  // with as e.g., u_{1,z} being the partial derivative with respect to z of
554  // u_{1}
555  if (_dimensions==3) {
556  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
557  // output component 0
558  // u_{2,y} - u_{1,z}
559  {
560  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
561  for (int j=0; j<target_NP; ++j) {
562  P_target_row(offset, j) = 0;
563  }
564  // role of input 0 on component 0 of curl
565  // (no contribution)
566 
567  offset = getTargetOffsetIndexDevice(i, 1, 0, 0);
568  for (int j=0; j<target_NP; ++j) {
569  P_target_row(offset, j) = 0;
570  }
571  // role of input 1 on component 0 of curl
572  // -u_{1,z}
573  P_target_row(offset, 3) = -std::pow(_epsilons(target_index), -1);
574 
575  offset = getTargetOffsetIndexDevice(i, 2, 0, 0);
576  for (int j=0; j<target_NP; ++j) {
577  P_target_row(offset, j) = 0;
578  }
579  // role of input 2 on component 0 of curl
580  // u_{2,y}
581  P_target_row(offset, 2) = std::pow(_epsilons(target_index), -1);
582  }
583 
584  // output component 1
585  // -u_{2,x} + u_{0,z}
586  {
587  int offset = getTargetOffsetIndexDevice(i, 0, 1, 0);
588  for (int j=0; j<target_NP; ++j) {
589  P_target_row(offset, j) = 0;
590  }
591  // role of input 0 on component 1 of curl
592  // u_{0,z}
593  P_target_row(offset, 3) = std::pow(_epsilons(target_index), -1);
594 
595  offset = getTargetOffsetIndexDevice(i, 1, 1, 0);
596  for (int j=0; j<target_NP; ++j) {
597  P_target_row(offset, j) = 0;
598  }
599  // role of input 1 on component 1 of curl
600  // (no contribution)
601 
602  offset = getTargetOffsetIndexDevice(i, 2, 1, 0);
603  for (int j=0; j<target_NP; ++j) {
604  P_target_row(offset, j) = 0;
605  }
606  // role of input 2 on component 1 of curl
607  // -u_{2,x}
608  P_target_row(offset, 1) = -std::pow(_epsilons(target_index), -1);
609  }
610 
611  // output component 2
612  // u_{1,x} - u_{0,y}
613  {
614  int offset = getTargetOffsetIndexDevice(i, 0, 2, 0);
615  for (int j=0; j<target_NP; ++j) {
616  P_target_row(offset, j) = 0;
617  }
618  // role of input 0 on component 1 of curl
619  // -u_{0,y}
620  P_target_row(offset, 2) = -std::pow(_epsilons(target_index), -1);
621 
622  offset = getTargetOffsetIndexDevice(i, 1, 2, 0);
623  for (int j=0; j<target_NP; ++j) {
624  P_target_row(offset, j) = 0;
625  }
626  // role of input 1 on component 1 of curl
627  // u_{1,x}
628  P_target_row(offset, 1) = std::pow(_epsilons(target_index), -1);
629 
630  offset = getTargetOffsetIndexDevice(i, 2, 2, 0);
631  for (int j=0; j<target_NP; ++j) {
632  P_target_row(offset, j) = 0;
633  }
634  // role of input 2 on component 1 of curl
635  // (no contribution)
636  }
637  });
638  } else if (_dimensions==2) {
639  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
640  // output component 0
641  // u_{1,y}
642  {
643  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
644  for (int j=0; j<target_NP; ++j) {
645  P_target_row(offset, j) = 0;
646  }
647  // role of input 0 on component 0 of curl
648  // (no contribution)
649 
650  offset = getTargetOffsetIndexDevice(i, 1, 0, 0);
651  for (int j=0; j<target_NP; ++j) {
652  P_target_row(offset, j) = 0;
653  }
654  // role of input 1 on component 0 of curl
655  // -u_{1,z}
656  P_target_row(offset, 2) = std::pow(_epsilons(target_index), -1);
657  }
658 
659  // output component 1
660  // -u_{0,x}
661  {
662  int offset = getTargetOffsetIndexDevice(i, 0, 1, 0);
663  for (int j=0; j<target_NP; ++j) {
664  P_target_row(offset, j) = 0;
665  }
666  // role of input 0 on component 1 of curl
667  // u_{0,z}
668  P_target_row(offset, 1) = -std::pow(_epsilons(target_index), -1);
669 
670  offset = getTargetOffsetIndexDevice(i, 1, 1, 0);
671  for (int j=0; j<target_NP; ++j) {
672  P_target_row(offset, j) = 0;
673  }
674  // role of input 1 on component 1 of curl
675  // (no contribution)
676  }
677  });
678  }
679  } else {
680  compadre_kernel_assert_release((false) && "Functionality not yet available.");
681  }
682 
683  /*
684  * End of VectorOfScalarClonesTaylorPolynomial basis
685  */
686 
687  } else if (_reconstruction_space == ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial) {
688 
689  /*
690  * Beginning of DivergenceFreeVectorTaylorPolynomial basis
691  */
692 
693  if (_operations(i) == TargetOperation::VectorPointEvaluation) {
694  // copied from VectorTaylorPolynomial
695  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int e) {
696  for (int m0=0; m0<_sampling_multiplier; ++m0) {
697  for (int m1=0; m1<_sampling_multiplier; ++m1) {
698 
699  this->calcPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(m1+1) /* target is neighbor, but also which component */, 1 /*alpha*/, _dimensions, _poly_order, false /*bool on only specific order*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e /* evaluate at target */);
700 
701  int offset = getTargetOffsetIndexDevice(i, m0 /*in*/, m1 /*out*/, e /*no additional*/);
702  for (int j=0; j<target_NP; ++j) {
703  P_target_row(offset, j) = delta(j);
704  }
705  }
706  }
707  });
708  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
709  } else if (_operations(i) == TargetOperation::CurlOfVectorPointEvaluation) {
710  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int e) {
711  for (int m0=0; m0<_sampling_multiplier; ++m0) { // input components
712  for (int m1=0; m1<_sampling_multiplier; ++m1) { // output components
713  int offset = getTargetOffsetIndexDevice(i, m0 /*in*/, m1 /*out*/, e /*no additional*/);
714  if (_dimensions==3) {
715  switch (m1) {
716  case 0:
717  // output component 0
718  this->calcGradientPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
719  // u2y
720  for (int j=0; j<target_NP; ++j) {
721  P_target_row(offset, j) = delta(j);
722  }
723  this->calcGradientPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 2 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
724  // -u1z
725  for (int j=0; j<target_NP; ++j) {
726  P_target_row(offset, j) -= delta(j);
727  }
728  // u2y - u1z
729  break;
730  case 1:
731  // output component 1
732  this->calcGradientPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
733  // -u2x
734  for (int j=0; j<target_NP; ++j) {
735  P_target_row(offset, j) = -delta(j);
736  }
737  this->calcGradientPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 2 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
738  // u0z
739  for (int j=0; j<target_NP; ++j) {
740  P_target_row(offset, j) += delta(j);
741  }
742  // -u2x + u0z
743  break;
744  default:
745  // output component 2
746  this->calcGradientPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
747  // u1x
748  for (int j=0; j<target_NP; ++j) {
749  P_target_row(offset, j) = delta(j);
750  }
751  this->calcGradientPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
752  // -u0y
753  for (int j=0; j<target_NP; ++j) {
754  P_target_row(offset, j) -= delta(j);
755  }
756  // u1x - u0y
757  break;
758  }
759  } else {
760  if (m1==0) {
761  // curl results in 1D output
762  P_target_row(offset, 2) = -std::pow(_epsilons(target_index), -1);
763  P_target_row(offset, 3) = std::pow(_epsilons(target_index), -1);
764 
765  this->calcGradientPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
766  // u1x
767  for (int j=0; j<target_NP; ++j) {
768  P_target_row(offset, j) = delta(j);
769  }
770  this->calcGradientPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
771  // -u0y
772  for (int j=0; j<target_NP; ++j) {
773  P_target_row(offset, j) -= delta(j);
774  }
775  // u1x - u0y
776  }
777  }
778  }
779  }
780  });
781  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
782  } else if (_operations(i) == TargetOperation::CurlCurlOfVectorPointEvaluation) {
783  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int e) {
784  for (int m0=0; m0<_sampling_multiplier; ++m0) { // input components
785  for (int m1=0; m1<_sampling_multiplier; ++m1) { // output components
786  int offset = getTargetOffsetIndexDevice(i, m0 /*in*/, m1 /*out*/, e /*no additional*/);
787  if (_dimensions == 3) {
788  switch (m1) {
789  // manually compute the output components
790  case 0:
791  // output component 0
792  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 1 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
793  // u1xy
794  for (int j=0; j<target_NP; ++j) {
795  P_target_row(offset, j) = delta(j);
796  }
797  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 1 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
798  // -u0yy
799  for (int j=0; j<target_NP; ++j) {
800  P_target_row(offset, j) -= delta(j);
801  }
802  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 2 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
803  // u2xz
804  for (int j=0; j<target_NP; ++j) {
805  P_target_row(offset, j) += delta(j);
806  }
807  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 2 /*partial_direction_1*/, 2 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
808  // -u0zz
809  for (int j=0; j<target_NP; ++j) {
810  P_target_row(offset, j) -= delta(j);
811  }
812  // u1xy - u0yy + u2xz - u0zz
813  break;
814  case 1:
815  // output component 1
816  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 0 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
817  // -u1xx
818  for (int j=0; j<target_NP; ++j) {
819  P_target_row(offset, j) = -delta(j);
820  }
821  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 0 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
822  // u0yx
823  for (int j=0; j<target_NP; ++j) {
824  P_target_row(offset, j) += delta(j);
825  }
826  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 2 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
827  // u2yz
828  for (int j=0; j<target_NP; ++j) {
829  P_target_row(offset, j) += delta(j);
830  }
831  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 2 /*partial_direction_1*/, 2 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
832  // -u1zz
833  for (int j=0; j<target_NP; ++j) {
834  P_target_row(offset, j) -= delta(j);
835  }
836  // -u1xx + u0yx + u2yz - u1zz
837  break;
838  default:
839  // output component 2
840  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 0 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
841  // -u2xx
842  for (int j=0; j<target_NP; ++j) {
843  P_target_row(offset, j) = -delta(j);
844  }
845  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 2 /*partial_direction_1*/, 0 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
846  // u0zx
847  for (int j=0; j<target_NP; ++j) {
848  P_target_row(offset, j) += delta(j);
849  }
850  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(2+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 1 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
851  // -u2yy
852  for (int j=0; j<target_NP; ++j) {
853  P_target_row(offset, j) -= delta(j);
854  }
855  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 2 /*partial_direction_1*/, 1 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
856  // u1zy
857  for (int j=0; j<target_NP; ++j) {
858  P_target_row(offset, j) += delta(j);
859  }
860  // -u2xx + u0zx - u2yy + u1zy
861  break;
862  }
863  }
864  if (_dimensions == 2) {
865  // uses curl curl u = grad ( div (u) ) - laplace (u)
866  switch (m1) {
867  case 0:
868  // output component 0
869  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 0 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
870  // u0xx
871  for (int j=0; j<target_NP; ++j) {
872  P_target_row(offset, j) = delta(j);
873  }
874  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 0 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
875  // u1yx
876  for (int j=0; j<target_NP; ++j) {
877  P_target_row(offset, j) += delta(j);
878  }
879  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 0 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
880  // -u0xx
881  for (int j=0; j<target_NP; ++j) {
882  P_target_row(offset, j) -= delta(j);
883  }
884  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 1 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
885  // -u0yy
886  for (int j=0; j<target_NP; ++j) {
887  P_target_row(offset, j) -= delta(j);
888  }
889  // u0xx + u1yx - u0xx - u0yy
890  break;
891  case 1:
892  // output component 1
893  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(0+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 1 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
894  // u0xy
895  for (int j=0; j<target_NP; ++j) {
896  P_target_row(offset, j) = delta(j);
897  }
898  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 1 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
899  // u1yy
900  for (int j=0; j<target_NP; ++j) {
901  P_target_row(offset, j) += delta(j);
902  }
903  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 0 /*partial_direction_1*/, 0 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
904  // -u1xx
905  for (int j=0; j<target_NP; ++j) {
906  P_target_row(offset, j) -= delta(j);
907  }
908  this->calcHessianPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(1+1) /* -(component+1) */, 1 /*alpha*/, 1 /*partial_direction_1*/, 1 /*partial_direction_2*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
909  // -u1yy
910  for (int j=0; j<target_NP; ++j) {
911  P_target_row(offset, j) -= delta(j);
912  }
913  // u0xy + u1yy - u1xx - u1yy
914  break;
915  }
916  }
917  }
918  }
919  });
920  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
921  } else if (_operations(i) == TargetOperation::GradientOfVectorPointEvaluation) {
922  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int e) {
923  for (int m0=0; m0<_sampling_multiplier; ++m0) { // input components
924  for (int m1=0; m1<_sampling_multiplier; ++m1) { // output components 1
925  for (int m2=0; m2<_sampling_multiplier; ++m2) { // output components 2
926  int offset = getTargetOffsetIndexDevice(i, m0 /*in*/, m1*_sampling_multiplier+m2 /*out*/, e);
927  this->calcGradientPij(teamMember, delta.data(), thread_workspace.data(), target_index, -(m1+1) /* target is neighbor */, 1 /*alpha*/, m2 /*partial_direction*/, _dimensions, _poly_order, false /*specific order only*/, NULL /*&V*/, ReconstructionSpace::DivergenceFreeVectorTaylorPolynomial, VectorPointSample, e);
928  for (int j=0; j<target_NP; ++j) {
929  P_target_row(offset, j) = delta(j);
930  }
931  }
932  }
933  }
934  });
935  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
936  }
937  } else {
938  compadre_kernel_assert_release((false) && "Functionality not yet available.");
939  }
940 
941  /*
942  * End of DivergenceFreeVectorTaylorPolynomial basis
943  */
944  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.");
945  } // !operation_handled
946  teamMember.team_barrier();
947  }
948 }
949 
950 KOKKOS_INLINE_FUNCTION
951 void GMLS::computeCurvatureFunctionals(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) const {
952 
953  compadre_kernel_assert_release(((int)thread_workspace.extent(0)>=(_curvature_poly_order+1)*_local_dimensions) && "Workspace thread_workspace not large enough.");
954 
955  const int target_index = _initial_index_for_batch + teamMember.league_rank();
956 
957  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, P_target_row.extent(0)), [&] (const int j) {
958  Kokkos::parallel_for(Kokkos::ThreadVectorRange(teamMember, P_target_row.extent(1)),
959  [=] (const int& k) {
960  P_target_row(j,k) = 0;
961  });
962  });
963  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, delta.extent(0)), [&] (const int j) { delta(j) = 0; });
964  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, thread_workspace.extent(0)), [&] (const int j) { thread_workspace(j) = 0; });
965  teamMember.team_barrier();
966 
967  const int manifold_NP = this->getNP(_curvature_poly_order, _dimensions-1, ReconstructionSpace::ScalarTaylorPolynomial);
968  for (size_t i=0; i<_curvature_support_operations.size(); ++i) {
969  if (_curvature_support_operations(i) == TargetOperation::ScalarPointEvaluation) {
970  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
971  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
972  this->calcPij(teamMember, delta.data(), thread_workspace.data(), target_index, local_neighbor_index, 0 /*alpha*/, _dimensions-1, _curvature_poly_order, false /*bool on only specific order*/, V, ReconstructionSpace::ScalarTaylorPolynomial, PointSample);
973  for (int j=0; j<manifold_NP; ++j) {
974  P_target_row(offset, j) = delta(j);
975  }
976  });
977  } else if (_curvature_support_operations(i) == TargetOperation::GradientOfScalarPointEvaluation) {
978  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
979  //int offset = i*manifold_NP;
980  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
981  this->calcGradientPij(teamMember, delta.data(), thread_workspace.data(), target_index, local_neighbor_index, 0 /*alpha*/, 0 /*partial_direction*/, _dimensions-1, _curvature_poly_order, false /*specific order only*/, V, ReconstructionSpace::ScalarTaylorPolynomial, PointSample);
982  for (int j=0; j<manifold_NP; ++j) {
983  P_target_row(offset, j) = delta(j);
984  }
985  if (_dimensions>2) { // _dimensions-1 > 1
986  //offset = (i+1)*manifold_NP;
987  offset = getTargetOffsetIndexDevice(i, 0, 1, 0);
988  this->calcGradientPij(teamMember, delta.data(), thread_workspace.data(), target_index, local_neighbor_index, 0 /*alpha*/, 1 /*partial_direction*/, _dimensions-1, _curvature_poly_order, false /*specific order only*/, V, ReconstructionSpace::ScalarTaylorPolynomial, PointSample);
989  for (int j=0; j<manifold_NP; ++j) {
990  P_target_row(offset, j) = delta(j);
991  }
992  }
993  });
994  } else {
995  compadre_kernel_assert_release((false) && "Functionality not yet available.");
996  }
997  }
998 }
999 
1000 KOKKOS_INLINE_FUNCTION
1001 void GMLS::computeTargetFunctionalsOnManifold(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_matrix_right_type G_inv, scratch_vector_type curvature_coefficients, scratch_vector_type curvature_gradients) const {
1002 
1003  compadre_kernel_assert_release(((int)thread_workspace.extent(0)>=(_poly_order+1)*_local_dimensions) && "Workspace thread_workspace not large enough.");
1004 
1005  // only designed for 2D manifold embedded in 3D space
1006  const int target_index = _initial_index_for_batch + teamMember.league_rank();
1007  const int target_NP = this->getNP(_poly_order, _dimensions-1, _reconstruction_space);
1008 
1009  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, P_target_row.extent(0)), [&] (const int j) {
1010  Kokkos::parallel_for(Kokkos::ThreadVectorRange(teamMember, P_target_row.extent(1)),
1011  [=] (const int& k) {
1012  P_target_row(j,k) = 0;
1013  });
1014  });
1015  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, delta.extent(0)), [&] (const int j) { delta(j) = 0; });
1016  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, thread_workspace.extent(0)), [&] (const int j) { thread_workspace(j) = 0; });
1017  teamMember.team_barrier();
1018 
1019  // determine if additional evaluation sites are requested by user and handled by target operations
1020  bool additional_evaluation_sites_need_handled =
1021  (_additional_evaluation_coordinates.extent(0) > 0) ? true : false; // additional evaluation sites are specified
1022 
1023  const int num_evaluation_sites = getNEvaluationSitesPerTarget(target_index);
1024 
1025  for (size_t i=0; i<_operations.size(); ++i) {
1026 
1027  bool additional_evaluation_sites_handled = false; // target operations that can handle these sites should flip this flag
1028 
1029  bool operation_handled = true;
1030 
1031  // USER defined targets on the manifold should be added to this file
1032  // if the USER defined targets don't catch this operation, then operation_handled will be false
1033  #include "Compadre_USER_ManifoldTargetFunctionals.hpp"
1034 
1035  // if the user didn't handle the operation, we pass it along to the toolkit's targets
1036  if (!operation_handled) {
1037  if (_dimensions>2) {
1038  if (_operations(i) == TargetOperation::ScalarPointEvaluation) {
1039  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int j) {
1040  this->calcPij(teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, _dimensions-1, _poly_order, false /*bool on only specific order*/, &V, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, j);
1041  int offset = getTargetOffsetIndexDevice(i, 0, 0, j);
1042  for (int k=0; k<target_NP; ++k) {
1043  P_target_row(offset, k) = delta(k);
1044  }
1045  });
1046  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
1047  } else if (_operations(i) == TargetOperation::VectorPointEvaluation) {
1048  // vector basis
1049  if (_reconstruction_space_rank == 1) {
1050  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int k) {
1051  // output component 0
1052  this->calcPij(teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, _dimensions-1, _poly_order, false /*bool on only specific order*/, &V, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, k);
1053  int offset = getTargetOffsetIndexDevice(i, 0, 0, k);
1054  for (int j=0; j<target_NP; ++j) {
1055  P_target_row(offset, j) = delta(j);
1056  P_target_row(offset, target_NP + j) = 0;
1057  }
1058  offset = getTargetOffsetIndexDevice(i, 1, 0, k);
1059  for (int j=0; j<target_NP; ++j) {
1060  P_target_row(offset, j) = 0;
1061  P_target_row(offset, target_NP + j) = 0;
1062  }
1063 
1064  // output component 1
1065  offset = getTargetOffsetIndexDevice(i, 0, 1, k);
1066  for (int j=0; j<target_NP; ++j) {
1067  P_target_row(offset, j) = 0;
1068  P_target_row(offset, target_NP + j) = 0;
1069  }
1070  offset = getTargetOffsetIndexDevice(i, 1, 1, k);
1071  for (int j=0; j<target_NP; ++j) {
1072  P_target_row(offset, j) = 0;
1073  P_target_row(offset, target_NP + j) = delta(j);
1074  }
1075  });
1076  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
1077  // scalar basis times number of components in the vector
1078  } else if (_reconstruction_space_rank == 0) {
1079  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int k) {
1080  // output component 0
1081  this->calcPij(teamMember, delta.data(), thread_workspace.data(), target_index, -1 /* target is neighbor */, 1 /*alpha*/, _dimensions-1, _poly_order, false /*bool on only specific order*/, &V, ReconstructionSpace::ScalarTaylorPolynomial, PointSample, k);
1082  int offset = getTargetOffsetIndexDevice(i, 0, 0, k);
1083  for (int j=0; j<target_NP; ++j) {
1084  P_target_row(offset, j) = delta(j);
1085  }
1086  offset = getTargetOffsetIndexDevice(i, 1, 0, k);
1087  for (int j=0; j<target_NP; ++j) {
1088  P_target_row(offset, j) = 0;
1089  }
1090 
1091  // output component 1
1092  offset = getTargetOffsetIndexDevice(i, 0, 1, k);
1093  for (int j=0; j<target_NP; ++j) {
1094  P_target_row(offset, j) = 0;
1095  }
1096  offset = getTargetOffsetIndexDevice(i, 1, 1, k);
1097  for (int j=0; j<target_NP; ++j) {
1098  P_target_row(offset, j) = delta(j);
1099  }
1100  });
1101  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
1102  }
1103 
1104  } else if (_operations(i) == TargetOperation::LaplacianOfScalarPointEvaluation) {
1105  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1106 
1107  double h = _epsilons(target_index);
1108  double a1=0, a2=0, a3=0, a4=0, a5=0;
1109  if (_curvature_poly_order > 0) {
1110  a1 = curvature_coefficients(1);
1111  a2 = curvature_coefficients(2);
1112  }
1113  if (_curvature_poly_order > 1) {
1114  a3 = curvature_coefficients(3);
1115  a4 = curvature_coefficients(4);
1116  a5 = curvature_coefficients(5);
1117  }
1118  double den = (h*h + a1*a1 + a2*a2);
1119 
1120  // Gaussian Curvature sanity check
1121  //double K_curvature = ( - a4*a4 + a3*a5) / den / den;
1122  //std::cout << "Gaussian curvature is: " << K_curvature << std::endl;
1123 
1124 
1125  const int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
1126  for (int j=0; j<target_NP; ++j) {
1127  P_target_row(offset, j) = 0;
1128  }
1129  // scaled
1130  if (_poly_order > 0 && _curvature_poly_order > 1) {
1131  P_target_row(offset, 1) = (-a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5))/den/den/(h*h);
1132  P_target_row(offset, 2) = (-a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5))/den/den/(h*h);
1133  }
1134  if (_poly_order > 1 && _curvature_poly_order > 0) {
1135  P_target_row(offset, 3) = (h*h+a2*a2)/den/(h*h);
1136  P_target_row(offset, 4) = -2*a1*a2/den/(h*h);
1137  P_target_row(offset, 5) = (h*h+a1*a1)/den/(h*h);
1138  }
1139 
1140  });
1141  } else if (_operations(i) == TargetOperation::ChainedStaggeredLaplacianOfScalarPointEvaluation) {
1142  if (_reconstruction_space == ReconstructionSpace::VectorTaylorPolynomial) {
1143  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1144 
1145  double h = _epsilons(target_index);
1146  double a1=0, a2=0, a3=0, a4=0, a5=0;
1147  if (_curvature_poly_order > 0) {
1148  a1 = curvature_coefficients(1);
1149  a2 = curvature_coefficients(2);
1150  }
1151  if (_curvature_poly_order > 1) {
1152  a3 = curvature_coefficients(3);
1153  a4 = curvature_coefficients(4);
1154  a5 = curvature_coefficients(5);
1155  }
1156  double den = (h*h + a1*a1 + a2*a2);
1157 
1158  double c0a = -a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5)/den/den/h;
1159  double c1a = (h*h+a2*a2)/den/h;
1160  double c2a = -a1*a2/den/h;
1161 
1162  double c0b = -a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5)/den/den/h;
1163  double c1b = -a1*a2/den/h;
1164  double c2b = (h*h+a1*a1)/den/h;
1165 
1166  // 1st input component
1167  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
1168  for (int j=0; j<target_NP; ++j) {
1169  P_target_row(offset, j) = 0;
1170  P_target_row(offset, target_NP + j) = 0;
1171  }
1172  P_target_row(offset, 0) = c0a;
1173  P_target_row(offset, 1) = c1a;
1174  P_target_row(offset, 2) = c2a;
1175  P_target_row(offset, target_NP + 0) = c0b;
1176  P_target_row(offset, target_NP + 1) = c1b;
1177  P_target_row(offset, target_NP + 2) = c2b;
1178  });
1179  } else if (_reconstruction_space == ReconstructionSpace::ScalarTaylorPolynomial) {
1180  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1181 
1182  double h = _epsilons(target_index);
1183  double a1=0, a2=0, a3=0, a4=0, a5=0;
1184  if (_curvature_poly_order > 0) {
1185  a1 = curvature_coefficients(1);
1186  a2 = curvature_coefficients(2);
1187  }
1188  if (_curvature_poly_order > 1) {
1189  a3 = curvature_coefficients(3);
1190  a4 = curvature_coefficients(4);
1191  a5 = curvature_coefficients(5);
1192  }
1193  double den = (h*h + a1*a1 + a2*a2);
1194 
1195  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
1196  for (int j=0; j<target_NP; ++j) {
1197  P_target_row(offset, j) = 0;
1198  }
1199 
1200  // verified
1201  if (_poly_order > 0 && _curvature_poly_order > 1) {
1202  P_target_row(offset, 1) = (-a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5))/den/den/(h*h);
1203  P_target_row(offset, 2) = (-a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5))/den/den/(h*h);
1204  }
1205  if (_poly_order > 1 && _curvature_poly_order > 0) {
1206  P_target_row(offset, 3) = (h*h+a2*a2)/den/(h*h);
1207  P_target_row(offset, 4) = -2*a1*a2/den/(h*h);
1208  P_target_row(offset, 5) = (h*h+a1*a1)/den/(h*h);
1209  }
1210 
1211  });
1212  }
1213  } else if (_operations(i) == TargetOperation::VectorLaplacianPointEvaluation) {
1214  // vector basis
1215  if (_reconstruction_space_rank == 1) {
1216  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1217 
1218  double h = _epsilons(target_index);
1219  double a1=0, a2=0, a3=0, a4=0, a5=0;
1220  if (_curvature_poly_order > 0) {
1221  a1 = curvature_coefficients(1);
1222  a2 = curvature_coefficients(2);
1223  }
1224  if (_curvature_poly_order > 1) {
1225  a3 = curvature_coefficients(3);
1226  a4 = curvature_coefficients(4);
1227  a5 = curvature_coefficients(5);
1228  }
1229  double den = (h*h + a1*a1 + a2*a2);
1230 
1231  for (int j=0; j<target_NP; ++j) {
1232  delta(j) = 0;
1233  }
1234 
1235  // 1/Sqrt[Det[G[r, s]]])*Div[Sqrt[Det[G[r, s]]]*Inv[G]*P
1236  if (_poly_order > 0 && _curvature_poly_order > 1) {
1237  delta(1) = (-a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4+(h*h+a1*a1)*a5))/den/den/(h*h);
1238  delta(2) = (-a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4+(h*h+a1*a1)*a5))/den/den/(h*h);
1239  }
1240  if (_poly_order > 1 && _curvature_poly_order > 0) {
1241  delta(3) = (h*h+a2*a2)/den/(h*h);
1242  delta(4) = -2*a1*a2/den/(h*h);
1243  delta(5) = (h*h+a1*a1)/den/(h*h);
1244  }
1245 
1246  // output component 0
1247  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
1248  for (int j=0; j<target_NP; ++j) {
1249  P_target_row(offset, j) = delta(j);
1250  P_target_row(offset, target_NP + j) = 0;
1251  }
1252  offset = getTargetOffsetIndexDevice(i, 1, 0, 0);
1253  for (int j=0; j<target_NP; ++j) {
1254  P_target_row(offset, j) = 0;
1255  P_target_row(offset, target_NP + j) = 0;
1256  }
1257 
1258  // output component 1
1259  offset = getTargetOffsetIndexDevice(i, 0, 1, 0);
1260  for (int j=0; j<target_NP; ++j) {
1261  P_target_row(offset, j) = 0;
1262  P_target_row(offset, target_NP + j) = 0;
1263  }
1264  offset = getTargetOffsetIndexDevice(i, 1, 1, 0);
1265  for (int j=0; j<target_NP; ++j) {
1266  P_target_row(offset, j) = 0;
1267  P_target_row(offset, target_NP + j) = delta(j);
1268  }
1269 
1270  });
1271  // scalar basis times number of components in the vector
1272  } else if (_reconstruction_space_rank == 0) {
1273  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1274 
1275  double h = _epsilons(target_index);
1276  double a1=0, a2=0, a3=0, a4=0, a5=0;
1277  if (_curvature_poly_order > 0) {
1278  a1 = curvature_coefficients(1);
1279  a2 = curvature_coefficients(2);
1280  }
1281  if (_curvature_poly_order > 1) {
1282  a3 = curvature_coefficients(3);
1283  a4 = curvature_coefficients(4);
1284  a5 = curvature_coefficients(5);
1285  }
1286  double den = (h*h + a1*a1 + a2*a2);
1287 
1288  for (int j=0; j<target_NP; ++j) {
1289  delta(j) = 0;
1290  }
1291 
1292  // 1/Sqrt[Det[G[r, s]]])*Div[Sqrt[Det[G[r, s]]]*Inv[G]*P
1293  if (_poly_order > 0 && _curvature_poly_order > 1) {
1294  delta(1) = (-a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4+(h*h+a1*a1)*a5))/den/den/(h*h);
1295  delta(2) = (-a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4+(h*h+a1*a1)*a5))/den/den/(h*h);
1296  }
1297  if (_poly_order > 1 && _curvature_poly_order > 0) {
1298  delta(3) = (h*h+a2*a2)/den/(h*h);
1299  delta(4) = -2*a1*a2/den/(h*h);
1300  delta(5) = (h*h+a1*a1)/den/(h*h);
1301  }
1302 
1303  // output component 0
1304  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
1305  for (int j=0; j<target_NP; ++j) {
1306  P_target_row(offset, j) = delta(j);
1307  }
1308  offset = getTargetOffsetIndexDevice(i, 1, 0, 0);
1309  for (int j=0; j<target_NP; ++j) {
1310  P_target_row(offset, j) = 0;
1311  }
1312 
1313  // output component 1
1314  offset = getTargetOffsetIndexDevice(i, 0, 1, 0);
1315  for (int j=0; j<target_NP; ++j) {
1316  P_target_row(offset, j) = 0;
1317  }
1318  offset = getTargetOffsetIndexDevice(i, 1, 1, 0);
1319  for (int j=0; j<target_NP; ++j) {
1320  P_target_row(offset, j) = delta(j);
1321  }
1322  });
1323  }
1324  } else if (_operations(i) == TargetOperation::GradientOfScalarPointEvaluation) {
1325  if (_reconstruction_space_rank == 0
1326  && (_polynomial_sampling_functional == PointSample
1327  || _polynomial_sampling_functional == ManifoldVectorPointSample)) {
1328  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1329 
1330  double h = _epsilons(target_index);
1331  double a1 = curvature_coefficients(1);
1332  double a2 = curvature_coefficients(2);
1333 
1334  double q1 = (h*h + a2*a2)/(h*h*h + h*a1*a1 + h*a2*a2);
1335  double q2 = -(a1*a2)/(h*h*h + h*a1*a1 + h*a2*a2);
1336  double q3 = (h*h + a1*a1)/(h*h*h + h*a1*a1 + h*a2*a2);
1337 
1338  double t1a = q1*1 + q2*0;
1339  double t2a = q1*0 + q2*1;
1340 
1341  double t1b = q2*1 + q3*0;
1342  double t2b = q2*0 + q3*1;
1343 
1344  // scaled
1345  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
1346  for (int j=0; j<target_NP; ++j) {
1347  P_target_row(offset, j) = 0;
1348  }
1349  if (_poly_order > 0 && _curvature_poly_order > 0) {
1350  P_target_row(offset, 1) = t1a + t2a;
1351  P_target_row(offset, 2) = 0;
1352  }
1353 
1354  offset = getTargetOffsetIndexDevice(i, 0, 1, 0);
1355  for (int j=0; j<target_NP; ++j) {
1356  P_target_row(offset, j) = 0;
1357  }
1358  if (_poly_order > 0 && _curvature_poly_order > 0) {
1359  P_target_row(offset, 1) = 0;
1360  P_target_row(offset, 2) = t1b + t2b;
1361  }
1362 
1363  });
1364  // staggered gradient w/ edge integrals performed by numerical integration
1365  // with a vector basis
1366  } else if (_reconstruction_space_rank == 1
1367  && _polynomial_sampling_functional
1368  == StaggeredEdgeIntegralSample) {
1369  compadre_kernel_assert_release((false) && "Functionality not yet available.");
1370 
1371  // staggered gradient w/ edge integrals known analytically, using a basis
1372  // of potentials
1373  } else if (_reconstruction_space_rank == 0
1374  && _polynomial_sampling_functional
1375  == StaggeredEdgeAnalyticGradientIntegralSample) {
1376  compadre_kernel_assert_release((false) && "Functionality not yet available.");
1377 
1378  } else {
1379  compadre_kernel_assert_release((false) && "Functionality not yet available.");
1380  }
1381  } else if (_operations(i) == TargetOperation::DivergenceOfVectorPointEvaluation) {
1382  // vector basis
1383  if (_reconstruction_space_rank == 1
1384  && _polynomial_sampling_functional == ManifoldVectorPointSample) {
1385  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1386 
1387  double h = _epsilons(target_index);
1388  double a1=0, a2=0, a3=0, a4=0, a5=0;
1389  if (_curvature_poly_order > 0) {
1390  a1 = curvature_coefficients(1);
1391  a2 = curvature_coefficients(2);
1392  }
1393  if (_curvature_poly_order > 1) {
1394  a3 = curvature_coefficients(3);
1395  a4 = curvature_coefficients(4);
1396  a5 = curvature_coefficients(5);
1397  }
1398  double den = (h*h + a1*a1 + a2*a2);
1399 
1400  // 1/Sqrt[Det[G[r, s]]])*Div[Sqrt[Det[G[r, s]]]*P
1401  // i.e. P recovers G^{-1}*grad of scalar
1402  double c0a = (a1*a3+a2*a4)/(h*den);
1403  double c1a = 1./h;
1404  double c2a = 0;
1405 
1406  double c0b = (a1*a4+a2*a5)/(h*den);
1407  double c1b = 0;
1408  double c2b = 1./h;
1409 
1410  // 1st input component
1411  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
1412  for (int j=0; j<target_NP; ++j) {
1413  P_target_row(offset, j) = 0;
1414  P_target_row(offset, target_NP + j) = 0;
1415  }
1416  P_target_row(offset, 0) = c0a;
1417  P_target_row(offset, 1) = c1a;
1418  P_target_row(offset, 2) = c2a;
1419 
1420  // 2nd input component
1421  offset = getTargetOffsetIndexDevice(i, 1, 0, 0);
1422  for (int j=0; j<target_NP; ++j) {
1423  P_target_row(offset, j) = 0;
1424  P_target_row(offset, target_NP + j) = 0;
1425  }
1426  P_target_row(offset, target_NP + 0) = c0b;
1427  P_target_row(offset, target_NP + 1) = c1b;
1428  P_target_row(offset, target_NP + 2) = c2b;
1429  });
1430  // scalar basis times number of components in the vector
1431  } else if (_reconstruction_space_rank == 0
1432  && _polynomial_sampling_functional == ManifoldVectorPointSample) {
1433  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1434 
1435  double h = _epsilons(target_index);
1436  double a1=0, a2=0, a3=0, a4=0, a5=0;
1437  if (_curvature_poly_order > 0) {
1438  a1 = curvature_coefficients(1);
1439  a2 = curvature_coefficients(2);
1440  }
1441  if (_curvature_poly_order > 1) {
1442  a3 = curvature_coefficients(3);
1443  a4 = curvature_coefficients(4);
1444  a5 = curvature_coefficients(5);
1445  }
1446  double den = (h*h + a1*a1 + a2*a2);
1447 
1448  // 1/Sqrt[Det[G[r, s]]])*Div[Sqrt[Det[G[r, s]]]*P
1449  // i.e. P recovers G^{-1}*grad of scalar
1450  double c0a = (a1*a3+a2*a4)/(h*den);
1451  double c1a = 1./h;
1452  double c2a = 0;
1453 
1454  double c0b = (a1*a4+a2*a5)/(h*den);
1455  double c1b = 0;
1456  double c2b = 1./h;
1457 
1458  // 1st input component
1459  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
1460  for (int j=0; j<target_NP; ++j) {
1461  P_target_row(offset, j) = 0;
1462  }
1463  P_target_row(offset, 0) = c0a;
1464  P_target_row(offset, 1) = c1a;
1465  P_target_row(offset, 2) = c2a;
1466 
1467  // 2nd input component
1468  offset = getTargetOffsetIndexDevice(i, 1, 0, 0);
1469  for (int j=0; j<target_NP; ++j) {
1470  P_target_row(offset, j) = 0;
1471  }
1472  P_target_row(offset, 0) = c0b;
1473  P_target_row(offset, 1) = c1b;
1474  P_target_row(offset, 2) = c2b;
1475  });
1476  // staggered divergence acting on vector polynomial space
1477  } else if (_reconstruction_space_rank == 1
1478  && _polynomial_sampling_functional
1479  == StaggeredEdgeIntegralSample) {
1480 
1481  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1482 
1483  double h = _epsilons(target_index);
1484  double a1=0, a2=0, a3=0, a4=0, a5=0;
1485  if (_curvature_poly_order > 0) {
1486  a1 = curvature_coefficients(1);
1487  a2 = curvature_coefficients(2);
1488  }
1489  if (_curvature_poly_order > 1) {
1490  a3 = curvature_coefficients(3);
1491  a4 = curvature_coefficients(4);
1492  a5 = curvature_coefficients(5);
1493  }
1494  double den = (h*h + a1*a1 + a2*a2);
1495 
1496  // 1/Sqrt[Det[G[r, s]]])*Div[Sqrt[Det[G[r, s]]]*Inv[G].P
1497  // i.e. P recovers grad of scalar
1498  double c0a = -a1*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5)/den/den/h;
1499  double c1a = (h*h+a2*a2)/den/h;
1500  double c2a = -a1*a2/den/h;
1501 
1502  double c0b = -a2*((h*h+a2*a2)*a3 - 2*a1*a2*a4 + (h*h+a1*a1)*a5)/den/den/h;
1503  double c1b = -a1*a2/den/h;
1504  double c2b = (h*h+a1*a1)/den/h;
1505 
1506  // 1st input component
1507  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
1508  for (int j=0; j<target_NP; ++j) {
1509  P_target_row(offset, j) = 0;
1510  P_target_row(offset, target_NP + j) = 0;
1511  }
1512  P_target_row(offset, 0) = c0a;
1513  P_target_row(offset, 1) = c1a;
1514  P_target_row(offset, 2) = c2a;
1515  P_target_row(offset, target_NP + 0) = c0b;
1516  P_target_row(offset, target_NP + 1) = c1b;
1517  P_target_row(offset, target_NP + 2) = c2b;
1518 
1519  });
1520  } else {
1521  compadre_kernel_assert_release((false) && "Functionality not yet available.");
1522  }
1523  } else if (_operations(i) == TargetOperation::GaussianCurvaturePointEvaluation) {
1524  double h = _epsilons(target_index);
1525 
1526  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int k) {
1527  XYZ relative_coord;
1528  if (k > 0) {
1529  for (int d=0; d<_dimensions-1; ++d) {
1530  relative_coord[d] = getTargetAuxiliaryCoordinate(target_index, k, d, &V);
1531  relative_coord[d] -= getTargetCoordinate(target_index, d, &V);
1532  }
1533  } else {
1534  for (int i=0; i<3; ++i) relative_coord[i] = 0;
1535  }
1536 
1537  int offset = getTargetOffsetIndexDevice(i, 0, 0, k);
1538  P_target_row(offset, 0) = GaussianCurvature(curvature_coefficients, h, relative_coord.x, relative_coord.y);
1539  });
1540  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
1541  } else if (_operations(i) == TargetOperation::CurlOfVectorPointEvaluation) {
1542  double h = _epsilons(target_index);
1543 
1544  Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember, num_evaluation_sites), [=] (const int k) {
1545  int alphax, alphay;
1546  XYZ relative_coord;
1547  if (k > 0) {
1548  for (int d=0; d<_dimensions-1; ++d) {
1549  relative_coord[d] = getTargetAuxiliaryCoordinate(target_index, k, d, &V);
1550  relative_coord[d] -= getTargetCoordinate(target_index, d, &V);
1551  }
1552  } else {
1553  for (int i=0; i<3; ++i) relative_coord[i] = 0;
1554  }
1555 
1556  int offset = getTargetOffsetIndexDevice(i, 0, 0, k);
1557  int index = 0;
1558  for (int n = 0; n <= _poly_order; n++){
1559  for (alphay = 0; alphay <= n; alphay++){
1560  alphax = n - alphay;
1561  P_target_row(offset, index) = SurfaceCurlOfScalar(curvature_coefficients, h, relative_coord.x, relative_coord.y, alphax, alphay, 0);
1562  index++;
1563  }
1564  }
1565 
1566  offset = getTargetOffsetIndexDevice(i, 0, 1, k);
1567  index = 0;
1568  for (int n = 0; n <= _poly_order; n++){
1569  for (alphay = 0; alphay <= n; alphay++){
1570  alphax = n - alphay;
1571  P_target_row(offset, index) = SurfaceCurlOfScalar(curvature_coefficients, h, relative_coord.x, relative_coord.y, alphax, alphay, 1);
1572  index++;
1573  }
1574  }
1575  });
1576  additional_evaluation_sites_handled = true; // additional non-target site evaluations handled
1577  } else if (_operations(i) == TargetOperation::ScalarFaceAverageEvaluation) {
1578  //printf("ScalarFaceAverageEvaluation\n");
1579  const double factorial[15] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
1580  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
1581 
1582 
1583  // Calculate basis matrix for NON MANIFOLD problems
1584  double cutoff_p = _epsilons(target_index);
1585  int alphax, alphay;
1586  double alphaf;
1587 
1588  // global dimension cannot be determined in a constexpr way, so we use a largest case scenario
1589  // of dimensions 3 for _global_dimension
1590  double triangle_coords[3/*_global_dimensions*/*3];
1591  for (int i=0; i<_global_dimensions*3; ++i) triangle_coords[i] = 0;
1592  // 3 is for # vertices in sub-triangle
1593  scratch_matrix_right_type triangle_coords_matrix(triangle_coords, _global_dimensions, 3);
1594 
1595  scratch_vector_type midpoint(delta.data(), _global_dimensions);
1596  getMidpointFromCellVertices(teamMember, midpoint, _target_extra_data, target_index, _global_dimensions /*dim*/);
1597  for (int j=0; j<_global_dimensions; ++j) {
1598  triangle_coords_matrix(j, 0) = midpoint(j);
1599  }
1600 
1601  size_t num_vertices = _target_extra_data.extent(1) / _global_dimensions;
1602  double reference_cell_area = 0.5;
1603  double entire_cell_area = 0.0;
1604  auto T=triangle_coords_matrix;
1605  for (size_t v=0; v<num_vertices; ++v) {
1606  int v1 = v;
1607  int v2 = (v+1) % num_vertices;
1608  for (int j=0; j<_global_dimensions; ++j) {
1609  triangle_coords_matrix(j,1) = _target_extra_data(target_index, v1*_global_dimensions+j) - triangle_coords_matrix(j,0);
1610  triangle_coords_matrix(j,2) = _target_extra_data(target_index, v2*_global_dimensions+j) - triangle_coords_matrix(j,0);
1611  }
1612  entire_cell_area += 0.5 * getAreaFromVectors(teamMember,
1613  Kokkos::subview(T, Kokkos::ALL(), 1), Kokkos::subview(T, Kokkos::ALL(), 2));
1614  }
1615 
1616 
1617  XYZ relative_coord;
1618  // loop over each two vertices
1619  for (size_t v=0; v<num_vertices; ++v) {
1620  int v1 = v;
1621  int v2 = (v+1) % num_vertices;
1622 
1623  for (int j=0; j<_global_dimensions; ++j) {
1624  triangle_coords_matrix(j,1) = _target_extra_data(target_index, v1*_global_dimensions+j) - triangle_coords_matrix(j,0);
1625  triangle_coords_matrix(j,2) = _target_extra_data(target_index, v2*_global_dimensions+j) - triangle_coords_matrix(j,0);
1626  }
1627  // triangle_coords now has:
1628  // (midpoint_x, midpoint_y, midpoint_z,
1629  // v1_x-midpoint_x, v1_y-midpoint_y, v1_z-midpoint_z,
1630  // v2_x-midpoint_x, v2_y-midpoint_y, v2_z-midpoint_z);
1631  for (int quadrature = 0; quadrature<_qm.getNumberOfQuadraturePoints(); ++quadrature) {
1632  double transformed_qp[3] = {0,0,0};
1633  for (int j=0; j<_global_dimensions; ++j) {
1634  for (int k=1; k<3; ++k) { // 3 is for # vertices in subtriangle
1635  transformed_qp[j] += T(j,k)*_qm.getSite(quadrature, k-1);
1636  }
1637  transformed_qp[j] += T(j,0);
1638  }
1639  // half the norm of the cross-product is the area of the triangle
1640  // so scaling is area / reference area (0.5) = the norm of the cross-product
1641  double sub_cell_area = 0.5 * getAreaFromVectors(teamMember,
1642  Kokkos::subview(T, Kokkos::ALL(), 1), Kokkos::subview(T, Kokkos::ALL(), 2));
1643  double scaling_factor = sub_cell_area / reference_cell_area;
1644 
1645  XYZ qp = XYZ(transformed_qp[0], transformed_qp[1], transformed_qp[2]);
1646  for (int j=0; j<2; ++j) {
1647  relative_coord[j] = convertGlobalToLocalCoordinate(qp,j,&V) - getTargetCoordinate(target_index,j,&V); // shift quadrature point by target site
1648  relative_coord[2] = 0;
1649  }
1650 
1651  int k = 0;
1652  const int start_index = 0;
1653  for (int n = start_index; n <= _poly_order; n++){
1654  for (alphay = 0; alphay <= n; alphay++){
1655  alphax = n - alphay;
1656  alphaf = factorial[alphax]*factorial[alphay];
1657  double val_to_sum = (scaling_factor * _qm.getWeight(quadrature)
1658  * std::pow(relative_coord.x/cutoff_p,alphax)
1659  *std::pow(relative_coord.y/cutoff_p,alphay)/alphaf) / entire_cell_area;
1660  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1661  if (quadrature==0 && v==0) P_target_row(offset, k) = val_to_sum;
1662  else P_target_row(offset, k) += val_to_sum;
1663  });
1664  k++;
1665  }
1666  }
1667  }
1668  }
1669 
1670  } else {
1671  compadre_kernel_assert_release((false) && "Functionality not yet available.");
1672  }
1673  } else if (_dimensions==2) { // 1D manifold in 2D problem
1674  if (_operations(i) == TargetOperation::GradientOfScalarPointEvaluation) {
1675  Kokkos::single(Kokkos::PerTeam(teamMember), [&] () {
1676  int offset = getTargetOffsetIndexDevice(i, 0, 0, 0);
1677  for (int j=0; j<target_NP; ++j) {
1678  P_target_row(offset, j) = 0;
1679  delta(j) = (j == 1) ? std::pow(_epsilons(target_index), -1) : 0;
1680  }
1681  for (int j=0; j<target_NP; ++j) {
1682  double v1 = delta(j)*G_inv(0,0);
1683  P_target_row(offset, j) = v1;
1684  }
1685  });
1686  } else {
1687  compadre_kernel_assert_release((false) && "Functionality not yet available.");
1688  }
1689  }
1690  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.");
1691  } // !operation_handled
1692  }
1693 }
1694 
1695 } // Compadre
1696 #endif
Compadre::DivergenceFreeVectorTaylorPolynomial
Divergence-free vector polynomial basis.
Definition: Compadre_Operators.hpp:82
Compadre::CurlCurlOfVectorPointEvaluation
Point evaluation of the curl of a curl of a vector (results in a vector)
Definition: Compadre_Operators.hpp:30
Compadre::GMLS::_operations
Kokkos::View< TargetOperation * > _operations
vector containing target functionals to be applied for reconstruction problem (device) ...
Definition: Compadre_GMLS.hpp:162
Compadre::CurlOfVectorPointEvaluation
Point evaluation of the curl of a vector (results in a vector)
Definition: Compadre_Operators.hpp:28
Compadre::getMidpointFromCellVertices
KOKKOS_INLINE_FUNCTION void getMidpointFromCellVertices(const member_type &teamMember, scratch_vector_type midpoint_storage, scratch_matrix_right_type cell_coordinates, const int cell_num, const int dim=3)
Definition: Compadre_Utilities.hpp:10
Compadre::VectorLaplacianPointEvaluation
Point evaluation of the laplacian of each component of a vector.
Definition: Compadre_Operators.hpp:20
Compadre::PartialYOfScalarPointEvaluation
Point evaluation of the partial with respect to y of a scalar.
Definition: Compadre_Operators.hpp:34
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:129
Compadre::ManifoldVectorPointSample
constexpr SamplingFunctional ManifoldVectorPointSample
Point evaluations of the entire vector source function (but on a manifold, so it includes a transform...
Definition: Compadre_Operators.hpp:148
Compadre::VectorPointEvaluation
Point evaluation of a vector (reconstructs entire vector at once, requiring a ReconstructionSpace hav...
Definition: Compadre_Operators.hpp:16
member_type
team_policy::member_type member_type
Definition: Compadre_Typedefs.hpp:48
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:50
Compadre::GMLS::_epsilons
Kokkos::View< double * > _epsilons
h supports determined through neighbor search (device)
Definition: Compadre_GMLS.hpp:86
Compadre::GMLS::getTargetOffsetIndexDevice
KOKKOS_INLINE_FUNCTION int getTargetOffsetIndexDevice(const int lro_num, const int input_component, const int output_component, const int additional_evaluation_local_index=0) const
Handles offset from operation input/output + extra evaluation sites.
Definition: Compadre_GMLS.hpp:562
Compadre_GMLS.hpp
Compadre::GMLS::calcGradientPij
KOKKOS_INLINE_FUNCTION void calcGradientPij(const member_type &teamMember, double *delta, double *thread_workspace, const int target_index, int neighbor_index, const double alpha, const int partial_direction, const int dimension, const int poly_order, bool specific_order_only, const scratch_matrix_right_type *V, const ReconstructionSpace reconstruction_space, const SamplingFunctional sampling_strategy, const int additional_evaluation_local_index=0) const
Evaluates the gradient of a polynomial basis under the Dirac Delta (pointwise) sampling function...
Definition: Compadre_GMLS_Basis.hpp:468
Compadre::GMLS::_sampling_multiplier
int _sampling_multiplier
actual dimension of the sampling functional e.g.
Definition: Compadre_GMLS.hpp:188
Compadre::GMLS::_reconstruction_space_rank
int _reconstruction_space_rank
actual rank of reconstruction basis
Definition: Compadre_GMLS.hpp:140
Compadre::GMLS::_curvature_poly_order
int _curvature_poly_order
order of basis for curvature reconstruction
Definition: Compadre_GMLS.hpp:122
Compadre::getAreaFromVectors
KOKKOS_INLINE_FUNCTION double getAreaFromVectors(const member_type &teamMember, view_type_1 v1, view_type_2 v2)
Definition: Compadre_Utilities.hpp:24
Compadre::GMLS::calcPij
KOKKOS_INLINE_FUNCTION void calcPij(const member_type &teamMember, double *delta, double *thread_workspace, const int target_index, int neighbor_index, const double alpha, const int dimension, const int poly_order, bool specific_order_only=false, const scratch_matrix_right_type *V=NULL, const ReconstructionSpace reconstruction_space=ReconstructionSpace::ScalarTaylorPolynomial, const SamplingFunctional sampling_strategy=PointSample, const int additional_evaluation_local_index=0) const
Evaluates the polynomial basis under a particular sampling function. Generally used to fill a row of ...
Definition: Compadre_GMLS_Basis.hpp:9
Compadre::GMLS::_data_sampling_multiplier
int _data_sampling_multiplier
effective dimension of the data sampling functional e.g.
Definition: Compadre_GMLS.hpp:192
Compadre_USER_StandardTargetFunctionals.hpp
Compadre::GMLS::_qm
Quadrature _qm
manages and calculates quadrature
Definition: Compadre_GMLS.hpp:281
Compadre::Quadrature::getWeight
KOKKOS_INLINE_FUNCTION double getWeight(const int index) const
Definition: Compadre_Quadrature.hpp:3276
Compadre::GMLS::convertGlobalToLocalCoordinate
KOKKOS_INLINE_FUNCTION double convertGlobalToLocalCoordinate(const XYZ global_coord, const int dim, const scratch_matrix_right_type *V) const
Returns a component of the local coordinate after transformation from global to local under the ortho...
Definition: Compadre_GMLS.hpp:528
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:151
Compadre::ScalarTaylorPolynomial
Scalar polynomial basis centered at the target site and scaled by sum of basis powers e...
Definition: Compadre_Operators.hpp:75
Compadre::Quadrature::getSite
KOKKOS_INLINE_FUNCTION double getSite(const int index, const int component) const
Definition: Compadre_Quadrature.hpp:3281
Compadre::GaussianCurvaturePointEvaluation
Point evaluation of Gaussian curvature.
Definition: Compadre_Operators.hpp:41
Compadre::GMLS::getNEvaluationSitesPerTarget
KOKKOS_INLINE_FUNCTION int getNEvaluationSitesPerTarget(const int target_index) const
(OPTIONAL) Returns number of additional evaluation sites for a particular target
Definition: Compadre_GMLS.hpp:454
Compadre::ScalarPointEvaluation
Point evaluation of a scalar.
Definition: Compadre_Operators.hpp:13
Compadre::GMLS::_initial_index_for_batch
int _initial_index_for_batch
initial index for current batch
Definition: Compadre_GMLS.hpp:217
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:821
Compadre::GMLS::getTargetCoordinate
KOKKOS_INLINE_FUNCTION double getTargetCoordinate(const int target_index, const int dim, const scratch_matrix_right_type *V=NULL) const
Returns one component of the target coordinate for a particular target.
Definition: Compadre_GMLS.hpp:471
Compadre::GMLS::computeTargetFunctionalsOnManifold
KOKKOS_INLINE_FUNCTION void computeTargetFunctionalsOnManifold(const member_type &teamMember, scratch_vector_type t1, scratch_vector_type t2, scratch_matrix_right_type P_target_row, scratch_matrix_right_type V, scratch_matrix_right_type G_inv, scratch_vector_type curvature_coefficients, scratch_vector_type curvature_gradients) const
Evaluates a polynomial basis with a target functional applied, using information from the manifold cu...
Definition: Compadre_GMLS_Targets.hpp:1001
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:154
Compadre::VectorOfScalarClonesTaylorPolynomial
Scalar basis reused as many times as there are components in the vector resulting in a much cheaper p...
Definition: Compadre_Operators.hpp:80
Compadre::GMLS::_polynomial_sampling_functional
const SamplingFunctional _polynomial_sampling_functional
polynomial sampling functional used to construct P matrix, set at GMLS class instantiation ...
Definition: Compadre_GMLS.hpp:152
scratch_matrix_right_type
Kokkos::View< double **, layout_right, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_matrix_right_type
Definition: Compadre_Typedefs.hpp:59
scratch_vector_type
Kokkos::View< double *, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_vector_type
Definition: Compadre_Typedefs.hpp:63
Compadre::GMLS::_target_extra_data
Kokkos::View< double **, layout_right > _target_extra_data
Extra data available to target operations (optional)
Definition: Compadre_GMLS.hpp:71
Compadre::GMLS::computeTargetFunctionals
KOKKOS_INLINE_FUNCTION void computeTargetFunctionals(const member_type &teamMember, scratch_vector_type t1, scratch_vector_type t2, scratch_matrix_right_type P_target_row) const
Evaluates a polynomial basis with a target functional applied to each member of the basis...
Definition: Compadre_GMLS_Targets.hpp:10
Compadre::LaplacianOfScalarPointEvaluation
Point evaluation of the laplacian of a scalar (could be on a manifold or not)
Definition: Compadre_Operators.hpp:18
Compadre::ChainedStaggeredLaplacianOfScalarPointEvaluation
Point evaluation of the chained staggered Laplacian acting on VectorTaylorPolynomial basis + Staggere...
Definition: Compadre_Operators.hpp:39
Compadre::DivergenceOfVectorPointEvaluation
Point evaluation of the divergence of a vector (results in a scalar)
Definition: Compadre_Operators.hpp:26
Compadre::GradientOfScalarPointEvaluation
Point evaluation of the gradient of a scalar.
Definition: Compadre_Operators.hpp:22
Compadre::GMLS::calcHessianPij
KOKKOS_INLINE_FUNCTION void calcHessianPij(const member_type &teamMember, double *delta, double *thread_workspace, const int target_index, int neighbor_index, const double alpha, const int partial_direction_1, const int partial_direction_2, const int dimension, const int poly_order, bool specific_order_only, const scratch_matrix_right_type *V, const ReconstructionSpace reconstruction_space, const SamplingFunctional sampling_strategy, const int additional_evaluation_local_index=0) const
Evaluates the Hessian of a polynomial basis under the Dirac Delta (pointwise) sampling function...
Definition: Compadre_GMLS_Basis.hpp:531
Compadre::GMLS::_curvature_support_operations
Kokkos::View< TargetOperation * > _curvature_support_operations
vector containing target functionals to be applied for curvature
Definition: Compadre_GMLS.hpp:159
Compadre::GMLS::computeCurvatureFunctionals
KOKKOS_INLINE_FUNCTION void computeCurvatureFunctionals(const member_type &teamMember, scratch_vector_type t1, scratch_vector_type t2, scratch_matrix_right_type P_target_row, const scratch_matrix_right_type *V, const local_index_type local_neighbor_index=-1) const
Evaluates a polynomial basis for the curvature with a gradient target functional applied.
Definition: Compadre_GMLS_Targets.hpp:951
d
import subprocess import os import re import math import sys import argparse d d d(20 *num_target_sites *pow(4, grid_num))
Compadre::GMLS::_basis_multiplier
int _basis_multiplier
dimension of the reconstructed function e.g.
Definition: Compadre_GMLS.hpp:183
Compadre::GMLS::_global_dimensions
int _global_dimensions
spatial dimension of the points, set at class instantiation only
Definition: Compadre_GMLS.hpp:128
Compadre::GMLS::_additional_evaluation_coordinates
Kokkos::View< double **, layout_right > _additional_evaluation_coordinates
(OPTIONAL) user provided additional coordinates for target operation evaluation (device) ...
Definition: Compadre_GMLS.hpp:106
Compadre::PartialZOfScalarPointEvaluation
Point evaluation of the partial with respect to z of a scalar.
Definition: Compadre_Operators.hpp:36
Compadre::PartialXOfScalarPointEvaluation
Point evaluation of the partial with respect to x of a scalar.
Definition: Compadre_Operators.hpp:32
Compadre::Quadrature::getNumberOfQuadraturePoints
KOKKOS_INLINE_FUNCTION int getNumberOfQuadraturePoints() const
Definition: Compadre_Quadrature.hpp:3248
Compadre::GMLS::_reconstruction_space
ReconstructionSpace _reconstruction_space
reconstruction space for GMLS problems, set at GMLS class instantiation
Definition: Compadre_GMLS.hpp:137
Compadre::GMLS::getTargetAuxiliaryCoordinate
KOKKOS_INLINE_FUNCTION double getTargetAuxiliaryCoordinate(const int target_index, const int additional_list_num, const int dim, const scratch_matrix_right_type *V=NULL) const
(OPTIONAL) Returns one component of the additional evaluation coordinates.
Definition: Compadre_GMLS.hpp:487
Compadre::PointSample
constexpr SamplingFunctional PointSample
Available sampling functionals.
Definition: Compadre_Operators.hpp:141
Compadre::GMLS::_dimensions
int _dimensions
dimension of the problem, set at class instantiation only
Definition: Compadre_GMLS.hpp:134
Compadre::GradientOfVectorPointEvaluation
Point evaluation of the gradient of a vector (results in a matrix, NOT CURRENTLY IMPLEMENTED) ...
Definition: Compadre_Operators.hpp:24
local_index_type
int local_index_type
Definition: Compadre_Typedefs.hpp:26
Compadre::GMLS::_local_dimensions
int _local_dimensions
dimension of the problem, set at class instantiation only. For manifolds, generally _global_dimension...
Definition: Compadre_GMLS.hpp:131
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:96
Compadre_Manifold_Functions.hpp
Compadre::VectorTaylorPolynomial
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...
Definition: Compadre_Operators.hpp:77
Compadre::VectorPointSample
constexpr SamplingFunctional VectorPointSample
Point evaluations of the entire vector source function.
Definition: Compadre_Operators.hpp:144
Compadre_USER_ManifoldTargetFunctionals.hpp
Compadre::GMLS::_poly_order
int _poly_order
order of basis for polynomial reconstruction
Definition: Compadre_GMLS.hpp:119
compadre_kernel_assert_debug
#define compadre_kernel_assert_debug(condition)
Definition: Compadre_Typedefs.hpp:151
Compadre::ScalarFaceAverageEvaluation
Average of values in a face of a cell using quadrature 2D in 3D problem, 1D in 2D problem...
Definition: Compadre_Operators.hpp:44
Generated on Thu Sep 26 2024 09:27:55 for Compadre by   doxygen 1.8.5