Compadre  1.3.3
GMLS_Manifold_Multiple_Evaluation_Sites.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 //! [Ambient to Local Back To Ambient Helper Function]
26 void AmbientLocalAmbient(XYZ& u, double* T_data, double* P_data) {
27  // reconstructions with vector output on a manifold are defined
28  // in their local chart (2D). Next, they are mapped to 3D using
29  // the tangent vector calculated at the target site.
30  //
31  // We are interested in a mapping using the tangent vector at
32  // additional evaluation site instead, so we map back to the
33  // local chart using the tangent vector defined at the target
34  // site (T), then mapping from the local chart to ambient space
35  // using the tangent vector calculated at the extra site (P).
36 
37 
38  host_scratch_matrix_right_type T(T_data, 3, 3);
39  host_scratch_matrix_right_type P(P_data, 3, 3);
40 
41  // first get back to local chart
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];
46  }
47  // second go to ambient space using tangent for first neighbor
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];
52  }
53 
54 }
55 
56 //! [Ambient to Local Back To Ambient Helper Function]
57 
58 // called from command line
59 int main (int argc, char* args[]) {
60 
61 // initializes MPI (if available) with command line arguments given
62 #ifdef COMPADRE_USE_MPI
63 MPI_Init(&argc, &args);
64 #endif
65 
66 // initializes Kokkos with command line arguments given
67 Kokkos::initialize(argc, args);
68 
69 // code block to reduce scope for all Kokkos View allocations
70 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
71 {
72 
73  CommandLineProcessor clp(argc, args);
74  auto order = clp.order;
75  auto dimension = clp.dimension;
76  auto number_target_coords = clp.number_target_coords;
77  auto constraint_name = clp.constraint_name;
78  auto solver_name = clp.solver_name;
79  auto problem_name = clp.problem_name;
80  int N_pts_on_sphere = (clp.number_source_coords>=0) ? clp.number_source_coords : 1000;
81 
82  // minimum neighbors for unisolvency is the same as the size of the polynomial basis
83  // dimension has one subtracted because it is a D-1 manifold represented in D dimensions
84  const int min_neighbors = Compadre::GMLS::getNP(order, dimension-1);
85 
86  //! [Parse Command Line Arguments]
87  Kokkos::Timer timer;
88  Kokkos::Profiling::pushRegion("Setup Point Data");
89  //! [Setting Up The Point Cloud]
90 
91 
92  // coordinates of source sites, bigger than needed then resized later
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);
96 
97  double r = 1.0;
98 
99  // set number of source coordinates from what was calculated
100  int number_source_coords;
101 
102  {
103  // fill source coordinates from surface of a sphere with quasiuniform points
104  // algorithm described at https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf
105  int N_count = 0;
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;
118  N_count++;
119  }
120  }
121 
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);
128  //printf("%f %f %f\n", source_coords(i,0), source_coords(i,1), source_coords(i,2));
129  }
130  number_source_coords = N_count;
131  }
132 
133  // resize source_coords to the size actually needed
134  Kokkos::resize(source_coords, number_source_coords, 3);
135  Kokkos::resize(source_coords_device, number_source_coords, 3);
136 
137  // coordinates of target sites
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);
141 
142  {
143  bool enough_pts_found = false;
144  // fill target coordinates from surface of a sphere with quasiuniform points
145  // stop when enough points are found
146  int N_count = 0;
147  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)
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;
159  N_count++;
160  if (N_count == number_target_coords) {
161  enough_pts_found = true;
162  break;
163  }
164  }
165  if (enough_pts_found) break;
166  }
167 
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);
174  //printf("%f %f %f\n", target_coords(i,0), target_coords(i,1), target_coords(i,2));
175  }
176  }
177 
178 
179  //! [Setting Up The Point Cloud]
180 
181  Kokkos::Profiling::popRegion();
182  Kokkos::fence();
183  Kokkos::Profiling::pushRegion("Creating Data");
184 
185  //! [Creating The Data]
186 
187 
188  // source coordinates need copied to device before using to construct sampling data
189  Kokkos::deep_copy(source_coords_device, source_coords);
190 
191  // target coordinates copied next, because it is a convenient time to send them to device
192  Kokkos::deep_copy(target_coords_device, target_coords);
193 
194  // ensure that source coordinates are sent to device before evaluating sampling data based on them
195  Kokkos::fence();
196 
197 
198  // need Kokkos View storing true solution (for samples)
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);
204 
205  // need Kokkos View storing true vector solution (for samples)
206  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> sampling_vector_data_device("samples of vector true solution",
207  source_coords_device.extent(0), 3);
208 
209  Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
210  (0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
211 
212  // coordinates of source site 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;
216 
217  // data for targets with scalar input
218  sampling_data_device(i) = sphere_harmonic54(xval, yval, zval);
219 
220  for (int j=0; j<3; ++j) {
221  double gradient[3] = {0,0,0};
222  gradient_sphereHarmonic54_ambient(gradient, xval, yval, zval);
223  sampling_vector_data_device(i,j) = gradient[j];
224  }
225  });
226 
227 
228  //! [Creating The Data]
229 
230  Kokkos::Profiling::popRegion();
231  Kokkos::Profiling::pushRegion("Neighbor Search");
232 
233  //! [Performing Neighbor Search]
234 
235 
236  // Point cloud construction for neighbor search
237  // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
238  auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
239 
240  // each row is a neighbor list for a target site, with the first column of each row containing
241  // the number of neighbors for that rows corresponding target site
242  double epsilon_multiplier = 1.7;
243  int estimated_upper_bound_number_neighbors =
244  point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
245 
246  Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device("neighbor lists",
247  number_target_coords, estimated_upper_bound_number_neighbors); // first column is # of neighbors
248  Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
249 
250  // each target site has a window size
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);
253 
254  // query the point cloud to generate the neighbor lists using a kdtree to produce the n nearest neighbor
255  // to each target site, adding (epsilon_multiplier-1)*100% to whatever the distance away the further neighbor used is from
256  // each target to the view for epsilon
257  point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/, target_coords, neighbor_lists,
258  epsilon, min_neighbors, epsilon_multiplier);
259 
260  //! [Performing Neighbor Search]
261 
262  Kokkos::Profiling::popRegion();
263  Kokkos::fence(); // let call to build neighbor lists complete before copying back to device
264  timer.reset();
265 
266  //! [Setting Up The GMLS Object]
267 
268 
269  // Copy data back to device (they were filled on the host)
270  // We could have filled Kokkos Views with memory space on the host
271  // and used these instead, and then the copying of data to the device
272  // would be performed in the GMLS class
273  Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
274  Kokkos::deep_copy(epsilon_device, epsilon);
275 
276  // initialize an instance of the GMLS class for problems with a scalar basis and
277  // traditional point sampling as the sampling functional
278  GMLS my_GMLS_scalar(order, dimension,
279  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
280  order /*manifold order*/);
281 
282  // pass in neighbor lists, source coordinates, target coordinates, and window sizes
283  //
284  // neighbor lists have the format:
285  // dimensions: (# number of target sites) X (# maximum number of neighbors for any given target + 1)
286  // the first column contains the number of neighbors for that rows corresponding target index
287  //
288  // source coordinates have the format:
289  // dimensions: (# number of source sites) X (dimension)
290  // entries in the neighbor lists (integers) correspond to rows of this 2D array
291  //
292  // target coordinates have the format:
293  // dimensions: (# number of target sites) X (dimension)
294  // # of target sites is same as # of rows of neighbor lists
295  //
296  my_GMLS_scalar.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
297 
298  // set up additional sites to evaluate target operators
299  // these sites will be the neighbors of the target site
300  my_GMLS_scalar.setAdditionalEvaluationSitesData(neighbor_lists_device, source_coords_device);
301 
302  // set a reference outward normal direction, causing a surface orientation after
303  // the GMLS instance computes an approximate tangent bundle
304  // on a sphere, the ambient coordinates are the outward normal direction
305  my_GMLS_scalar.setReferenceOutwardNormalDirection(target_coords_device, true /* use to orient surface */);
306 
307  // create a vector of target operations
308  std::vector<TargetOperation> lro_scalar(3);
309  lro_scalar[0] = ScalarPointEvaluation;
310  lro_scalar[1] = GaussianCurvaturePointEvaluation;
311  lro_scalar[2] = CurlOfVectorPointEvaluation;
312 
313  // and then pass them to the GMLS class
314  my_GMLS_scalar.addTargets(lro_scalar);
315 
316  // sets the weighting kernel function from WeightingFunctionType for curvature
317  my_GMLS_scalar.setCurvatureWeightingType(WeightingFunctionType::Power);
318 
319  // power to use in the weighting kernel function for curvature coefficients
320  my_GMLS_scalar.setCurvatureWeightingPower(2);
321 
322  // sets the weighting kernel function from WeightingFunctionType
323  my_GMLS_scalar.setWeightingType(WeightingFunctionType::Power);
324 
325  // power to use in that weighting kernel function
326  my_GMLS_scalar.setWeightingPower(2);
327 
328  // generate the alphas that to be combined with data for each target operation requested in lro
329  my_GMLS_scalar.generateAlphas();
330 
331  Kokkos::Profiling::pushRegion("Full Polynomial Basis GMLS Solution");
332  // initialize another instance of the GMLS class for problems with a vector basis on a manifold and point
333  // evaluation of that vector as the sampling functional
334  // VectorTaylorPolynomial indicates that the basis will be a polynomial with as many components as the
335  // dimension of the manifold. This differs from another possibility, which follows this class.
337  order, dimension,
338  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
339  order /*manifold order*/);
340 
341  my_GMLS_vector.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
342  my_GMLS_vector.setAdditionalEvaluationSitesData(neighbor_lists_device, source_coords_device);
343  std::vector<TargetOperation> lro_vector(2);
344  lro_vector[0] = VectorPointEvaluation;
345  //lro_vector[1] = DivergenceOfVectorPointEvaluation;
346  my_GMLS_vector.addTargets(lro_vector);
347  my_GMLS_vector.setCurvatureWeightingType(WeightingFunctionType::Power);
348  my_GMLS_vector.setCurvatureWeightingPower(2);
349  my_GMLS_vector.setWeightingType(WeightingFunctionType::Power);
350  my_GMLS_vector.setWeightingPower(2);
351  my_GMLS_vector.generateAlphas();
352  Kokkos::Profiling::popRegion();
353 
354  Kokkos::Profiling::pushRegion("Scalar Polynomial Basis Repeated to Form a Vector GMLS Solution");
355  // initialize another instance of the GMLS class for problems with a vector basis on a manifold and point
356  // evaluation of that vector as the sampling functional
357  // VectorOfScalarClonesTaylorPolynomial indicates a scalar polynomial will be solved for, since
358  // each component of the reconstructed vector are independent. This results in solving a smaller system
359  // for each GMLS problem, and is the suggested way to do vector reconstructions when sampling functionals
360  // acting on the basis would not create non-zero offdiagonal blocks.
361  //
362  // As an example, consider a 2D manifold in 3D ambient space. The GMLS problem is posed in the local chart,
363  // meaning that the problem being solved looks like
364  //
365  // [P_0 0 | where P_0 has dimension #number of neighbors for a target X #dimension of a scalar basis
366  // | 0 P_1]
367  //
368  // P_1 is similarly shaped, and for sampling functional that is a point evaluation, P_0 and P_1 are
369  // identical and their degrees of freedom in this system are disjoint, allowing us to solve for the
370  // degrees of freedom for either block independently. Additionally, the will produce the exact
371  // same polynomial coefficients for the point sampling functional, therefore it makes sense to use
372  // VectorOfScalarClonesTaylorPolynomial.
373  //
374  // In the print-out for this program, we include the timings and errors on this and VectorTaylorPolynomial
375  // in order to demonstrate that they produce exactly the same answer, but that one is much more efficient.
376  //
378  order, dimension,
379  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
380  order /*manifold order*/);
381 
382  my_GMLS_vector_of_scalar_clones.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
383  my_GMLS_vector_of_scalar_clones.setAdditionalEvaluationSitesData(neighbor_lists_device, source_coords_device);
384  std::vector<TargetOperation> lro_vector_of_scalar_clones(2);
385  lro_vector_of_scalar_clones[0] = VectorPointEvaluation;
386  //lro_vector_of_scalar_clones[1] = DivergenceOfVectorPointEvaluation;
387  my_GMLS_vector_of_scalar_clones.addTargets(lro_vector_of_scalar_clones);
388  my_GMLS_vector_of_scalar_clones.setCurvatureWeightingType(WeightingFunctionType::Power);
389  my_GMLS_vector_of_scalar_clones.setCurvatureWeightingPower(2);
390  my_GMLS_vector_of_scalar_clones.setWeightingType(WeightingFunctionType::Power);
391  my_GMLS_vector_of_scalar_clones.setWeightingPower(2);
392  my_GMLS_vector_of_scalar_clones.generateAlphas();
393  Kokkos::Profiling::popRegion();
394 
395 
396  //! [Setting Up The GMLS Object]
397 
398  double instantiation_time = timer.seconds();
399  std::cout << "Took " << instantiation_time << "s to complete alphas generation." << std::endl;
400  Kokkos::fence(); // let generateAlphas finish up before using alphas
401  Kokkos::Profiling::pushRegion("Apply Alphas to Data");
402 
403  //! [Apply GMLS Alphas To Data]
404 
405 
406  // it is important to note that if you expect to use the data as a 1D view, then you should use double*
407  // however, if you know that the target operation will result in a 2D view (vector or matrix output),
408  // then you should template with double** as this is something that can not be infered from the input data
409  // or the target operator at compile time. Additionally, a template argument is required indicating either
410  // Kokkos::HostSpace or Kokkos::DefaultExecutionSpace::memory_space()
411 
412  // The Evaluator class takes care of handling input data views as well as the output data views.
413  // It uses information from the GMLS class to determine how many components are in the input
414  // as well as output for any choice of target functionals and then performs the contactions
415  // on the data using the alpha coefficients generated by the GMLS class, all on the device.
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);
419 
420  auto output_value = scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
421  (sampling_data_device, ScalarPointEvaluation, PointSample,
422  true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
423 
424  auto output_gaussian_curvature = scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
425  (ones_data_device, GaussianCurvaturePointEvaluation, PointSample,
426  true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
427 
428  auto output_curl = scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
429  (sampling_data_device, CurlOfVectorPointEvaluation, PointSample,
430  true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
431 
432  //auto output_laplacian = scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
433  // (sampling_data_device, LaplacianOfScalarPointEvaluation, PointSample,
434  // true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
435 
436  //auto output_gradient = scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
437  // (sampling_data_device, GradientOfScalarPointEvaluation, PointSample,
438  // true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
439 
440  auto output_vector = vector_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
441  (sampling_vector_data_device, VectorPointEvaluation, VaryingManifoldVectorPointSample,
442  true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
443 
444  //auto output_divergence = vector_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
445  // (sampling_vector_data_device, DivergenceOfVectorPointEvaluation, ManifoldVectorPointSample,
446  // true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
447 
448  auto output_vector_of_scalar_clones =
449  vector_gmls_evaluator_of_scalar_clones.applyAlphasToDataAllComponentsAllTargetSites<double**,
450  Kokkos::HostSpace>(sampling_vector_data_device, VectorPointEvaluation, VaryingManifoldVectorPointSample,
451  true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
452 
453  //auto output_divergence_of_scalar_clones =
454  // vector_gmls_evaluator_of_scalar_clones.applyAlphasToDataAllComponentsAllTargetSites<double*,
455  // Kokkos::HostSpace>(sampling_vector_data_device, DivergenceOfVectorPointEvaluation, ManifoldVectorPointSample,
456  // true /*scalar_as_vector_if_needed*/, 1 /*evaluation site index*/);
457 
458 
459  // Kokkos::fence(); // let application of alphas to data finish before using results
460  //
461  //// move gradient data to device so that it can be transformed into velocity
462  //auto output_gradient_device_mirror = Kokkos::create_mirror(Kokkos::DefaultExecutionSpace::memory_space(), output_gradient);
463  //Kokkos::deep_copy(output_gradient_device_mirror, output_gradient);
464  //Kokkos::parallel_for("Create Velocity From Surface Gradient", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
465  // (0,target_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
466  //
467  // // coordinates of target site i
468  // double xval = target_coords_device(i,0);
469  // double yval = (dimension>1) ? target_coords_device(i,1) : 0;
470  // double zval = (dimension>2) ? target_coords_device(i,2) : 0;
471 
472  // double gradx = output_gradient_device_mirror(i,0);
473  // double grady = output_gradient_device_mirror(i,1);
474  // double gradz = output_gradient_device_mirror(i,2);
475  //
476  // // overwrites gradient with velocity
477  // output_gradient_device_mirror(i,0) = (grady*zval - yval*gradz);
478  // output_gradient_device_mirror(i,1) = (-(gradx*zval - xval*gradz));
479  // output_gradient_device_mirror(i,2) = (gradx*yval - xval*grady);
480  //
481  //});
482  //Kokkos::deep_copy(output_gradient, output_gradient_device_mirror);
483 
484 
485  //! [Apply GMLS Alphas To Data]
486 
487  Kokkos::fence(); // let application of alphas to data finish before using results
488  Kokkos::Profiling::popRegion();
489  // times the Comparison in Kokkos
490  Kokkos::Profiling::pushRegion("Comparison");
491 
492  //! [Check That Solutions Are Correct]
493 
494  double tangent_bundle_error = 0;
495  double tangent_bundle_norm = 0;
496  double values_error = 0;
497  double values_norm = 0;
498  double gc_error = 0;
499  double gc_norm = 0;
500  double curl_ambient_error = 0;
501  double curl_ambient_norm = 0;
502  //double laplacian_error = 0;
503  //double laplacian_norm = 0;
504  //double gradient_ambient_error = 0;
505  //double gradient_ambient_norm = 0;
506  double vector_ambient_error = 0;
507  double vector_ambient_norm = 0;
508  //double divergence_ambient_error = 0;
509  //double divergence_ambient_norm = 0;
510  double vector_of_scalar_clones_ambient_error = 0;
511  double vector_of_scalar_clones_ambient_norm = 0;
512  //double divergence_of_scalar_clones_ambient_error = 0;
513  //double divergence_of_scalar_clones_ambient_norm = 0;
514 
515  // tangent vectors for each source coordinate are stored here
516  auto d_prestencil_weights = my_GMLS_vector_of_scalar_clones.getPrestencilWeights();
517  auto prestencil_weights = Kokkos::create_mirror_view(d_prestencil_weights);
518  Kokkos::deep_copy(prestencil_weights, d_prestencil_weights);
519 
520  // tangent vector at target sites are stored here
521  auto d_tangent_directions = my_GMLS_vector_of_scalar_clones.getTangentDirections();
522  auto tangent_directions = Kokkos::create_mirror_view(d_tangent_directions);
523  Kokkos::deep_copy(tangent_directions, d_tangent_directions);
524 
525  // loop through the target sites
526  for (int i=0; i<number_target_coords; i++) {
527 
529  (tangent_directions.data() + TO_GLOBAL(i)*TO_GLOBAL(3)*TO_GLOBAL(3), 3, 3);
530  XYZ u;
531 
532 
533  // load value from output
534  double GMLS_value = output_value(i);
535  //printf("GMLS val: %f, %d\n", GMLS_value, i);
536 
537  double GMLS_gc = output_gaussian_curvature(i);
538  //printf("GMLS gc: %f, %d\n", GMLS_gc, i);
539 
540  // // load laplacian from output
541  // double GMLS_Laplacian = output_laplacian(i);
542 
543  // target site i's coordinate
544  //double xval = target_coords(i,0);
545  //double yval = (dimension>1) ? target_coords(i,1) : 0;
546  //double zval = (dimension>2) ? target_coords(i,2) : 0;
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};
551 
552  // get tangent vector and see if orthgonal to coordinate (it should be on a sphere)
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 /* local neighbor index */, j, k);
557  }
558  tangent_bundle_error += tangent_inner_prod * tangent_inner_prod;
559  }
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);
563  }
564  // inner product could be plus or minus 1 (depends on tangent direction ordering)
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;
568 
569  // evaluation of various exact solutions
570  double actual_value = sphere_harmonic54(xval, yval, zval);
571  double actual_gc = 1.0; // Gaussian curvature is constant
572  // double actual_Laplacian = laplace_beltrami_sphere_harmonic54(xval, yval, zval);
573  double actual_Gradient_ambient[3] = {0,0,0}; // initialized for 3, but only filled up to dimension
574  gradient_sphereHarmonic54_ambient(actual_Gradient_ambient, xval, yval, zval);
575  // //velocity_sphereHarmonic54_ambient(actual_Gradient_ambient, xval, yval, zval);
576  double actual_Curl_ambient[3] = {0,0,0};
577  curl_sphere_harmonic54(actual_Curl_ambient, xval, yval, zval);
578 
579  values_error += (GMLS_value - actual_value)*(GMLS_value - actual_value);
580  values_norm += actual_value*actual_value;
581 
582  gc_error += (GMLS_gc - actual_gc)*(GMLS_gc - actual_gc);
583  gc_norm += 1;
584 
585  for (int j=0; j<dimension; ++j) u[j] = output_curl(i,j);
586  AmbientLocalAmbient(u, T.data(), &prestencil_weights(0, i, 0 /* local neighbor index */, 0, 0));
587  for (int j=0; j<dimension; ++j) output_curl(i,j) = u[j];
588 
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];
592  }
593 
594  // laplacian_error += (GMLS_Laplacian - actual_Laplacian)*(GMLS_Laplacian - actual_Laplacian);
595  // laplacian_norm += actual_Laplacian*actual_Laplacian;
596 
597  // for (int j=0; j<dimension; ++j) {
598  // gradient_ambient_error += (output_gradient(i,j) - actual_Gradient_ambient[j])*(output_gradient(i,j) - actual_Gradient_ambient[j]);
599  // gradient_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
600  // }
601 
602  for (int j=0; j<dimension; ++j) u[j] = output_vector(i,j);
603  AmbientLocalAmbient(u, T.data(), &prestencil_weights(0, i, 0 /* local neighbor index */, 0, 0));
604  for (int j=0; j<dimension; ++j) output_vector(i,j) = u[j];
605 
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];
609  }
610 
611  // divergence_ambient_error += (output_divergence(i) - actual_Laplacian)*(output_divergence(i) - actual_Laplacian);
612  // divergence_ambient_norm += actual_Laplacian*actual_Laplacian;
613 
614 
615  for (int j=0; j<dimension; ++j) u[j] = output_vector_of_scalar_clones(i,j);
616  AmbientLocalAmbient(u, T.data(), &prestencil_weights(0, i, 0 /* local neighbor index */, 0, 0));
617  for (int j=0; j<dimension; ++j) output_vector_of_scalar_clones(i,j) = u[j];
618  //// first get back to local chart
619  //double local_vec[3] = {0,0};
620  //for (int j=0; j<dimension; ++j) {
621  // local_vec[0] += T(0,j) * output_vector_of_scalar_clones(i,j);
622  // local_vec[1] += T(1,j) * output_vector_of_scalar_clones(i,j);
623  //}
624  //// second go to ambient space using tangent for first neighbor
625  //for (int j=0; j<dimension; ++j) output_vector_of_scalar_clones(i,j) = 0;
626  //for (int j=0; j<dimension; ++j) {
627  // output_vector_of_scalar_clones(i,j) += prestencil_weights(0, i, 0 /* local neighbor index */, 0, j) * local_vec[0];
628  // output_vector_of_scalar_clones(i,j) += prestencil_weights(0, i, 0 /* local neighbor index */, 1, j) * local_vec[1];
629  //}
630 
631 
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];
635  }
636 
637  // divergence_of_scalar_clones_ambient_error += (output_divergence_of_scalar_clones(i) - actual_Laplacian)*(output_divergence_of_scalar_clones(i) - actual_Laplacian);
638  // divergence_of_scalar_clones_ambient_norm += actual_Laplacian*actual_Laplacian;
639 
640  }
641 
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);
646 
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);
651 
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);
656 
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);
661 
662  //laplacian_error /= number_target_coords;
663  //laplacian_error = std::sqrt(laplacian_error);
664  //laplacian_norm /= number_target_coords;
665  //laplacian_norm = std::sqrt(laplacian_norm);
666 
667  //gradient_ambient_error /= number_target_coords;
668  //gradient_ambient_error = std::sqrt(gradient_ambient_error);
669  //gradient_ambient_norm /= number_target_coords;
670  //gradient_ambient_norm = std::sqrt(gradient_ambient_norm);
671 
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);
676 
677  //divergence_ambient_error /= number_target_coords;
678  //divergence_ambient_error = std::sqrt(divergence_ambient_error);
679  //divergence_ambient_norm /= number_target_coords;
680  //divergence_ambient_norm = std::sqrt(divergence_ambient_norm);
681 
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);
686 
687  //divergence_of_scalar_clones_ambient_error /= number_target_coords;
688  //divergence_of_scalar_clones_ambient_error = std::sqrt(divergence_of_scalar_clones_ambient_error);
689  //divergence_of_scalar_clones_ambient_norm /= number_target_coords;
690  //divergence_of_scalar_clones_ambient_norm = std::sqrt(divergence_of_scalar_clones_ambient_norm);
691 
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);
696  //printf("Laplace-Beltrami Error: %g\n", laplacian_error / laplacian_norm);
697  //printf("Surface Gradient (Ambient) Error: %g\n", gradient_ambient_error / gradient_ambient_norm);
698  printf("Surface Vector (VectorBasis) Error: %g\n", vector_ambient_error / vector_ambient_norm);
699  //printf("Surface Divergence (VectorBasis) Error: %g\n", divergence_ambient_error / divergence_ambient_norm);
700  printf("Surface Vector (ScalarClones) Error: %g\n",
701  vector_of_scalar_clones_ambient_error / vector_of_scalar_clones_ambient_norm);
702  //printf("Surface Divergence (ScalarClones) Error: %g\n",
703  // divergence_of_scalar_clones_ambient_error / divergence_of_scalar_clones_ambient_norm);
704  //! [Check That Solutions Are Correct]
705  // popRegion hidden from tutorial
706  // stop timing comparison loop
707  Kokkos::Profiling::popRegion();
708  //! [Finalize Program]
709 
710 
711 } // end of code block to reduce scope, causing Kokkos View de-allocations
712 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
713 
714 // finalize Kokkos and MPI (if available)
715 Kokkos::finalize();
716 #ifdef COMPADRE_USE_MPI
717 MPI_Finalize();
718 #endif
719 
720 return 0;
721 
722 } // main
723 
724 
725 //! [Finalize Program]
KOKKOS_INLINE_FUNCTION double sphere_harmonic54(double x, double y, double z)
Point evaluation of the curl of a vector (results in a vector)
Lightweight Evaluator Helper This class is a lightweight wrapper for extracting and applying all rele...
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...
constexpr SamplingFunctional VaryingManifoldVectorPointSample
For integrating polynomial dotted with normal over an edge.
#define PI
KOKKOS_INLINE_FUNCTION void curl_sphere_harmonic54(double *curl, double x, double y, double z)
Point evaluation of Gaussian curvature.
void AmbientLocalAmbient(XYZ &u, double *T_data, double *P_data)
[Ambient to Local Back To Ambient Helper Function]
#define TO_GLOBAL(variable)
Point evaluation of a scalar.
Kokkos::View< double **, layout_right, host_execution_space, Kokkos::MemoryTraits< Kokkos::Unmanaged > > host_scratch_matrix_right_type
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)
int main(int argc, char *args[])
[Ambient to Local Back To Ambient Helper Function]
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) ...
constexpr SamplingFunctional PointSample
Available sampling functionals.
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...