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