9#include <Compadre_Config.h>
17#ifdef COMPADRE_USE_MPI
21#include <Kokkos_Timer.hpp>
22#include <Kokkos_Core.hpp>
42 double local_vec[3] = {0,0};
43 for (
int j=0; j<3; ++j) {
44 local_vec[0] += T(0,j) * u[j];
45 local_vec[1] += T(1,j) * u[j];
48 for (
int j=0; j<3; ++j) u[j] = 0;
49 for (
int j=0; j<3; ++j) {
50 u[j] += P(0, j) * local_vec[0];
51 u[j] += P(1, j) * local_vec[1];
59int main (
int argc,
char* args[]) {
62#ifdef COMPADRE_USE_MPI
63MPI_Init(&argc, &args);
67Kokkos::initialize(argc, args);
74 auto order = clp.
order;
88 Kokkos::Profiling::pushRegion(
"Setup Point Data");
93 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device(
"source coordinates",
94 1.25*N_pts_on_sphere, 3);
95 Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
100 int number_source_coords;
106 double a = 4*
PI*r*r/N_pts_on_sphere;
107 double d = std::sqrt(a);
108 int M_theta = std::round(
PI/d);
109 double d_theta =
PI/M_theta;
110 double d_phi = a/d_theta;
111 for (
int i=0; i<M_theta; ++i) {
112 double theta =
PI*(i + 0.5)/M_theta;
113 int M_phi = std::round(2*
PI*std::sin(theta)/d_phi);
114 for (
int j=0; j<M_phi; ++j) {
115 double phi = 2*
PI*j/M_phi;
116 source_coords(N_count, 0) = theta;
117 source_coords(N_count, 1) = phi;
122 for (
int i=0; i<N_count; ++i) {
123 double theta = source_coords(i,0);
124 double phi = source_coords(i,1);
125 source_coords(i,0) = r*std::sin(theta)*std::cos(phi);
126 source_coords(i,1) = r*std::sin(theta)*std::sin(phi);
127 source_coords(i,2) = r*cos(theta);
130 number_source_coords = N_count;
134 Kokkos::resize(source_coords, number_source_coords, 3);
135 Kokkos::resize(source_coords_device, number_source_coords, 3);
138 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device (
"target coordinates",
139 number_target_coords, 3);
140 Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
143 bool enough_pts_found =
false;
147 double a = 4.0*
PI*r*r/((double)(5.0*number_target_coords));
148 double d = std::sqrt(a);
149 int M_theta = std::round(
PI/d);
150 double d_theta =
PI/((double)M_theta);
151 double d_phi = a/d_theta;
152 for (
int i=0; i<M_theta; ++i) {
153 double theta =
PI*(i + 0.5)/M_theta;
154 int M_phi = std::round(2*
PI*std::sin(theta)/d_phi);
155 for (
int j=0; j<M_phi; ++j) {
156 double phi = 2*
PI*j/M_phi;
157 target_coords(N_count, 0) = theta;
158 target_coords(N_count, 1) = phi;
160 if (N_count == number_target_coords) {
161 enough_pts_found =
true;
165 if (enough_pts_found)
break;
168 for (
int i=0; i<N_count; ++i) {
169 double theta = target_coords(i,0);
170 double phi = target_coords(i,1);
171 target_coords(i,0) = r*std::sin(theta)*std::cos(phi);
172 target_coords(i,1) = r*std::sin(theta)*std::sin(phi);
173 target_coords(i,2) = r*cos(theta);
181 Kokkos::Profiling::popRegion();
183 Kokkos::Profiling::pushRegion(
"Creating Data");
189 Kokkos::deep_copy(source_coords_device, source_coords);
192 Kokkos::deep_copy(target_coords_device, target_coords);
199 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device(
"samples of true solution",
200 source_coords_device.extent(0));
201 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> ones_data_device(
"samples of ones",
202 source_coords_device.extent(0));
203 Kokkos::deep_copy(ones_data_device, 1.0);
206 Kokkos::View<double**, Kokkos::DefaultExecutionSpace> sampling_vector_data_device(
"samples of vector true solution",
207 source_coords_device.extent(0), 3);
209 Kokkos::parallel_for(
"Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
210 (0,source_coords.extent(0)), KOKKOS_LAMBDA(
const int i) {
213 double xval = source_coords_device(i,0);
214 double yval = (dimension>1) ? source_coords_device(i,1) : 0;
215 double zval = (dimension>2) ? source_coords_device(i,2) : 0;
220 for (
int j=0; j<3; ++j) {
221 double gradient[3] = {0,0,0};
223 sampling_vector_data_device(i,j) = gradient[j];
230 Kokkos::Profiling::popRegion();
231 Kokkos::Profiling::pushRegion(
"Neighbor Search");
242 double epsilon_multiplier = 1.7;
243 int estimated_upper_bound_number_neighbors =
244 point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
246 Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device(
"neighbor lists",
247 number_target_coords, estimated_upper_bound_number_neighbors);
248 Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
251 Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device(
"h supports", number_target_coords);
252 Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
257 point_cloud_search.generate2DNeighborListsFromKNNSearch(
false , target_coords, neighbor_lists,
258 epsilon, min_neighbors, epsilon_multiplier);
262 Kokkos::Profiling::popRegion();
273 Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
274 Kokkos::deep_copy(epsilon_device, epsilon);
278 GMLS my_GMLS_scalar(order, dimension,
279 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
296 my_GMLS_scalar.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
308 std::vector<TargetOperation> lro_scalar(3);
331 Kokkos::Profiling::pushRegion(
"Full Polynomial Basis GMLS Solution");
338 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
341 my_GMLS_vector.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
343 std::vector<TargetOperation> lro_vector(2);
352 Kokkos::Profiling::popRegion();
354 Kokkos::Profiling::pushRegion(
"Scalar Polynomial Basis Repeated to Form a Vector GMLS Solution");
379 solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
382 my_GMLS_vector_of_scalar_clones.
setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
384 std::vector<TargetOperation> lro_vector_of_scalar_clones(2);
387 my_GMLS_vector_of_scalar_clones.
addTargets(lro_vector_of_scalar_clones);
390 my_GMLS_vector_of_scalar_clones.
setWeightingType(WeightingFunctionType::Power);
393 Kokkos::Profiling::popRegion();
398 double instantiation_time = timer.seconds();
399 std::cout <<
"Took " << instantiation_time <<
"s to complete alphas generation." << std::endl;
401 Kokkos::Profiling::pushRegion(
"Apply Alphas to Data");
416 Evaluator scalar_gmls_evaluator(&my_GMLS_scalar);
417 Evaluator vector_gmls_evaluator(&my_GMLS_vector);
418 Evaluator vector_gmls_evaluator_of_scalar_clones(&my_GMLS_vector_of_scalar_clones);
448 auto output_vector_of_scalar_clones =
488 Kokkos::Profiling::popRegion();
490 Kokkos::Profiling::pushRegion(
"Comparison");
494 double tangent_bundle_error = 0;
495 double tangent_bundle_norm = 0;
496 double values_error = 0;
497 double values_norm = 0;
500 double curl_ambient_error = 0;
501 double curl_ambient_norm = 0;
506 double vector_ambient_error = 0;
507 double vector_ambient_norm = 0;
510 double vector_of_scalar_clones_ambient_error = 0;
511 double vector_of_scalar_clones_ambient_norm = 0;
517 auto prestencil_weights = Kokkos::create_mirror_view(d_prestencil_weights);
518 Kokkos::deep_copy(prestencil_weights, d_prestencil_weights);
522 auto tangent_directions = Kokkos::create_mirror_view(d_tangent_directions);
523 Kokkos::deep_copy(tangent_directions, d_tangent_directions);
526 for (
int i=0; i<number_target_coords; i++) {
534 double GMLS_value = output_value(i);
537 double GMLS_gc = output_gaussian_curvature(i);
547 double xval = source_coords(neighbor_lists(i,1),0);
548 double yval = (dimension>1) ? source_coords(neighbor_lists(i,1),1) : 0;
549 double zval = (dimension>2) ? source_coords(neighbor_lists(i,1),2) : 0;
550 double coord[3] = {xval, yval, zval};
553 for (
int j=0; j<dimension-1; ++j) {
554 double tangent_inner_prod = 0;
555 for (
int k=0; k<std::min(dimension,3); ++k) {
556 tangent_inner_prod += coord[k] * prestencil_weights(0, i, 0 , j, k);
558 tangent_bundle_error += tangent_inner_prod * tangent_inner_prod;
560 double normal_inner_prod = 0;
561 for (
int k=0; k<dimension; ++k) {
562 normal_inner_prod += coord[k] * my_GMLS_scalar.
getTangentBundle(i, dimension-1, k);
565 double normal_error_to_sum = (normal_inner_prod > 0) ? normal_inner_prod - 1 : normal_inner_prod + 1;
566 tangent_bundle_error += normal_error_to_sum * normal_error_to_sum;
567 tangent_bundle_norm += 1;
571 double actual_gc = 1.0;
573 double actual_Gradient_ambient[3] = {0,0,0};
576 double actual_Curl_ambient[3] = {0,0,0};
579 values_error += (GMLS_value - actual_value)*(GMLS_value - actual_value);
580 values_norm += actual_value*actual_value;
582 gc_error += (GMLS_gc - actual_gc)*(GMLS_gc - actual_gc);
585 for (
int j=0; j<dimension; ++j) u[j] = output_curl(i,j);
587 for (
int j=0; j<dimension; ++j) output_curl(i,j) = u[j];
589 for (
int j=0; j<dimension; ++j) {
590 curl_ambient_error += (output_curl(i,j) - actual_Curl_ambient[j])*(output_curl(i,j) - actual_Curl_ambient[j]);
591 curl_ambient_norm += actual_Curl_ambient[j]*actual_Curl_ambient[j];
602 for (
int j=0; j<dimension; ++j) u[j] = output_vector(i,j);
604 for (
int j=0; j<dimension; ++j) output_vector(i,j) = u[j];
606 for (
int j=0; j<dimension; ++j) {
607 vector_ambient_error += (output_vector(i,j) - actual_Gradient_ambient[j])*(output_vector(i,j) - actual_Gradient_ambient[j]);
608 vector_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
615 for (
int j=0; j<dimension; ++j) u[j] = output_vector_of_scalar_clones(i,j);
617 for (
int j=0; j<dimension; ++j) output_vector_of_scalar_clones(i,j) = u[j];
632 for (
int j=0; j<dimension; ++j) {
633 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]);
634 vector_of_scalar_clones_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
642 tangent_bundle_error /= number_target_coords;
643 tangent_bundle_error = std::sqrt(tangent_bundle_error);
644 tangent_bundle_norm /= number_target_coords;
645 tangent_bundle_norm = std::sqrt(tangent_bundle_norm);
647 values_error /= number_target_coords;
648 values_error = std::sqrt(values_error);
649 values_norm /= number_target_coords;
650 values_norm = std::sqrt(values_norm);
652 gc_error /= number_target_coords;
653 gc_error = std::sqrt(gc_error);
654 gc_norm /= number_target_coords;
655 gc_norm = std::sqrt(gc_norm);
657 curl_ambient_error /= number_target_coords;
658 curl_ambient_error = std::sqrt(curl_ambient_error);
659 curl_ambient_norm /= number_target_coords;
660 curl_ambient_norm = std::sqrt(curl_ambient_norm);
672 vector_ambient_error /= number_target_coords;
673 vector_ambient_error = std::sqrt(vector_ambient_error);
674 vector_ambient_norm /= number_target_coords;
675 vector_ambient_norm = std::sqrt(vector_ambient_norm);
682 vector_of_scalar_clones_ambient_error /= number_target_coords;
683 vector_of_scalar_clones_ambient_error = std::sqrt(vector_of_scalar_clones_ambient_error);
684 vector_of_scalar_clones_ambient_norm /= number_target_coords;
685 vector_of_scalar_clones_ambient_norm = std::sqrt(vector_of_scalar_clones_ambient_norm);
692 printf(
"Tangent Bundle Error: %g\n", tangent_bundle_error / tangent_bundle_norm);
693 printf(
"Point Value Error: %g\n", values_error / values_norm);
694 printf(
"Gaussian Curvature Error: %g\n", gc_error / gc_norm);
695 printf(
"Surface Curl (Ambient) Error: %g\n", curl_ambient_error / curl_ambient_norm);
698 printf(
"Surface Vector (VectorBasis) Error: %g\n", vector_ambient_error / vector_ambient_norm);
700 printf(
"Surface Vector (ScalarClones) Error: %g\n",
701 vector_of_scalar_clones_ambient_error / vector_of_scalar_clones_ambient_norm);
707 Kokkos::Profiling::popRegion();
716#ifdef COMPADRE_USE_MPI
#define TO_GLOBAL(variable)
Kokkos::View< double **, layout_right, host_execution_space, Kokkos::MemoryTraits< Kokkos::Unmanaged > > host_scratch_matrix_right_type
KOKKOS_INLINE_FUNCTION double sphere_harmonic54(double x, double y, double z)
KOKKOS_INLINE_FUNCTION void gradient_sphereHarmonic54_ambient(double *gradient, double x, double y, double z)
KOKKOS_INLINE_FUNCTION void curl_sphere_harmonic54(double *curl, double x, double y, double z)
void AmbientLocalAmbient(XYZ &u, double *T_data, double *P_data)
[Ambient to Local Back To Ambient Helper Function]
int main(int argc, char *args[])
[Ambient to Local Back To Ambient Helper Function]
Lightweight Evaluator Helper This class is a lightweight wrapper for extracting and applying all rele...
Kokkos::View< output_data_type, output_array_layout, output_memory_space > applyAlphasToDataAllComponentsAllTargetSites(view_type_input_data sampling_data, TargetOperation lro, const SamplingFunctional sro_in=PointSample, bool scalar_as_vector_if_needed=true, const int evaluation_site_local_index=0) const
Transformation of data under GMLS (allocates memory for output)
Generalized Moving Least Squares (GMLS)
void addTargets(TargetOperation lro)
Adds a target to the vector of target functional to be applied to the reconstruction.
void setReferenceOutwardNormalDirection(view_type outward_normal_directions, bool use_to_orient_surface=true)
(OPTIONAL) Sets outward normal direction.
void setCurvatureWeightingParameter(int wp, int index=0)
Parameter for weighting kernel for curvature index = 0 sets p paramater for weighting kernel index = ...
void setAdditionalEvaluationSitesData(view_type_1 additional_evaluation_indices, view_type_2 additional_evaluation_coordinates)
(OPTIONAL) Sets additional evaluation sites for each target site
void setWeightingParameter(int wp, int index=0)
Parameter for weighting kernel for GMLS problem index = 0 sets p paramater for weighting kernel index...
void generateAlphas(const int number_of_batches=1, const bool keep_coefficients=false, const bool clear_cache=true)
Meant to calculate target operations and apply the evaluations to the previously constructed polynomi...
void setProblemData(view_type_1 neighbor_lists, view_type_2 source_coordinates, view_type_3 target_coordinates, view_type_4 epsilons)
Sets basic problem data (neighbor lists, source coordinates, and target coordinates)
void setWeightingType(const std::string &wt)
Type for weighting kernel for GMLS problem.
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....
decltype(_T) * getTangentDirections()
Get a view (device) of all tangent direction bundles.
void setCurvatureWeightingType(const std::string &wt)
Type for weighting kernel for curvature.
double getTangentBundle(const int target_index, const int direction, const int component) const
Get component of tangent or normal directions for manifold problems.
decltype(_prestencil_weights) getPrestencilWeights() const
Get a view (device) of all rank 2 preprocessing tensors This is a rank 5 tensor that is able to provi...
PointCloudSearch< view_type > CreatePointCloudSearch(view_type src_view, const local_index_type dimensions=-1, const local_index_type max_leaf=-1)
CreatePointCloudSearch allows for the construction of an object of type PointCloudSearch with templat...
constexpr SamplingFunctional VaryingManifoldVectorPointSample
For integrating polynomial dotted with normal over an edge.
constexpr SamplingFunctional PointSample
Available sampling functionals.
@ CurlOfVectorPointEvaluation
Point evaluation of the curl of a vector (results in a vector)
@ GaussianCurvaturePointEvaluation
Point evaluation of Gaussian curvature.
@ VectorPointEvaluation
Point evaluation of a vector (reconstructs entire vector at once, requiring a ReconstructionSpace hav...
@ ScalarPointEvaluation
Point evaluation of a scalar.
std::string constraint_name