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) {
31 const int my_num_neighbors = data._pc._nla.getNumberOfNeighborsDevice(target_index);
34 const double factorial[15] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
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) {
44 component = -(neighbor_index+1);
48 if (neighbor_index > -1) {
50 for (
int i=0; i<dimension; ++i) {
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);
55 }
else if (evaluation_site_local_index > 0) {
57 for (
int i=0; i<dimension; ++i) {
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);
65 for (
int i=0; i<3; ++i) relative_coord[i] = 0;
69 if ((polynomial_sampling_functional ==
PointSample ||
75 double cutoff_p = data._epsilons(target_index);
76 const int start_index = specific_order_only ? poly_order : 0;
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;
90 for (
int d=0; d<dimension; ++d) {
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);
94 for (
int n=0; n<dimension_offset; ++n) {
95 *(delta+d*dimension_offset+n) = 0;
102 double cutoff_p = data._epsilons(target_index);
107 double cutoff_p = data._epsilons(target_index);
113 double cutoff_p = data._epsilons(target_index);
114 const int start_index = specific_order_only ? poly_order : 0;
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);
122 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
124 double cutoff_p = data._epsilons(target_index);
127 const int start_index = specific_order_only ? poly_order : 0;
129 for (int quadrature = 0; quadrature<data._qm.getNumberOfQuadraturePoints(); ++quadrature) {
133 XYZ tangent_quadrature_coord_2d;
134 XYZ quadrature_coord_2d;
135 for (int j=0; j<dimension; ++j) {
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);
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++){
146 alphaf = factorial[alphax]*factorial[alphay];
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;
155 double dot_product = tangent_quadrature_coord_2d[0]*v0 + tangent_quadrature_coord_2d[1]*v1;
159 *(delta+i) = dot_product * data._qm.getWeight(quadrature);
161 *(delta+i) += dot_product * data._qm.getWeight(quadrature);
170 double cutoff_p = data._epsilons(target_index);
171 int alphax, alphay, alphaz;
173 const int start_index = specific_order_only ? poly_order : 0;
175 for (
int quadrature = 0; quadrature<data._qm.getNumberOfQuadraturePoints(); ++quadrature) {
179 XYZ quadrature_coord_3d;
180 XYZ tangent_quadrature_coord_3d;
181 for (
int j=0; j<dimension; ++j) {
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);
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++){
193 for (alphay = 0; alphay <= s; alphay++){
195 alphaf = factorial[alphax]*factorial[alphay]*factorial[alphaz];
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;
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;
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;
223 double dot_product = tangent_quadrature_coord_3d[0]*v0 + tangent_quadrature_coord_3d[1]*v1 + tangent_quadrature_coord_3d[2]*v2;
226 if (quadrature == 0) {
227 *(delta+i) = dot_product * data._qm.getWeight(quadrature);
229 *(delta+i) += dot_product * data._qm.getWeight(quadrature);
234 }
else if (dimension == 2) {
235 for (alphay = 0; alphay <= n; alphay++){
237 alphaf = factorial[alphax]*factorial[alphay];
244 v1 = std::pow(quadrature_coord_3d.
x/cutoff_p,alphax)
245 *std::pow(quadrature_coord_3d.
y/cutoff_p,alphay)/alphaf;
248 v0 = std::pow(quadrature_coord_3d.
x/cutoff_p,alphax)
249 *std::pow(quadrature_coord_3d.
y/cutoff_p,alphay)/alphaf;
254 double dot_product = tangent_quadrature_coord_3d[0]*v0 + tangent_quadrature_coord_3d[1]*v1;
257 if (quadrature == 0) {
258 *(delta+i) = dot_product * data._qm.getWeight(quadrature);
260 *(delta+i) += dot_product * data._qm.getWeight(quadrature);
277 "FaceNormalIntegralSample, EdgeTangentIntegralSample, FaceNormalAverageSample, \
278 and EdgeTangentAverageSample only support 2d or 3d with 2d manifold");
281 &&
"Only 1d quadrature may be specified for edge integrals");
284 &&
"Quadrature points not generated");
289 "Sampling functional does not support specific_order_only");
291 double cutoff_p = data._epsilons(target_index);
292 auto global_neighbor_index = data._pc.getNeighborIndex(target_index, neighbor_index);
303 int quadrature_point_loop = data._qm.getNumberOfQuadraturePoints();
306 double G_data[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;
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);
325 radius = std::sqrt(radius);
329 auto E = edge_coords_matrix;
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];
338 theta = std::atan(norm_a_cross_b / a_dot_b);
342 double entire_edge_length = 0.0;
343 for (
int quadrature = 0; quadrature<quadrature_point_loop; ++quadrature) {
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);
351 unscaled_transformed_qp[j] += E(j,0);
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];
364 transformed_qp_norm = std::sqrt(transformed_qp_norm);
366 for (
int j=0; j<data._global_dimensions; ++j) {
367 scaled_transformed_qp[j] = unscaled_transformed_qp[j] * radius / transformed_qp_norm;
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);
377 relative_coord[2] = 0;
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);
386 relative_coord[2] = 0;
393 for (
int j=0; j<data._global_dimensions; ++j) {
395 direction[j] = data._source_extra_data(global_neighbor_index, 2*data._global_dimensions + j);
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)};
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;
412 for (
int j=0; j<data._global_dimensions; ++j) {
414 direction[j] = data._source_extra_data(global_neighbor_index, 3*data._global_dimensions + j);
422 for (
int j=0; j<data._basis_multiplier; ++j) {
426 local_direction[j] = data._pc.convertGlobalToLocalCoordinate(direction,j,*V);
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++){
435 alphaf = factorial[alphax]*factorial[alphay];
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;
445 double dot_product = 0.0;
448 dot_product = local_direction[0]*v0 + local_direction[1]*v1;
450 dot_product = direction[0]*v0 + direction[1]*v1;
455 *(delta+i) = dot_product * data._qm.getWeight(quadrature) * G_determinant;
457 *(delta+i) += dot_product * data._qm.getWeight(quadrature) * G_determinant;
463 entire_edge_length += G_determinant * data._qm.getWeight(quadrature);
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;
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);
488 double triangle_coords[3*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;
501 for (
int j=0; j<data._global_dimensions; ++j) {
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);
506 radius = std::sqrt(radius);
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);
519 num_vertices = num_vertices / data._global_dimensions;
520 auto T = triangle_coords_matrix;
524 double entire_cell_area = 0.0;
525 for (
int v=0; v<num_vertices; ++v) {
527 int v2 = (v+1) % num_vertices;
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);
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) {
549 unscaled_transformed_qp[j] += T(j,k)*data._qm.getSite(quadrature, k-1);
552 unscaled_transformed_qp[j] += T(j,0);
556 double G_determinant = 1.0;
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];
566 transformed_qp_norm = std::sqrt(transformed_qp_norm);
568 for (
int j=0; j<data._global_dimensions; ++j) {
569 scaled_transformed_qp[j] = unscaled_transformed_qp[j] * radius / transformed_qp_norm;
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));
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);
610 relative_coord[2] = 0;
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);
619 relative_coord[2] = 0;
624 "CellAverageSample does not support specific_order_only");
625 for (
int n = 0; n <= poly_order; n++){
626 for (alphay = 0; alphay <= 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;
637 entire_cell_area += G_determinant * data._qm.getWeight(quadrature);
642 for (
int n = 0; n <= poly_order; n++){
643 for (alphay = 0; alphay <= n; alphay++){
644 *(delta+k) /= entire_cell_area;