29int main (
int argc,
char* args[]) {
32#ifdef COMPADRE_USE_MPI
33MPI_Init(&argc, &args);
37Kokkos::initialize(argc, args);
44 auto order = clp.
order;
58 Kokkos::Profiling::pushRegion(
"Setup Point Data");
63 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device(
"source coordinates",
64 1.25*N_pts_on_sphere, 3);
65 Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
70 int number_source_coords;
76 double a = 4*
PI*r*r/N_pts_on_sphere;
77 double d = std::sqrt(a);
78 int M_theta = std::round(
PI/d);
79 double d_theta =
PI/M_theta;
80 double d_phi = a/d_theta;
81 for (
int i=0; i<M_theta; ++i) {
82 double theta =
PI*(i + 0.5)/M_theta;
83 int M_phi = std::ceil(2*
PI*std::sin(theta)/d_phi);
84 for (
int j=0; j<M_phi; ++j) {
85 double phi = 2*
PI*j/M_phi;
86 source_coords(N_count, 0) = theta;
87 source_coords(N_count, 1) = phi;
92 for (
int i=0; i<N_count; ++i) {
93 double theta = source_coords(i,0);
94 double phi = source_coords(i,1);
95 source_coords(i,0) = r*std::sin(theta)*std::cos(phi);
96 source_coords(i,1) = r*std::sin(theta)*std::sin(phi);
97 source_coords(i,2) = r*cos(theta);
100 number_source_coords = N_count;
104 Kokkos::resize(source_coords, number_source_coords, 3);
105 Kokkos::resize(source_coords_device, number_source_coords, 3);
108 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device (
"target coordinates",
109 number_target_coords, 3);
110 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
113 bool enough_pts_found =
false;
117 double a = 4.0*
PI*r*r/((double)(1.0*number_target_coords));
118 double d = std::sqrt(a);
119 int M_theta = std::round(
PI/d);
120 double d_theta =
PI/((double)M_theta);
121 double d_phi = a/d_theta;
122 for (
int i=0; i<M_theta; ++i) {
123 double theta =
PI*(i + 0.5)/M_theta;
124 int M_phi = std::ceil(2*
PI*std::sin(theta)/d_phi);
125 for (
int j=0; j<M_phi; ++j) {
126 double phi = 2*
PI*j/M_phi;
127 target_coords(N_count, 0) = theta;
128 target_coords(N_count, 1) = phi;
130 if (N_count == number_target_coords) {
131 enough_pts_found =
true;
135 if (enough_pts_found)
break;
138 for (
int i=0; i<N_count; ++i) {
139 double theta = target_coords(i,0);
140 double phi = target_coords(i,1);
141 target_coords(i,0) = r*std::sin(theta)*std::cos(phi);
142 target_coords(i,1) = r*std::sin(theta)*std::sin(phi);
143 target_coords(i,2) = r*cos(theta);
154 Kokkos::Profiling::popRegion();
156 Kokkos::Profiling::pushRegion(
"Creating Data");
162 Kokkos::deep_copy(source_coords_device, source_coords);
165 Kokkos::deep_copy(target_coords_device, target_coords);
172 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
173 source_coords_device.extent(0));
175 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> ones_data_device(
"samples of ones",
176 source_coords_device.extent(0));
177 Kokkos::deep_copy(ones_data_device, 1.0);
180 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> sampling_vector_data_device(
"samples of vector true solution",
181 source_coords_device.extent(0), 3);
183 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
184 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
187 double xval = source_coords_device(i,0);
188 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
189 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
194 for (
int j=0; j<3; ++j) {
195 double gradient[3] = {0,0,0};
197 sampling_vector_data_device(i,j) = gradient[j];
204 Kokkos::Profiling::popRegion();
205 Kokkos::Profiling::pushRegion(
"Neighbor Search");
209 double epsilon_multiplier = 1.9;
215 Kokkos::View<int*> neighbor_lists_device(
"neighbor lists",
217 Kokkos::View<int*>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
220 Kokkos::View<int*> number_of_neighbors_list_device(
"number of neighbor lists",
221 number_target_coords);
222 Kokkos::View<int*>::HostMirror number_of_neighbors_list = Kokkos::create_mirror_view(number_of_neighbors_list_device);
225 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
226 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
228 size_t storage_size = point_cloud_search.generateCRNeighborListsFromKNNSearch(
true , target_coords, neighbor_lists,
229 number_of_neighbors_list, epsilon, min_neighbors, epsilon_multiplier);
232 Kokkos::resize(neighbor_lists_device, storage_size);
233 neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
236 point_cloud_search.generateCRNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
237 number_of_neighbors_list, epsilon, min_neighbors, epsilon_multiplier);
251 Kokkos::Profiling::popRegion();
262 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
263 Kokkos::deep_copy(number_of_neighbors_list_device, number_of_neighbors_list);
264 Kokkos::deep_copy(epsilon_device, epsilon);
268 GMLS my_GMLS_scalar(order, dimension,
269 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
286 my_GMLS_scalar.
setProblemData(neighbor_lists_device, number_of_neighbors_list_device, source_coords_device, target_coords_device, epsilon_device);
294 std::vector<TargetOperation> lro_scalar(4);
318 Kokkos::Profiling::pushRegion(
"Full Polynomial Basis GMLS Solution");
325 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
328 my_GMLS_vector.
setProblemData(neighbor_lists_device, number_of_neighbors_list_device, source_coords_device, target_coords_device, epsilon_device);
329 std::vector<TargetOperation> lro_vector(2);
338 Kokkos::Profiling::popRegion();
340 Kokkos::Profiling::pushRegion(
"Scalar Polynomial Basis Repeated to Form a Vector GMLS Solution");
365 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
368 my_GMLS_vector_of_scalar_clones.
setProblemData(neighbor_lists_device, number_of_neighbors_list_device, source_coords_device, target_coords_device, epsilon_device);
369 std::vector<TargetOperation> lro_vector_of_scalar_clones(2);
372 my_GMLS_vector_of_scalar_clones.
addTargets(lro_vector_of_scalar_clones);
375 my_GMLS_vector_of_scalar_clones.
setWeightingType(WeightingFunctionType::Power);
378 Kokkos::Profiling::popRegion();
383 double instantiation_time = timer.seconds();
384 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
386 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
401 Evaluator scalar_gmls_evaluator(&my_GMLS_scalar);
402 Evaluator vector_gmls_evaluator(&my_GMLS_vector);
403 Evaluator vector_gmls_evaluator_of_scalar_clones(&my_GMLS_vector_of_scalar_clones);
423 auto output_vector_of_scalar_clones =
427 auto output_divergence_of_scalar_clones =
461 Kokkos::Profiling::popRegion();
463 Kokkos::Profiling::pushRegion(
"Comparison");
467 double tangent_bundle_error = 0;
468 double tangent_bundle_norm = 0;
469 double values_error = 0;
470 double values_norm = 0;
471 double laplacian_error = 0;
472 double laplacian_norm = 0;
473 double gradient_ambient_error = 0;
474 double gradient_ambient_norm = 0;
477 double vector_ambient_error = 0;
478 double vector_ambient_norm = 0;
479 double divergence_ambient_error = 0;
480 double divergence_ambient_norm = 0;
481 double vector_of_scalar_clones_ambient_error = 0;
482 double vector_of_scalar_clones_ambient_norm = 0;
483 double divergence_of_scalar_clones_ambient_error = 0;
484 double divergence_of_scalar_clones_ambient_norm = 0;
487 for (
int i=0; i<number_target_coords; i++) {
490 double GMLS_value = output_value(i);
491 double GMLS_gc = output_gc(i);
494 double GMLS_Laplacian = output_laplacian(i);
497 double xval = target_coords(i,0);
498 double yval = (dimension>1) ? target_coords(i,1) : 0;
499 double zval = (dimension>2) ? target_coords(i,2) : 0;
500 double coord[3] = {xval, yval, zval};
503 for (
unsigned int j=0; j<static_cast<unsigned int>(dimension-1); ++j) {
506 double tangent_inner_prod = 0;
507 for (
int k=0; k<dimension; ++k) {
510 tangent_bundle_error += tangent_inner_prod * tangent_inner_prod;
512 double normal_inner_prod = 0;
513 for (
int k=0; k<dimension; ++k) {
514 normal_inner_prod += coord[k] * my_GMLS_scalar.
getTangentBundle(i, dimension-1, k);
517 double normal_error_to_sum = (normal_inner_prod > 0) ? normal_inner_prod - 1 : normal_inner_prod + 1;
518 tangent_bundle_error += normal_error_to_sum * normal_error_to_sum;
519 tangent_bundle_norm += 1;
524 double actual_Gradient_ambient[3] = {0,0,0};
528 values_error += (GMLS_value - actual_value)*(GMLS_value - actual_value);
529 values_norm += actual_value*actual_value;
531 laplacian_error += (GMLS_Laplacian - actual_Laplacian)*(GMLS_Laplacian - actual_Laplacian);
532 laplacian_norm += actual_Laplacian*actual_Laplacian;
534 double actual_gc = 1.0;
535 gc_error += (GMLS_gc - actual_gc)*(GMLS_gc - actual_gc);
536 gc_norm += actual_gc*actual_gc;
538 for (
int j=0; j<dimension; ++j) {
539 gradient_ambient_error += (output_gradient(i,j) - actual_Gradient_ambient[j])*(output_gradient(i,j) - actual_Gradient_ambient[j]);
540 gradient_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
543 for (
int j=0; j<dimension; ++j) {
544 vector_ambient_error += (output_vector(i,j) - actual_Gradient_ambient[j])*(output_vector(i,j) - actual_Gradient_ambient[j]);
545 vector_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
548 divergence_ambient_error += (output_divergence(i) - actual_Laplacian)*(output_divergence(i) - actual_Laplacian);
549 divergence_ambient_norm += actual_Laplacian*actual_Laplacian;
551 for (
int j=0; j<dimension; ++j) {
552 vector_of_scalar_clones_ambient_error += (output_vector_of_scalar_clones(i,j) - actual_Gradient_ambient[j])*(output_vector_of_scalar_clones(i,j) - actual_Gradient_ambient[j]);
553 vector_of_scalar_clones_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
556 divergence_of_scalar_clones_ambient_error += (output_divergence_of_scalar_clones(i) - actual_Laplacian)*(output_divergence_of_scalar_clones(i) - actual_Laplacian);
557 divergence_of_scalar_clones_ambient_norm += actual_Laplacian*actual_Laplacian;
561 tangent_bundle_error /= number_target_coords;
562 tangent_bundle_error = std::sqrt(tangent_bundle_error);
563 tangent_bundle_norm /= number_target_coords;
564 tangent_bundle_norm = std::sqrt(tangent_bundle_norm);
566 values_error /= number_target_coords;
567 values_error = std::sqrt(values_error);
568 values_norm /= number_target_coords;
569 values_norm = std::sqrt(values_norm);
571 laplacian_error /= number_target_coords;
572 laplacian_error = std::sqrt(laplacian_error);
573 laplacian_norm /= number_target_coords;
574 laplacian_norm = std::sqrt(laplacian_norm);
576 gradient_ambient_error /= number_target_coords;
577 gradient_ambient_error = std::sqrt(gradient_ambient_error);
578 gradient_ambient_norm /= number_target_coords;
579 gradient_ambient_norm = std::sqrt(gradient_ambient_norm);
581 gc_error /= number_target_coords;
582 gc_error = std::sqrt(gc_error);
583 gc_norm /= number_target_coords;
584 gc_norm = std::sqrt(gc_norm);
586 vector_ambient_error /= number_target_coords;
587 vector_ambient_error = std::sqrt(vector_ambient_error);
588 vector_ambient_norm /= number_target_coords;
589 vector_ambient_norm = std::sqrt(vector_ambient_norm);
591 divergence_ambient_error /= number_target_coords;
592 divergence_ambient_error = std::sqrt(divergence_ambient_error);
593 divergence_ambient_norm /= number_target_coords;
594 divergence_ambient_norm = std::sqrt(divergence_ambient_norm);
596 vector_of_scalar_clones_ambient_error /= number_target_coords;
597 vector_of_scalar_clones_ambient_error = std::sqrt(vector_of_scalar_clones_ambient_error);
598 vector_of_scalar_clones_ambient_norm /= number_target_coords;
599 vector_of_scalar_clones_ambient_norm = std::sqrt(vector_of_scalar_clones_ambient_norm);
601 divergence_of_scalar_clones_ambient_error /= number_target_coords;
602 divergence_of_scalar_clones_ambient_error = std::sqrt(divergence_of_scalar_clones_ambient_error);
603 divergence_of_scalar_clones_ambient_norm /= number_target_coords;
604 divergence_of_scalar_clones_ambient_norm = std::sqrt(divergence_of_scalar_clones_ambient_norm);
606 printf(
"Tangent Bundle Error: %g\n", tangent_bundle_error / tangent_bundle_norm);
607 printf(
"Point Value Error: %g\n", values_error / values_norm);
608 printf(
"Laplace-Beltrami Error: %g\n", laplacian_error / laplacian_norm);
609 printf(
"Gaussian Curvature Error: %g\n", gc_error / gc_norm);
610 printf(
"Surface Gradient (Ambient) Error: %g\n", gradient_ambient_error / gradient_ambient_norm);
611 printf(
"Surface Vector (VectorBasis) Error: %g\n", vector_ambient_error / vector_ambient_norm);
612 printf(
"Surface Divergence (VectorBasis) Error: %g\n", divergence_ambient_error / divergence_ambient_norm);
613 printf(
"Surface Vector (ScalarClones) Error: %g\n",
614 vector_of_scalar_clones_ambient_error / vector_of_scalar_clones_ambient_norm);
615 printf(
"Surface Divergence (ScalarClones) Error: %g\n",
616 divergence_of_scalar_clones_ambient_error / divergence_of_scalar_clones_ambient_norm);
620 Kokkos::Profiling::popRegion();
629#ifdef COMPADRE_USE_MPI