1 #ifndef _COMPADRE_BASIS_HPP_
2 #define _COMPADRE_BASIS_HPP_
24 template <
typename BasisData>
25 KOKKOS_INLINE_FUNCTION
26 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) {
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;
671 template <
typename BasisData>
672 KOKKOS_INLINE_FUNCTION
673 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) {
679 const int my_num_neighbors = data._pc._nla.getNumberOfNeighborsDevice(target_index);
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) {
689 component = -(neighbor_index+1);
697 if (neighbor_index > -1) {
698 for (
int i=0; i<dimension; ++i) {
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);
703 }
else if (evaluation_site_local_index > 0) {
704 for (
int i=0; i<dimension; ++i) {
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);
711 for (
int i=0; i<3; ++i) relative_coord[i] = 0;
714 double cutoff_p = data._epsilons(target_index);
715 const int start_index = specific_order_only ? poly_order : 0;
717 if ((polynomial_sampling_functional ==
PointSample ||
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);
729 DivergenceFreePolynomialBasis::evaluatePartialDerivative(teamMember, delta, thread_workspace, dimension, poly_order, component, partial_direction, cutoff_p, relative_coord.
x, relative_coord.
y, relative_coord.
z);
734 BernsteinPolynomialBasis::evaluatePartialDerivative(teamMember, delta, thread_workspace, dimension, poly_order, component, partial_direction, cutoff_p, relative_coord.
x, relative_coord.
y, relative_coord.
z);
759 template <
typename BasisData>
760 KOKKOS_INLINE_FUNCTION
761 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) {
767 const int my_num_neighbors = data._pc._nla.getNumberOfNeighborsDevice(target_index);
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) {
777 component = -(neighbor_index+1);
785 if (neighbor_index > -1) {
786 for (
int i=0; i<dimension; ++i) {
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);
791 }
else if (evaluation_site_local_index > 0) {
792 for (
int i=0; i<dimension; ++i) {
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);
799 for (
int i=0; i<3; ++i) relative_coord[i] = 0;
802 double cutoff_p = data._epsilons(target_index);
804 if ((polynomial_sampling_functional ==
PointSample ||
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);
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);
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);
840 template <
typename BasisData>
841 KOKKOS_INLINE_FUNCTION
842 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) {
846 const int target_index = data._initial_index_for_batch + teamMember.league_rank();
855 int my_num_neighbors = data._pc._nla.getNumberOfNeighborsDevice(target_index);
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) {
863 for (
int j = 0; j < thread_workspace.extent_int(0); ++j) {
864 thread_workspace(j) = 0;
866 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,data._pc._nla.getNumberOfNeighborsDevice(target_index)),
869 for (int d=0; d<data._sampling_multiplier; ++d) {
876 double alpha_weight = 1;
877 if (data._polynomial_sampling_functional==StaggeredEdgeIntegralSample
878 || data._polynomial_sampling_functional==StaggeredEdgeAnalyticGradientIntegralSample) {
884 r = data._pc.EuclideanVectorLength(data._pc.getRelativeCoord(target_index, i, dimension) * alpha_weight, dimension);
886 r = data._pc.EuclideanVectorLength(data._pc.getRelativeCoord(target_index, i, dimension, V) * alpha_weight, dimension);
890 w(i+my_num_neighbors*d) = GMLS::Wab(r, data._epsilons(target_index), data._weighting_type, data._weighting_p, data._weighting_n);
892 calcPij<BasisData>(data, teamMember, delta.data(), thread_workspace.data(), target_index, i + d*my_num_neighbors, 0 , dimension, polynomial_order, false , V, reconstruction_space, polynomial_sampling_functional);
895 for (int j = 0; j < storage_size; ++j) {
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.");
902 for (int j = 0; j < storage_size; ++j) {
904 P(i+my_num_neighbors*d, j) = delta[j];
906 compadre_kernel_assert_extreme_debug(delta[j]==delta[j] &&
"NaN in P matrix.");
911 teamMember.team_barrier();
937 template <
typename BasisData>
938 KOKKOS_INLINE_FUNCTION
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) {
951 for (
int j = 0; j < thread_workspace.extent_int(0); ++j) {
952 thread_workspace(j) = 0;
954 Kokkos::parallel_for(Kokkos::TeamThreadRange(teamMember,data._pc._nla.getNumberOfNeighborsDevice(target_index)),
962 r = data._pc.EuclideanVectorLength(data._pc.getRelativeCoord(target_index, i, dimension), dimension);
964 r = data._pc.EuclideanVectorLength(data._pc.getRelativeCoord(target_index, i, dimension, V), dimension);
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 , dimension, 1,
true );
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 , dimension, data._curvature_poly_order,
false , V);
976 for (
int j = 0; j < storage_size; ++j) {
977 P(i, j) = delta[j] * std::sqrt(w(i));
981 teamMember.team_barrier();
#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 evaluatePartialDerivative(const member_type &teamMember, double *delta, double *workspace, const int dimension, const int max_degree, const int component, const int partial_direction, 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 first partial derivatives of the Bernstein polynomial basis delta[j] = weight_of_origin...
KOKKOS_INLINE_FUNCTION void evaluateSecondPartialDerivative(const member_type &teamMember, double *delta, double *workspace, const int dimension, const int max_degree, const int component, const int partial_direction_1, const int partial_direction_2, 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 second partial derivatives of the Bernstein polynomial basis delta[j] = weight_of_origi...
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)