Compadre 1.5.5
Loading...
Searching...
No Matches
Compadre_Basis.hpp
Go to the documentation of this file.
1#ifndef _COMPADRE_BASIS_HPP_
2#define _COMPADRE_BASIS_HPP_
3
4#include "Compadre_GMLS.hpp"
5
6namespace Compadre {
7
8/*! \brief Evaluates the polynomial basis under a particular sampling function. Generally used to fill a row of P.
9 \param data [in] - GMLSBasisData struct
10 \param teamMember [in] - Kokkos::TeamPolicy member type (created by parallel_for)
11 \param delta [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large as the _basis_multipler*the dimension of the polynomial basis.
12 \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 _poly_order*the spatial dimension of the polynomial basis.
13 \param target_index [in] - target number
14 \param neighbor_index [in] - index of neighbor for this target with respect to local numbering [0,...,number of neighbors for target]
15 \param alpha [in] - double to determine convex combination of target and neighbor site at which to evaluate polynomials. (1-alpha)*neighbor + alpha*target
16 \param dimension [in] - spatial dimension of basis to evaluate. e.g. dimension two basis of order one is 1, x, y, whereas for dimension 3 it is 1, x, y, z
17 \param poly_order [in] - polynomial basis degree
18 \param specific_order_only [in] - boolean for only evaluating one degree of polynomial when true
19 \param V [in] - orthonormal basis matrix size _dimensions * _dimensions whose first _dimensions-1 columns are an approximation of the tangent plane
20 \param reconstruction_space [in] - space of polynomial that a sampling functional is to evaluate
21 \param sampling_strategy [in] - sampling functional specification
22 \param evaluation_site_local_index [in] - local index for evaluation sites (0 is target site)
23*/
24template <typename BasisData>
25KOKKOS_INLINE_FUNCTION
26void calcPij(const BasisData& data, const member_type& teamMember, double* delta, double* thread_workspace, const int target_index, int neighbor_index, const double alpha, const int dimension, const int poly_order, bool specific_order_only = false, const scratch_matrix_right_type* V = NULL, const ReconstructionSpace reconstruction_space = ReconstructionSpace::ScalarTaylorPolynomial, const SamplingFunctional polynomial_sampling_functional = PointSample, const int evaluation_site_local_index = 0) {
27/*
28 * This class is under two levels of hierarchical parallelism, so we
29 * do not put in any finer grain parallelism in this function
30 */
31 const int my_num_neighbors = data._pc._nla.getNumberOfNeighborsDevice(target_index);
32
33 // store precalculated factorials for speedup
34 const double factorial[15] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
35
36 int component = 0;
37 if (neighbor_index >= my_num_neighbors) {
38 component = neighbor_index / my_num_neighbors;
39 neighbor_index = neighbor_index % my_num_neighbors;
40 } else if (neighbor_index < 0) {
41 // -1 maps to 0 component
42 // -2 maps to 1 component
43 // -3 maps to 2 component
44 component = -(neighbor_index+1);
45 }
46
47 XYZ relative_coord;
48 if (neighbor_index > -1) {
49 // Evaluate at neighbor site
50 for (int i=0; i<dimension; ++i) {
51 // calculates (alpha*target+(1-alpha)*neighbor)-1*target = (alpha-1)*target + (1-alpha)*neighbor
52 relative_coord[i] = (alpha-1)*data._pc.getTargetCoordinate(target_index, i, V);
53 relative_coord[i] += (1-alpha)*data._pc.getNeighborCoordinate(target_index, neighbor_index, i, V);
54 }
55 } else if (evaluation_site_local_index > 0) {
56 // Extra evaluation site
57 for (int i=0; i<dimension; ++i) {
58 // evaluation_site_local_index is for evaluation site, which includes target site
59 // the -1 converts to the local index for ADDITIONAL evaluation sites
60 relative_coord[i] = data._additional_pc.getNeighborCoordinate(target_index, evaluation_site_local_index-1, i, V);
61 relative_coord[i] -= data._pc.getTargetCoordinate(target_index, i, V);
62 }
63 } else {
64 // Evaluate at the target site
65 for (int i=0; i<3; ++i) relative_coord[i] = 0;
66 }
67
68 // basis ActualReconstructionSpaceRank is 0 (evaluated like a scalar) and sampling functional is traditional
69 if ((polynomial_sampling_functional == PointSample ||
70 polynomial_sampling_functional == VectorPointSample ||
71 polynomial_sampling_functional == ManifoldVectorPointSample ||
72 polynomial_sampling_functional == VaryingManifoldVectorPointSample)&&
73 (reconstruction_space == ScalarTaylorPolynomial || reconstruction_space == VectorOfScalarClonesTaylorPolynomial)) {
74
75 double cutoff_p = data._epsilons(target_index);
76 const int start_index = specific_order_only ? poly_order : 0; // only compute specified order if requested
77
78 ScalarTaylorPolynomialBasis::evaluate(teamMember, delta, thread_workspace, dimension, poly_order, cutoff_p, relative_coord.x, relative_coord.y, relative_coord.z, start_index);
79
80 // basis ActualReconstructionSpaceRank is 1 (is a true vector basis) and sampling functional is traditional
81 } else if ((polynomial_sampling_functional == VectorPointSample ||
82 polynomial_sampling_functional == ManifoldVectorPointSample ||
83 polynomial_sampling_functional == VaryingManifoldVectorPointSample) &&
84 (reconstruction_space == VectorTaylorPolynomial)) {
85
86 const int dimension_offset = GMLS::getNP(data._poly_order, dimension, reconstruction_space);
87 double cutoff_p = data._epsilons(target_index);
88 const int start_index = specific_order_only ? poly_order : 0; // only compute specified order if requested
89
90 for (int d=0; d<dimension; ++d) {
91 if (d==component) {
92 ScalarTaylorPolynomialBasis::evaluate(teamMember, delta+component*dimension_offset, thread_workspace, dimension, poly_order, cutoff_p, relative_coord.x, relative_coord.y, relative_coord.z, start_index);
93 } else {
94 for (int n=0; n<dimension_offset; ++n) {
95 *(delta+d*dimension_offset+n) = 0;
96 }
97 }
98 }
99 } else if ((polynomial_sampling_functional == VectorPointSample) &&
100 (reconstruction_space == DivergenceFreeVectorTaylorPolynomial)) {
101 // Divergence free vector polynomial basis
102 double cutoff_p = data._epsilons(target_index);
103
104 DivergenceFreePolynomialBasis::evaluate(teamMember, delta, thread_workspace, dimension, poly_order, component, cutoff_p, relative_coord.x, relative_coord.y, relative_coord.z);
105 } else if (reconstruction_space == BernsteinPolynomial) {
106 // Bernstein vector polynomial basis
107 double cutoff_p = data._epsilons(target_index);
108
109 BernsteinPolynomialBasis::evaluate(teamMember, delta, thread_workspace, dimension, poly_order, component, cutoff_p, relative_coord.x, relative_coord.y, relative_coord.z);
110
111 } else if ((polynomial_sampling_functional == StaggeredEdgeAnalyticGradientIntegralSample) &&
112 (reconstruction_space == ScalarTaylorPolynomial)) {
113 double cutoff_p = data._epsilons(target_index);
114 const int start_index = specific_order_only ? poly_order : 0; // only compute specified order if requested
115 // basis is actually scalar with staggered sampling functional
116 ScalarTaylorPolynomialBasis::evaluate(teamMember, delta, thread_workspace, dimension, poly_order, cutoff_p, relative_coord.x, relative_coord.y, relative_coord.z, start_index, 0.0, -1.0);
117 relative_coord.x = 0;
118 relative_coord.y = 0;
119 relative_coord.z = 0;
120 ScalarTaylorPolynomialBasis::evaluate(teamMember, delta, thread_workspace, dimension, poly_order, cutoff_p, relative_coord.x, relative_coord.y, relative_coord.z, start_index, 1.0, 1.0);
121 } else if (polynomial_sampling_functional == StaggeredEdgeIntegralSample) {
122 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
123 if (data._problem_type == ProblemType::MANIFOLD) {
124 double cutoff_p = data._epsilons(target_index);
125 int alphax, alphay;
126 double alphaf;
127 const int start_index = specific_order_only ? poly_order : 0; // only compute specified order if requested
128
129 for (int quadrature = 0; quadrature<data._qm.getNumberOfQuadraturePoints(); ++quadrature) {
130
131 int i = 0;
132
133 XYZ tangent_quadrature_coord_2d;
134 XYZ quadrature_coord_2d;
135 for (int j=0; j<dimension; ++j) {
136 // calculates (alpha*target+(1-alpha)*neighbor)-1*target = (alpha-1)*target + (1-alpha)*neighbor
137 quadrature_coord_2d[j] = (data._qm.getSite(quadrature,0)-1)*data._pc.getTargetCoordinate(target_index, j, V);
138 quadrature_coord_2d[j] += (1-data._qm.getSite(quadrature,0))*data._pc.getNeighborCoordinate(target_index, neighbor_index, j, V);
139 tangent_quadrature_coord_2d[j] = data._pc.getTargetCoordinate(target_index, j, V);
140 tangent_quadrature_coord_2d[j] -= data._pc.getNeighborCoordinate(target_index, neighbor_index, j, V);
141 }
142 for (int j=0; j<data._basis_multiplier; ++j) {
143 for (int n = start_index; n <= poly_order; n++){
144 for (alphay = 0; alphay <= n; alphay++){
145 alphax = n - alphay;
146 alphaf = factorial[alphax]*factorial[alphay];
147
148 // local evaluation of vector [0,p] or [p,0]
149 double v0, v1;
150 v0 = (j==0) ? std::pow(quadrature_coord_2d.x/cutoff_p,alphax)
151 *std::pow(quadrature_coord_2d.y/cutoff_p,alphay)/alphaf : 0;
152 v1 = (j==0) ? 0 : std::pow(quadrature_coord_2d.x/cutoff_p,alphax)
153 *std::pow(quadrature_coord_2d.y/cutoff_p,alphay)/alphaf;
154
155 double dot_product = tangent_quadrature_coord_2d[0]*v0 + tangent_quadrature_coord_2d[1]*v1;
156
157 // multiply by quadrature weight
158 if (quadrature==0) {
159 *(delta+i) = dot_product * data._qm.getWeight(quadrature);
160 } else {
161 *(delta+i) += dot_product * data._qm.getWeight(quadrature);
162 }
163 i++;
164 }
165 }
166 }
167 }
168 } else {
169 // Calculate basis matrix for NON MANIFOLD problems
170 double cutoff_p = data._epsilons(target_index);
171 int alphax, alphay, alphaz;
172 double alphaf;
173 const int start_index = specific_order_only ? poly_order : 0; // only compute specified order if requested
174
175 for (int quadrature = 0; quadrature<data._qm.getNumberOfQuadraturePoints(); ++quadrature) {
176
177 int i = 0;
178
179 XYZ quadrature_coord_3d;
180 XYZ tangent_quadrature_coord_3d;
181 for (int j=0; j<dimension; ++j) {
182 // calculates (alpha*target+(1-alpha)*neighbor)-1*target = (alpha-1)*target + (1-alpha)*neighbor
183 quadrature_coord_3d[j] = (data._qm.getSite(quadrature,0)-1)*data._pc.getTargetCoordinate(target_index, j, NULL);
184 quadrature_coord_3d[j] += (1-data._qm.getSite(quadrature,0))*data._pc.getNeighborCoordinate(target_index, neighbor_index, j, NULL);
185 tangent_quadrature_coord_3d[j] = data._pc.getTargetCoordinate(target_index, j, NULL);
186 tangent_quadrature_coord_3d[j] -= data._pc.getNeighborCoordinate(target_index, neighbor_index, j, NULL);
187 }
188 for (int j=0; j<data._basis_multiplier; ++j) {
189 for (int n = start_index; n <= poly_order; n++) {
190 if (dimension == 3) {
191 for (alphaz = 0; alphaz <= n; alphaz++){
192 int s = n - alphaz;
193 for (alphay = 0; alphay <= s; alphay++){
194 alphax = s - alphay;
195 alphaf = factorial[alphax]*factorial[alphay]*factorial[alphaz];
196
197 // local evaluation of vector [p, 0, 0], [0, p, 0] or [0, 0, p]
198 double v0, v1, v2;
199 switch(j) {
200 case 1:
201 v0 = 0.0;
202 v1 = std::pow(quadrature_coord_3d.x/cutoff_p,alphax)
203 *std::pow(quadrature_coord_3d.y/cutoff_p,alphay)
204 *std::pow(quadrature_coord_3d.z/cutoff_p,alphaz)/alphaf;
205 v2 = 0.0;
206 break;
207 case 2:
208 v0 = 0.0;
209 v1 = 0.0;
210 v2 = std::pow(quadrature_coord_3d.x/cutoff_p,alphax)
211 *std::pow(quadrature_coord_3d.y/cutoff_p,alphay)
212 *std::pow(quadrature_coord_3d.z/cutoff_p,alphaz)/alphaf;
213 break;
214 default:
215 v0 = std::pow(quadrature_coord_3d.x/cutoff_p,alphax)
216 *std::pow(quadrature_coord_3d.y/cutoff_p,alphay)
217 *std::pow(quadrature_coord_3d.z/cutoff_p,alphaz)/alphaf;
218 v1 = 0.0;
219 v2 = 0.0;
220 break;
221 }
222
223 double dot_product = tangent_quadrature_coord_3d[0]*v0 + tangent_quadrature_coord_3d[1]*v1 + tangent_quadrature_coord_3d[2]*v2;
224
225 // multiply by quadrature weight
226 if (quadrature == 0) {
227 *(delta+i) = dot_product * data._qm.getWeight(quadrature);
228 } else {
229 *(delta+i) += dot_product * data._qm.getWeight(quadrature);
230 }
231 i++;
232 }
233 }
234 } else if (dimension == 2) {
235 for (alphay = 0; alphay <= n; alphay++){
236 alphax = n - alphay;
237 alphaf = factorial[alphax]*factorial[alphay];
238
239 // local evaluation of vector [p, 0] or [0, p]
240 double v0, v1;
241 switch(j) {
242 case 1:
243 v0 = 0.0;
244 v1 = std::pow(quadrature_coord_3d.x/cutoff_p,alphax)
245 *std::pow(quadrature_coord_3d.y/cutoff_p,alphay)/alphaf;
246 break;
247 default:
248 v0 = std::pow(quadrature_coord_3d.x/cutoff_p,alphax)
249 *std::pow(quadrature_coord_3d.y/cutoff_p,alphay)/alphaf;
250 v1 = 0.0;
251 break;
252 }
253
254 double dot_product = tangent_quadrature_coord_3d[0]*v0 + tangent_quadrature_coord_3d[1]*v1;
255
256 // multiply by quadrature weight
257 if (quadrature == 0) {
258 *(delta+i) = dot_product * data._qm.getWeight(quadrature);
259 } else {
260 *(delta+i) += dot_product * data._qm.getWeight(quadrature);
261 }
262 i++;
263 }
264 }
265 }
266 }
267 }
268 } // NON MANIFOLD PROBLEMS
269 });
270 } else if ((polynomial_sampling_functional == FaceNormalIntegralSample ||
271 polynomial_sampling_functional == EdgeTangentIntegralSample ||
272 polynomial_sampling_functional == FaceNormalAverageSample ||
273 polynomial_sampling_functional == EdgeTangentAverageSample) &&
274 reconstruction_space == VectorTaylorPolynomial) {
275
276 compadre_kernel_assert_debug(data._local_dimensions==2 &&
277 "FaceNormalIntegralSample, EdgeTangentIntegralSample, FaceNormalAverageSample, \
278 and EdgeTangentAverageSample only support 2d or 3d with 2d manifold");
279
280 compadre_kernel_assert_debug(data._qm.getDimensionOfQuadraturePoints()==1 \
281 && "Only 1d quadrature may be specified for edge integrals");
282
283 compadre_kernel_assert_debug(data._qm.getNumberOfQuadraturePoints()>=1 \
284 && "Quadrature points not generated");
285
286 compadre_kernel_assert_debug(data._source_extra_data.extent(0)>0 && "Extra data used but not set.");
287
288 compadre_kernel_assert_debug(!specific_order_only &&
289 "Sampling functional does not support specific_order_only");
290
291 double cutoff_p = data._epsilons(target_index);
292 auto global_neighbor_index = data._pc.getNeighborIndex(target_index, neighbor_index);
293 int alphax, alphay;
294 double alphaf;
295
296 /*
297 * requires quadrature points defined on an edge, not a target/source edge (spoke)
298 *
299 * data._source_extra_data will contain the endpoints (2 for 2D, 3 for 3D) and then the unit normals
300 * (e0_x, e0_y, e1_x, e1_y, n_x, n_y, t_x, t_y)
301 */
302
303 int quadrature_point_loop = data._qm.getNumberOfQuadraturePoints();
304
305 const int TWO = 2; // used because of # of vertices on an edge
306 double G_data[3*TWO]; // max(2,3)*TWO
307 double edge_coords[3*TWO];
308 for (int i=0; i<data._global_dimensions*TWO; ++i) G_data[i] = 0;
309 for (int i=0; i<data._global_dimensions*TWO; ++i) edge_coords[i] = 0;
310 // 2 is for # vertices on an edge
311 scratch_matrix_right_type G(G_data, data._global_dimensions, TWO);
312 scratch_matrix_right_type edge_coords_matrix(edge_coords, data._global_dimensions, TWO);
313
314 // neighbor coordinate is assumed to be midpoint
315 // could be calculated, but is correct for sphere
316 // and also for non-manifold problems
317 // uses given midpoint, rather than computing the midpoint from vertices
318 double radius = 0.0;
319 // this midpoint should lie on the sphere, or this will be the wrong radius
320 for (int j=0; j<data._global_dimensions; ++j) {
321 edge_coords_matrix(j, 0) = data._source_extra_data(global_neighbor_index, j);
322 edge_coords_matrix(j, 1) = data._source_extra_data(global_neighbor_index, data._global_dimensions + j) - edge_coords_matrix(j, 0);
323 radius += edge_coords_matrix(j, 0)*edge_coords_matrix(j, 0);
324 }
325 radius = std::sqrt(radius);
326
327 // edge_coords now has:
328 // (v0_x, v0_y, v1_x-v0_x, v1_y-v0_y)
329 auto E = edge_coords_matrix;
330
331 // get arc length of edge on manifold
332 double theta = 0.0;
333 if (data._problem_type == ProblemType::MANIFOLD) {
334 XYZ a = {E(0,0), E(1,0), E(2,0)};
335 XYZ b = {E(0,1)+E(0,0), E(1,1)+E(1,0), E(2,1)+E(2,0)};
336 double a_dot_b = a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
337 double norm_a_cross_b = getAreaFromVectors(teamMember, a, b);
338 theta = std::atan(norm_a_cross_b / a_dot_b);
339 }
340
341 // loop
342 double entire_edge_length = 0.0;
343 for (int quadrature = 0; quadrature<quadrature_point_loop; ++quadrature) {
344
345 double G_determinant = 1.0;
346 double scaled_transformed_qp[3] = {0,0,0};
347 double unscaled_transformed_qp[3] = {0,0,0};
348 for (int j=0; j<data._global_dimensions; ++j) {
349 unscaled_transformed_qp[j] += E(j,1)*data._qm.getSite(quadrature, 0);
350 // adds back on shift by endpoint
351 unscaled_transformed_qp[j] += E(j,0);
352 }
353
354 // project onto the sphere
355 if (data._problem_type == ProblemType::MANIFOLD) {
356 // unscaled_transformed_qp now lives on cell edge, but if on manifold,
357 // not directly on the sphere, just near by
358
359 // normalize to project back onto sphere
360 double transformed_qp_norm = 0;
361 for (int j=0; j<data._global_dimensions; ++j) {
362 transformed_qp_norm += unscaled_transformed_qp[j]*unscaled_transformed_qp[j];
363 }
364 transformed_qp_norm = std::sqrt(transformed_qp_norm);
365 // transformed_qp made unit length
366 for (int j=0; j<data._global_dimensions; ++j) {
367 scaled_transformed_qp[j] = unscaled_transformed_qp[j] * radius / transformed_qp_norm;
368 }
369
370 G_determinant = radius * theta;
371 XYZ qp = XYZ(scaled_transformed_qp[0], scaled_transformed_qp[1], scaled_transformed_qp[2]);
372 for (int j=0; j<data._local_dimensions; ++j) {
373 relative_coord[j] = data._pc.convertGlobalToLocalCoordinate(qp,j,*V)
374 - data._pc.getTargetCoordinate(target_index,j,V);
375 // shift quadrature point by target site
376 }
377 relative_coord[2] = 0;
378 } else { // not on a manifold, but still integrated
379 XYZ endpoints_difference = {E(0,1), E(1,1), 0};
380 G_determinant = data._pc.EuclideanVectorLength(endpoints_difference, 2);
381 for (int j=0; j<data._local_dimensions; ++j) {
382 relative_coord[j] = unscaled_transformed_qp[j]
383 - data._pc.getTargetCoordinate(target_index,j,V);
384 // shift quadrature point by target site
385 }
386 relative_coord[2] = 0;
387 }
388
389 // get normal or tangent direction (ambient)
390 XYZ direction;
391 if (polynomial_sampling_functional == FaceNormalIntegralSample
392 || polynomial_sampling_functional == FaceNormalAverageSample) {
393 for (int j=0; j<data._global_dimensions; ++j) {
394 // normal direction
395 direction[j] = data._source_extra_data(global_neighbor_index, 2*data._global_dimensions + j);
396 }
397 } else {
398 if (data._problem_type == ProblemType::MANIFOLD) {
399 // generate tangent from outward normal direction of the sphere and edge normal
400 XYZ k = {scaled_transformed_qp[0], scaled_transformed_qp[1], scaled_transformed_qp[2]};
401 double k_norm = std::sqrt(k[0]*k[0]+k[1]*k[1]+k[2]*k[2]);
402 k[0] = k[0]/k_norm; k[1] = k[1]/k_norm; k[2] = k[2]/k_norm;
403 XYZ n = {data._source_extra_data(global_neighbor_index, 2*data._global_dimensions + 0),
404 data._source_extra_data(global_neighbor_index, 2*data._global_dimensions + 1),
405 data._source_extra_data(global_neighbor_index, 2*data._global_dimensions + 2)};
406
407 double norm_k_cross_n = getAreaFromVectors(teamMember, k, n);
408 direction[0] = (k[1]*n[2] - k[2]*n[1]) / norm_k_cross_n;
409 direction[1] = (k[2]*n[0] - k[0]*n[2]) / norm_k_cross_n;
410 direction[2] = (k[0]*n[1] - k[1]*n[0]) / norm_k_cross_n;
411 } else {
412 for (int j=0; j<data._global_dimensions; ++j) {
413 // tangent direction
414 direction[j] = data._source_extra_data(global_neighbor_index, 3*data._global_dimensions + j);
415 }
416 }
417 }
418
419 // convert direction to local chart (for manifolds)
420 XYZ local_direction;
421 if (data._problem_type == ProblemType::MANIFOLD) {
422 for (int j=0; j<data._basis_multiplier; ++j) {
423 // Project ambient normal direction onto local chart basis as a local direction.
424 // Using V alone to provide vectors only gives tangent vectors at
425 // the target site. This could result in accuracy < 3rd order.
426 local_direction[j] = data._pc.convertGlobalToLocalCoordinate(direction,j,*V);
427 }
428 }
429
430 int i = 0;
431 for (int j=0; j<data._basis_multiplier; ++j) {
432 for (int n = 0; n <= poly_order; n++){
433 for (alphay = 0; alphay <= n; alphay++){
434 alphax = n - alphay;
435 alphaf = factorial[alphax]*factorial[alphay];
436
437 // local evaluation of vector [0,p] or [p,0]
438 double v0, v1;
439 v0 = (j==0) ? std::pow(relative_coord.x/cutoff_p,alphax)
440 *std::pow(relative_coord.y/cutoff_p,alphay)/alphaf : 0;
441 v1 = (j==1) ? std::pow(relative_coord.x/cutoff_p,alphax)
442 *std::pow(relative_coord.y/cutoff_p,alphay)/alphaf : 0;
443
444 // either n*v or t*v
445 double dot_product = 0.0;
446 if (data._problem_type == ProblemType::MANIFOLD) {
447 // alternate option for projection
448 dot_product = local_direction[0]*v0 + local_direction[1]*v1;
449 } else {
450 dot_product = direction[0]*v0 + direction[1]*v1;
451 }
452
453 // multiply by quadrature weight
454 if (quadrature==0) {
455 *(delta+i) = dot_product * data._qm.getWeight(quadrature) * G_determinant;
456 } else {
457 *(delta+i) += dot_product * data._qm.getWeight(quadrature) * G_determinant;
458 }
459 i++;
460 }
461 }
462 }
463 entire_edge_length += G_determinant * data._qm.getWeight(quadrature);
464 }
465 if (polynomial_sampling_functional == FaceNormalAverageSample ||
466 polynomial_sampling_functional == EdgeTangentAverageSample) {
467 int k = 0;
468 for (int j=0; j<data._basis_multiplier; ++j) {
469 for (int n = 0; n <= poly_order; n++){
470 for (alphay = 0; alphay <= n; alphay++){
471 *(delta+k) /= entire_edge_length;
472 k++;
473 }
474 }
475 }
476 }
477 } else if (polynomial_sampling_functional == CellAverageSample ||
478 polynomial_sampling_functional == CellIntegralSample) {
479
480 compadre_kernel_assert_debug(data._local_dimensions==2 &&
481 "CellAverageSample only supports 2d or 3d with 2d manifold");
482 auto global_neighbor_index = data._pc.getNeighborIndex(target_index, neighbor_index);
483 double cutoff_p = data._epsilons(target_index);
484 int alphax, alphay;
485 double alphaf;
486
487 double G_data[3*3]; //data._global_dimensions*3
488 double triangle_coords[3*3]; //data._global_dimensions*3
489 for (int i=0; i<data._global_dimensions*3; ++i) G_data[i] = 0;
490 for (int i=0; i<data._global_dimensions*3; ++i) triangle_coords[i] = 0;
491 // 3 is for # vertices in sub-triangle
492 scratch_matrix_right_type G(G_data, data._global_dimensions, 3);
493 scratch_matrix_right_type triangle_coords_matrix(triangle_coords, data._global_dimensions, 3);
494
495 // neighbor coordinate is assumed to be midpoint
496 // could be calculated, but is correct for sphere
497 // and also for non-manifold problems
498 // uses given midpoint, rather than computing the midpoint from vertices
499 double radius = 0.0;
500 // this midpoint should lie on the sphere, or this will be the wrong radius
501 for (int j=0; j<data._global_dimensions; ++j) {
502 // midpoint
503 triangle_coords_matrix(j, 0) = data._pc.getNeighborCoordinate(target_index, neighbor_index, j);
504 radius += triangle_coords_matrix(j, 0)*triangle_coords_matrix(j, 0);
505 }
506 radius = std::sqrt(radius);
507
508 // NaN in entry (data._global_dimensions) is a convention for indicating fewer vertices
509 // for this cell and NaN is checked by entry!=entry
510 int num_vertices = 0;
511 for (int j=0; j<data._source_extra_data.extent_int(1); ++j) {
512 auto val = data._source_extra_data(global_neighbor_index, j);
513 if (val != val) {
514 break;
515 } else {
516 num_vertices++;
517 }
518 }
519 num_vertices = num_vertices / data._global_dimensions;
520 auto T = triangle_coords_matrix;
521
522 // loop over each two vertices
523 // made for flat surfaces (either dim=2 or on a manifold)
524 double entire_cell_area = 0.0;
525 for (int v=0; v<num_vertices; ++v) {
526 int v1 = v;
527 int v2 = (v+1) % num_vertices;
528
529 for (int j=0; j<data._global_dimensions; ++j) {
530 triangle_coords_matrix(j,1) = data._source_extra_data(global_neighbor_index,
531 v1*data._global_dimensions+j)
532 - triangle_coords_matrix(j,0);
533 triangle_coords_matrix(j,2) = data._source_extra_data(global_neighbor_index,
534 v2*data._global_dimensions+j)
535 - triangle_coords_matrix(j,0);
536 }
537
538 // triangle_coords now has:
539 // (midpoint_x, midpoint_y, midpoint_z,
540 // v1_x-midpoint_x, v1_y-midpoint_y, v1_z-midpoint_z,
541 // v2_x-midpoint_x, v2_y-midpoint_y, v2_z-midpoint_z);
542 for (int quadrature = 0; quadrature<data._qm.getNumberOfQuadraturePoints(); ++quadrature) {
543 double unscaled_transformed_qp[3] = {0,0,0};
544 double scaled_transformed_qp[3] = {0,0,0};
545 for (int j=0; j<data._global_dimensions; ++j) {
546 for (int k=1; k<3; ++k) { // 3 is for # vertices in subtriangle
547 // uses vertex-midpoint as one direction
548 // and other vertex-midpoint as other direction
549 unscaled_transformed_qp[j] += T(j,k)*data._qm.getSite(quadrature, k-1);
550 }
551 // adds back on shift by midpoint
552 unscaled_transformed_qp[j] += T(j,0);
553 }
554
555 // project onto the sphere
556 double G_determinant = 1.0;
557 if (data._problem_type == ProblemType::MANIFOLD) {
558 // unscaled_transformed_qp now lives on cell, but if on manifold,
559 // not directly on the sphere, just near by
560
561 // normalize to project back onto sphere
562 double transformed_qp_norm = 0;
563 for (int j=0; j<data._global_dimensions; ++j) {
564 transformed_qp_norm += unscaled_transformed_qp[j]*unscaled_transformed_qp[j];
565 }
566 transformed_qp_norm = std::sqrt(transformed_qp_norm);
567 // transformed_qp made unit length
568 for (int j=0; j<data._global_dimensions; ++j) {
569 scaled_transformed_qp[j] = unscaled_transformed_qp[j] * radius / transformed_qp_norm;
570 }
571
572
573 // u_qp = midpoint + r_qp[1]*(v_1-midpoint) + r_qp[2]*(v_2-midpoint)
574 // s_qp = u_qp * radius / norm(u_qp) = radius * u_qp / norm(u_qp)
575 //
576 // so G(:,i) is \partial{s_qp}/ \partial{r_qp[i]}
577 // where r_qp is reference quadrature point (R^2 in 2D manifold in R^3)
578 //
579 // G(:,i) = radius * ( \partial{u_qp}/\partial{r_qp[i]} * (\sum_m u_qp[k]^2)^{-1/2}
580 // + u_qp * \partial{(\sum_m u_qp[k]^2)^{-1/2}}/\partial{r_qp[i]} )
581 //
582 // = radius * ( T(:,i)/norm(u_qp) + u_qp*(-1/2)*(\sum_m u_qp[k]^2)^{-3/2}
583 // *2*(\sum_k u_qp[k]*\partial{u_qp[k]}/\partial{r_qp[i]}) )
584 //
585 // = radius * ( T(:,i)/norm(u_qp) + u_qp*(-1/2)*(\sum_m u_qp[k]^2)^{-3/2}
586 // *2*(\sum_k u_qp[k]*T(k,i)) )
587 //
588 // NOTE: we do not multiply G by radius before determining area from vectors,
589 // so we multiply by radius**2 after calculation
590 double qp_norm_sq = transformed_qp_norm*transformed_qp_norm;
591 for (int j=0; j<data._global_dimensions; ++j) {
592 G(j,1) = T(j,1)/transformed_qp_norm;
593 G(j,2) = T(j,2)/transformed_qp_norm;
594 for (int k=0; k<data._global_dimensions; ++k) {
595 G(j,1) += unscaled_transformed_qp[j]*(-0.5)*std::pow(qp_norm_sq,-1.5)
596 *2*(unscaled_transformed_qp[k]*T(k,1));
597 G(j,2) += unscaled_transformed_qp[j]*(-0.5)*std::pow(qp_norm_sq,-1.5)
598 *2*(unscaled_transformed_qp[k]*T(k,2));
599 }
600 }
601 G_determinant = getAreaFromVectors(teamMember,
602 Kokkos::subview(G, Kokkos::ALL(), 1), Kokkos::subview(G, Kokkos::ALL(), 2));
603 G_determinant *= radius*radius;
604 XYZ qp = XYZ(scaled_transformed_qp[0], scaled_transformed_qp[1], scaled_transformed_qp[2]);
605 for (int j=0; j<data._local_dimensions; ++j) {
606 relative_coord[j] = data._pc.convertGlobalToLocalCoordinate(qp,j,*V)
607 - data._pc.getTargetCoordinate(target_index,j,V);
608 // shift quadrature point by target site
609 }
610 relative_coord[2] = 0;
611 } else {
612 G_determinant = getAreaFromVectors(teamMember,
613 Kokkos::subview(T, Kokkos::ALL(), 1), Kokkos::subview(T, Kokkos::ALL(), 2));
614 for (int j=0; j<data._local_dimensions; ++j) {
615 relative_coord[j] = unscaled_transformed_qp[j]
616 - data._pc.getTargetCoordinate(target_index,j,V);
617 // shift quadrature point by target site
618 }
619 relative_coord[2] = 0;
620 }
621
622 int k = 0;
623 compadre_kernel_assert_debug(!specific_order_only &&
624 "CellAverageSample does not support specific_order_only");
625 for (int n = 0; n <= poly_order; n++){
626 for (alphay = 0; alphay <= n; alphay++){
627 alphax = n - alphay;
628 alphaf = factorial[alphax]*factorial[alphay];
629 double val_to_sum = G_determinant * (data._qm.getWeight(quadrature)
630 * std::pow(relative_coord.x/cutoff_p,alphax)
631 * std::pow(relative_coord.y/cutoff_p,alphay) / alphaf);
632 if (quadrature==0 && v==0) *(delta+k) = val_to_sum;
633 else *(delta+k) += val_to_sum;
634 k++;
635 }
636 }
637 entire_cell_area += G_determinant * data._qm.getWeight(quadrature);
638 }
639 }
640 if (polynomial_sampling_functional == CellAverageSample) {
641 int k = 0;
642 for (int n = 0; n <= poly_order; n++){
643 for (alphay = 0; alphay <= n; alphay++){
644 *(delta+k) /= entire_cell_area;
645 k++;
646 }
647 }
648 }
649 } else {
650 compadre_kernel_assert_release((false) && "Sampling and basis space combination not defined.");
651 }
652}
653
654/*! \brief Evaluates the gradient of a polynomial basis under the Dirac Delta (pointwise) sampling function.
655 \param data [in] - GMLSBasisData struct
656 \param teamMember [in] - Kokkos::TeamPolicy member type (created by parallel_for)
657 \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.
658 \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 _poly_order*the spatial dimension of the polynomial basis.
659 \param target_index [in] - target number
660 \param neighbor_index [in] - index of neighbor for this target with respect to local numbering [0,...,number of neighbors for target]
661 \param alpha [in] - double to determine convex combination of target and neighbor site at which to evaluate polynomials. (1-alpha)*neighbor + alpha*target
662 \param partial_direction [in] - direction that partial is taken with respect to, e.g. 0 is x direction, 1 is y direction
663 \param dimension [in] - spatial dimension of basis to evaluate. e.g. dimension two basis of order one is 1, x, y, whereas for dimension 3 it is 1, x, y, z
664 \param poly_order [in] - polynomial basis degree
665 \param specific_order_only [in] - boolean for only evaluating one degree of polynomial when true
666 \param V [in] - orthonormal basis matrix size _dimensions * _dimensions whose first _dimensions-1 columns are an approximation of the tangent plane
667 \param reconstruction_space [in] - space of polynomial that a sampling functional is to evaluate
668 \param sampling_strategy [in] - sampling functional specification
669 \param evaluation_site_local_index [in] - local index for evaluation sites (0 is target site)
670*/
671template <typename BasisData>
672KOKKOS_INLINE_FUNCTION
673void calcGradientPij(const BasisData& data, const member_type& teamMember, double* delta, double* thread_workspace, const int target_index, int neighbor_index, const double alpha, const int partial_direction, const int dimension, const int poly_order, bool specific_order_only, const scratch_matrix_right_type* V, const ReconstructionSpace reconstruction_space, const SamplingFunctional polynomial_sampling_functional, const int evaluation_site_local_index = 0) {
674/*
675 * This class is under two levels of hierarchical parallelism, so we
676 * do not put in any finer grain parallelism in this function
677 */
678
679 const int my_num_neighbors = data._pc._nla.getNumberOfNeighborsDevice(target_index);
680
681 int component = 0;
682 if (neighbor_index >= my_num_neighbors) {
683 component = neighbor_index / my_num_neighbors;
684 neighbor_index = neighbor_index % my_num_neighbors;
685 } else if (neighbor_index < 0) {
686 // -1 maps to 0 component
687 // -2 maps to 1 component
688 // -3 maps to 2 component
689 component = -(neighbor_index+1);
690 }
691
692 // alpha corresponds to the linear combination of target_index and neighbor_index coordinates
693 // coordinate to evaluate = alpha*(target_index's coordinate) + (1-alpha)*(neighbor_index's coordinate)
694
695 // partial_direction - 0=x, 1=y, 2=z
696 XYZ relative_coord;
697 if (neighbor_index > -1) {
698 for (int i=0; i<dimension; ++i) {
699 // calculates (alpha*target+(1-alpha)*neighbor)-1*target = (alpha-1)*target + (1-alpha)*neighbor
700 relative_coord[i] = (alpha-1)*data._pc.getTargetCoordinate(target_index, i, V);
701 relative_coord[i] += (1-alpha)*data._pc.getNeighborCoordinate(target_index, neighbor_index, i, V);
702 }
703 } else if (evaluation_site_local_index > 0) {
704 for (int i=0; i<dimension; ++i) {
705 // evaluation_site_local_index is for evaluation site, which includes target site
706 // the -1 converts to the local index for ADDITIONAL evaluation sites
707 relative_coord[i] = data._additional_pc.getNeighborCoordinate(target_index, evaluation_site_local_index-1, i, V);
708 relative_coord[i] -= data._pc.getTargetCoordinate(target_index, i, V);
709 }
710 } else {
711 for (int i=0; i<3; ++i) relative_coord[i] = 0;
712 }
713
714 double cutoff_p = data._epsilons(target_index);
715 const int start_index = specific_order_only ? poly_order : 0; // only compute specified order if requested
716
717 if ((polynomial_sampling_functional == PointSample ||
718 polynomial_sampling_functional == VectorPointSample ||
719 polynomial_sampling_functional == ManifoldVectorPointSample ||
720 polynomial_sampling_functional == VaryingManifoldVectorPointSample) &&
721 (reconstruction_space == ScalarTaylorPolynomial || reconstruction_space == VectorOfScalarClonesTaylorPolynomial)) {
722
723 ScalarTaylorPolynomialBasis::evaluatePartialDerivative(teamMember, delta, thread_workspace, dimension, poly_order, partial_direction, cutoff_p, relative_coord.x, relative_coord.y, relative_coord.z, start_index);
724
725 } else if ((polynomial_sampling_functional == VectorPointSample) &&
726 (reconstruction_space == DivergenceFreeVectorTaylorPolynomial)) {
727
728 // Divergence free vector polynomial basis
729 DivergenceFreePolynomialBasis::evaluatePartialDerivative(teamMember, delta, thread_workspace, dimension, poly_order, component, partial_direction, cutoff_p, relative_coord.x, relative_coord.y, relative_coord.z);
730
731 } else if (reconstruction_space == BernsteinPolynomial) {
732
733 // Bernstein vector polynomial basis
734 BernsteinPolynomialBasis::evaluatePartialDerivative(teamMember, delta, thread_workspace, dimension, poly_order, component, partial_direction, cutoff_p, relative_coord.x, relative_coord.y, relative_coord.z);
735
736 } else {
737 compadre_kernel_assert_release((false) && "Sampling and basis space combination not defined.");
738 }
739}
740
741/*! \brief Evaluates the Hessian of a polynomial basis under the Dirac Delta (pointwise) sampling function.
742 \param data [in] - GMLSBasisData struct
743 \param teamMember [in] - Kokkos::TeamPolicy member type (created by parallel_for)
744 \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.
745 \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 _poly_order*the spatial dimension of the polynomial basis.
746 \param target_index [in] - target number
747 \param neighbor_index [in] - index of neighbor for this target with respect to local numbering [0,...,number of neighbors for target]
748 \param alpha [in] - double to determine convex combination of target and neighbor site at which to evaluate polynomials. (1-alpha)*neighbor + alpha*target
749 \param partial_direction_1 [in] - first direction that partial is taken with respect to, e.g. 0 is x direction, 1 is y direction
750 \param partial_direction_2 [in] - second direction that partial is taken with respect to, e.g. 0 is x direction, 1 is y direction
751 \param dimension [in] - spatial dimension of basis to evaluate. e.g. dimension two basis of order one is 1, x, y, whereas for dimension 3 it is 1, x, y, z
752 \param poly_order [in] - polynomial basis degree
753 \param specific_order_only [in] - boolean for only evaluating one degree of polynomial when true
754 \param V [in] - orthonormal basis matrix size _dimensions * _dimensions whose first _dimensions-1 columns are an approximation of the tangent plane
755 \param reconstruction_space [in] - space of polynomial that a sampling functional is to evaluate
756 \param sampling_strategy [in] - sampling functional specification
757 \param evaluation_site_local_index [in] - local index for evaluation sites (0 is target site)
758*/
759template <typename BasisData>
760KOKKOS_INLINE_FUNCTION
761void calcHessianPij(const BasisData& data, const member_type& teamMember, double* delta, double* thread_workspace, const int target_index, int neighbor_index, const double alpha, const int partial_direction_1, const int partial_direction_2, const int dimension, const int poly_order, bool specific_order_only, const scratch_matrix_right_type* V, const ReconstructionSpace reconstruction_space, const SamplingFunctional polynomial_sampling_functional, const int evaluation_site_local_index = 0) {
762/*
763 * This class is under two levels of hierarchical parallelism, so we
764 * do not put in any finer grain parallelism in this function
765 */
766
767 const int my_num_neighbors = data._pc._nla.getNumberOfNeighborsDevice(target_index);
768
769 int component = 0;
770 if (neighbor_index >= my_num_neighbors) {
771 component = neighbor_index / my_num_neighbors;
772 neighbor_index = neighbor_index % my_num_neighbors;
773 } else if (neighbor_index < 0) {
774 // -1 maps to 0 component
775 // -2 maps to 1 component
776 // -3 maps to 2 component
777 component = -(neighbor_index+1);
778 }
779
780 // alpha corresponds to the linear combination of target_index and neighbor_index coordinates
781 // coordinate to evaluate = alpha*(target_index's coordinate) + (1-alpha)*(neighbor_index's coordinate)
782
783 // partial_direction - 0=x, 1=y, 2=z
784 XYZ relative_coord;
785 if (neighbor_index > -1) {
786 for (int i=0; i<dimension; ++i) {
787 // calculates (alpha*target+(1-alpha)*neighbor)-1*target = (alpha-1)*target + (1-alpha)*neighbor
788 relative_coord[i] = (alpha-1)*data._pc.getTargetCoordinate(target_index, i, V);
789 relative_coord[i] += (1-alpha)*data._pc.getNeighborCoordinate(target_index, neighbor_index, i, V);
790 }
791 } else if (evaluation_site_local_index > 0) {
792 for (int i=0; i<dimension; ++i) {
793 // evaluation_site_local_index is for evaluation site, which includes target site
794 // the -1 converts to the local index for ADDITIONAL evaluation sites
795 relative_coord[i] = data._additional_pc.getNeighborCoordinate(target_index, evaluation_site_local_index-1, i, V);
796 relative_coord[i] -= data._pc.getTargetCoordinate(target_index, i, V);
797 }
798 } else {
799 for (int i=0; i<3; ++i) relative_coord[i] = 0;
800 }
801
802 double cutoff_p = data._epsilons(target_index);
803
804 if ((polynomial_sampling_functional == PointSample ||
805 polynomial_sampling_functional == VectorPointSample ||
806 polynomial_sampling_functional == ManifoldVectorPointSample ||
807 polynomial_sampling_functional == VaryingManifoldVectorPointSample) &&
808 (reconstruction_space == ScalarTaylorPolynomial || reconstruction_space == VectorOfScalarClonesTaylorPolynomial)) {
809
810 ScalarTaylorPolynomialBasis::evaluateSecondPartialDerivative(teamMember, delta, thread_workspace, dimension, poly_order, partial_direction_1, partial_direction_2, cutoff_p, relative_coord.x, relative_coord.y, relative_coord.z);
811
812 } else if ((polynomial_sampling_functional == VectorPointSample) &&
813 (reconstruction_space == DivergenceFreeVectorTaylorPolynomial)) {
814
815 DivergenceFreePolynomialBasis::evaluateSecondPartialDerivative(teamMember, delta, thread_workspace, dimension, poly_order, component, partial_direction_1, partial_direction_2, cutoff_p, relative_coord.x, relative_coord.y, relative_coord.z);
816
817 } else if (reconstruction_space == BernsteinPolynomial) {
818
819 BernsteinPolynomialBasis::evaluateSecondPartialDerivative(teamMember, delta, thread_workspace, dimension, poly_order, component, partial_direction_1, partial_direction_2, cutoff_p, relative_coord.x, relative_coord.y, relative_coord.z);
820
821 } else {
822 compadre_kernel_assert_release((false) && "Sampling and basis space combination not defined.");
823 }
824}
825
826/*! \brief Fills the _P matrix with either P or P*sqrt(w)
827 \param data [in] - GMLSBasisData struct
828 \param teamMember [in] - Kokkos::TeamPolicy member type (created by parallel_for)
829 \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.
830 \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 _poly_order*the spatial dimension of the polynomial basis.
831 \param P [out] - 2D Kokkos View which will contain evaluation of sampling functional on polynomial basis for each neighbor the target has (stored column major)
832 \param w [out] - 1D Kokkos View which will contain weighting kernel values for the target with each neighbor if weight_p = true
833 \param dimension [in] - spatial dimension of basis to evaluate. e.g. dimension two basis of order one is 1, x, y, whereas for dimension 3 it is 1, x, y, z
834 \param polynomial_order [in] - polynomial basis degree
835 \param weight_p [in] - boolean whether to fill w with kernel weights
836 \param V [in] - orthonormal basis matrix size _dimensions * _dimensions whose first _dimensions-1 columns are an approximation of the tangent plane
837 \param reconstruction_space [in] - space of polynomial that a sampling functional is to evaluate
838 \param sampling_strategy [in] - sampling functional specification
839*/
840template <typename BasisData>
841KOKKOS_INLINE_FUNCTION
842void createWeightsAndP(const BasisData& data, const member_type& teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P, scratch_vector_type w, const int dimension, int polynomial_order, bool weight_p = false, scratch_matrix_right_type* V = NULL, const ReconstructionSpace reconstruction_space = ReconstructionSpace::ScalarTaylorPolynomial, const SamplingFunctional polynomial_sampling_functional = PointSample) {
843 /*
844 * Creates sqrt(W)*P
845 */
846 const int target_index = data._initial_index_for_batch + teamMember.league_rank();
847// printf("specific order: %d\n", specific_order);
848// {
849// const int storage_size = (specific_order > 0) ? GMLS::getNP(specific_order, dimension)-GMLS::getNP(specific_order-1, dimension) : GMLS::getNP(data._poly_order, dimension);
850// printf("storage size: %d\n", storage_size);
851// }
852// printf("weight_p: %d\n", weight_p);
853
854 // not const b.c. of gcc 7.2 issue
855 int my_num_neighbors = data._pc._nla.getNumberOfNeighborsDevice(target_index);
856
857 // storage_size needs to change based on the size of the basis
858 int storage_size = GMLS::getNP(polynomial_order, dimension, reconstruction_space);
859 storage_size *= data._basis_multiplier;
860 for (int j = 0; j < delta.extent_int(0); ++j) {
861 delta(j) = 0;
862 }
863 for (int j = 0; j < thread_workspace.extent_int(0); ++j) {
864 thread_workspace(j) = 0;
865 }
866 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,data._pc._nla.getNumberOfNeighborsDevice(target_index)),
867 [&] (const int i) {
868
869 for (int d=0; d<data._sampling_multiplier; ++d) {
870 // in 2d case would use distance between SVD coordinates
871
872 // ignores V when calculating weights from a point, i.e. uses actual point values
873 double r;
874
875 // coefficient muliplied by relative distance (allows for mid-edge weighting if applicable)
876 double alpha_weight = 1;
877 if (data._polynomial_sampling_functional==StaggeredEdgeIntegralSample
878 || data._polynomial_sampling_functional==StaggeredEdgeAnalyticGradientIntegralSample) {
879 alpha_weight = 0.5;
880 }
881
882 // get Euchlidean distance of scaled relative coordinate from the origin
883 if (V==NULL) {
884 r = data._pc.EuclideanVectorLength(data._pc.getRelativeCoord(target_index, i, dimension) * alpha_weight, dimension);
885 } else {
886 r = data._pc.EuclideanVectorLength(data._pc.getRelativeCoord(target_index, i, dimension, V) * alpha_weight, dimension);
887 }
888
889 // generate weight vector from distances and window sizes
890 w(i+my_num_neighbors*d) = GMLS::Wab(r, data._epsilons(target_index), data._weighting_type, data._weighting_p, data._weighting_n);
891
892 calcPij<BasisData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, i + d*my_num_neighbors, 0 /*alpha*/, dimension, polynomial_order, false /*bool on only specific order*/, V, reconstruction_space, polynomial_sampling_functional);
893
894 if (weight_p) {
895 for (int j = 0; j < storage_size; ++j) {
896 // no need to convert offsets to global indices because the sum will never be large
897 P(i+my_num_neighbors*d, j) = delta[j] * std::sqrt(w(i+my_num_neighbors*d));
898 compadre_kernel_assert_extreme_debug(delta[j]==delta[j] && "NaN in sqrt(W)*P matrix.");
899 }
900
901 } else {
902 for (int j = 0; j < storage_size; ++j) {
903 // no need to convert offsets to global indices because the sum will never be large
904 P(i+my_num_neighbors*d, j) = delta[j];
905
906 compadre_kernel_assert_extreme_debug(delta[j]==delta[j] && "NaN in P matrix.");
907 }
908 }
909 }
910 });
911 teamMember.team_barrier();
912// Kokkos::single(Kokkos::PerTeam(teamMember), [=] () {
913// for (int k=0; k<data._pc._nla.getNumberOfNeighborsDevice(target_index); k++) {
914// for (int l=0; l<_NP; l++) {
915// printf("%i %i %0.16f\n", k, l, P(k,l) );// << " " << l << " " << R(k,l) << std::endl;
916// }
917// }
918// });
919}
920
921/*! \brief Fills the _P matrix with P*sqrt(w) for use in solving for curvature
922
923 Uses _curvature_poly_order as the polynomial order of the basis
924
925 \param data [in] - GMLSBasisData struct
926 \param teamMember [in] - Kokkos::TeamPolicy member type (created by parallel_for)
927 \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
928s_multipler*the dimension of the polynomial basis.
929 \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
930_order*the spatial dimension of the polynomial basis.
931 \param P [out] - 2D Kokkos View which will contain evaluation of sampling functional on polynomial basis for each neighbor the target has (stored column major)
932 \param w [out] - 1D Kokkos View which will contain weighting kernel values for the target with each neighbor if weight_p = true
933 \param dimension [in] - spatial dimension of basis to evaluate. e.g. dimension two basis of order one is 1, x, y, whereas for dimension 3 it is 1, x, y, z
934 \param only_specific_order [in] - boolean for only evaluating one degree of polynomial when true
935 \param V [in] - orthonormal basis matrix size _dimensions * _dimensions whose first _dimensions-1 columns are an approximation of the tangent plane
936*/
937template <typename BasisData>
938KOKKOS_INLINE_FUNCTION
939void createWeightsAndPForCurvature(const BasisData& data, const member_type& teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P, scratch_vector_type w, const int dimension, bool only_specific_order, scratch_matrix_right_type* V = NULL) {
940/*
941 * This function has two purposes
942 * 1.) Used to calculate specifically for 1st order polynomials, from which we can reconstruct a tangent plane
943 * 2.) Used to calculate a polynomial of data._curvature_poly_order, which we use to calculate curvature of the manifold
944 */
945
946 const int target_index = data._initial_index_for_batch + teamMember.league_rank();
947 int storage_size = only_specific_order ? GMLS::getNP(1, dimension)-GMLS::getNP(0, dimension) : GMLS::getNP(data._curvature_poly_order, dimension);
948 for (int j = 0; j < delta.extent_int(0); ++j) {
949 delta(j) = 0;
950 }
951 for (int j = 0; j < thread_workspace.extent_int(0); ++j) {
952 thread_workspace(j) = 0;
953 }
954 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,data._pc._nla.getNumberOfNeighborsDevice(target_index)),
955 [&] (const int i) {
956
957 // ignores V when calculating weights from a point, i.e. uses actual point values
958 double r;
959
960 // get Euclidean distance of scaled relative coordinate from the origin
961 if (V==NULL) {
962 r = data._pc.EuclideanVectorLength(data._pc.getRelativeCoord(target_index, i, dimension), dimension);
963 } else {
964 r = data._pc.EuclideanVectorLength(data._pc.getRelativeCoord(target_index, i, dimension, V), dimension);
965 }
966
967 // generate weight vector from distances and window sizes
968 if (only_specific_order) {
969 w(i) = GMLS::Wab(r, data._epsilons(target_index), data._curvature_weighting_type, data._curvature_weighting_p, data._curvature_weighting_n);
970 calcPij<BasisData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, i, 0 /*alpha*/, dimension, 1, true /*bool on only specific order*/);
971 } else {
972 w(i) = GMLS::Wab(r, data._epsilons(target_index), data._curvature_weighting_type, data._curvature_weighting_p, data._curvature_weighting_n);
973 calcPij<BasisData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, i, 0 /*alpha*/, dimension, data._curvature_poly_order, false /*bool on only specific order*/, V);
974 }
975
976 for (int j = 0; j < storage_size; ++j) {
977 P(i, j) = delta[j] * std::sqrt(w(i));
978 }
979
980 });
981 teamMember.team_barrier();
982}
983
984} // Compadre
985#endif
#define compadre_kernel_assert_debug(condition)
Kokkos::View< double *, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_vector_type
#define compadre_kernel_assert_release(condition)
compadre_kernel_assert_release is similar to compadre_assert_release, but is a call on the device,...
team_policy::member_type member_type
Kokkos::View< double **, layout_right, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_matrix_right_type
if(some_conditions_for_a_user_defined_operation)
static KOKKOS_INLINE_FUNCTION int getNP(const int m, const int dimension=3, const ReconstructionSpace r_space=ReconstructionSpace::ScalarTaylorPolynomial)
Returns size of the basis for a given polynomial order and dimension General to dimension 1....
KOKKOS_INLINE_FUNCTION void evaluate(const member_type &teamMember, double *delta, double *workspace, const int dimension, const int max_degree, const int component, const double h, const double x, const double y, const double z, const int starting_order=0, const double weight_of_original_value=0.0, const double weight_of_new_value=1.0)
Evaluates the Bernstein polynomial basis delta[j] = weight_of_original_value * delta[j] + weight_of_n...
KOKKOS_INLINE_FUNCTION void evaluate(const member_type &teamMember, double *delta, double *workspace, const int dimension, const int max_degree, const int component, const double h, const double x, const double y, const double z, const int starting_order=0, const double weight_of_original_value=0.0, const double weight_of_new_value=1.0)
Evaluates the divergence-free polynomial basis delta[j] = weight_of_original_value * delta[j] + weigh...
KOKKOS_INLINE_FUNCTION void evaluate(const member_type &teamMember, double *delta, double *workspace, const int dimension, const int max_degree, const double h, const double x, const double y, const double z, const int starting_order=0, const double weight_of_original_value=0.0, const double weight_of_new_value=1.0)
Evaluates the scalar Taylor polynomial basis delta[j] = weight_of_original_value * delta[j] + weight_...
@ MANIFOLD
Solve GMLS problem on a manifold (will use QR or SVD to solve the resultant GMLS problem dependent on...
constexpr SamplingFunctional VaryingManifoldVectorPointSample
For integrating polynomial dotted with normal over an edge.
constexpr SamplingFunctional FaceNormalIntegralSample
For integrating polynomial dotted with normal over an edge.
constexpr SamplingFunctional StaggeredEdgeIntegralSample
Samples consist of the result of integrals of a vector dotted with the tangent along edges between ne...
KOKKOS_INLINE_FUNCTION void createWeightsAndP(const BasisData &data, const member_type &teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P, scratch_vector_type w, const int dimension, int polynomial_order, bool weight_p=false, scratch_matrix_right_type *V=NULL, const ReconstructionSpace reconstruction_space=ReconstructionSpace::ScalarTaylorPolynomial, const SamplingFunctional polynomial_sampling_functional=PointSample)
Fills the _P matrix with either P or P*sqrt(w)
constexpr SamplingFunctional PointSample
Available sampling functionals.
constexpr SamplingFunctional StaggeredEdgeAnalyticGradientIntegralSample
Analytical integral of a gradient source vector is just a difference of the scalar source at neighbor...
constexpr SamplingFunctional CellAverageSample
For polynomial integrated on cells.
KOKKOS_INLINE_FUNCTION void createWeightsAndPForCurvature(const BasisData &data, const member_type &teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P, scratch_vector_type w, const int dimension, bool only_specific_order, scratch_matrix_right_type *V=NULL)
Fills the _P matrix with P*sqrt(w) for use in solving for curvature.
constexpr SamplingFunctional CellIntegralSample
For polynomial integrated on cells.
constexpr SamplingFunctional FaceNormalAverageSample
For polynomial dotted with normal on edge.
KOKKOS_INLINE_FUNCTION void calcPij(const BasisData &data, const member_type &teamMember, double *delta, double *thread_workspace, const int target_index, int neighbor_index, const double alpha, const int dimension, const int poly_order, bool specific_order_only=false, const scratch_matrix_right_type *V=NULL, const ReconstructionSpace reconstruction_space=ReconstructionSpace::ScalarTaylorPolynomial, const SamplingFunctional polynomial_sampling_functional=PointSample, const int evaluation_site_local_index=0)
Evaluates the polynomial basis under a particular sampling function. Generally used to fill a row of ...
constexpr SamplingFunctional ManifoldVectorPointSample
Point evaluations of the entire vector source function (but on a manifold, so it includes a transform...
constexpr SamplingFunctional VectorPointSample
Point evaluations of the entire vector source function.
KOKKOS_INLINE_FUNCTION void calcHessianPij(const BasisData &data, const member_type &teamMember, double *delta, double *thread_workspace, const int target_index, int neighbor_index, const double alpha, const int partial_direction_1, const int partial_direction_2, const int dimension, const int poly_order, bool specific_order_only, const scratch_matrix_right_type *V, const ReconstructionSpace reconstruction_space, const SamplingFunctional polynomial_sampling_functional, const int evaluation_site_local_index=0)
Evaluates the Hessian of a polynomial basis under the Dirac Delta (pointwise) sampling function.
constexpr SamplingFunctional EdgeTangentAverageSample
For polynomial dotted with tangent.
KOKKOS_INLINE_FUNCTION void calcGradientPij(const BasisData &data, const member_type &teamMember, double *delta, double *thread_workspace, const int target_index, int neighbor_index, const double alpha, const int partial_direction, const int dimension, const int poly_order, bool specific_order_only, const scratch_matrix_right_type *V, const ReconstructionSpace reconstruction_space, const SamplingFunctional polynomial_sampling_functional, const int evaluation_site_local_index=0)
Evaluates the gradient of a polynomial basis under the Dirac Delta (pointwise) sampling function.
constexpr SamplingFunctional EdgeTangentIntegralSample
For integrating polynomial dotted with tangent over an edge.
ReconstructionSpace
Space in which to reconstruct polynomial.
@ DivergenceFreeVectorTaylorPolynomial
Divergence-free vector polynomial basis.
@ BernsteinPolynomial
Bernstein polynomial basis.
@ VectorTaylorPolynomial
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...
@ ScalarTaylorPolynomial
Scalar polynomial basis centered at the target site and scaled by sum of basis powers e....
@ VectorOfScalarClonesTaylorPolynomial
Scalar basis reused as many times as there are components in the vector resulting in a much cheaper p...
KOKKOS_INLINE_FUNCTION double getAreaFromVectors(const member_type &teamMember, view_type_1 v1, view_type_2 v2)