Compadre  1.2.0
GMLS_Staggered_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 
16 #ifdef COMPADRE_USE_MPI
17 #include <mpi.h>
18 #endif
19 
20 #include <Kokkos_Timer.hpp>
21 #include <Kokkos_Core.hpp>
22 
23 using namespace Compadre;
24 
25 //! [Parse Command Line Arguments]
26 
27 // called from command line
28 int main (int argc, char* args[]) {
29 
30 // initializes MPI (if available) with command line arguments given
31 #ifdef COMPADRE_USE_MPI
32 MPI_Init(&argc, &args);
33 #endif
34 
35 // initializes Kokkos with command line arguments given
36 Kokkos::initialize(argc, args);
37 
38 // code block to reduce scope for all Kokkos View allocations
39 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
40 {
41  // check if 8 arguments are given from the command line, the first being the program name
42  // constraint_type used in solving each GMLS problem:
43  // 0 - No constraints used in solving each GMLS problem
44  // 1 - Neumann Gradient Scalar used in solving each GMLS problem
45  int constraint_type = 0; // No constraints by default
46  if (argc >= 8) {
47  int arg8toi = atoi(args[7]);
48  if (arg8toi > 0) {
49  constraint_type = arg8toi;
50  }
51  }
52 
53  // check if 7 arguments are given from the command line, the first being the program name
54  // problem_type used in solving each GMLS problem:
55  // 0 - Standard GMLS problem
56  // 1 - Manifold GMLS problem
57  int problem_type = 1; // Manifold for this example
58  if (argc >= 7) {
59  int arg7toi = atoi(args[6]);
60  if (arg7toi > 0) {
61  problem_type = arg7toi;
62  }
63  }
64 
65  // check if 6 arguments are given from the command line, the first being the program name
66  // solver_type used for factorization in solving each GMLS problem:
67  // 0 - SVD used for factorization in solving each GMLS problem
68  // 1 - QR used for factorization in solving each GMLS problem
69  // 2 - LU used for factorization in solving each GMLS problem
70  int solver_type = 1; // QR by default
71  if (argc >= 6) {
72  int arg6toi = atoi(args[5]);
73  if (arg6toi >= 0) {
74  solver_type = arg6toi;
75  }
76  }
77 
78  // check if 5 arguments are given from the command line, the first being the program name
79  // N_pts_on_sphere used to determine spatial resolution
80  int N_pts_on_sphere = 1000; // 1000 points by default
81  if (argc >= 5) {
82  int arg5toi = atoi(args[4]);
83  if (arg5toi > 0) {
84  N_pts_on_sphere = arg5toi;
85  }
86  }
87 
88  // check if 4 arguments are given from the command line
89  // dimension for the coordinates and the solution
90  int dimension = 3; // dimension 3 by default
91  if (argc >= 4) {
92  int arg4toi = atoi(args[3]);
93  if (arg4toi > 0) {
94  dimension = arg4toi;
95  }
96  }
97 
98  // check if 3 arguments are given from the command line
99  // set the number of target sites where we will reconstruct the target functionals at
100  int number_target_coords = 200; // 200 target sites by default
101  if (argc >= 3) {
102  int arg3toi = atoi(args[2]);
103  if (arg3toi > 0) {
104  number_target_coords = arg3toi;
105  }
106  }
107 
108  // check if 2 arguments are given from the command line
109  // set the number of target sites where we will reconstruct the target functionals at
110  int order = 3; // 3rd degree polynomial basis by default
111  if (argc >= 2) {
112  int arg2toi = atoi(args[1]);
113  if (arg2toi > 0) {
114  order = arg2toi;
115  }
116  }
117 
118  // minimum neighbors for unisolvency is the same as the size of the polynomial basis
119  // dimension has one subtracted because it is a D-1 manifold represented in D dimensions
120  const int min_neighbors = Compadre::GMLS::getNP(order, dimension-1);
121 
122  //! [Parse Command Line Arguments]
123  Kokkos::Timer timer;
124  Kokkos::Profiling::pushRegion("Setup Point Data");
125  //! [Setting Up The Point Cloud]
126 
127 
128  // coordinates of source sites, bigger than needed then resized later
129  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device("source coordinates",
130  1.25*N_pts_on_sphere, 3);
131  Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
132 
133  double r = 1.0;
134 
135  // set number of source coordinates from what was calculated
136  int number_source_coords;
137 
138  { // fill source coordinates from surface of a sphere with quasiuniform points
139  // algorithm described at https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf
140  int N_count = 0;
141  double a = 4*PI*r*r/N_pts_on_sphere;
142  double d = std::sqrt(a);
143  int M_theta = std::round(PI/d);
144  double d_theta = PI/M_theta;
145  double d_phi = a/d_theta;
146  for (int i=0; i<M_theta; ++i) {
147  double theta = PI*(i + 0.5)/M_theta;
148  int M_phi = std::round(2*PI*std::sin(theta)/d_phi);
149  for (int j=0; j<M_phi; ++j) {
150  double phi = 2*PI*j/M_phi;
151  source_coords(N_count, 0) = theta;
152  source_coords(N_count, 1) = phi;
153  N_count++;
154  }
155  }
156 
157  for (int i=0; i<N_count; ++i) {
158  double theta = source_coords(i,0);
159  double phi = source_coords(i,1);
160  source_coords(i,0) = r*std::sin(theta)*std::cos(phi);
161  source_coords(i,1) = r*std::sin(theta)*std::sin(phi);
162  source_coords(i,2) = r*cos(theta);
163  //printf("%f %f %f\n", source_coords(i,0), source_coords(i,1), source_coords(i,2));
164  }
165  number_source_coords = N_count;
166  }
167 
168  // resize source_coords to the size actually needed
169  Kokkos::resize(source_coords, number_source_coords, 3);
170  Kokkos::resize(source_coords_device, number_source_coords, 3);
171 
172  // coordinates of target sites
173  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device("target coordinates",
174  number_target_coords, 3);
175  Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
176 
177  // seed random number generator
178  std::mt19937 rng(50);
179 
180  // generate random integers in [0...number_source_coords-1] (used to pick target sites)
181  std::uniform_int_distribution<int> gen_num_neighbors(0, number_source_coords-1); // uniform, unbiased
182 
183  // fill target sites with random selections from source sites
184  for (int i=0; i<number_target_coords; ++i) {
185  const int source_site_to_copy = gen_num_neighbors(rng);
186  for (int j=0; j<3; ++j) {
187  target_coords(i,j) = source_coords(source_site_to_copy,j);
188  }
189  }
190 
191 
192  //! [Setting Up The Point Cloud]
193 
194  Kokkos::Profiling::popRegion();
195  Kokkos::fence();
196  Kokkos::Profiling::pushRegion("Creating Data");
197 
198  //! [Creating The Data]
199 
200 
201  // source coordinates need copied to device before using to construct sampling data
202  Kokkos::deep_copy(source_coords_device, source_coords);
203  Kokkos::deep_copy(target_coords_device, target_coords);
204 
205  // ensure that source coordinates are sent to device before evaluating sampling data based on them
206  Kokkos::fence();
207 
208 
209  // need Kokkos View storing true solution (for samples)
210  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device("samples of true solution",
211  source_coords_device.extent(0));
212 
213  // need Kokkos View storing true vector solution (for samples)
214  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> sampling_vector_data_device("samples of vector true solution",
215  source_coords_device.extent(0), 3);
216 
217  Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
218  (0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
219 
220  // coordinates of source site i
221  double xval = source_coords_device(i,0);
222  double yval = (dimension>1) ? source_coords_device(i,1) : 0;
223  double zval = (dimension>2) ? source_coords_device(i,2) : 0;
224 
225  // data for targets with scalar input
226  sampling_data_device(i) = sphere_harmonic54(xval, yval, zval);
227  //printf("%f\n", sampling_data_device(i));
228 
229  for (int j=0; j<3; ++j) {
230  double gradient[3] = {0,0,0};
231  gradient_sphereHarmonic54_ambient(gradient, xval, yval, zval);
232  sampling_vector_data_device(i,j) = gradient[j];
233  }
234  //printf("%f %f %f\n", sampling_vector_data_device(i,0), sampling_vector_data_device(i,1), sampling_vector_data_device(i,2));
235  });
236 
237 
238  //! [Creating The Data]
239 
240  Kokkos::Profiling::popRegion();
241  Kokkos::Profiling::pushRegion("Neighbor Search");
242 
243  //! [Performing Neighbor Search]
244 
245 
246  // Point cloud construction for neighbor search
247  // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
248  auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
249 
250  // each row is a neighbor list for a target site, with the first column of each row containing
251  // the number of neighbors for that rows corresponding target site
252  double epsilon_multiplier = 1.5;
253  int estimated_upper_bound_number_neighbors =
254  point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
255 
256  Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device("neighbor lists",
257  number_target_coords, estimated_upper_bound_number_neighbors); // first column is # of neighbors
258  Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
259 
260  // each target site has a window size
261  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device("h supports", number_target_coords);
262  Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
263 
264  // query the point cloud to generate the neighbor lists using a kdtree to produce the n nearest neighbor
265  // to each target site, adding (epsilon_multiplier-1)*100% to whatever the distance away the further neighbor used is from
266  // each target to the view for epsilon
267  point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/, target_coords, neighbor_lists,
268  epsilon, min_neighbors, epsilon_multiplier);
269 
270 
271  //! [Performing Neighbor Search]
272 
273  Kokkos::Profiling::popRegion();
274  Kokkos::fence(); // let call to build neighbor lists complete before copying back to device
275  timer.reset();
276 
277  //! [Setting Up The GMLS Object]
278 
279 
280  // Copy data back to device (they were filled on the host)
281  // We could have filled Kokkos Views with memory space on the host
282  // and used these instead, and then the copying of data to the device
283  // would be performed in the GMLS class
284  Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
285  Kokkos::deep_copy(epsilon_device, epsilon);
286 
287  // solver name for passing into the GMLS class
288  std::string solver_name;
289  if (solver_type == 0) { // SVD
290  solver_name = "SVD";
291  } else if (solver_type == 1) { // QR
292  solver_name = "QR";
293  } else if (solver_type == 2) { // LU
294  solver_name = "LU";
295  }
296 
297  // problem name for passing into the GMLS class
298  std::string problem_name;
299  if (problem_type == 0) { // Standard
300  problem_name = "STANDARD";
301  } else if (problem_type == 1) { // Manifold
302  problem_name = "MANIFOLD";
303  }
304 
305  // boundary name for passing into the GMLS class
306  std::string constraint_name;
307  if (constraint_type == 0) { // No constraints
308  constraint_name = "NO_CONSTRAINT";
309  } else if (constraint_type == 1) { // Neumann Gradient Scalar
310  constraint_name = "NEUMANN_GRAD_SCALAR";
311  }
312 
313  // initialize an instance of the GMLS class
316  order, dimension,
317  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
318  order /*manifold order*/);
319 
320  // pass in neighbor lists, source coordinates, target coordinates, and window sizes
321  //
322  // neighbor lists have the format:
323  // dimensions: (# number of target sites) X (# maximum number of neighbors for any given target + 1)
324  // the first column contains the number of neighbors for that rows corresponding target index
325  //
326  // source coordinates have the format:
327  // dimensions: (# number of source sites) X (dimension)
328  // entries in the neighbor lists (integers) correspond to rows of this 2D array
329  //
330  // target coordinates have the format:
331  // dimensions: (# number of target sites) X (dimension)
332  // # of target sites is same as # of rows of neighbor lists
333  //
334  my_GMLS_vector_1.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
335 
336  // create a vector of target operations
337  std::vector<TargetOperation> lro_vector_1(1);
338  lro_vector_1[0] = DivergenceOfVectorPointEvaluation;
339 
340  // and then pass them to the GMLS class
341  my_GMLS_vector_1.addTargets(lro_vector_1);
342 
343  // sets the weighting kernel function from WeightingFunctionType for curvature
344  my_GMLS_vector_1.setCurvatureWeightingType(WeightingFunctionType::Power);
345 
346  // power to use in the weighting kernel function for curvature coefficients
347  my_GMLS_vector_1.setCurvatureWeightingPower(2);
348 
349  // sets the weighting kernel function from WeightingFunctionType
350  my_GMLS_vector_1.setWeightingType(WeightingFunctionType::Power);
351 
352  // power to use in that weighting kernel function
353  my_GMLS_vector_1.setWeightingPower(2);
354 
355  // setup quadrature for StaggeredEdgeIntegralSample
356  my_GMLS_vector_1.setOrderOfQuadraturePoints(2);
357  my_GMLS_vector_1.setDimensionOfQuadraturePoints(1);
358  my_GMLS_vector_1.setQuadratureType("LINE");
359 
360  // generate the alphas that to be combined with data for each target operation requested in lro
361  my_GMLS_vector_1.generateAlphas();
362 
363  // initialize another instance of the GMLS class
367  order, dimension,
368  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
369  order /*manifold order*/);
370 
371  my_GMLS_vector_2.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
372  std::vector<TargetOperation> lro_vector_2(2);
374  lro_vector_2[1] = DivergenceOfVectorPointEvaluation;
375  //lro_vector_2[2] = GradientOfScalarPointEvaluation;
376  my_GMLS_vector_2.addTargets(lro_vector_2);
377  my_GMLS_vector_2.setCurvatureWeightingType(WeightingFunctionType::Power);
378  my_GMLS_vector_2.setCurvatureWeightingPower(2);
379  my_GMLS_vector_2.setWeightingType(WeightingFunctionType::Power);
380  my_GMLS_vector_2.setWeightingPower(2);
381  my_GMLS_vector_2.setOrderOfQuadraturePoints(2);
382  my_GMLS_vector_2.setDimensionOfQuadraturePoints(1);
383  my_GMLS_vector_2.setQuadratureType("LINE");
384  my_GMLS_vector_2.generateAlphas();
385 
386  // initialize another instance of the GMLS class
389  order, dimension,
390  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
391  order /*manifold order*/);
392 
393  my_GMLS_scalar.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
394 
395  std::vector<TargetOperation> lro_scalar(1);
397  //lro_scalar[1] = GradientOfScalarPointEvaluation;
398  my_GMLS_scalar.addTargets(lro_scalar);
399  my_GMLS_scalar.setCurvatureWeightingType(WeightingFunctionType::Power);
400  my_GMLS_scalar.setCurvatureWeightingPower(2);
401  my_GMLS_scalar.setWeightingType(WeightingFunctionType::Power);
402  my_GMLS_scalar.setWeightingPower(2);
403  my_GMLS_scalar.generateAlphas();
404 
405 
406  //! [Setting Up The GMLS Object]
407 
408  double instantiation_time = timer.seconds();
409  std::cout << "Took " << instantiation_time << "s to complete alphas generation." << std::endl;
410  Kokkos::fence(); // let generateAlphas finish up before using alphas
411  Kokkos::Profiling::pushRegion("Apply Alphas to Data");
412 
413  //! [Apply GMLS Alphas To Data]
414 
415 
416  // it is important to note that if you expect to use the data as a 1D view, then you should use double*
417  // however, if you know that the target operation will result in a 2D view (vector or matrix output),
418  // then you should template with double** as this is something that can not be infered from the input data
419  // or the target operator at compile time. Additionally, a template argument is required indicating either
420  // Kokkos::HostSpace or Kokkos::DefaultExecutionSpace::memory_space()
421 
422  // The Evaluator class takes care of handling input data views as well as the output data views.
423  // It uses information from the GMLS class to determine how many components are in the input
424  // as well as output for any choice of target functionals and then performs the contactions
425  // on the data using the alpha coefficients generated by the GMLS class, all on the device.
426  Evaluator vector_1_gmls_evaluator(&my_GMLS_vector_1);
427  Evaluator vector_2_gmls_evaluator(&my_GMLS_vector_2);
428  Evaluator scalar_gmls_evaluator(&my_GMLS_scalar);
429 
430 
431  //auto output_gradient_vectorbasis =
432  // vector_2_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
433  // (sampling_data_device, GradientOfScalarPointEvaluation, StaggeredEdgeAnalyticGradientIntegralSample);
434 
435  //auto output_gradient_scalarbasis =
436  // scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
437  // (sampling_data_device, GradientOfScalarPointEvaluation, StaggeredEdgeAnalyticGradientIntegralSample);
438 
439  auto output_divergence_vectorsamples =
440  vector_1_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
441  (sampling_vector_data_device, DivergenceOfVectorPointEvaluation, StaggeredEdgeIntegralSample);
442 
443  auto output_divergence_scalarsamples =
444  vector_2_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
446 
447  auto output_laplacian_vectorbasis =
448  vector_2_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
450 
451  auto output_laplacian_scalarbasis =
452  scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
454 
455 
456  //! [Apply GMLS Alphas To Data]
457 
458  Kokkos::fence(); // let application of alphas to data finish before using results
459  Kokkos::Profiling::popRegion();
460  // times the Comparison in Kokkos
461  Kokkos::Profiling::pushRegion("Comparison");
462 
463  //! [Check That Solutions Are Correct]
464 
465  double laplacian_vectorbasis_error = 0;
466  double laplacian_vectorbasis_norm = 0;
467 
468  double laplacian_scalarbasis_error = 0;
469  double laplacian_scalarbasis_norm = 0;
470 
471  double gradient_vectorbasis_ambient_error = 0;
472  double gradient_vectorbasis_ambient_norm = 0;
473 
474  double gradient_scalarbasis_ambient_error = 0;
475  double gradient_scalarbasis_ambient_norm = 0;
476 
477  double divergence_vectorsamples_ambient_error = 0;
478  double divergence_vectorsamples_ambient_norm = 0;
479 
480  double divergence_scalarsamples_ambient_error = 0;
481  double divergence_scalarsamples_ambient_norm = 0;
482 
483  // loop through the target sites
484  for (int i=0; i<number_target_coords; i++) {
485 
486  // target site i's coordinate
487  double xval = target_coords(i,0);
488  double yval = (dimension>1) ? target_coords(i,1) : 0;
489  double zval = (dimension>2) ? target_coords(i,2) : 0;
490 
491  // evaluation of various exact solutions
492  double actual_Laplacian = laplace_beltrami_sphere_harmonic54(xval, yval, zval);
493  double actual_Gradient_ambient[3] = {0,0,0}; // initialized for 3, but only filled up to dimension
494  gradient_sphereHarmonic54_ambient(actual_Gradient_ambient, xval, yval, zval);
495 
496  laplacian_vectorbasis_error += (output_laplacian_vectorbasis(i) - actual_Laplacian)*(output_laplacian_vectorbasis(i) - actual_Laplacian);
497  laplacian_vectorbasis_norm += actual_Laplacian*actual_Laplacian;
498 
499  //printf("Error of %f, %f vs %f\n", (output_laplacian_scalarbasis(i) - actual_Laplacian), output_laplacian_scalarbasis(i), actual_Laplacian);
500  laplacian_scalarbasis_error += (output_laplacian_scalarbasis(i) - actual_Laplacian)*(output_laplacian_scalarbasis(i) - actual_Laplacian);
501  laplacian_scalarbasis_norm += actual_Laplacian*actual_Laplacian;
502 
503  //for (int j=0; j<dimension; ++j) {
504  // //printf("VectorBasis Error of %f, %f vs %f\n", (output_gradient_vectorbasis(i,j) - actual_Gradient_ambient[j]), output_gradient_vectorbasis(i,j), actual_Gradient_ambient[j]);
505  // gradient_vectorbasis_ambient_error += (output_gradient_vectorbasis(i,j) - actual_Gradient_ambient[j])*(output_gradient_vectorbasis(i,j) - actual_Gradient_ambient[j]);
506  // gradient_vectorbasis_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
507  //}
508 
509  //for (int j=0; j<dimension; ++j) {
510  // //printf("ScalarBasis Error of %f, %f vs %f\n", (output_gradient_scalarbasis(i,j) - actual_Gradient_ambient[j]), output_gradient_scalarbasis(i,j), actual_Gradient_ambient[j]);
511  // gradient_scalarbasis_ambient_error += (output_gradient_scalarbasis(i,j) - actual_Gradient_ambient[j])*(output_gradient_scalarbasis(i,j) - actual_Gradient_ambient[j]);
512  // gradient_scalarbasis_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
513  //}
514 
515  //printf("Error of %f, %f vs %f\n", (output_divergence(i) - actual_Laplacian), output_divergence(i), actual_Laplacian);
516  divergence_vectorsamples_ambient_error += (output_divergence_vectorsamples(i) - actual_Laplacian)*(output_divergence_vectorsamples(i) - actual_Laplacian);
517  divergence_vectorsamples_ambient_norm += actual_Laplacian*actual_Laplacian;
518 
519  divergence_scalarsamples_ambient_error += (output_divergence_scalarsamples(i) - actual_Laplacian)*(output_divergence_scalarsamples(i) - actual_Laplacian);
520  divergence_scalarsamples_ambient_norm += actual_Laplacian*actual_Laplacian;
521 
522  }
523 
524  laplacian_vectorbasis_error /= number_target_coords;
525  laplacian_vectorbasis_error = std::sqrt(laplacian_vectorbasis_error);
526  laplacian_vectorbasis_norm /= number_target_coords;
527  laplacian_vectorbasis_norm = std::sqrt(laplacian_vectorbasis_norm);
528 
529  laplacian_scalarbasis_error /= number_target_coords;
530  laplacian_scalarbasis_error = std::sqrt(laplacian_scalarbasis_error);
531  laplacian_scalarbasis_norm /= number_target_coords;
532  laplacian_scalarbasis_norm = std::sqrt(laplacian_scalarbasis_norm);
533 
534  gradient_vectorbasis_ambient_error /= number_target_coords;
535  gradient_vectorbasis_ambient_error = std::sqrt(gradient_vectorbasis_ambient_error);
536  gradient_vectorbasis_ambient_norm /= number_target_coords;
537  gradient_vectorbasis_ambient_norm = std::sqrt(gradient_vectorbasis_ambient_norm);
538 
539  gradient_scalarbasis_ambient_error /= number_target_coords;
540  gradient_scalarbasis_ambient_error = std::sqrt(gradient_scalarbasis_ambient_error);
541  gradient_scalarbasis_ambient_norm /= number_target_coords;
542  gradient_scalarbasis_ambient_norm = std::sqrt(gradient_scalarbasis_ambient_norm);
543 
544  divergence_vectorsamples_ambient_error /= number_target_coords;
545  divergence_vectorsamples_ambient_error = std::sqrt(divergence_vectorsamples_ambient_error);
546  divergence_vectorsamples_ambient_norm /= number_target_coords;
547  divergence_vectorsamples_ambient_norm = std::sqrt(divergence_vectorsamples_ambient_norm);
548 
549  divergence_scalarsamples_ambient_error /= number_target_coords;
550  divergence_scalarsamples_ambient_error = std::sqrt(divergence_scalarsamples_ambient_error);
551  divergence_scalarsamples_ambient_norm /= number_target_coords;
552  divergence_scalarsamples_ambient_norm = std::sqrt(divergence_scalarsamples_ambient_norm);
553 
554  printf("Staggered Laplace-Beltrami (VectorBasis) Error: %g\n", laplacian_vectorbasis_error / laplacian_vectorbasis_norm);
555  printf("Staggered Laplace-Beltrami (ScalarBasis) Error: %g\n", laplacian_scalarbasis_error / laplacian_scalarbasis_norm);
556  printf("Surface Staggered Gradient (VectorBasis) Error: %g\n", gradient_vectorbasis_ambient_error / gradient_vectorbasis_ambient_norm);
557  printf("Surface Staggered Gradient (ScalarBasis) Error: %g\n", gradient_scalarbasis_ambient_error / gradient_scalarbasis_ambient_norm);
558  printf("Surface Staggered Divergence (VectorSamples) Error: %g\n", divergence_vectorsamples_ambient_error / divergence_vectorsamples_ambient_norm);
559  printf("Surface Staggered Divergence (ScalarSamples) Error: %g\n", divergence_scalarsamples_ambient_error / divergence_scalarsamples_ambient_norm);
560  //! [Check That Solutions Are Correct]
561  // popRegion hidden from tutorial
562  // stop timing comparison loop
563  Kokkos::Profiling::popRegion();
564  //! [Finalize Program]
565 
566 
567 } // end of code block to reduce scope, causing Kokkos View de-allocations
568 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
569 
570 // finalize Kokkos and MPI (if available)
571 Kokkos::finalize();
572 #ifdef COMPADRE_USE_MPI
573 MPI_Finalize();
574 #endif
575 
576 return 0;
577 
578 } // main
579 
580 
581 //! [Finalize Program]
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...
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
constexpr SamplingFunctional StaggeredEdgeAnalyticGradientIntegralSample
Analytical integral of a gradient source vector is just a difference of the scalar source at neighbor...
Scalar polynomial basis centered at the target site and scaled by sum of basis powers e...
int main(int argc, char *args[])
[Parse Command Line Arguments]
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...
constexpr SamplingFunctional StaggeredEdgeIntegralSample
Samples consist of the result of integrals of a vector dotted with the tangent along edges between ne...
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.
Point evaluation of the chained staggered Laplacian acting on VectorTaylorPolynomial basis + Staggere...
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)
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...