Compadre  1.3.3
GMLS_Staggered.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_Tutorial.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 // becomes false if the computed solution not within the failure_threshold of the actual solution
40 bool all_passed = true;
41 
42 // code block to reduce scope for all Kokkos View allocations
43 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
44 {
45 
46  CommandLineProcessor clp(argc, args);
47  auto order = clp.order;
48  auto dimension = clp.dimension;
49  auto number_target_coords = clp.number_target_coords;
50  auto constraint_name = clp.constraint_name;
51  auto solver_name = clp.solver_name;
52  auto problem_name = clp.problem_name;
53 
54  // the functions we will be seeking to reconstruct are in the span of the basis
55  // of the reconstruction space we choose for GMLS, so the error should be very small
56  const double failure_tolerance = 9e-8;
57 
58  // minimum neighbors for unisolvency is the same as the size of the polynomial basis
59  const int min_neighbors = Compadre::GMLS::getNP(order+1, dimension);
60 
61  //! [Parse Command Line Arguments]
62  Kokkos::Timer timer;
63  Kokkos::Profiling::pushRegion("Setup Point Data");
64  //! [Setting Up The Point Cloud]
65 
66  // approximate spacing of source sites
67  double h_spacing = 0.05;
68  int n_neg1_to_1 = 2*(1/h_spacing) + 1; // always odd
69 
70  // number of source coordinate sites that will fill a box of [-1,1]x[-1,1]x[-1,1] with a spacing approximately h
71  const int number_source_coords = std::pow(n_neg1_to_1, dimension);
72 
73  // coordinates of source sites
74  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> source_coords_device("source coordinates",
75  number_source_coords, 3);
76  Kokkos::View<double**>::HostMirror source_coords = Kokkos::create_mirror_view(source_coords_device);
77 
78  // coordinates of target sites
79  Kokkos::View<double**, Kokkos::DefaultExecutionSpace> target_coords_device ("target coordinates", number_target_coords, 3);
80  Kokkos::View<double**>::HostMirror target_coords = Kokkos::create_mirror_view(target_coords_device);
81 
82  // tangent bundle for each target sites
83  Kokkos::View<double***, Kokkos::DefaultExecutionSpace> tangent_bundles_device ("tangent bundles", number_target_coords, 3, 3);
84  Kokkos::View<double***>::HostMirror tangent_bundles = Kokkos::create_mirror_view(tangent_bundles_device);
85 
86  // fill source coordinates with a uniform grid
87  int source_index = 0;
88  double this_coord[3] = {0,0,0};
89  for (int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
90  this_coord[0] = i*h_spacing;
91  for (int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
92  this_coord[1] = j*h_spacing;
93  for (int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
94  this_coord[2] = k*h_spacing;
95  if (dimension==3) {
96  source_coords(source_index,0) = this_coord[0];
97  source_coords(source_index,1) = this_coord[1];
98  source_coords(source_index,2) = this_coord[2];
99  source_index++;
100  }
101  }
102  if (dimension==2) {
103  source_coords(source_index,0) = this_coord[0];
104  source_coords(source_index,1) = this_coord[1];
105  source_coords(source_index,2) = 0;
106  source_index++;
107  }
108  }
109  if (dimension==1) {
110  source_coords(source_index,0) = this_coord[0];
111  source_coords(source_index,1) = 0;
112  source_coords(source_index,2) = 0;
113  source_index++;
114  }
115  }
116 
117  // Generate target points - these are random permutation from available source points
118  // Note that this is assuming that the number of targets in this test will not exceed
119  // the number of source coords, which is 41^3 = 68921
120  // seed random number generator
121  std::mt19937 rng(50);
122  // generate random integers in [0..number_source_coords-1] (used to pick target sites)
123  std::uniform_int_distribution<int> gen_num_neighbours(0, source_index);
124  // fill target sites with random selections from source sites
125  for (int i=0; i<number_target_coords; ++i) {
126  const int source_site_to_copy = gen_num_neighbours(rng);
127  for (int j=0; j<3; j++) {
128  target_coords(i, j) = source_coords(source_site_to_copy, j);
129  }
130  if (constraint_name == "NEUMANN_GRAD_SCALAR") { // create bundles of normal vectors
131  if (dimension == 3) {
132  tangent_bundles(i, 0, 0) = 0.0;
133  tangent_bundles(i, 0, 1) = 0.0;
134  tangent_bundles(i, 0, 2) = 0.0;
135  tangent_bundles(i, 1, 0) = 0.0;
136  tangent_bundles(i, 1, 1) = 0.0;
137  tangent_bundles(i, 1, 2) = 0.0;
138  tangent_bundles(i, 2, 0) = 1.0/(sqrt(3.0));
139  tangent_bundles(i, 2, 1) = 1.0/(sqrt(3.0));
140  tangent_bundles(i, 2, 2) = 1.0/(sqrt(3.0));
141  }
142  if (dimension == 2) {
143  tangent_bundles(i, 0, 0) = 0.0;
144  tangent_bundles(i, 0, 1) = 0.0;
145  tangent_bundles(i, 1, 0) = 1.0/(sqrt(2.0));
146  tangent_bundles(i, 1, 1) = 1.0/(sqrt(2.0));
147  }
148  }
149  }
150 
151  //! [Setting Up The Point Cloud]
152 
153  Kokkos::Profiling::popRegion();
154  Kokkos::Profiling::pushRegion("Creating Data");
155 
156  //! [Creating The Data]
157 
158 
159  // source coordinates need copied to device before using to construct sampling data
160  Kokkos::deep_copy(source_coords_device, source_coords);
161 
162  // target coordinates copied next, because it is a convenient time to send them to device
163  Kokkos::deep_copy(target_coords_device, target_coords);
164 
165  // need Kokkos View storing true solution
166  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> sampling_data_device("samples of true solution",
167  source_coords_device.extent(0));
168 
169  Kokkos::parallel_for("Sampling Manufactured Solutions", Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>
170  (0,source_coords.extent(0)), KOKKOS_LAMBDA(const int i) {
171  // coordinates of source site i
172  double xval = source_coords_device(i,0);
173  double yval = (dimension>1) ? source_coords_device(i,1) : 0;
174  double zval = (dimension>2) ? source_coords_device(i,2) : 0;
175 
176  // data for targets with scalar input
177  sampling_data_device(i) = trueSolution(xval, yval, zval, order, dimension);
178  });
179 
180 
181  //! [Creating The Data]
182 
183  Kokkos::Profiling::popRegion();
184  Kokkos::Profiling::pushRegion("Neighbor Search");
185 
186  //! [Performing Neighbor Search]
187 
188 
189  // Point cloud construction for neighbor search
190  // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
191  auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
192 
193  // each row is a neighbor list for a target site, with the first column of each row containing
194  // the number of neighbors for that rows corresponding target site
195  // for the default values in this test, the multiplier is suggested to be 2.2
196  double epsilon_multiplier = 2.2;
197  int estimated_upper_bound_number_neighbors =
198  point_cloud_search.getEstimatedNumberNeighborsUpperBound(min_neighbors, dimension, epsilon_multiplier);
199 
200  Kokkos::View<int**, Kokkos::DefaultExecutionSpace> neighbor_lists_device("neighbor lists",
201  number_target_coords, estimated_upper_bound_number_neighbors); // first column is # of neighbors
202  Kokkos::View<int**>::HostMirror neighbor_lists = Kokkos::create_mirror_view(neighbor_lists_device);
203 
204  // each target site has a window size
205  Kokkos::View<double*, Kokkos::DefaultExecutionSpace> epsilon_device("h supports", number_target_coords);
206  Kokkos::View<double*>::HostMirror epsilon = Kokkos::create_mirror_view(epsilon_device);
207 
208  // query the point cloud to generate the neighbor lists using a kdtree to produce the n nearest neighbor
209  // to each target site, adding (epsilon_multiplier-1)*100% to whatever the distance away the further neighbor used is from
210  // each target to the view for epsilon
211  point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/, target_coords, neighbor_lists,
212  epsilon, min_neighbors, epsilon_multiplier);
213 
214  //! [Performing Neighbor Search]
215 
216  Kokkos::Profiling::popRegion();
217  Kokkos::fence(); // let call to build neighbor lists complete before copying back to device
218  timer.reset();
219 
220  //! [Setting Up The GMLS Object]
221 
222 
223  // Copy data back to device (they were filled on the host)
224  // We could have filled Kokkos Views with memory space on the host
225  // and used these instead, and then the copying of data to the device
226  // would be performed in the GMLS class
227  Kokkos::deep_copy(neighbor_lists_device, neighbor_lists);
228  Kokkos::deep_copy(epsilon_device, epsilon);
229 
230  // initialize an instance of the GMLS class
231  // NULL in manifold order indicates non-manifold case
232  // First, analytica gradient on scalar polynomial basis
233  GMLS scalar_basis_gmls(ScalarTaylorPolynomial,
235  order, dimension,
236  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
237  0 /*manifold order*/);
238 
239  // Another class performing Gaussian quadrature integration on vector polynomial basis
240  GMLS vector_basis_gmls(VectorTaylorPolynomial,
243  order, dimension,
244  solver_name.c_str(), problem_name.c_str(), constraint_name.c_str(),
245  0 /*manifold order*/);
246 
247  // pass in neighbor lists, source coordinates, target coordinates, and window sizes
248  //
249  // neighbor lists have the format:
250  // dimensions: (# number of target sites) X (# maximum number of neighbors for any given target + 1)
251  // the first column contains the number of neighbors for that rows corresponding target index
252  //
253  // source coordinates have the format:
254  // dimensions: (# number of source sites) X (dimension)
255  // entries in the neighbor lists (integers) correspond to rows of this 2D array
256  //
257  // target coordinates have the format:
258  // dimensions: (# number of target sites) X (dimension)
259  // # of target sites is same as # of rows of neighbor lists
260  //
261  scalar_basis_gmls.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
262  vector_basis_gmls.setProblemData(neighbor_lists_device, source_coords_device, target_coords_device, epsilon_device);
263  if (constraint_name == "NEUMANN_GRAD_SCALAR") {
264  scalar_basis_gmls.setTangentBundle(tangent_bundles_device);
265  vector_basis_gmls.setTangentBundle(tangent_bundles_device);
266  }
267 
268  // create a vector of target operations
269  std::vector<TargetOperation> lro(2);
272 
273  // and then pass them to the GMLS class
274  scalar_basis_gmls.addTargets(lro);
275  vector_basis_gmls.addTargets(lro);
276 
277  // sets the weighting kernel function from WeightingFunctionType
278  scalar_basis_gmls.setWeightingType(WeightingFunctionType::Power);
279  vector_basis_gmls.setWeightingType(WeightingFunctionType::Power);
280 
281  // power to use in that weighting kernel function
282  scalar_basis_gmls.setWeightingPower(2);
283  vector_basis_gmls.setWeightingPower(2);
284 
285  // setup quadrature for StaggeredEdgeIntegralSample
286  vector_basis_gmls.setOrderOfQuadraturePoints(order);
287  vector_basis_gmls.setDimensionOfQuadraturePoints(1);
288  vector_basis_gmls.setQuadratureType("LINE");
289 
290  // generate the alphas that to be combined with data for each target operation requested in lro
291  scalar_basis_gmls.generateAlphas();
292  vector_basis_gmls.generateAlphas();
293 
294  //! [Setting Up The GMLS Object]
295 
296  double instantiation_time = timer.seconds();
297  std::cout << "Took " << instantiation_time << "s to complete alphas generation." << std::endl;
298  Kokkos::fence(); // let generateAlphas finish up before using alphas
299  Kokkos::Profiling::pushRegion("Apply Alphas to Data");
300 
301  //! [Apply GMLS Alphas To Data]
302 
303  // it is important to note that if you expect to use the data as a 1D view, then you should use double*
304  // however, if you know that the target operation will result in a 2D view (vector or matrix output),
305  // then you should template with double** as this is something that can not be infered from the input data
306  // or the target operator at compile time. Additionally, a template argument is required indicating either
307  // Kokkos::HostSpace or Kokkos::DefaultExecutionSpace::memory_space()
308 
309  // The Evaluator class takes care of handling input data views as well as the output data views.
310  // It uses information from the GMLS class to determine how many components are in the input
311  // as well as output for any choice of target functionals and then performs the contactions
312  // on the data using the alpha coefficients generated by the GMLS class, all on the device.
313  Evaluator gmls_evaluator_scalar(&scalar_basis_gmls);
314  Evaluator gmls_evaluator_vector(&vector_basis_gmls);
315 
316  auto output_divergence_scalar = gmls_evaluator_scalar.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
318 
319  auto output_divergence_vector = gmls_evaluator_vector.applyAlphasToDataAllComponentsAllTargetSites<double*, Kokkos::HostSpace>
321 
322  // evaluating the gradient
323  auto output_gradient_scalar = gmls_evaluator_scalar.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
325 
326  auto output_gradient_vector = gmls_evaluator_vector.applyAlphasToDataAllComponentsAllTargetSites<double**, Kokkos::HostSpace>
328 
329  //! [Apply GMLS Alphas To Data]
330 
331  Kokkos::fence(); // let application of alphas to data finish before using results
332  Kokkos::Profiling::popRegion();
333  // times the Comparison in Kokkos
334  Kokkos::Profiling::pushRegion("Comparison");
335 
336  //! [Check That Solutions Are Correct]
337 
338  // loop through the target sites
339  for (int i=0; i<number_target_coords; i++) {
340 
341  // target site i's coordinate
342  double xval = target_coords(i,0);
343  double yval = (dimension>1) ? target_coords(i,1) : 0;
344  double zval = (dimension>2) ? target_coords(i,2) : 0;
345 
346  // evaluation of various exact solutions
347  double actual_Divergence;
348  actual_Divergence = trueLaplacian(xval, yval, zval, order, dimension);
349  double actual_Gradient[3] = {0,0,0}; // initialized for 3, but only filled up to dimension
350  trueGradient(actual_Gradient, xval, yval, zval, order, dimension);
351 
352  // calculate correction factor
353  double g = (dimension == 3) ? (1.0/sqrt(3.0))*(actual_Gradient[0] + actual_Gradient[1] + actual_Gradient[2])
354  : (1.0/sqrt(2.0))*(actual_Gradient[0] + actual_Gradient[1]);
355  // obtain number of neighbor for each target
356  // in order to exploit the index where the value for the Lagrange multiplier is stored
357  int num_neigh_i = neighbor_lists(i, 0);
358 
359  // load divergence from output
360  double GMLS_Divergence_Scalar = output_divergence_scalar(i);
361  double GMLS_Divergence_Vector = output_divergence_vector(i);
362 
363  // obtain adjusted value for divergence
364  if (constraint_name == "NEUMANN_GRAD_SCALAR") {
365  double b_i_scalar = scalar_basis_gmls.getAlpha0TensorTo0Tensor(DivergenceOfVectorPointEvaluation, i, num_neigh_i);
366  GMLS_Divergence_Scalar = GMLS_Divergence_Scalar + b_i_scalar*g;
367 
368  double b_i_vector = vector_basis_gmls.getAlpha0TensorTo0Tensor(DivergenceOfVectorPointEvaluation, i, num_neigh_i);
369  GMLS_Divergence_Vector = GMLS_Divergence_Vector + b_i_vector*g;
370  }
371 
372  // check divergence - scalar basis
373  if(std::abs(actual_Divergence - GMLS_Divergence_Scalar) > failure_tolerance) {
374  all_passed = false;
375  std::cout << i << " Failed Divergence on SCALAR basis by: " << std::abs(actual_Divergence - GMLS_Divergence_Scalar) << std::endl;
376  std::cout << i << " GMLS " << GMLS_Divergence_Scalar << " adjusted " << GMLS_Divergence_Scalar << " actual " << actual_Divergence << std::endl;
377  }
378 
379  // check divergence - vector basis
380  if(std::abs(actual_Divergence - GMLS_Divergence_Vector) > failure_tolerance) {
381  all_passed = false;
382  std::cout << i << " Failed Divergence on VECTOR basis by: " << std::abs(actual_Divergence - GMLS_Divergence_Vector) << std::endl;
383  std::cout << i << " GMLS " << GMLS_Divergence_Vector << " adjusted " << GMLS_Divergence_Vector << " actual " << actual_Divergence << std::endl;
384  }
385 
386  // load gradient from output
387  double Scalar_GMLS_GradX = (dimension>1) ? output_gradient_scalar(i,0) : 0;
388  double Scalar_GMLS_GradY = (dimension>1) ? output_gradient_scalar(i,1) : 0;
389  double Scalar_GMLS_GradZ = (dimension>2) ? output_gradient_scalar(i,2) : 0;
390 
391  double Vector_GMLS_GradX = (dimension>1) ? output_gradient_vector(i,0) : 0;
392  double Vector_GMLS_GradY = (dimension>1) ? output_gradient_vector(i,1) : 0;
393  double Vector_GMLS_GradZ = (dimension>2) ? output_gradient_vector(i,2) : 0;
394 
395  // Obtain adjusted value
396  if (constraint_name == "NEUMANN_GRAD_SCALAR") {
397  double bx_i_scalar = scalar_basis_gmls.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 0, num_neigh_i);
398  double by_i_scalar = scalar_basis_gmls.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 1, num_neigh_i);
399  double bz_i_scalar = scalar_basis_gmls.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 2, num_neigh_i);
400  Scalar_GMLS_GradX = Scalar_GMLS_GradX + bx_i_scalar*g;
401  Scalar_GMLS_GradY = Scalar_GMLS_GradY + by_i_scalar*g;
402  Scalar_GMLS_GradZ = Scalar_GMLS_GradZ + bz_i_scalar*g;
403 
404  double bx_i_vector = vector_basis_gmls.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 0, num_neigh_i);
405  double by_i_vector = vector_basis_gmls.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 1, num_neigh_i);
406  double bz_i_vector = vector_basis_gmls.getAlpha0TensorTo1Tensor(GradientOfScalarPointEvaluation, i, 2, num_neigh_i);
407  Vector_GMLS_GradX = Vector_GMLS_GradX + bx_i_vector*g;
408  Vector_GMLS_GradY = Vector_GMLS_GradY + by_i_vector*g;
409  Vector_GMLS_GradZ = Vector_GMLS_GradZ + bz_i_vector*g;
410  }
411 
412  // check gradient - scalar gmls
413  if(std::abs(actual_Gradient[0] - Scalar_GMLS_GradX) > failure_tolerance) {
414  all_passed = false;
415  std::cout << i << " Failed Scalar GradX by: " << std::abs(actual_Gradient[0] - Scalar_GMLS_GradX) << std::endl;
416  if (dimension>1) {
417  if(std::abs(actual_Gradient[1] - Scalar_GMLS_GradY) > failure_tolerance) {
418  all_passed = false;
419  std::cout << i << " Failed Scalar GradY by: " << std::abs(actual_Gradient[1] - Scalar_GMLS_GradY) << std::endl;
420  }
421  }
422  if (dimension>2) {
423  if(std::abs(actual_Gradient[2] - Scalar_GMLS_GradZ) > failure_tolerance) {
424  all_passed = false;
425  std::cout << i << " Failed Scalar GradZ by: " << std::abs(actual_Gradient[2] - Scalar_GMLS_GradZ) << std::endl;
426  }
427  }
428  }
429 
430  // check gradient - vector gmls
431  if(std::abs(actual_Gradient[0] - Vector_GMLS_GradX) > failure_tolerance) {
432  all_passed = false;
433  std::cout << i << " Failed Vector GradX by: " << std::abs(actual_Gradient[0] - Vector_GMLS_GradX) << std::endl;
434  if (dimension>1) {
435  if(std::abs(actual_Gradient[1] - Vector_GMLS_GradY) > failure_tolerance) {
436  all_passed = false;
437  std::cout << i << " Failed Vector GradY by: " << std::abs(actual_Gradient[1] - Vector_GMLS_GradY) << std::endl;
438  }
439  }
440  if (dimension>2) {
441  if(std::abs(actual_Gradient[2] - Vector_GMLS_GradZ) > failure_tolerance) {
442  all_passed = false;
443  std::cout << i << " Failed Vector GradZ by: " << std::abs(actual_Gradient[2] - Vector_GMLS_GradZ) << std::endl;
444  }
445  }
446  }
447  }
448 
449 
450  //! [Check That Solutions Are Correct]
451  // popRegion hidden from tutorial
452  // stop timing comparison loop
453  Kokkos::Profiling::popRegion();
454  //! [Finalize Program]
455 
456 } // end of code block to reduce scope, causing Kokkos View de-allocations
457 // otherwise, Views may be deallocating when we call Kokkos::finalize() later
458 
459 // finalize Kokkos and MPI (if available)
460 Kokkos::finalize();
461 #ifdef COMPADRE_USE_MPI
462 MPI_Finalize();
463 #endif
464 
465 // output to user that test passed or failed
466 if(all_passed) {
467  fprintf(stdout, "Passed test \n");
468  return 0;
469 } else {
470  fprintf(stdout, "Failed test \n");
471  return -1;
472 }
473 
474 } // main
475 
476 
477 //! [Finalize Program]
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...
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...
KOKKOS_INLINE_FUNCTION void trueGradient(double *ans, double x, double y, double z, int order, int dimension)
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 divergence of a vector (results in a scalar)
Point evaluation of the gradient of a scalar.
KOKKOS_INLINE_FUNCTION double trueLaplacian(double x, double y, double z, int order, int dimension)
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 trueSolution(double x, double y, double z, int order, int dimension)
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...