Compadre  1.3.3
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 #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  { // fill source coordinates from surface of a sphere with quasiuniform points
73  // algorithm described at https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf
74  int N_count = 0;
75  double a = 4*PI*r*r/N_pts_on_sphere;
76  double d = std::sqrt(a);
77  int M_theta = std::round(PI/d);
78  double d_theta = PI/M_theta;
79  double d_phi = a/d_theta;
80  for (int i=0; i<M_theta; ++i) {
81  double theta = PI*(i + 0.5)/M_theta;
82  int M_phi = std::round(2*PI*std::sin(theta)/d_phi);
83  for (int j=0; j<M_phi; ++j) {
84  double phi = 2*PI*j/M_phi;
85  source_coords(N_count, 0) = theta;
86  source_coords(N_count, 1) = phi;
87  N_count++;
88  }
89  }
90 
91  for (int i=0; i<N_count; ++i) {
92  double theta = source_coords(i,0);
93  double phi = source_coords(i,1);
94  source_coords(i,0) = r*std::sin(theta)*std::cos(phi);
95  source_coords(i,1) = r*std::sin(theta)*std::sin(phi);
96  source_coords(i,2) = r*cos(theta);
97  //printf("%f %f %f\n", source_coords(i,0), source_coords(i,1), source_coords(i,2));
98  }
99  number_source_coords = N_count;
100  }
101 
102  // resize source_coords to the size actually needed
103  Kokkos::resize(source_coords, number_source_coords, 3);
104  Kokkos::resize(source_coords_device, number_source_coords, 3);
105 
106  // coordinates of target sites
107  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device("target coordinates",
108  number_target_coords, 3);
109  Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
110 
111  // seed random number generator
112  std::mt19937 rng(50);
113 
114  // generate random integers in [0...number_source_coords-1] (used to pick target sites)
115  std::uniform_int_distribution<int> gen_num_neighbors(0, number_source_coords-1); // uniform, unbiased
116 
117  // fill target sites with random selections from source sites
118  for (int i=0; i<number_target_coords; ++i) {
119  const int source_site_to_copy = gen_num_neighbors(rng);
120  for (int j=0; j<3; ++j) {
121  target_coords(i,j) = source_coords(source_site_to_copy,j);
122  }
123  }
124 
125 
126  //! [Setting Up The Point Cloud]
127 
128  Kokkos::Profiling::popRegion();
129  Kokkos::fence();
130  Kokkos::Profiling::pushRegion("Creating Data");
131 
132  //! [Creating The Data]
133 
134 
135  // source coordinates need copied to device before using to construct sampling data
136  Kokkos::deep_copy(source_coords_device, source_coords);
137  Kokkos::deep_copy(target_coords_device, target_coords);
138 
139  // ensure that source coordinates are sent to device before evaluating sampling data based on them
140  Kokkos::fence();
141 
142 
143  // need Kokkos View storing true solution (for samples)
144  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device("samples of true solution",
145  source_coords_device.extent(0));
146 
147  // need Kokkos View storing true vector solution (for samples)
148  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> sampling_vector_data_device("samples of vector true solution",
149  source_coords_device.extent(0), 3);
150 
151  Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
152  (0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
153 
154  // coordinates of source site i
155  double xval = source_coords_device(i,0);
156  double yval = (dimension>1) ? source_coords_device(i,1) : 0;
157  double zval = (dimension>2) ? source_coords_device(i,2) : 0;
158 
159  // data for targets with scalar input
160  sampling_data_device(i) = sphere_harmonic54(xval, yval, zval);
161  //printf("%f\n", sampling_data_device(i));
162 
163  for (int j=0; j<3; ++j) {
164  double gradient[3] = {0,0,0};
165  gradient_sphereHarmonic54_ambient(gradient, xval, yval, zval);
166  sampling_vector_data_device(i,j) = gradient[j];
167  }
168  //printf("%f %f %f\n", sampling_vector_data_device(i,0), sampling_vector_data_device(i,1), sampling_vector_data_device(i,2));
169  });
170 
171 
172  //! [Creating The Data]
173 
174  Kokkos::Profiling::popRegion();
175  Kokkos::Profiling::pushRegion("Neighbor Search");
176 
177  //! [Performing Neighbor Search]
178 
179 
180  // Point cloud construction for neighbor search
181  // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
182  auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
183 
184  // each row is a neighbor list for a target site, with the first column of each row containing
185  // the number of neighbors for that rows corresponding target site
186  double epsilon_multiplier = 1.5;
187  int estimated_upper_bound_number_neighbors =
188  point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
189 
190  Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device("neighbor lists",
191  number_target_coords, estimated_upper_bound_number_neighbors); // first column is # of neighbors
192  Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
193 
194  // each target site has a window size
195  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device("h supports", number_target_coords);
196  Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
197 
198  // query the point cloud to generate the neighbor lists using a kdtree to produce the n nearest neighbor
199  // to each target site, adding (epsilon_multiplier-1)*100% to whatever the distance away the further neighbor used is from
200  // each target to the view for epsilon
201  point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/, target_coords, neighbor_lists,
202  epsilon, min_neighbors, epsilon_multiplier);
203 
204 
205  //! [Performing Neighbor Search]
206 
207  Kokkos::Profiling::popRegion();
208  Kokkos::fence(); // let call to build neighbor lists complete before copying back to device
209  timer.reset();
210 
211  //! [Setting Up The GMLS Object]
212 
213 
214  // Copy data back to device (they were filled on the host)
215  // We could have filled Kokkos Views with memory space on the host
216  // and used these instead, and then the copying of data to the device
217  // would be performed in the GMLS class
218  Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
219  Kokkos::deep_copy(epsilon_device, epsilon);
220 
221  // initialize an instance of the GMLS class
224  order, dimension,
225  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
226  order /*manifold order*/);
227 
228  // pass in neighbor lists, source coordinates, target coordinates, and window sizes
229  //
230  // neighbor lists have the format:
231  // dimensions: (# number of target sites) X (# maximum number of neighbors for any given target + 1)
232  // the first column contains the number of neighbors for that rows corresponding target index
233  //
234  // source coordinates have the format:
235  // dimensions: (# number of source sites) X (dimension)
236  // entries in the neighbor lists (integers) correspond to rows of this 2D array
237  //
238  // target coordinates have the format:
239  // dimensions: (# number of target sites) X (dimension)
240  // # of target sites is same as # of rows of neighbor lists
241  //
242  my_GMLS_vector_1.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
243 
244  // create a vector of target operations
245  std::vector<TargetOperation> lro_vector_1(1);
246  lro_vector_1[0] = DivergenceOfVectorPointEvaluation;
247 
248  // and then pass them to the GMLS class
249  my_GMLS_vector_1.addTargets(lro_vector_1);
250 
251  // sets the weighting kernel function from WeightingFunctionType for curvature
252  my_GMLS_vector_1.setCurvatureWeightingType(WeightingFunctionType::Power);
253 
254  // power to use in the weighting kernel function for curvature coefficients
255  my_GMLS_vector_1.setCurvatureWeightingPower(2);
256 
257  // sets the weighting kernel function from WeightingFunctionType
258  my_GMLS_vector_1.setWeightingType(WeightingFunctionType::Power);
259 
260  // power to use in that weighting kernel function
261  my_GMLS_vector_1.setWeightingPower(2);
262 
263  // setup quadrature for StaggeredEdgeIntegralSample
264  my_GMLS_vector_1.setOrderOfQuadraturePoints(2);
265  my_GMLS_vector_1.setDimensionOfQuadraturePoints(1);
266  my_GMLS_vector_1.setQuadratureType("LINE");
267 
268  // generate the alphas that to be combined with data for each target operation requested in lro
269  my_GMLS_vector_1.generateAlphas();
270 
271  // initialize another instance of the GMLS class
275  order, dimension,
276  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
277  order /*manifold order*/);
278 
279  my_GMLS_vector_2.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
280  std::vector<TargetOperation> lro_vector_2(2);
282  lro_vector_2[1] = DivergenceOfVectorPointEvaluation;
283  //lro_vector_2[2] = GradientOfScalarPointEvaluation;
284  my_GMLS_vector_2.addTargets(lro_vector_2);
285  my_GMLS_vector_2.setCurvatureWeightingType(WeightingFunctionType::Power);
286  my_GMLS_vector_2.setCurvatureWeightingPower(2);
287  my_GMLS_vector_2.setWeightingType(WeightingFunctionType::Power);
288  my_GMLS_vector_2.setWeightingPower(2);
289  my_GMLS_vector_2.setOrderOfQuadraturePoints(2);
290  my_GMLS_vector_2.setDimensionOfQuadraturePoints(1);
291  my_GMLS_vector_2.setQuadratureType("LINE");
292  my_GMLS_vector_2.generateAlphas();
293 
294  // initialize another instance of the GMLS class
297  order, dimension,
298  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
299  order /*manifold order*/);
300 
301  my_GMLS_scalar.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
302 
303  std::vector<TargetOperation> lro_scalar(1);
305  //lro_scalar[1] = GradientOfScalarPointEvaluation;
306  my_GMLS_scalar.addTargets(lro_scalar);
307  my_GMLS_scalar.setCurvatureWeightingType(WeightingFunctionType::Power);
308  my_GMLS_scalar.setCurvatureWeightingPower(2);
309  my_GMLS_scalar.setWeightingType(WeightingFunctionType::Power);
310  my_GMLS_scalar.setWeightingPower(2);
311  my_GMLS_scalar.generateAlphas();
312 
313 
314  //! [Setting Up The GMLS Object]
315 
316  double instantiation_time = timer.seconds();
317  std::cout << "Took " << instantiation_time << "s to complete alphas generation." << std::endl;
318  Kokkos::fence(); // let generateAlphas finish up before using alphas
319  Kokkos::Profiling::pushRegion("Apply Alphas to Data");
320 
321  //! [Apply GMLS Alphas To Data]
322 
323 
324  // it is important to note that if you expect to use the data as a 1D view, then you should use double*
325  // however, if you know that the target operation will result in a 2D view (vector or matrix output),
326  // then you should template with double** as this is something that can not be infered from the input data
327  // or the target operator at compile time. Additionally, a template argument is required indicating either
328  // Kokkos::HostSpace or Kokkos::DefaultExecutionSpace::memory_space()
329 
330  // The Evaluator class takes care of handling input data views as well as the output data views.
331  // It uses information from the GMLS class to determine how many components are in the input
332  // as well as output for any choice of target functionals and then performs the contactions
333  // on the data using the alpha coefficients generated by the GMLS class, all on the device.
334  Evaluator vector_1_gmls_evaluator(&my_GMLS_vector_1);
335  Evaluator vector_2_gmls_evaluator(&my_GMLS_vector_2);
336  Evaluator scalar_gmls_evaluator(&my_GMLS_scalar);
337 
338 
339  //auto output_gradient_vectorbasis =
340  // vector_2_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
341  // (sampling_data_device, GradientOfScalarPointEvaluation, StaggeredEdgeAnalyticGradientIntegralSample);
342 
343  //auto output_gradient_scalarbasis =
344  // scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
345  // (sampling_data_device, GradientOfScalarPointEvaluation, StaggeredEdgeAnalyticGradientIntegralSample);
346 
347  auto output_divergence_vectorsamples =
348  vector_1_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
349  (sampling_vector_data_device, DivergenceOfVectorPointEvaluation, StaggeredEdgeIntegralSample);
350 
351  auto output_divergence_scalarsamples =
352  vector_2_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
354 
355  auto output_laplacian_vectorbasis =
356  vector_2_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
358 
359  auto output_laplacian_scalarbasis =
360  scalar_gmls_evaluator.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
362 
363 
364  //! [Apply GMLS Alphas To Data]
365 
366  Kokkos::fence(); // let application of alphas to data finish before using results
367  Kokkos::Profiling::popRegion();
368  // times the Comparison in Kokkos
369  Kokkos::Profiling::pushRegion("Comparison");
370 
371  //! [Check That Solutions Are Correct]
372 
373  double laplacian_vectorbasis_error = 0;
374  double laplacian_vectorbasis_norm = 0;
375 
376  double laplacian_scalarbasis_error = 0;
377  double laplacian_scalarbasis_norm = 0;
378 
379  double gradient_vectorbasis_ambient_error = 0;
380  double gradient_vectorbasis_ambient_norm = 0;
381 
382  double gradient_scalarbasis_ambient_error = 0;
383  double gradient_scalarbasis_ambient_norm = 0;
384 
385  double divergence_vectorsamples_ambient_error = 0;
386  double divergence_vectorsamples_ambient_norm = 0;
387 
388  double divergence_scalarsamples_ambient_error = 0;
389  double divergence_scalarsamples_ambient_norm = 0;
390 
391  // loop through the target sites
392  for (int i=0; i<number_target_coords; i++) {
393 
394  // target site i's coordinate
395  double xval = target_coords(i,0);
396  double yval = (dimension>1) ? target_coords(i,1) : 0;
397  double zval = (dimension>2) ? target_coords(i,2) : 0;
398 
399  // evaluation of various exact solutions
400  double actual_Laplacian = laplace_beltrami_sphere_harmonic54(xval, yval, zval);
401  double actual_Gradient_ambient[3] = {0,0,0}; // initialized for 3, but only filled up to dimension
402  gradient_sphereHarmonic54_ambient(actual_Gradient_ambient, xval, yval, zval);
403 
404  laplacian_vectorbasis_error += (output_laplacian_vectorbasis(i) - actual_Laplacian)*(output_laplacian_vectorbasis(i) - actual_Laplacian);
405  laplacian_vectorbasis_norm += actual_Laplacian*actual_Laplacian;
406 
407  //printf("Error of %f, %f vs %f\n", (output_laplacian_scalarbasis(i) - actual_Laplacian), output_laplacian_scalarbasis(i), actual_Laplacian);
408  laplacian_scalarbasis_error += (output_laplacian_scalarbasis(i) - actual_Laplacian)*(output_laplacian_scalarbasis(i) - actual_Laplacian);
409  laplacian_scalarbasis_norm += actual_Laplacian*actual_Laplacian;
410 
411  //for (int j=0; j<dimension; ++j) {
412  // //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]);
413  // gradient_vectorbasis_ambient_error += (output_gradient_vectorbasis(i,j) - actual_Gradient_ambient[j])*(output_gradient_vectorbasis(i,j) - actual_Gradient_ambient[j]);
414  // gradient_vectorbasis_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
415  //}
416 
417  //for (int j=0; j<dimension; ++j) {
418  // //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]);
419  // gradient_scalarbasis_ambient_error += (output_gradient_scalarbasis(i,j) - actual_Gradient_ambient[j])*(output_gradient_scalarbasis(i,j) - actual_Gradient_ambient[j]);
420  // gradient_scalarbasis_ambient_norm += actual_Gradient_ambient[j]*actual_Gradient_ambient[j];
421  //}
422 
423  //printf("Error of %f, %f vs %f\n", (output_divergence(i) - actual_Laplacian), output_divergence(i), actual_Laplacian);
424  divergence_vectorsamples_ambient_error += (output_divergence_vectorsamples(i) - actual_Laplacian)*(output_divergence_vectorsamples(i) - actual_Laplacian);
425  divergence_vectorsamples_ambient_norm += actual_Laplacian*actual_Laplacian;
426 
427  divergence_scalarsamples_ambient_error += (output_divergence_scalarsamples(i) - actual_Laplacian)*(output_divergence_scalarsamples(i) - actual_Laplacian);
428  divergence_scalarsamples_ambient_norm += actual_Laplacian*actual_Laplacian;
429 
430  }
431 
432  laplacian_vectorbasis_error /= number_target_coords;
433  laplacian_vectorbasis_error = std::sqrt(laplacian_vectorbasis_error);
434  laplacian_vectorbasis_norm /= number_target_coords;
435  laplacian_vectorbasis_norm = std::sqrt(laplacian_vectorbasis_norm);
436 
437  laplacian_scalarbasis_error /= number_target_coords;
438  laplacian_scalarbasis_error = std::sqrt(laplacian_scalarbasis_error);
439  laplacian_scalarbasis_norm /= number_target_coords;
440  laplacian_scalarbasis_norm = std::sqrt(laplacian_scalarbasis_norm);
441 
442  gradient_vectorbasis_ambient_error /= number_target_coords;
443  gradient_vectorbasis_ambient_error = std::sqrt(gradient_vectorbasis_ambient_error);
444  gradient_vectorbasis_ambient_norm /= number_target_coords;
445  gradient_vectorbasis_ambient_norm = std::sqrt(gradient_vectorbasis_ambient_norm);
446 
447  gradient_scalarbasis_ambient_error /= number_target_coords;
448  gradient_scalarbasis_ambient_error = std::sqrt(gradient_scalarbasis_ambient_error);
449  gradient_scalarbasis_ambient_norm /= number_target_coords;
450  gradient_scalarbasis_ambient_norm = std::sqrt(gradient_scalarbasis_ambient_norm);
451 
452  divergence_vectorsamples_ambient_error /= number_target_coords;
453  divergence_vectorsamples_ambient_error = std::sqrt(divergence_vectorsamples_ambient_error);
454  divergence_vectorsamples_ambient_norm /= number_target_coords;
455  divergence_vectorsamples_ambient_norm = std::sqrt(divergence_vectorsamples_ambient_norm);
456 
457  divergence_scalarsamples_ambient_error /= number_target_coords;
458  divergence_scalarsamples_ambient_error = std::sqrt(divergence_scalarsamples_ambient_error);
459  divergence_scalarsamples_ambient_norm /= number_target_coords;
460  divergence_scalarsamples_ambient_norm = std::sqrt(divergence_scalarsamples_ambient_norm);
461 
462  printf("Staggered Laplace-Beltrami (VectorBasis) Error: %g\n", laplacian_vectorbasis_error / laplacian_vectorbasis_norm);
463  printf("Staggered Laplace-Beltrami (ScalarBasis) Error: %g\n", laplacian_scalarbasis_error / laplacian_scalarbasis_norm);
464  printf("Surface Staggered Gradient (VectorBasis) Error: %g\n", gradient_vectorbasis_ambient_error / gradient_vectorbasis_ambient_norm);
465  printf("Surface Staggered Gradient (ScalarBasis) Error: %g\n", gradient_scalarbasis_ambient_error / gradient_scalarbasis_ambient_norm);
466  printf("Surface Staggered Divergence (VectorSamples) Error: %g\n", divergence_vectorsamples_ambient_error / divergence_vectorsamples_ambient_norm);
467  printf("Surface Staggered Divergence (ScalarSamples) Error: %g\n", divergence_scalarsamples_ambient_error / divergence_scalarsamples_ambient_norm);
468  //! [Check That Solutions Are Correct]
469  // popRegion hidden from tutorial
470  // stop timing comparison loop
471  Kokkos::Profiling::popRegion();
472  //! [Finalize Program]
473 
474 
475 } // end of code block to reduce scope, causing Kokkos View de-allocations
476 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
477 
478 // finalize Kokkos and MPI (if available)
479 Kokkos::finalize();
480 #ifdef COMPADRE_USE_MPI
481 MPI_Finalize();
482 #endif
483 
484 return 0;
485 
486 } // main
487 
488 
489 //! [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 (allocates memory for output)
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...