Compadre  1.3.3
GMLS_Manifold.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <string>
3 #include <vector>
4 #include <map>
5 #include <stdlib.h>
6 #include <cstdio>
7 #include <random>
8 
9 #include <Compadre_Config.h>
10 #include <Compadre_GMLS.hpp>
11 #include <Compadre_Evaluator.hpp>
13 
14 #include "GMLS_Manifold.hpp"
15 #include "CommandLineProcessor.hpp"
16 
17 #ifdef COMPADRE_USE_MPI
18 #include <mpi.h>
19 #endif
20 
21 #include <Kokkos_Timer.hpp>
22 #include <Kokkos_Core.hpp>
23 
24 using namespace Compadre;
25 
26 //! [Parse Command Line Arguments]
27 
28 // called from command line
29 int main (int argc, char* args[]) {
30 
31 // initializes MPI (if available) with command line arguments given
32 #ifdef COMPADRE_USE_MPI
33 MPI_Init(&argc, &args);
34 #endif
35 
36 // initializes Kokkos with command line arguments given
37 Kokkos::initialize(argc, args);
38 
39 // code block to reduce scope for all Kokkos View allocations
40 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
41 {
42 
43  CommandLineProcessor clp(argc, args);
44  auto order = clp.order;
45  auto dimension = clp.dimension;
46  auto number_target_coords = clp.number_target_coords;
47  auto constraint_name = clp.constraint_name;
48  auto solver_name = clp.solver_name;
49  auto problem_name = clp.problem_name;
50  int N_pts_on_sphere = (clp.number_source_coords>=0) ? clp.number_source_coords : 1000;
51 
52  // minimum neighbors for unisolvency is the same as the size of the polynomial basis
53  // dimension has one subtracted because it is a D-1 manifold represented in D dimensions
54  const int min_neighbors = Compadre::GMLS::getNP(order, dimension-1);
55 
56  //! [Parse Command Line Arguments]
57  Kokkos::Timer timer;
58  Kokkos::Profiling::pushRegion("Setup Point Data");
59  //! [Setting Up The Point Cloud]
60 
61 
62  // coordinates of source sites, bigger than needed then resized later
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);
66 
67  double r = 1.0;
68 
69  // set number of source coordinates from what was calculated
70  int number_source_coords;
71 
72  {
73  // fill source coordinates from surface of a sphere with quasiuniform points
74  // algorithm described at https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf
75  int N_count = 0;
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::round(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;
88  N_count++;
89  }
90  }
91 
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);
98  //printf("%f %f %f\n", source_coords(i,0), source_coords(i,1), source_coords(i,2));
99  }
100  number_source_coords = N_count;
101  }
102 
103  // resize source_coords to the size actually needed
104  Kokkos::resize(source_coords, number_source_coords, 3);
105  Kokkos::resize(source_coords_device, number_source_coords, 3);
106 
107  // coordinates of target sites
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);
111 
112  {
113  bool enough_pts_found = false;
114  // fill target coordinates from surface of a sphere with quasiuniform points
115  // stop when enough points are found
116  int N_count = 0;
117  double a = 4.0*PI*r*r/((double)(5.0*number_target_coords)); // 5.0 is in case number_target_coords is set to something very small (like 1)
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::round(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;
129  N_count++;
130  if (N_count == number_target_coords) {
131  enough_pts_found = true;
132  break;
133  }
134  }
135  if (enough_pts_found) break;
136  }
137 
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);
144  //printf("%f %f %f\n", target_coords(i,0), target_coords(i,1), target_coords(i,2));
145  }
146  }
147 
148 
149  //! [Setting Up The Point Cloud]
150 
151  Kokkos::Profiling::popRegion();
152  Kokkos::fence();
153  Kokkos::Profiling::pushRegion("Creating Data");
154 
155  //! [Creating The Data]
156 
157 
158  // source coordinates need copied to device before using to construct sampling data
159  Kokkos::deep_copy(source_coords_device, source_coords);
160 
161  // target coordinates copied next, because it is a convenient time to send them to device
162  Kokkos::deep_copy(target_coords_device, target_coords);
163 
164  // ensure that source coordinates are sent to device before evaluating sampling data based on them
165  Kokkos::fence();
166 
167 
168  // need Kokkos View storing true solution (for samples)
169  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device("samples of true solution",
170  source_coords_device.extent(0));
171 
172  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> ones_data_device("samples of ones",
173  source_coords_device.extent(0));
174  Kokkos::deep_copy(ones_data_device, 1.0);
175 
176  // need Kokkos View storing true vector solution (for samples)
177  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> sampling_vector_data_device("samples of vector true solution",
178  source_coords_device.extent(0), 3);
179 
180  Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
181  (0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
182 
183  // coordinates of source site i
184  double xval = source_coords_device(i,0);
185  double yval = (dimension>1) ? source_coords_device(i,1) : 0;
186  double zval = (dimension>2) ? source_coords_device(i,2) : 0;
187 
188  // data for targets with scalar input
189  sampling_data_device(i) = sphere_harmonic54(xval, yval, zval);
190 
191  for (int j=0; j<3; ++j) {
192  double gradient[3] = {0,0,0};
193  gradient_sphereHarmonic54_ambient(gradient, xval, yval, zval);
194  sampling_vector_data_device(i,j) = gradient[j];
195  }
196  });
197 
198 
199  //! [Creating The Data]
200 
201  Kokkos::Profiling::popRegion();
202  Kokkos::Profiling::pushRegion("Neighbor Search");
203 
204  //! [Performing Neighbor Search]
205 
206 
207  // Point cloud construction for neighbor search
208  // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
209  auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
210 
211  // each row is a neighbor list for a target site, with the first column of each row containing
212  // the number of neighbors for that rows corresponding target site
213  double epsilon_multiplier = 1.9;
214  int estimated_upper_bound_number_neighbors =
215  point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
216 
217  Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device("neighbor lists",
218  number_target_coords, estimated_upper_bound_number_neighbors); // first column is # of neighbors
219  Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
220 
221  // each target site has a window size
222  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device("h supports", number_target_coords);
223  Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
224 
225  // query the point cloud to generate the neighbor lists using a kdtree to produce the n nearest neighbor
226  // to each target site, adding (epsilon_multiplier-1)*100% to whatever the distance away the further neighbor used is from
227  // each target to the view for epsilon
228  point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/, target_coords, neighbor_lists,
229  epsilon, min_neighbors, epsilon_multiplier);
230 
231  //! [Performing Neighbor Search]
232 
233  Kokkos::Profiling::popRegion();
234  Kokkos::fence(); // let call to build neighbor lists complete before copying back to device
235  timer.reset();
236 
237  //! [Setting Up The GMLS Object]
238 
239 
240  // Copy data back to device (they were filled on the host)
241  // We could have filled Kokkos Views with memory space on the host
242  // and used these instead, and then the copying of data to the device
243  // would be performed in the GMLS class
244  Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
245  Kokkos::deep_copy(epsilon_device, epsilon);
246 
247  // initialize an instance of the GMLS class for problems with a scalar basis and
248  // traditional point sampling as the sampling functional
249  GMLS my_GMLS_scalar(order, dimension,
250  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
251  order /*manifold order*/);
252 
253  // pass in neighbor lists, source coordinates, target coordinates, and window sizes
254  //
255  // neighbor lists have the format:
256  // dimensions: (# number of target sites) X (# maximum number of neighbors for any given target + 1)
257  // the first column contains the number of neighbors for that rows corresponding target index
258  //
259  // source coordinates have the format:
260  // dimensions: (# number of source sites) X (dimension)
261  // entries in the neighbor lists (integers) correspond to rows of this 2D array
262  //
263  // target coordinates have the format:
264  // dimensions: (# number of target sites) X (dimension)
265  // # of target sites is same as # of rows of neighbor lists
266  //
267  my_GMLS_scalar.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
268 
269  // set a reference outward normal direction, causing a surface orientation after
270  // the GMLS instance computes an approximate tangent bundle
271  // on a sphere, the ambient coordinates are the outward normal direction
272  my_GMLS_scalar.setReferenceOutwardNormalDirection(target_coords_device, true /* use to orient surface */);
273 
274  // create a vector of target operations
275  std::vector<TargetOperation> lro_scalar(4);
276  lro_scalar[0] = ScalarPointEvaluation;
277  lro_scalar[1] = LaplacianOfScalarPointEvaluation;
278  lro_scalar[2] = GradientOfScalarPointEvaluation;
279  lro_scalar[3] = GaussianCurvaturePointEvaluation;
280 
281  // and then pass them to the GMLS class
282  my_GMLS_scalar.addTargets(lro_scalar);
283 
284  // sets the weighting kernel function from WeightingFunctionType for curvature
285  my_GMLS_scalar.setCurvatureWeightingType(WeightingFunctionType::Power);
286 
287  // power to use in the weighting kernel function for curvature coefficients
288  my_GMLS_scalar.setCurvatureWeightingPower(2);
289 
290  // sets the weighting kernel function from WeightingFunctionType
291  my_GMLS_scalar.setWeightingType(WeightingFunctionType::Power);
292 
293  // power to use in that weighting kernel function
294  my_GMLS_scalar.setWeightingPower(2);
295 
296  // generate the alphas that to be combined with data for each target operation requested in lro
297  my_GMLS_scalar.generateAlphas();
298 
299  Kokkos::Profiling::pushRegion("Full Polynomial Basis GMLS Solution");
300  // initialize another instance of the GMLS class for problems with a vector basis on a manifold and point
301  // evaluation of that vector as the sampling functional
302  // VectorTaylorPolynomial indicates that the basis will be a polynomial with as many components as the
303  // dimension of the manifold. This differs from another possibility, which follows this class.
305  order, dimension,
306  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
307  order /*manifold order*/);
308 
309  my_GMLS_vector.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
310  std::vector<TargetOperation> lro_vector(2);
311  lro_vector[0] = VectorPointEvaluation;
312  lro_vector[1] = DivergenceOfVectorPointEvaluation;
313  my_GMLS_vector.addTargets(lro_vector);
314  my_GMLS_vector.setCurvatureWeightingType(WeightingFunctionType::Power);
315  my_GMLS_vector.setCurvatureWeightingPower(2);
316  my_GMLS_vector.setWeightingType(WeightingFunctionType::Power);
317  my_GMLS_vector.setWeightingPower(2);
318  my_GMLS_vector.generateAlphas();
319  Kokkos::Profiling::popRegion();
320 
321  Kokkos::Profiling::pushRegion("Scalar Polynomial Basis Repeated to Form a Vector GMLS Solution");
322  // initialize another instance of the GMLS class for problems with a vector basis on a manifold and point
323  // evaluation of that vector as the sampling functional
324  // VectorOfScalarClonesTaylorPolynomial indicates a scalar polynomial will be solved for, since
325  // each component of the reconstructed vector are independent. This results in solving a smaller system
326  // for each GMLS problem, and is the suggested way to do vector reconstructions when sampling functionals
327  // acting on the basis would not create non-zero offdiagonal blocks.
328  //
329  // As an example, consider a 2D manifold in 3D ambient space. The GMLS problem is posed in the local chart,
330  // meaning that the problem being solved looks like
331  //
332  // [P_0 0 | where P_0 has dimension #number of neighbors for a target X #dimension of a scalar basis
333  // | 0 P_1]
334  //
335  // P_1 is similarly shaped, and for sampling functional that is a point evaluation, P_0 and P_1 are
336  // identical and their degrees of freedom in this system are disjoint, allowing us to solve for the
337  // degrees of freedom for either block independently. Additionally, the will produce the exact
338  // same polynomial coefficients for the point sampling functional, therefore it makes sense to use
339  // VectorOfScalarClonesTaylorPolynomial.
340  //
341  // In the print-out for this program, we include the timings and errors on this and VectorTaylorPolynomial
342  // in order to demonstrate that they produce exactly the same answer, but that one is much more efficient.
343  //
345  order, dimension,
346  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
347  order /*manifold order*/);
348 
349  my_GMLS_vector_of_scalar_clones.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
350  std::vector<TargetOperation> lro_vector_of_scalar_clones(2);
351  lro_vector_of_scalar_clones[0] = VectorPointEvaluation;
352  lro_vector_of_scalar_clones[1] = DivergenceOfVectorPointEvaluation;
353  my_GMLS_vector_of_scalar_clones.addTargets(lro_vector_of_scalar_clones);
354  my_GMLS_vector_of_scalar_clones.setCurvatureWeightingType(WeightingFunctionType::Power);
355  my_GMLS_vector_of_scalar_clones.setCurvatureWeightingPower(2);
356  my_GMLS_vector_of_scalar_clones.setWeightingType(WeightingFunctionType::Power);
357  my_GMLS_vector_of_scalar_clones.setWeightingPower(2);
358  my_GMLS_vector_of_scalar_clones.generateAlphas();
359  Kokkos::Profiling::popRegion();
360 
361 
362  //! [Setting Up The GMLS Object]
363 
364  double instantiation_time = timer.seconds();
365  std::cout << "Took " << instantiation_time << "s to complete alphas generation." << std::endl;
366  Kokkos::fence(); // let generateAlphas finish up before using alphas
367  Kokkos::Profiling::pushRegion("Apply Alphas to Data");
368 
369  //! [Apply GMLS Alphas To Data]
370 
371 
372  // it is important to note that if you expect to use the data as a 1D view, then you should use double*
373  // however, if you know that the target operation will result in a 2D view (vector or matrix output),
374  // then you should template with double** as this is something that can not be infered from the input data
375  // or the target operator at compile time. Additionally, a template argument is required indicating either
376  // Kokkos::HostSpace or Kokkos::DefaultExecutionSpace::memory_space()
377 
378  // The Evaluator class takes care of handling input data views as well as the output data views.
379  // It uses information from the GMLS class to determine how many components are in the input
380  // as well as output for any choice of target functionals and then performs the contactions
381  // on the data using the alpha coefficients generated by the GMLS class, all on the device.
382  Evaluator scalar_gmls_evaluator(&my_GMLS_scalar);
383  Evaluator vector_gmls_evaluator(&my_GMLS_vector);
384  Evaluator vector_gmls_evaluator_of_scalar_clones(&my_GMLS_vector_of_scalar_clones);
385 
386  auto output_value = scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
387  (sampling_data_device, ScalarPointEvaluation);
388 
389  auto output_laplacian = scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
390  (sampling_data_device, LaplacianOfScalarPointEvaluation);
391 
392  auto output_gradient = scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
393  (sampling_data_device, GradientOfScalarPointEvaluation);
394 
395  auto output_gc = scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
396  (ones_data_device, GaussianCurvaturePointEvaluation);
397 
398  auto output_vector = vector_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
399  (sampling_vector_data_device, VectorPointEvaluation, ManifoldVectorPointSample);
400 
401  auto output_divergence = vector_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
402  (sampling_vector_data_device, DivergenceOfVectorPointEvaluation, ManifoldVectorPointSample);
403 
404  auto output_vector_of_scalar_clones =
405  vector_gmls_evaluator_of_scalar_clones.applyAlphasToDataAllComponentsAllTargetSites<double**,
406  Kokkos::HostSpace>(sampling_vector_data_device, VectorPointEvaluation, ManifoldVectorPointSample);
407 
408  auto output_divergence_of_scalar_clones =
409  vector_gmls_evaluator_of_scalar_clones.applyAlphasToDataAllComponentsAllTargetSites<double*,
410  Kokkos::HostSpace>(sampling_vector_data_device, DivergenceOfVectorPointEvaluation, ManifoldVectorPointSample);
411 
412 
413  // Kokkos::fence(); // let application of alphas to data finish before using results
414  //
415  //// move gradient data to device so that it can be transformed into velocity
416  //auto output_gradient_device_mirror = Kokkos::create_mirror(Kokkos::DefaultExecutionSpace::memory_space(), output_gradient);
417  //Kokkos::deep_copy(output_gradient_device_mirror, output_gradient);
418  //Kokkos::parallel_for("Create Velocity From Surface Gradient", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
419  // (0,target_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
420  //
421  // // coordinates of target site i
422  // double xval = target_coords_device(i,0);
423  // double yval = (dimension>1) ? target_coords_device(i,1) : 0;
424  // double zval = (dimension>2) ? target_coords_device(i,2) : 0;
425 
426  // double gradx = output_gradient_device_mirror(i,0);
427  // double grady = output_gradient_device_mirror(i,1);
428  // double gradz = output_gradient_device_mirror(i,2);
429  //
430  // // overwrites gradient with velocity
431  // output_gradient_device_mirror(i,0) = (grady*zval - yval*gradz);
432  // output_gradient_device_mirror(i,1) = (-(gradx*zval - xval*gradz));
433  // output_gradient_device_mirror(i,2) = (gradx*yval - xval*grady);
434  //
435  //});
436  //Kokkos::deep_copy(output_gradient, output_gradient_device_mirror);
437 
438 
439  //! [Apply GMLS Alphas To Data]
440 
441  Kokkos::fence(); // let application of alphas to data finish before using results
442  Kokkos::Profiling::popRegion();
443  // times the Comparison in Kokkos
444  Kokkos::Profiling::pushRegion("Comparison");
445 
446  //! [Check That Solutions Are Correct]
447 
448  double tangent_bundle_error = 0;
449  double tangent_bundle_norm = 0;
450  double values_error = 0;
451  double values_norm = 0;
452  double laplacian_error = 0;
453  double laplacian_norm = 0;
454  double gradient_ambient_error = 0;
455  double gradient_ambient_norm = 0;
456  double gc_error = 0;
457  double gc_norm = 0;
458  double vector_ambient_error = 0;
459  double vector_ambient_norm = 0;
460  double divergence_ambient_error = 0;
461  double divergence_ambient_norm = 0;
462  double vector_of_scalar_clones_ambient_error = 0;
463  double vector_of_scalar_clones_ambient_norm = 0;
464  double divergence_of_scalar_clones_ambient_error = 0;
465  double divergence_of_scalar_clones_ambient_norm = 0;
466 
467  // loop through the target sites
468  for (int i=0; i<number_target_coords; i++) {
469 
470  // load value from output
471  double GMLS_value = output_value(i);
472  double GMLS_gc = output_gc(i);
473 
474  // load laplacian from output
475  double GMLS_Laplacian = output_laplacian(i);
476 
477  // target site i's coordinate
478  double xval = target_coords(i,0);
479  double yval = (dimension>1) ? target_coords(i,1) : 0;
480  double zval = (dimension>2) ? target_coords(i,2) : 0;
481  double coord[3] = {xval, yval, zval};
482 
483  // get tangent vector and see if orthgonal to coordinate (it should be on a sphere)
484  for (int j=0; j<dimension-1; ++j) {
485  double tangent_inner_prod = 0;
486  for (int k=0; k<dimension; ++k) {
487  tangent_inner_prod += coord[k] * my_GMLS_scalar.getTangentBundle(i, j, k);
488  }
489  tangent_bundle_error += tangent_inner_prod * tangent_inner_prod;
490  }
491  double normal_inner_prod = 0;
492  for (int k=0; k<dimension; ++k) {
493  normal_inner_prod += coord[k] * my_GMLS_scalar.getTangentBundle(i, dimension-1, k);
494  }
495  // inner product could be plus or minus 1 (depends on tangent direction ordering)
496  double normal_error_to_sum = (normal_inner_prod > 0) ? normal_inner_prod - 1 : normal_inner_prod + 1;
497  tangent_bundle_error += normal_error_to_sum * normal_error_to_sum;
498  tangent_bundle_norm += 1;
499 
500  // evaluation of various exact solutions
501  double actual_value = sphere_harmonic54(xval, yval, zval);
502  double actual_Laplacian = laplace_beltrami_sphere_harmonic54(xval, yval, zval);
503  double actual_Gradient_ambient[3] = {0,0,0}; // initialized for 3, but only filled up to dimension
504  gradient_sphereHarmonic54_ambient(actual_Gradient_ambient, xval, yval, zval);
505  //velocity_sphereHarmonic54_ambient(actual_Gradient_ambient, xval, yval, zval);
506 
507  values_error += (GMLS_value - actual_value)*(GMLS_value - actual_value);
508  values_norm += actual_value*actual_value;
509 
510  laplacian_error += (GMLS_Laplacian - actual_Laplacian)*(GMLS_Laplacian - actual_Laplacian);
511  laplacian_norm += actual_Laplacian*actual_Laplacian;
512 
513  double actual_gc = 1.0;
514  gc_error += (GMLS_gc - actual_gc)*(GMLS_gc - actual_gc);
515  gc_norm += actual_gc*actual_gc;
516 
517  for (int j=0; j<dimension; ++j) {
518  gradient_ambient_error += (output_gradient(i,j) - actual_Gradient_ambient[j])*(output_gradient(i,j) - actual_Gradient_ambient[j]);
519  gradient_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
520  }
521 
522  for (int j=0; j<dimension; ++j) {
523  vector_ambient_error += (output_vector(i,j) - actual_Gradient_ambient[j])*(output_vector(i,j) - actual_Gradient_ambient[j]);
524  vector_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
525  }
526 
527  divergence_ambient_error += (output_divergence(i) - actual_Laplacian)*(output_divergence(i) - actual_Laplacian);
528  divergence_ambient_norm += actual_Laplacian*actual_Laplacian;
529 
530  for (int j=0; j<dimension; ++j) {
531  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]);
532  vector_of_scalar_clones_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
533  }
534 
535  divergence_of_scalar_clones_ambient_error += (output_divergence_of_scalar_clones(i) - actual_Laplacian)*(output_divergence_of_scalar_clones(i) - actual_Laplacian);
536  divergence_of_scalar_clones_ambient_norm += actual_Laplacian*actual_Laplacian;
537 
538  }
539 
540  tangent_bundle_error /= number_target_coords;
541  tangent_bundle_error = std::sqrt(tangent_bundle_error);
542  tangent_bundle_norm /= number_target_coords;
543  tangent_bundle_norm = std::sqrt(tangent_bundle_norm);
544 
545  values_error /= number_target_coords;
546  values_error = std::sqrt(values_error);
547  values_norm /= number_target_coords;
548  values_norm = std::sqrt(values_norm);
549 
550  laplacian_error /= number_target_coords;
551  laplacian_error = std::sqrt(laplacian_error);
552  laplacian_norm /= number_target_coords;
553  laplacian_norm = std::sqrt(laplacian_norm);
554 
555  gradient_ambient_error /= number_target_coords;
556  gradient_ambient_error = std::sqrt(gradient_ambient_error);
557  gradient_ambient_norm /= number_target_coords;
558  gradient_ambient_norm = std::sqrt(gradient_ambient_norm);
559 
560  gc_error /= number_target_coords;
561  gc_error = std::sqrt(gc_error);
562  gc_norm /= number_target_coords;
563  gc_norm = std::sqrt(gc_norm);
564 
565  vector_ambient_error /= number_target_coords;
566  vector_ambient_error = std::sqrt(vector_ambient_error);
567  vector_ambient_norm /= number_target_coords;
568  vector_ambient_norm = std::sqrt(vector_ambient_norm);
569 
570  divergence_ambient_error /= number_target_coords;
571  divergence_ambient_error = std::sqrt(divergence_ambient_error);
572  divergence_ambient_norm /= number_target_coords;
573  divergence_ambient_norm = std::sqrt(divergence_ambient_norm);
574 
575  vector_of_scalar_clones_ambient_error /= number_target_coords;
576  vector_of_scalar_clones_ambient_error = std::sqrt(vector_of_scalar_clones_ambient_error);
577  vector_of_scalar_clones_ambient_norm /= number_target_coords;
578  vector_of_scalar_clones_ambient_norm = std::sqrt(vector_of_scalar_clones_ambient_norm);
579 
580  divergence_of_scalar_clones_ambient_error /= number_target_coords;
581  divergence_of_scalar_clones_ambient_error = std::sqrt(divergence_of_scalar_clones_ambient_error);
582  divergence_of_scalar_clones_ambient_norm /= number_target_coords;
583  divergence_of_scalar_clones_ambient_norm = std::sqrt(divergence_of_scalar_clones_ambient_norm);
584 
585  printf("Tangent Bundle Error: %g\n", tangent_bundle_error / tangent_bundle_norm);
586  printf("Point Value Error: %g\n", values_error / values_norm);
587  printf("Laplace-Beltrami Error: %g\n", laplacian_error / laplacian_norm);
588  printf("Gaussian Curvature Error: %g\n", gc_error / gc_norm);
589  printf("Surface Gradient (Ambient) Error: %g\n", gradient_ambient_error / gradient_ambient_norm);
590  printf("Surface Vector (VectorBasis) Error: %g\n", vector_ambient_error / vector_ambient_norm);
591  printf("Surface Divergence (VectorBasis) Error: %g\n", divergence_ambient_error / divergence_ambient_norm);
592  printf("Surface Vector (ScalarClones) Error: %g\n",
593  vector_of_scalar_clones_ambient_error / vector_of_scalar_clones_ambient_norm);
594  printf("Surface Divergence (ScalarClones) Error: %g\n",
595  divergence_of_scalar_clones_ambient_error / divergence_of_scalar_clones_ambient_norm);
596  //! [Check That Solutions Are Correct]
597  // popRegion hidden from tutorial
598  // stop timing comparison loop
599  Kokkos::Profiling::popRegion();
600  //! [Finalize Program]
601 
602 
603 } // end of code block to reduce scope, causing Kokkos View de-allocations
604 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
605 
606 // finalize Kokkos and MPI (if available)
607 Kokkos::finalize();
608 #ifdef COMPADRE_USE_MPI
609 MPI_Finalize();
610 #endif
611 
612 return 0;
613 
614 } // main
615 
616 
617 //! [Finalize Program]
int main(int argc, char *args[])
[Parse Command Line Arguments]
KOKKOS_INLINE_FUNCTION double sphere_harmonic54(double x, double y, double z)
Lightweight Evaluator Helper This class is a lightweight wrapper for extracting and applying all rele...
constexpr SamplingFunctional ManifoldVectorPointSample
Point evaluations of the entire vector source function (but on a manifold, so it includes a transform...
Point evaluation of a vector (reconstructs entire vector at once, requiring a ReconstructionSpace hav...
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...
#define PI
Point evaluation of Gaussian curvature.
Point evaluation of a scalar.
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...
Scalar basis reused as many times as there are components in the vector resulting in a much cheaper p...
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)
Point evaluation of the laplacian of a scalar (could be on a manifold or not)
Point evaluation of the divergence of a vector (results in a scalar)
KOKKOS_INLINE_FUNCTION void gradient_sphereHarmonic54_ambient(double *gradient, double x, double y, double z)
Point evaluation of the gradient of a scalar.
Generalized Moving Least Squares (GMLS)
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) ...
KOKKOS_INLINE_FUNCTION double laplace_beltrami_sphere_harmonic54(double x, double y, double z)
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...