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