Compadre 1.5.5
Loading...
Searching...
No Matches
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>
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//! [Ambient to Local Back To Ambient Helper Function]
26void 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
59int main (int argc, char* args[]) {
60
61// initializes MPI (if available) with command line arguments given
62#ifdef COMPADRE_USE_MPI
63MPI_Init(&argc, &args);
64#endif
65
66// initializes Kokkos with command line arguments given
67Kokkos::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.setCurvatureWeightingParameter(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.setWeightingParameter(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.
336 GMLS my_GMLS_vector(ReconstructionSpace::VectorTaylorPolynomial, VaryingManifoldVectorPointSample,
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.setCurvatureWeightingParameter(2);
349 my_GMLS_vector.setWeightingType(WeightingFunctionType::Power);
350 my_GMLS_vector.setWeightingParameter(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 //
377 GMLS my_GMLS_vector_of_scalar_clones(ReconstructionSpace::VectorOfScalarClonesTaylorPolynomial, VaryingManifoldVectorPointSample,
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.setCurvatureWeightingParameter(2);
390 my_GMLS_vector_of_scalar_clones.setWeightingType(WeightingFunctionType::Power);
391 my_GMLS_vector_of_scalar_clones.setWeightingParameter(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>
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)
715Kokkos::finalize();
716#ifdef COMPADRE_USE_MPI
717MPI_Finalize();
718#endif
719
720return 0;
721
722} // main
723
724
725//! [Finalize Program]
#define TO_GLOBAL(variable)
Kokkos::View< double **, layout_right, host_execution_space, Kokkos::MemoryTraits< Kokkos::Unmanaged > > host_scratch_matrix_right_type
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 void curl_sphere_harmonic54(double *curl, double x, double y, double z)
void AmbientLocalAmbient(XYZ &u, double *T_data, double *P_data)
[Ambient to Local Back To Ambient Helper Function]
int main(int argc, char *args[])
[Ambient to Local Back To Ambient Helper Function]
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 setAdditionalEvaluationSitesData(view_type_1 additional_evaluation_indices, view_type_2 additional_evaluation_coordinates)
(OPTIONAL) Sets additional evaluation sites for each target site
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....
decltype(_T) * getTangentDirections()
Get a view (device) of all tangent direction bundles.
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.
decltype(_prestencil_weights) getPrestencilWeights() const
Get a view (device) of all rank 2 preprocessing tensors This is a rank 5 tensor that is able to provi...
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.
constexpr SamplingFunctional PointSample
Available sampling functionals.
@ CurlOfVectorPointEvaluation
Point evaluation of the curl of a vector (results in a vector)
@ GaussianCurvaturePointEvaluation
Point evaluation of Gaussian curvature.
@ VectorPointEvaluation
Point evaluation of a vector (reconstructs entire vector at once, requiring a ReconstructionSpace hav...
@ ScalarPointEvaluation
Point evaluation of a scalar.