Compadre 1.5.5
Loading...
Searching...
No Matches
Compadre_GMLS.hpp
Go to the documentation of this file.
1#ifndef _COMPADRE_GMLS_HPP_
2#define _COMPADRE_GMLS_HPP_
3
4#include "Compadre_Config.h"
6
7#include "Compadre_Misc.hpp"
18
19namespace Compadre {
20
21class Evaluator;
22struct GMLSBasisData;
23struct GMLSSolutionData;
24
25//! Generalized Moving Least Squares (GMLS)
26/*!
27* This class sets up a batch of GMLS problems from a given set of neighbor lists, target sites, and source sites.
28* GMLS requires a target functional, reconstruction space, and sampling functional to be specified.
29* For a given choice of reconstruction space and sampling functional, multiple targets can be generated with very little
30* additional computation, which is why this class allows for multiple target functionals to be specified.
31*/
32class GMLS {
33
34friend class Evaluator;
35
36public:
37
39 Kokkos::View<double**, layout_right>,
42
44
45 typedef Kokkos::View<double**, layout_right> coordinates_type;
46
47private:
48
49 // matrices that may be needed for matrix factorization on the device
50 // supports batched matrix factorization dispatch
51
52 //! contains weights for all problems
53 Kokkos::View<double*> _w;
54
55 //! P*sqrt(w) matrix for all problems
56 Kokkos::View<double*> _P;
57
58 //! sqrt(w)*Identity matrix for all problems, later holds polynomial coefficients for all problems
59 Kokkos::View<double*> _RHS;
60
61 //! stores evaluations of targets applied to basis
62 Kokkos::View<double*> _Z;
63
64 //! Rank 3 tensor for high order approximation of tangent vectors for all problems. First rank is
65 //! for the target index, the second is for the local direction to the manifolds 0..(_dimensions-1)
66 //! are tangent, _dimensions is the normal, and the third is for the spatial dimension (_dimensions)
67 Kokkos::View<double*> _T;
68
69 //! Rank 2 tensor for high order approximation of tangent vectors for all problems. First rank is
70 //! for the target index, the second is for the spatial dimension (_dimensions)
71 Kokkos::View<double*> _ref_N;
72
73 //! tangent vectors information (host)
74 Kokkos::View<double*>::HostMirror _host_T;
75
76 //! reference outward normal vectors information (host)
77 Kokkos::View<double*>::HostMirror _host_ref_N;
78
79 //! metric tensor inverse for all problems
80 Kokkos::View<double*> _manifold_metric_tensor_inverse;
81
82 //! curvature polynomial coefficients for all problems
83 Kokkos::View<double*> _manifold_curvature_coefficients;
84
85 //! _dimension-1 gradient values for curvature for all problems
86 Kokkos::View<double*> _manifold_curvature_gradient;
87
88 //! Extra data available to basis functions (optional)
89 Kokkos::View<double**, layout_right> _source_extra_data;
90
91 //! Extra data available to target operations (optional)
92 Kokkos::View<double**, layout_right> _target_extra_data;
93
94 //! connections between points and neighbors
96
97 //! h supports determined through neighbor search (device)
98 Kokkos::View<double*> _epsilons;
99
100 //! generated weights for nontraditional samples required to transform data into expected sampling
101 //! functional form (device).
102 Kokkos::View<double*****, layout_right> _prestencil_weights;
103
104 //! generated weights for nontraditional samples required to transform data into expected sampling
105 //! functional form (host)
106 Kokkos::View<const double*****, layout_right>::HostMirror _host_prestencil_weights;
107
108 //! (OPTIONAL) connections between additional points and neighbors
110
111 //! Solution Set (contains all alpha values from solution and alpha layout methods)
112 // _h_ss is private so that getSolutionSetHost() must be called
113 // which ensures that the copy of the solution to device is necessary
116
117 //! order of basis for polynomial reconstruction
119
120 //! order of basis for curvature reconstruction
122
123 //! dimension of basis for polynomial reconstruction
124 int _NP;
125
126 //! spatial dimension of the points, set at class instantiation only
128
129 //! dimension of the problem, set at class instantiation only. For manifolds, generally _global_dimensions-1
131
132 //! dimension of the problem, set at class instantiation only
134
135 //! reconstruction space for GMLS problems, set at GMLS class instantiation
137
138 //! actual rank of reconstruction basis
140
141 //! solver type for GMLS problem - can be QR, SVD or LU
143
144 //! problem type for GMLS problem, can also be set to STANDARD for normal or MANIFOLD for manifold problems <br>
145 //! <b>NOTE: can only be set at object instantiation</b>
147
148 //! constraint type for GMLS problem
150
151 //! polynomial sampling functional used to construct P matrix, set at GMLS class instantiation <br>
152 //! <b>NOTE: can only be set at object instantiation</b>
154
155 //! generally the same as _polynomial_sampling_functional, but can differ if specified at
156 // can only be set at object instantiation
157 //! GMLS class instantiation
159
160 //! vector containing target functionals to be applied for curvature
161 Kokkos::View<TargetOperation*> _curvature_support_operations;
162
163 //! vector containing target functionals to be applied for reconstruction problem (device)
164 Kokkos::View<TargetOperation*> _operations;
165
166 //! vector containing target functionals to be applied for reconstruction problem (host)
167 Kokkos::View<TargetOperation*>::HostMirror _host_operations;
168
169 //! weighting kernel type for GMLS
171
172 //! weighting kernel type for curvature problem
174
175 //! first parameter to be used for weighting kernel
177
178 //! second parameter to be used for weighting kernel
180
181 //! first parameter to be used for weighting kernel for curvature
183
184 //! second parameter to be used for weighting kernel for curvature
186
187 //! dimension of the reconstructed function
188 //! e.g. reconstruction of vector on a 2D manifold in 3D would have _basis_multiplier of 2
190
191 //! actual dimension of the sampling functional
192 //! e.g. reconstruction of vector on a 2D manifold in 3D would have _basis_multiplier of 2
193 //! e.g. in 3D, a scalar will be 1, a vector will be 3, and a vector of reused scalars will be 1
195
196 //! effective dimension of the data sampling functional
197 //! e.g. in 3D, a scalar will be 1, a vector will be 3, and a vector of reused scalars will be 3
199
200 //! whether or not the orthonormal tangent directions were provided by the user. If they are not,
201 //! then for the case of calculations on manifolds, a GMLS approximation of the tangent space will
202 //! be made and stored for use.
204
205 //! whether or not the reference outward normal directions were provided by the user.
207
208 //! whether or not to use reference outward normal directions to orient the surface in a manifold problem.
210
211 //! whether entire calculation was computed at once
212 //! the alternative is that it was broken up over many smaller groups, in which case
213 //! this is false, and so the _RHS matrix can not be stored or requested
215
216 //! whether polynomial coefficients were requested to be stored (in a state not yet applied to data)
218
219 //! initial index for current batch
221
222 //! determines scratch level spaces and is used to call kernels
224
225 //! order of exact polynomial integration for quadrature rule
227
228 //! dimension of quadrature rule
230
231 //! quadrature rule type
232 std::string _quadrature_type;
233
234 //! manages and calculates quadrature
236
237private:
238
239/** @name Private Modifiers
240 * Private function because information lives on the device
241 */
242///@{
243
244 //! (OPTIONAL)
245 //! Sets additional points for evaluation of target operation on polynomial reconstruction.
246 //! If this is never called, then the target sites are the only locations where the target
247 //! operations will be evaluated and applied to polynomial reconstructions.
248 template <typename view_type>
249 void setAuxiliaryEvaluationCoordinates(view_type evaluation_coordinates) {
250 // allocate memory on device
251 auto additional_evaluation_coordinates = coordinates_type("device additional evaluation coordinates",
252 evaluation_coordinates.extent(0), evaluation_coordinates.extent(1));
253
254 typedef typename view_type::memory_space input_array_memory_space;
255 if (std::is_same<input_array_memory_space, device_memory_space>::value) {
256 // check if on the device, then copy directly
257 // if it is, then it doesn't match the internal layout we use
258 // then copy to the host mirror
259 // switches potential layout mismatches
260 Kokkos::deep_copy(additional_evaluation_coordinates, evaluation_coordinates);
261 } else {
262 // if is on the host, copy to the host mirror
263 // then copy to the device
264 // switches potential layout mismatches
265 auto host_additional_evaluation_coordinates = Kokkos::create_mirror_view(additional_evaluation_coordinates);
266 Kokkos::deep_copy(host_additional_evaluation_coordinates, evaluation_coordinates);
267 // switches memory spaces
268 Kokkos::deep_copy(additional_evaluation_coordinates, host_additional_evaluation_coordinates);
269 }
270 this->resetCoefficientData();
271 _additional_pc.setSourceCoordinates(additional_evaluation_coordinates);
272 }
273
274 //! (OPTIONAL)
275 //! Sets additional points for evaluation of target operation on polynomial reconstruction.
276 //! If this is never called, then the target sites are the only locations where the target
277 //! operations will be evaluated and applied to polynomial reconstructions. (device)
278 template <typename view_type>
280 this->resetCoefficientData();
281 _additional_pc.setSourceCoordinates(evaluation_coordinates);
282 }
283
284 //! (OPTIONAL)
285 //! Sets the additional target evaluation indices list information from compressed row format (if same view_type)
286 template <typename view_type>
287 typename std::enable_if<view_type::rank==1&&std::is_same<neighbor_lists_type::internal_view_type,view_type>::value==1, void>::type
288 setAuxiliaryEvaluationIndicesLists(view_type additional_evaluation_indices, view_type number_of_neighbors_list) {
289
290 auto additional_nla = NeighborLists<view_type>(additional_evaluation_indices, number_of_neighbors_list);
291 this->resetCoefficientData();
292 _additional_pc.setNeighborLists(additional_nla);
293
294 }
295
296 //! (OPTIONAL)
297 //! Sets the additional target evaluation indices list information from compressed row format (if different view_type)
298 template <typename view_type>
299 typename std::enable_if<view_type::rank==1&&std::is_same<neighbor_lists_type::internal_view_type,view_type>::value==0, void>::type
300 setAuxiliaryEvaluationIndicesLists(view_type additional_evaluation_indices, view_type number_of_neighbors_list) {
301
302 typedef neighbor_lists_type::internal_view_type gmls_view_type;
303 gmls_view_type d_additional_evaluation_indices("compressed row additional evaluation indices lists data", additional_evaluation_indices.extent(0));
304 gmls_view_type d_number_of_neighbors_list("number of additional evaluation indices", number_of_neighbors_list.extent(0));
305 Kokkos::deep_copy(d_additional_evaluation_indices, additional_evaluation_indices);
306 Kokkos::deep_copy(d_number_of_neighbors_list, number_of_neighbors_list);
307 Kokkos::fence();
308 auto additional_nla = NeighborLists<gmls_view_type>(d_additional_evaluation_indices, d_number_of_neighbors_list);
309 this->resetCoefficientData();
310 _additional_pc.setNeighborLists(additional_nla);
311
312 }
313
314 //! (OPTIONAL)
315 //! Sets the additional target evaluation indices list information. Should be # targets x maximum number of indices
316 //! evaluation indices for any target + 1. first entry in every row should be the number of indices for the corresponding target.
317 template <typename view_type>
318 typename std::enable_if<view_type::rank==2, void>::type setAuxiliaryEvaluationIndicesLists(view_type additional_evaluation_indices) {
319
320 auto additional_nla = Convert2DToCompressedRowNeighborLists<decltype(additional_evaluation_indices), Kokkos::View<int*> >(additional_evaluation_indices);
321 this->resetCoefficientData();
322 _additional_pc.setNeighborLists(additional_nla);
323
324 }
325
326///@}
327
328
329/** @name Private Utility
330 *
331 */
332///@{
333
334 //! Parses a string to determine solver type
335 static DenseSolverType parseSolverType(const std::string& dense_solver_type) {
336 std::string solver_type_to_lower = dense_solver_type;
337 transform(solver_type_to_lower.begin(), solver_type_to_lower.end(), solver_type_to_lower.begin(), ::tolower);
338 if (solver_type_to_lower == "lu") {
339 return DenseSolverType::LU;
340 } else {
341 return DenseSolverType::QR;
342 }
343 }
344
345 //! Parses a string to determine problem type
346 static ProblemType parseProblemType(const std::string& problem_type) {
347 std::string problem_type_to_lower = problem_type;
348 transform(problem_type_to_lower.begin(), problem_type_to_lower.end(), problem_type_to_lower.begin(), ::tolower);
349 if (problem_type_to_lower == "standard") {
351 } else if (problem_type_to_lower == "manifold") {
353 } else {
355 }
356 }
357
358 //! Parses a string to determine constraint type
359 static ConstraintType parseConstraintType(const std::string& constraint_type) {
360 std::string constraint_type_to_lower = constraint_type;
361 transform(constraint_type_to_lower.begin(), constraint_type_to_lower.end(), constraint_type_to_lower.begin(), ::tolower);
362 if (constraint_type_to_lower == "none") {
364 } else if (constraint_type_to_lower == "neumann_grad_scalar") {
366 } else {
368 }
369 }
370
371///@}
372
373public:
374
375/** @name Instantiation / Destruction
376 *
377 */
378///@{
379
380 //! Maximal constructor, not intended for users
381 GMLS(ReconstructionSpace reconstruction_space,
382 const SamplingFunctional polynomial_sampling_strategy,
383 const SamplingFunctional data_sampling_strategy,
384 const int poly_order,
385 const int dimensions,
386 const DenseSolverType dense_solver_type,
387 const ProblemType problem_type,
388 const ConstraintType constraint_type,
389 const int manifold_curvature_poly_order) :
390 _poly_order(poly_order),
391 _curvature_poly_order(manifold_curvature_poly_order),
392 _dimensions(dimensions),
393 _reconstruction_space(reconstruction_space),
394 _dense_solver_type(dense_solver_type),
395 _problem_type(problem_type),
396 _constraint_type(constraint_type),
398 && (polynomial_sampling_strategy == VectorPointSample)) ? ManifoldVectorPointSample : polynomial_sampling_strategy),
400 && (data_sampling_strategy == VectorPointSample)) ? ManifoldVectorPointSample : data_sampling_strategy)
401 {
402
403 compadre_assert_release(poly_order<11 && "Unsupported polynomial order (>=11).");
404 compadre_assert_release(manifold_curvature_poly_order<11 && "Unsupported curvature polynomial order (>=11).");
405 _NP = this->getNP(_poly_order, dimensions, _reconstruction_space);
406 Kokkos::fence();
407
408 // register curvature operations for manifold problems
410 _curvature_support_operations = Kokkos::View<TargetOperation*>
411 ("operations needed for manifold gradient reconstruction", 1);
412 auto curvature_support_operations_mirror =
413 Kokkos::create_mirror_view(_curvature_support_operations);
414 curvature_support_operations_mirror(0) =
416 Kokkos::deep_copy(_curvature_support_operations, curvature_support_operations_mirror);
417 }
418
419 // various initializations
420
423 _weighting_p = 2;
424 _weighting_n = 1;
427
429
432
437 _store_PTWP_inv_PTW = false;
438
440
441 _global_dimensions = dimensions;
443 _local_dimensions = dimensions-1;
444 } else {
445 _local_dimensions = dimensions;
446 }
447
450
457 }
458
459 //! Maximal constructor, but with string arguments
460 GMLS(ReconstructionSpace reconstruction_space,
461 const SamplingFunctional polynomial_sampling_strategy,
462 const SamplingFunctional data_sampling_strategy,
463 const int poly_order,
464 const int dimensions = 3,
465 const std::string dense_solver_type = std::string("QR"),
466 const std::string problem_type = std::string("STANDARD"),
467 const std::string constraint_type = std::string("NO_CONSTRAINT"),
468 const int manifold_curvature_poly_order = 2)
469 : GMLS(reconstruction_space, polynomial_sampling_strategy, data_sampling_strategy, poly_order, dimensions, parseSolverType(dense_solver_type), parseProblemType(problem_type), parseConstraintType(constraint_type), manifold_curvature_poly_order) {}
470
471 //! Constructor for the case when the data sampling functional does not match the polynomial
472 //! sampling functional. Only case anticipated is staggered Laplacian.
473 GMLS(const int poly_order,
474 const int dimensions = 3,
475 const std::string dense_solver_type = std::string("QR"),
476 const std::string problem_type = std::string("STANDARD"),
477 const std::string constraint_type = std::string("NO_CONSTRAINT"),
478 const int manifold_curvature_poly_order = 2)
479 : GMLS(ReconstructionSpace::VectorOfScalarClonesTaylorPolynomial, VectorPointSample, VectorPointSample, poly_order, dimensions, dense_solver_type, problem_type, constraint_type, manifold_curvature_poly_order) {}
480
481 //! Constructor for the case when nonstandard sampling functionals or reconstruction spaces
482 //! are to be used. Reconstruction space and sampling strategy can only be set at instantiation.
483 GMLS(ReconstructionSpace reconstruction_space,
484 SamplingFunctional dual_sampling_strategy,
485 const int poly_order,
486 const int dimensions = 3,
487 const std::string dense_solver_type = std::string("QR"),
488 const std::string problem_type = std::string("STANDARD"),
489 const std::string constraint_type = std::string("NO_CONSTRAINT"),
490 const int manifold_curvature_poly_order = 2)
491 : GMLS(reconstruction_space, dual_sampling_strategy, dual_sampling_strategy, poly_order, dimensions, dense_solver_type, problem_type, constraint_type, manifold_curvature_poly_order) {}
492
493///@}
494
495/** @name Public Utility
496 *
497 */
498///@{
499
500 //! Returns size of the basis for a given polynomial order and dimension
501 //! General to dimension 1..3 and polynomial order m
502 //! The divfree options will return the divergence-free basis if true
503 KOKKOS_INLINE_FUNCTION
504 static int getNP(const int m, const int dimension = 3, const ReconstructionSpace r_space = ReconstructionSpace::ScalarTaylorPolynomial) {
506 return DivergenceFreePolynomialBasis::getSize(m, dimension);
507 } else if (r_space == ReconstructionSpace::BernsteinPolynomial) {
508 return BernsteinPolynomialBasis::getSize(m, dimension);
509 } else {
510 return ScalarTaylorPolynomialBasis::getSize(m, dimension);
511 }
512 }
513
514 //! Returns number of neighbors needed for unisolvency for a given basis order and dimension
515 KOKKOS_INLINE_FUNCTION
516 static int getNN(const int m, const int dimension = 3, const ReconstructionSpace r_space = ReconstructionSpace::ScalarTaylorPolynomial) {
517 // may need div-free argument in the future
518 const int np = getNP(m, dimension, r_space);
519 int nn = np;
520 switch (dimension) {
521 case 3:
522 nn = np * (1.7 + m*0.1);
523 break;
524 case 2:
525 nn = np * (1.4 + m*0.03);
526 break;
527 case 1:
528 nn = np * 1.1;
529 }
530 return nn;
531 }
532
533 /*! \brief Evaluates the weighting kernel
534 \param r [in] - Euclidean distance of relative vector. Euclidean distance of (target - neighbor) in some basis.
535 \param h [in] - window size. Kernel is guaranteed to take on a value of zero if it exceeds h.
536 \param weighting_type [in] - weighting type to be evaluated as the kernel. e,g. power, Gaussian, etc..
537 \param p [in] - parameter to be given to the kernel (see Wab definition for context).
538 \param n [in] - parameter to be given to the kernel (see Wab definition for context).
539 */
540 KOKKOS_INLINE_FUNCTION
541 static double Wab(const double r, const double h, const WeightingFunctionType& weighting_type, const int p, const int n) {
542 if (weighting_type == WeightingFunctionType::Power) {
543 // compactly supported on [0,h]
544 // (1 - |r/h|^n)^p
545 // p=0,n=1 -> Uniform, boxcar
546 // p=1,n=1 -> triangular
547 // p=1,n=2 -> Epanechnikov, parabolic
548 // p=2,n=2 -> Quartic, biweight
549 // p=3,n=2 -> Triweight
550 // p=3,n=3 -> Tricube
551 double abs_r_over_h_to_n = std::abs(r/h);
552 if (n>1) abs_r_over_h_to_n = std::pow(abs_r_over_h_to_n, n);
553 return std::pow(1.0-abs_r_over_h_to_n, p) * double(1.0-abs_r_over_h_to_n>0.0);
554 } else if (weighting_type == WeightingFunctionType::CubicSpline) {
555 // compactly supported on [0,h]
556 // invariant to p and n
557 double x = std::abs(r/h);
558 return ((1-x)+x*(1-x)*(1-2*x)) * double(x<=1.0);
559 } else if (weighting_type == WeightingFunctionType::Cosine) {
560 // compactly supported on [0,h]
561 double pi = 3.14159265358979323846;
562 double abs_r_over_h_to_n = std::abs(r/h);
563 return std::cos(0.5*pi*r/h) * double(1.0-abs_r_over_h_to_n>0.0);
564 } else if (weighting_type == WeightingFunctionType::Gaussian) {
565 // NOT compactly supported on [0,h], but approximately 0 at h with >> p
566 // invariant to n, p is number of standard deviations at distance h
567 // 2.5066282746310002416124 = sqrt(2*pi)
568 double h_over_p = h/p;
569 double abs_r_over_h_to_n = std::abs(r/h);
570 return double(1.0-abs_r_over_h_to_n>0.0)/( h_over_p * 2.5066282746310002416124 ) * std::exp(-.5*r*r/(h_over_p*h_over_p));
571 } else if (weighting_type == WeightingFunctionType::Sigmoid) {
572 // NOT compactly supported on [0,h], but approximately 0 at h with >> p
573 // n=0 is sigmoid, n==2 is logistic, with larger p making Wab decay more quickly
574 double abs_r_over_h_to_n = std::abs(r/h);
575 return double(1.0-abs_r_over_h_to_n>0.0) / (std::exp(p*r) + std::exp(-p*r) + n);
576 } else { // unsupported type
577 compadre_kernel_assert_release(false && "Invalid WeightingFunctionType selected.");
578 return 0;
579 }
580 }
581
582
583///@}
584
585/** @name Accessors
586 * Retrieve member variables through public member functions
587 */
588///@{
589
590
591 //! Returns (size of the basis used in instance's polynomial reconstruction) x (data input dimension)
593 host_managed_local_index_type sizes("sizes", 2);
594 sizes(0) = _basis_multiplier*_NP;
596 return sizes;
597 }
598
599 //! Returns size of the basis used in instance's polynomial reconstruction
601 auto sizes = this->getPolynomialCoefficientsDomainRangeSize();
602 return sizes(0);
603 }
604
605 //! Returns 2D array size in memory on which coefficients are stored
607 auto M_by_N = this->getPolynomialCoefficientsDomainRangeSize();
608 compadre_assert_release(_entire_batch_computed_at_once
609 && "Entire batch not computed at once, so getFullPolynomialCoefficientsBasis() can not be called.");
611 && "generateAlphas() called with keep_coefficients set to false.");
612 host_managed_local_index_type sizes("sizes", 2);
614 getRHSDims(_dense_solver_type, _constraint_type, _reconstruction_space, _dimensions, M_by_N[1], M_by_N[0], sizes(0), sizes(1));
615 } else {
616 getPDims(_dense_solver_type, _constraint_type, _reconstruction_space, _dimensions, M_by_N[1], M_by_N[0], sizes(1), sizes(0));
617 }
618 return sizes;
619 }
620
621 //! Dimension of the GMLS problem, set only at class instantiation
622 int getDimensions() const { return _dimensions; }
623
624 //! Dimension of the GMLS problem's point data (spatial description of points in ambient space), set only at class instantiation
626
627 //! Local dimension of the GMLS problem (less than global dimension if on a manifold), set only at class instantiation
628 int getLocalDimensions() const { return _local_dimensions; }
629
630 //! Get dense solver type
632
633 //! Get problem type
635
636 //! Get constraint type
638
639 //! Get basis order used for reconstruction
640 int getPolynomialOrder() const { return _poly_order; }
641
642 //! Get basis order used for curvature reconstruction
644
645 //! Type for weighting kernel for GMLS problem
647
648 //! Type for weighting kernel for curvature
650
651 //! Get parameter for weighting kernel for GMLS problem
652 int getWeightingParameter(const int index = 0) const {
653 if (index==1) {
654 return _weighting_n;
655 } else {
656 return _weighting_p;
657 }
658 }
659
660 //! Get parameter for weighting kernel for curvature
661 int getManifoldWeightingParameter(const int index = 0) const {
662 if (index==1) {
664 } else {
666 }
667 }
668
669 //! Number of quadrature points
671
672 //! Order of quadrature points
674
675 //! Dimensions of quadrature points
677
678 //! Type of quadrature points
679 std::string getQuadratureType() const { return _quadrature_type; }
680
681 //! Get neighbor list accessor
683 return const_cast<neighbor_lists_type*>(&_pc._nla);
684 }
685
686 //! Get a view (device) of all point connection info
687 decltype(_pc)* getPointConnections() { return &_pc; }
688
689 //! (OPTIONAL) Get additional evaluation sites neighbor list-like accessor
691 return const_cast<neighbor_lists_type*>(&_additional_pc._nla);
692 }
693
694 //! (OPTIONAL) Get a view (device) of all additional evaluation point connection info
696
697 //! Get a view (device) of all window sizes
698 decltype(_epsilons)* getWindowSizes() { return &_epsilons; }
699
700 //! Get a view (device) of all tangent direction bundles.
701 decltype(_T)* getTangentDirections() { return &_T; }
702
703 //! Get a view (device) of all reference outward normal directions.
705
706 //! Get component of tangent or normal directions for manifold problems
707 double getTangentBundle(const int target_index, const int direction, const int component) const {
708 // Component index 0.._dimensions-2 will return tangent direction
709 // Component index _dimensions-1 will return the normal direction
710 scratch_matrix_right_type::HostMirror
711 T(_host_T.data() + target_index*_dimensions*_dimensions, _dimensions, _dimensions);
712 return T(direction, component);
713 }
714
715 //! Get component of tangent or normal directions for manifold problems
716 double getReferenceNormalDirection(const int target_index, const int component) const {
718 "getRefenceNormalDirection called, but reference outwrad normal directions were never provided.");
719 scratch_vector_type::HostMirror
720 ref_N(_host_ref_N.data() + target_index*_dimensions, _dimensions);
721 return ref_N(component);
722 }
723
724 //! Get a view (device) of all rank 2 preprocessing tensors
725 //! This is a rank 5 tensor that is able to provide data transformation
726 //! into a form that GMLS is able to operate on. The ranks are as follows:
727 //!
728 //! 1 - Either size 2 if it operates on the target site and neighbor site (staggered schemes)
729 //! or 1 if it operates only on the neighbor sites (almost every scheme)
730 //!
731 //! 2 - If the data transform varies with each target site (but could be the same for each neighbor of that target site), then this is the number of target sites
732 //!
733 //! 3 - If the data transform varies with each neighbor of each target site, then this is the number of neighbors for each respective target (max number of neighbors for all target sites is its uniform size)
734 //!
735 //! 4 - Data transform resulting in rank 1 data for the GMLS operator will have size _local_dimensions, otherwise 1
736 //!
737 //! 5 - Data transform taking in rank 1 data will have size _global_dimensions, otherwise 1
739 return _prestencil_weights;
740 }
741
742 //! Get a view (device) of all polynomial coefficients basis
745 && "Entire batch not computed at once, so getFullPolynomialCoefficientsBasis() can not be called.");
747 && "generateAlphas() called with keep_coefficients set to false.");
749 return _RHS;
750 } else {
751 return _P;
752 }
753 }
754
755 //! Get the polynomial sampling functional specified at instantiation
757
758 //! Get the data sampling functional specified at instantiation (often the same as the polynomial sampling functional)
760
761 //! Get the reconstruction space specified at instantiation
763
764 //! Returns a stencil to transform data from its existing state into the input expected
765 //! for some sampling functionals.
766 double getPreStencilWeight(SamplingFunctional sro, const int target_index, const int neighbor_index, bool for_target, const int output_component = 0, const int input_component = 0) const {
767 // for certain sampling strategies, linear combinations of the neighbor and target value are needed
768 // for the traditional PointSample, this value is 1 for the neighbor and 0 for the target
769 if (sro == PointSample ) {
770 if (for_target) return 0; else return 1;
771 }
772
773 // these check conditions on the sampling operator and change indexing on target and neighbors
774 // in order to reuse information, such as if the same data transformation is used, regardless
775 // of target site or neighbor site
776 const int target_index_in_weights =
779 target_index : 0;
780 const int neighbor_index_in_weights =
782 neighbor_index : 0;
783
784 return _host_prestencil_weights((int)for_target, target_index_in_weights, neighbor_index_in_weights,
785 output_component, input_component);
786 }
787
788 //! Get solution set on host
789 decltype(_h_ss)* getSolutionSetHost(bool alpha_validity_check=true) {
791 // solution solved for on device, but now solution
792 // requested on the host
794 }
795 compadre_assert_release((!alpha_validity_check || _h_ss._contains_valid_alphas) &&
796 "getSolutionSetHost() called with invalid alpha values.");
797 return &_h_ss;
798 }
799
800 //! Get solution set on device
801 decltype(_d_ss)* getSolutionSetDevice(bool alpha_validity_check=true) {
802 compadre_assert_release((!alpha_validity_check || _d_ss._contains_valid_alphas) &&
803 "getSolutionSetDevice() called with invalid alpha values.");
804 return &_d_ss;
805 }
806
807 //! Check if GMLS solution set contains valid alpha values (has generateAlphas been called)
808 bool containsValidAlphas() const { return this->_d_ss._contains_valid_alphas; }
809
810 //! Get GMLS solution data
812
813 //! Get GMLS basis data
814 const GMLSBasisData extractBasisData() const;
815
816///@}
817
818
819/** @name Modifiers
820 * Changed member variables through public member functions
821 */
822///@{
823
825 if (_RHS.extent(0) > 0)
826 _RHS = Kokkos::View<double*>("RHS",0);
829 }
830
831 //! Sets basic problem data (neighbor lists, source coordinates, and target coordinates)
832 template<typename view_type_1, typename view_type_2, typename view_type_3, typename view_type_4>
834 view_type_1 neighbor_lists,
835 view_type_2 source_coordinates,
836 view_type_3 target_coordinates,
837 view_type_4 epsilons) {
838 this->setNeighborLists<view_type_1>(neighbor_lists);
839 this->setSourceSites<view_type_2>(source_coordinates);
840 this->setTargetSites<view_type_3>(target_coordinates);
841 this->setWindowSizes<view_type_4>(epsilons);
842 }
843
844 //! Sets basic problem data (neighbor lists data, number of neighbors list, source coordinates, and target coordinates)
845 template<typename view_type_1, typename view_type_2, typename view_type_3, typename view_type_4>
847 view_type_1 cr_neighbor_lists,
848 view_type_1 number_of_neighbors_list,
849 view_type_2 source_coordinates,
850 view_type_3 target_coordinates,
851 view_type_4 epsilons) {
852 this->setNeighborLists<view_type_1>(cr_neighbor_lists, number_of_neighbors_list);
853 this->setSourceSites<view_type_2>(source_coordinates);
854 this->setTargetSites<view_type_3>(target_coordinates);
855 this->setWindowSizes<view_type_4>(epsilons);
856 }
857
858 //! (OPTIONAL) Sets additional evaluation sites for each target site
859 template<typename view_type_1, typename view_type_2>
861 view_type_1 additional_evaluation_indices,
862 view_type_2 additional_evaluation_coordinates) {
863 this->setAuxiliaryEvaluationIndicesLists<view_type_1>(additional_evaluation_indices);
864 this->setAuxiliaryEvaluationCoordinates<view_type_2>(additional_evaluation_coordinates);
865 }
866
867 //! (OPTIONAL) Sets additional evaluation sites for each target site
868 template<typename view_type_1, typename view_type_2>
870 view_type_1 cr_additional_evaluation_indices,
871 view_type_1 number_of_additional_evaluation_indices,
872 view_type_2 additional_evaluation_coordinates) {
873 this->setAuxiliaryEvaluationIndicesLists<view_type_1>(cr_additional_evaluation_indices,
874 number_of_additional_evaluation_indices);
875 this->setAuxiliaryEvaluationCoordinates<view_type_2>(additional_evaluation_coordinates);
876 }
877
878 //! Sets neighbor list information from compressed row neighborhood lists data (if same view_type).
879 template <typename view_type>
880 typename std::enable_if<view_type::rank==1&&std::is_same<neighbor_lists_type::internal_view_type,view_type>::value==1, void>::type
881 setNeighborLists(view_type neighbor_lists, view_type number_of_neighbors_list) {
882
883 auto nla = NeighborLists<view_type>(neighbor_lists, number_of_neighbors_list);
884 this->resetCoefficientData();
885 _pc.setNeighborLists(nla);
886 }
887
888 //! Sets neighbor list information from compressed row neighborhood lists data (if different view_type).
889 template <typename view_type>
890 typename std::enable_if<view_type::rank==1&&std::is_same<neighbor_lists_type::internal_view_type,view_type>::value==0, void>::type
891 setNeighborLists(view_type neighbor_lists, view_type number_of_neighbors_list) {
892
893 typedef neighbor_lists_type::internal_view_type gmls_view_type;
894 gmls_view_type d_neighbor_lists("compressed row neighbor lists data", neighbor_lists.extent(0));
895 gmls_view_type d_number_of_neighbors_list("number of neighbors list", number_of_neighbors_list.extent(0));
896 Kokkos::deep_copy(d_neighbor_lists, neighbor_lists);
897 Kokkos::deep_copy(d_number_of_neighbors_list, number_of_neighbors_list);
898 Kokkos::fence();
899 auto nla = NeighborLists<gmls_view_type>(d_neighbor_lists, d_number_of_neighbors_list);
900 this->resetCoefficientData();
901 _pc.setNeighborLists(nla);
902 }
903
904 //! Sets neighbor list information. Should be # targets x maximum number of neighbors for any target + 1.
905 //! first entry in ever row should be the number of neighbors for the corresponding target.
906 template <typename view_type>
907 typename std::enable_if<view_type::rank==2, void>::type setNeighborLists(view_type neighbor_lists) {
908
909 auto nla = Convert2DToCompressedRowNeighborLists<decltype(neighbor_lists), Kokkos::View<int*> >(neighbor_lists);
910 this->resetCoefficientData();
911 _pc.setNeighborLists(nla);
912 }
913
914 //! Sets source coordinate information. Rows of this 2D-array should correspond to neighbor IDs contained in the entries
915 //! of the neighbor lists 2D array.
916 template<typename view_type>
917 void setSourceSites(view_type source_coordinates) {
918
919 // allocate memory on device
920 auto sc = coordinates_type("device neighbor coordinates",
921 source_coordinates.extent(0), source_coordinates.extent(1));
922
923 typedef typename view_type::memory_space input_array_memory_space;
924 if (std::is_same<input_array_memory_space, device_memory_space>::value) {
925 // check if on the device, then copy directly
926 // if it is, then it doesn't match the internal layout we use
927 // then copy to the host mirror
928 // switches potential layout mismatches
929 Kokkos::deep_copy(sc, source_coordinates);
930 } else {
931 // if is on the host, copy to the host mirror
932 // then copy to the device
933 // switches potential layout mismatches
934 auto host_source_coordinates = Kokkos::create_mirror_view(sc);
935 Kokkos::deep_copy(host_source_coordinates, source_coordinates);
936 // switches memory spaces
937 Kokkos::deep_copy(sc, host_source_coordinates);
938 }
939 this->resetCoefficientData();
940 _pc.setSourceCoordinates(sc);
941 }
942
943 //! Sets source coordinate information. Rows of this 2D-array should correspond to neighbor IDs contained in the entries
944 //! of the neighbor lists 2D array.
945 template<typename view_type>
946 void setSourceSites(coordinates_type source_coordinates) {
947 // allocate memory on device
948 this->resetCoefficientData();
949 _pc.setSourceCoordinates(source_coordinates);
950 }
951
952 //! Sets target coordinate information. Rows of this 2D-array should correspond to rows of the neighbor lists.
953 template<typename view_type>
954 void setTargetSites(view_type target_coordinates) {
955 // allocate memory on device
956 auto tc = coordinates_type("device target coordinates",
957 target_coordinates.extent(0), target_coordinates.extent(1));
958
959 typedef typename view_type::memory_space input_array_memory_space;
960 if (std::is_same<input_array_memory_space, device_memory_space>::value) {
961 // check if on the device, then copy directly
962 // if it is, then it doesn't match the internal layout we use
963 // then copy to the host mirror
964 // switches potential layout mismatches
965 Kokkos::deep_copy(tc, target_coordinates);
966 } else {
967 // if is on the host, copy to the host mirror
968 // then copy to the device
969 // switches potential layout mismatches
970 auto host_target_coordinates = Kokkos::create_mirror_view(tc);
971 Kokkos::deep_copy(host_target_coordinates, target_coordinates);
972 // switches memory spaces
973 Kokkos::deep_copy(tc, host_target_coordinates);
974 }
975 if (this->getAdditionalEvaluationIndices()->getNumberOfTargets() != target_coordinates.extent_int(0)) {
977 Kokkos::View<int*>(),
978 Kokkos::View<int*>("number of additional evaluation indices",
979 target_coordinates.extent(0))
980 );
981 }
982 this->resetCoefficientData();
983 _pc.setTargetCoordinates(tc);
984 if (_additional_pc._target_coordinates.extent(0) != _pc._target_coordinates.extent(0)) {
986 }
987 }
988
989 //! Sets target coordinate information. Rows of this 2D-array should correspond to rows of the neighbor lists.
990 template<typename view_type>
991 void setTargetSites(coordinates_type target_coordinates) {
992 if (this->getAdditionalEvaluationIndices()->getNumberOfTargets() != target_coordinates.extent(0)) {
994 Kokkos::View<int*>(),
995 Kokkos::View<int*>("number of additional evaluation indices",
996 target_coordinates.extent(0))
997 );
998 }
999 this->resetCoefficientData();
1000 _pc.setTargetCoordinates(target_coordinates);
1001 if (_additional_pc._target_coordinates.extent(0) != _pc._target_coordinates.extent(0)) {
1002 _additional_pc.setTargetCoordinates(target_coordinates);
1003 }
1004 }
1005
1006 //! Sets window sizes, also called the support of the kernel
1007 template<typename view_type>
1008 void setWindowSizes(view_type epsilons) {
1009
1010 // allocate memory on device
1011 _epsilons = decltype(_epsilons)("device epsilons", epsilons.extent(0));
1012
1013 // copy data from host to device
1014 Kokkos::deep_copy(_epsilons, epsilons);
1015 this->resetCoefficientData();
1016 }
1017
1018 //! Sets window sizes, also called the support of the kernel (device)
1019 template<typename view_type>
1020 void setWindowSizes(decltype(_epsilons) epsilons) {
1021 // allocate memory on device
1022 _epsilons = epsilons;
1023 this->resetCoefficientData();
1024 }
1025
1026 //! (OPTIONAL)
1027 //! Sets orthonormal tangent directions for reconstruction on a manifold. The first rank of this 2D array
1028 //! corresponds to the target indices, i.e., rows of the neighbor lists 2D array. The second rank is the
1029 //! ordinal of the tangent direction (spatial dimensions-1 are tangent, last one is normal), and the third
1030 //! rank is indices into the spatial dimension.
1031 template<typename view_type>
1032 void setTangentBundle(view_type tangent_directions) {
1033 // accept input from user as a rank 3 tensor
1034 // but convert data to a rank 2 tensor with the last rank of dimension = _dimensions x _dimensions
1035 // this allows for nonstrided views on the device later
1036
1037 // allocate memory on device
1038 _T = decltype(_T)("device tangent directions", _pc._target_coordinates.extent(0)*_dimensions*_dimensions);
1039
1040 compadre_assert_release( (std::is_same<decltype(_T)::memory_space, typename view_type::memory_space>::value) &&
1041 "Memory space does not match between _T and tangent_directions");
1042
1043 auto this_dimensions = _dimensions;
1044 auto this_T = _T;
1045 // rearrange data on device from data given on host
1046 Kokkos::parallel_for("copy tangent vectors", Kokkos::RangePolicy<device_execution_space>(0, _pc._target_coordinates.extent(0)), KOKKOS_LAMBDA(const int i) {
1047 scratch_matrix_right_type T(this_T.data() + i*this_dimensions*this_dimensions, this_dimensions, this_dimensions);
1048 for (int j=0; j<this_dimensions; ++j) {
1049 for (int k=0; k<this_dimensions; ++k) {
1050 T(j,k) = tangent_directions(i, j, k);
1051 }
1052 }
1053 });
1055
1056 // copy data from device back to host in rearranged format
1057 _host_T = Kokkos::create_mirror_view(_T);
1058 Kokkos::deep_copy(_host_T, _T);
1059 this->resetCoefficientData();
1060 }
1061
1062 //! (OPTIONAL)
1063 //! Sets outward normal direction. For manifolds this may be used for orienting surface. It is also accessible
1064 //! for sampling operators that require a normal direction.
1065 template<typename view_type>
1066 void setReferenceOutwardNormalDirection(view_type outward_normal_directions, bool use_to_orient_surface = true) {
1067 // accept input from user as a rank 2 tensor
1068
1069 // allocate memory on device
1070 _ref_N = decltype(_ref_N)("device normal directions", _pc._target_coordinates.extent(0)*_dimensions);
1071 // to assist LAMBDA capture
1072 auto this_ref_N = this->_ref_N;
1073 auto this_dimensions = this->_dimensions;
1074
1075 // rearrange data on device from data given on host
1076 Kokkos::parallel_for("copy normal vectors", Kokkos::RangePolicy<device_execution_space>(0, _pc._target_coordinates.extent(0)), KOKKOS_LAMBDA(const int i) {
1077 for (int j=0; j<this_dimensions; ++j) {
1078 this_ref_N(i*this_dimensions + j) = outward_normal_directions(i, j);
1079 }
1080 });
1081 Kokkos::fence();
1084
1085 // copy data from device back to host in rearranged format
1086 _host_ref_N = Kokkos::create_mirror_view(_ref_N);
1087 Kokkos::deep_copy(_host_ref_N, _ref_N);
1088 this->resetCoefficientData();
1089 }
1090
1091 //! (OPTIONAL)
1092 //! Sets extra data to be used by sampling functionals in certain instances.
1093 template<typename view_type>
1094 void setSourceExtraData(view_type extra_data) {
1095
1096 // allocate memory on device
1097 _source_extra_data = decltype(_source_extra_data)("device source extra data", extra_data.extent(0), extra_data.extent(1));
1098
1099 auto host_extra_data = Kokkos::create_mirror_view(_source_extra_data);
1100 Kokkos::deep_copy(host_extra_data, extra_data);
1101 // copy data from host to device
1102 Kokkos::deep_copy(_source_extra_data, host_extra_data);
1103 this->resetCoefficientData();
1104 }
1105
1106 //! (OPTIONAL)
1107 //! Sets extra data to be used by sampling functionals in certain instances. (device)
1108 template<typename view_type>
1109 void setSourceExtraData(decltype(_source_extra_data) extra_data) {
1110 // allocate memory on device
1111 _source_extra_data = extra_data;
1112 this->resetCoefficientData();
1113 }
1114
1115 //! (OPTIONAL)
1116 //! Sets extra data to be used by target operations in certain instances.
1117 template<typename view_type>
1118 void setTargetExtraData(view_type extra_data) {
1119
1120 // allocate memory on device
1121 _target_extra_data = decltype(_target_extra_data)("device target extra data", extra_data.extent(0), extra_data.extent(1));
1122
1123 auto host_extra_data = Kokkos::create_mirror_view(_target_extra_data);
1124 Kokkos::deep_copy(host_extra_data, extra_data);
1125 // copy data from host to device
1126 Kokkos::deep_copy(_target_extra_data, host_extra_data);
1127 this->resetCoefficientData();
1128 }
1129
1130 //! (OPTIONAL)
1131 //! Sets extra data to be used by target operations in certain instances. (device)
1132 template<typename view_type>
1133 void setTargetExtraData(decltype(_target_extra_data) extra_data) {
1134 // allocate memory on device
1135 _target_extra_data = extra_data;
1136 this->resetCoefficientData();
1137 }
1138
1139 //! Set dense solver type type
1141 _dense_solver_type = dst;
1142 this->resetCoefficientData();
1143 }
1144
1145 //! Set constraint type
1147 _constraint_type = ct;
1148 this->resetCoefficientData();
1149 }
1150
1151 //! Type for weighting kernel for GMLS problem
1152 void setWeightingType( const std::string &wt) {
1153 std::string wt_to_lower = wt;
1154 transform(wt_to_lower.begin(), wt_to_lower.end(), wt_to_lower.begin(), ::tolower);
1155 if (wt_to_lower == "power") {
1157 } else if (wt_to_lower == "gaussian") {
1159 } else if (wt_to_lower == "cubicspline") {
1161 } else if (wt_to_lower == "cosine") {
1163 } else if (wt_to_lower == "sigmoid") {
1165 } else {
1166 // Power is default
1168 }
1169 this->resetCoefficientData();
1170 }
1171
1172 //! Type for weighting kernel for GMLS problem
1174 _weighting_type = wt;
1175 this->resetCoefficientData();
1176 }
1177
1178 //! Type for weighting kernel for curvature
1179 void setCurvatureWeightingType( const std::string &wt) {
1180 std::string wt_to_lower = wt;
1181 transform(wt_to_lower.begin(), wt_to_lower.end(), wt_to_lower.begin(), ::tolower);
1182 if (wt_to_lower == "power") {
1184 } else if (wt_to_lower == "gaussian") {
1186 } else if (wt_to_lower == "cubicspline") {
1188 } else if (wt_to_lower == "cosine") {
1190 } else if (wt_to_lower == "sigmoid") {
1192 } else {
1193 // Power is default
1195 }
1196 this->resetCoefficientData();
1197 }
1198
1199 //! Type for weighting kernel for curvature
1202 this->resetCoefficientData();
1203 }
1204
1205 //! Sets basis order to be used when reconstructing any function
1206 void setPolynomialOrder(const int poly_order) {
1207 compadre_assert_release(poly_order<11 && "Unsupported polynomial order (>=11).");
1208 _poly_order = poly_order;
1209 _NP = this->getNP(_poly_order, _dimensions, _reconstruction_space);
1210 this->resetCoefficientData();
1211 }
1212
1213 //! Sets basis order to be used when reconstruction curvature
1214 void setCurvaturePolynomialOrder(const int curvature_poly_order) {
1215 compadre_assert_release(curvature_poly_order<11 && "Unsupported curvature polynomial order (>=11).");
1216 _curvature_poly_order = curvature_poly_order;
1217 this->resetCoefficientData();
1218 }
1219
1220 //! Parameter for weighting kernel for GMLS problem
1221 //! index = 0 sets p paramater for weighting kernel
1222 //! index = 1 sets n paramater for weighting kernel
1223 void setWeightingParameter(int wp, int index = 0) {
1224 if (index==1) {
1225 _weighting_n = wp;
1226 } else {
1227 _weighting_p = wp;
1228 }
1229 this->resetCoefficientData();
1230 }
1231
1232 //! Parameter for weighting kernel for curvature
1233 //! index = 0 sets p paramater for weighting kernel
1234 //! index = 1 sets n paramater for weighting kernel
1235 void setCurvatureWeightingParameter(int wp, int index = 0) {
1236 if (index==1) {
1238 } else {
1240 }
1241 this->resetCoefficientData();
1242 }
1243
1244 //! Number quadrature points to use
1247 this->resetCoefficientData();
1248 }
1249
1250 //! Dimensions of quadrature points to use
1253 this->resetCoefficientData();
1254 }
1255
1256 //! Type of quadrature points
1257 void setQuadratureType(std::string quadrature_type) {
1258 _quadrature_type = quadrature_type;
1259 this->resetCoefficientData();
1260 }
1261
1262 //! Adds a target to the vector of target functional to be applied to the reconstruction
1264 _h_ss.addTargets(lro);
1265 this->resetCoefficientData();
1266 }
1267
1268 //! Adds a vector of target functionals to the vector of target functionals already to be applied to the reconstruction
1269 void addTargets(std::vector<TargetOperation> lro) {
1270 _h_ss.addTargets(lro);
1271 this->resetCoefficientData();
1272 }
1273
1274 //! Empties the vector of target functionals to apply to the reconstruction
1276 _h_ss.clearTargets();
1277 this->resetCoefficientData();
1278 }
1279
1280 //! Verify whether _pc is valid
1281 bool verifyPointConnections(bool assert_valid = false) {
1282 bool result = (_pc._target_coordinates.extent_int(0)==_pc._nla.getNumberOfTargets());
1283 compadre_assert_release((!assert_valid || result) &&
1284 "Target coordinates and neighbor lists have different size.");
1285
1286 result &= (_pc._source_coordinates.extent(1)==_pc._target_coordinates.extent(1));
1287 compadre_assert_release((!assert_valid || result) &&
1288 "Source coordinates and target coordinates have different dimensions.");
1289
1290 result &= (_pc._source_coordinates.extent(0)>0||_pc._target_coordinates.extent(0)==0);
1291 compadre_assert_release((!assert_valid || result) &&
1292 "Source coordinates not set in GMLS class before calling generatePolynomialCoefficients.");
1293 return result;
1294 }
1295
1296 //! Verify whether _additional_pc is valid
1297 bool verifyAdditionalPointConnections(bool assert_valid = false) {
1299 compadre_assert_release((!assert_valid || result) &&
1300 "Target coordinates and additional evaluation indices have different size.");
1301
1302 result &= (_pc._target_coordinates.extent(0)==_additional_pc._target_coordinates.extent(0));
1303 compadre_assert_release((!assert_valid || result) &&
1304 "Target coordinates and additional evaluation indices have different size.");
1305
1306 if (_additional_pc._source_coordinates.extent(0)>0) {
1308 compadre_assert_release((!assert_valid || result) &&
1309 "Target coordinates and additional evaluation coordinates have different dimensions.");
1310 }
1311 return result;
1312 }
1313
1314 /*! \brief Generates polynomial coefficients by setting up and solving least squares problems
1315 Sets up the batch of GMLS problems to be solved for. Provides alpha values
1316 that can later be contracted against data or degrees of freedom to form a
1317 global linear system.
1318
1319 \param number_of_batches [in] - how many batches to break up the total workload into (for storage)
1320 \param keep_coefficients [in] - whether to store (P^T W P)^-1 * P^T * W
1321 \param clear_cache [in] - clear whatever temporary memory was used for calculations that
1322 \a keep_coefficients hasn't indicated needs to be saved
1323
1324 \a clear_cache should be \a false to be used for debugging as it will leave all data structures used to
1325 compute the solution intact. If \a clear_cache is set \a true and \a keep_coefficients is \a true, it will
1326 allow the polynomial coefficients to persist after calculation.
1327 */
1328 void generatePolynomialCoefficients(const int number_of_batches = 1, const bool keep_coefficients = false,
1329 const bool clear_cache = true);
1330
1331 /*! \brief Meant to calculate target operations and apply the evaluations to the previously
1332 constructed polynomial coefficients. But now that is inside of generatePolynomialCoefficients because
1333 it must be to handle number_of_batches>1. Effectively, this just calls generatePolynomialCoefficients.
1334
1335 \param number_of_batches [in] - how many batches to break up the total workload into (for storage)
1336 \param keep_coefficients [in] - whether to store (P^T W P)^-1 * P^T * W
1337 \param clear_cache [in] - clear whatever temporary memory was used for calculations that
1338 \a keep_coefficients hasn't indicated needs to be saved
1339
1340 \a clear_cache should be \a false to be used for debugging as it will leave all data structures used to
1341 compute the solution intact. If \a clear_cache is set \a true and \a keep_coefficients is \a true, it will
1342 allow the polynomial coefficients to persist after calculation.
1343 */
1344 void generateAlphas(const int number_of_batches = 1, const bool keep_coefficients = false,
1345 const bool clear_cache = true);
1346
1347///@}
1348
1349
1350}; // GMLS Class
1351} // Compadre
1352
1353#endif
1354
1355
Kokkos::View< int *, host_execution_space > host_managed_local_index_type
#define compadre_kernel_assert_release(condition)
compadre_kernel_assert_release is similar to compadre_assert_release, but is a call on the device,...
#define compadre_assert_debug(condition)
compadre_assert_debug is used for assertions that are checked in loops, as these significantly impact...
Kokkos::View< double **, layout_right, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_matrix_right_type
#define compadre_assert_release(condition)
compadre_assert_release is used for assertions that should always be checked, but generally are not e...
Lightweight Evaluator Helper This class is a lightweight wrapper for extracting and applying all rele...
Generalized Moving Least Squares (GMLS)
SolutionSet< device_memory_space > _d_ss
int _sampling_multiplier
actual dimension of the sampling functional e.g.
neighbor_lists_type * getNeighborLists() const
Get neighbor list accessor.
std::enable_if< view_type::rank==2, void >::type setAuxiliaryEvaluationIndicesLists(view_type additional_evaluation_indices)
(OPTIONAL) Sets the additional target evaluation indices list information.
void setWindowSizes(decltype(_epsilons) epsilons)
Sets window sizes, also called the support of the kernel (device)
host_managed_local_index_type getPolynomialCoefficientsMemorySize() const
Returns 2D array size in memory on which coefficients are stored.
decltype(_RHS) getFullPolynomialCoefficientsBasis() const
Get a view (device) of all polynomial coefficients basis.
void addTargets(TargetOperation lro)
Adds a target to the vector of target functional to be applied to the reconstruction.
double getPreStencilWeight(SamplingFunctional sro, const int target_index, const int neighbor_index, bool for_target, const int output_component=0, const int input_component=0) const
Returns a stencil to transform data from its existing state into the input expected for some sampling...
void setDimensionOfQuadraturePoints(int dim)
Dimensions of quadrature points to use.
void setCurvatureWeightingType(const WeightingFunctionType wt)
Type for weighting kernel for curvature.
int _weighting_p
first parameter to be used for weighting kernel
std::string _quadrature_type
quadrature rule type
void setSourceExtraData(decltype(_source_extra_data) extra_data)
(OPTIONAL) Sets extra data to be used by sampling functionals in certain instances.
Kokkos::View< double * > _T
Rank 3 tensor for high order approximation of tangent vectors for all problems.
point_connections_type _additional_pc
(OPTIONAL) connections between additional points and neighbors
void setWeightingType(const WeightingFunctionType wt)
Type for weighting kernel for GMLS problem.
std::enable_if< view_type::rank==1 &&std::is_same< neighbor_lists_type::internal_view_type, view_type >::value==1, void >::type setNeighborLists(view_type neighbor_lists, view_type number_of_neighbors_list)
Sets neighbor list information from compressed row neighborhood lists data (if same view_type).
int getLocalDimensions() const
Local dimension of the GMLS problem (less than global dimension if on a manifold),...
void addTargets(std::vector< TargetOperation > lro)
Adds a vector of target functionals to the vector of target functionals already to be applied to the ...
bool _reference_outward_normal_direction_provided
whether or not the reference outward normal directions were provided by the user.
int _basis_multiplier
dimension of the reconstructed function e.g.
ReconstructionSpace _reconstruction_space
reconstruction space for GMLS problems, set at GMLS class instantiation
int _reconstruction_space_rank
actual rank of reconstruction basis
void setProblemData(view_type_1 cr_neighbor_lists, view_type_1 number_of_neighbors_list, view_type_2 source_coordinates, view_type_3 target_coordinates, view_type_4 epsilons)
Sets basic problem data (neighbor lists data, number of neighbors list, source coordinates,...
Kokkos::View< double * > _manifold_metric_tensor_inverse
metric tensor inverse for all problems
int getOrderOfQuadraturePoints() const
Order of quadrature points.
const GMLSSolutionData extractSolutionData() const
Get GMLS solution data.
void setReferenceOutwardNormalDirection(view_type outward_normal_directions, bool use_to_orient_surface=true)
(OPTIONAL) Sets outward normal direction.
int _order_of_quadrature_points
order of exact polynomial integration for quadrature rule
void setCurvatureWeightingParameter(int wp, int index=0)
Parameter for weighting kernel for curvature index = 0 sets p paramater for weighting kernel index = ...
int getDimensions() const
Dimension of the GMLS problem, set only at class instantiation.
int _global_dimensions
spatial dimension of the points, set at class instantiation only
void setSourceExtraData(view_type extra_data)
(OPTIONAL) Sets extra data to be used by sampling functionals in certain instances.
void generatePolynomialCoefficients(const int number_of_batches=1, const bool keep_coefficients=false, const bool clear_cache=true)
Generates polynomial coefficients by setting up and solving least squares problems Sets up the batch ...
ReconstructionSpace getReconstructionSpace() const
Get the reconstruction space specified at instantiation.
Kokkos::View< double * > _epsilons
h supports determined through neighbor search (device)
PointConnections< Kokkos::View< double **, layout_right >, Kokkos::View< double **, layout_right >, NeighborLists< Kokkos::View< int * > > > point_connections_type
GMLS(ReconstructionSpace reconstruction_space, SamplingFunctional dual_sampling_strategy, const int poly_order, const int dimensions=3, const std::string dense_solver_type=std::string("QR"), const std::string problem_type=std::string("STANDARD"), const std::string constraint_type=std::string("NO_CONSTRAINT"), const int manifold_curvature_poly_order=2)
Constructor for the case when nonstandard sampling functionals or reconstruction spaces are to be use...
WeightingFunctionType _weighting_type
weighting kernel type for GMLS
ProblemType _problem_type
problem type for GMLS problem, can also be set to STANDARD for normal or MANIFOLD for manifold proble...
int getPolynomialCoefficientsSize() const
Returns size of the basis used in instance's polynomial reconstruction.
neighbor_lists_type * getAdditionalEvaluationIndices() const
(OPTIONAL) Get additional evaluation sites neighbor list-like accessor
int getDimensionOfQuadraturePoints() const
Dimensions of quadrature points.
int _data_sampling_multiplier
effective dimension of the data sampling functional e.g.
std::string getQuadratureType() const
Type of quadrature points.
bool containsValidAlphas() const
Check if GMLS solution set contains valid alpha values (has generateAlphas been called)
int getPolynomialOrder() const
Get basis order used for reconstruction.
void setSourceSites(coordinates_type source_coordinates)
Sets source coordinate information.
void setTangentBundle(view_type tangent_directions)
(OPTIONAL) Sets orthonormal tangent directions for reconstruction on a manifold.
Kokkos::View< double **, layout_right > _source_extra_data
Extra data available to basis functions (optional)
WeightingFunctionType getManifoldWeightingType() const
Type for weighting kernel for curvature.
DenseSolverType _dense_solver_type
solver type for GMLS problem - can be QR, SVD or LU
Kokkos::View< TargetOperation * > _curvature_support_operations
vector containing target functionals to be applied for curvature
void setAdditionalEvaluationSitesData(view_type_1 additional_evaluation_indices, view_type_2 additional_evaluation_coordinates)
(OPTIONAL) Sets additional evaluation sites for each target site
DenseSolverType getDenseSolverType() const
Get dense solver type.
int getGlobalDimensions() const
Dimension of the GMLS problem's point data (spatial description of points in ambient space),...
void setDenseSolverType(const DenseSolverType dst)
Set dense solver type type.
Kokkos::View< double * > _Z
stores evaluations of targets applied to basis
int _initial_index_for_batch
initial index for current batch
void clearTargets()
Empties the vector of target functionals to apply to the reconstruction.
decltype(_ref_N) * getReferenceNormalDirections()
Get a view (device) of all reference outward normal directions.
Kokkos::View< TargetOperation * >::HostMirror _host_operations
vector containing target functionals to be applied for reconstruction problem (host)
void setAdditionalEvaluationSitesData(view_type_1 cr_additional_evaluation_indices, view_type_1 number_of_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...
Kokkos::View< constdouble *****, layout_right >::HostMirror _host_prestencil_weights
generated weights for nontraditional samples required to transform data into expected sampling functi...
int _NP
dimension of basis for polynomial reconstruction
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...
decltype(_pc) * getPointConnections()
Get a view (device) of all point connection info.
void setTargetExtraData(decltype(_target_extra_data) extra_data)
(OPTIONAL) Sets extra data to be used by target operations in certain instances.
void setCurvaturePolynomialOrder(const int curvature_poly_order)
Sets basis order to be used when reconstruction curvature.
int getWeightingParameter(const int index=0) const
Get parameter for weighting kernel for GMLS problem.
int _poly_order
order of basis for polynomial reconstruction
void setAuxiliaryEvaluationCoordinates(view_type evaluation_coordinates)
(OPTIONAL) Sets additional points for evaluation of target operation on polynomial reconstruction.
void setOrderOfQuadraturePoints(int order)
Number quadrature points to use.
bool verifyPointConnections(bool assert_valid=false)
Verify whether _pc is valid.
decltype(_d_ss) * getSolutionSetDevice(bool alpha_validity_check=true)
Get solution set on device.
Kokkos::View< double * > _manifold_curvature_coefficients
curvature polynomial coefficients for all problems
SamplingFunctional _data_sampling_functional
generally the same as _polynomial_sampling_functional, but can differ if specified at
SolutionSet< host_memory_space > _h_ss
Solution Set (contains all alpha values from solution and alpha layout methods)
ConstraintType _constraint_type
constraint type for GMLS problem
Kokkos::View< double **, layout_right > _target_extra_data
Extra data available to target operations (optional)
WeightingFunctionType getWeightingType() const
Type for weighting kernel for GMLS problem.
double getReferenceNormalDirection(const int target_index, const int component) const
Get component of tangent or normal directions for manifold problems.
decltype(_h_ss) * getSolutionSetHost(bool alpha_validity_check=true)
Get solution set on host.
static ProblemType parseProblemType(const std::string &problem_type)
Parses a string to determine problem type.
static ConstraintType parseConstraintType(const std::string &constraint_type)
Parses a string to determine constraint type.
int _dimensions
dimension of the problem, set at class instantiation only
int _curvature_poly_order
order of basis for curvature reconstruction
GMLS(const int poly_order, const int dimensions=3, const std::string dense_solver_type=std::string("QR"), const std::string problem_type=std::string("STANDARD"), const std::string constraint_type=std::string("NO_CONSTRAINT"), const int manifold_curvature_poly_order=2)
Constructor for the case when the data sampling functional does not match the polynomial sampling fun...
int _curvature_weighting_p
first parameter to be used for weighting kernel for curvature
decltype(_epsilons) * getWindowSizes()
Get a view (device) of all window sizes.
bool _use_reference_outward_normal_direction_provided_to_orient_surface
whether or not to use reference outward normal directions to orient the surface in a manifold problem...
void setQuadratureType(std::string quadrature_type)
Type of quadrature points.
std::enable_if< view_type::rank==2, void >::type setNeighborLists(view_type neighbor_lists)
Sets neighbor list information.
int getManifoldWeightingParameter(const int index=0) const
Get parameter for weighting kernel for curvature.
std::enable_if< view_type::rank==1 &&std::is_same< neighbor_lists_type::internal_view_type, view_type >::value==0, void >::type setNeighborLists(view_type neighbor_lists, view_type number_of_neighbors_list)
Sets neighbor list information from compressed row neighborhood lists data (if different view_type).
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)
ParallelManager _pm
determines scratch level spaces and is used to call kernels
WeightingFunctionType _curvature_weighting_type
weighting kernel type for curvature problem
void setWeightingType(const std::string &wt)
Type for weighting kernel for GMLS problem.
int _weighting_n
second parameter to be used for weighting kernel
int getCurvaturePolynomialOrder() const
Get basis order used for curvature reconstruction.
void setConstraintType(const ConstraintType ct)
Set constraint type.
Kokkos::View< TargetOperation * > _operations
vector containing target functionals to be applied for reconstruction problem (device)
bool _orthonormal_tangent_space_provided
whether or not the orthonormal tangent directions were provided by the user.
Kokkos::View< double * > _RHS
sqrt(w)*Identity matrix for all problems, later holds polynomial coefficients for all problems
SamplingFunctional _polynomial_sampling_functional
polynomial sampling functional used to construct P matrix, set at GMLS class instantiation NOTE: ca...
std::enable_if< view_type::rank==1 &&std::is_same< neighbor_lists_type::internal_view_type, view_type >::value==0, void >::type setAuxiliaryEvaluationIndicesLists(view_type additional_evaluation_indices, view_type number_of_neighbors_list)
(OPTIONAL) Sets the additional target evaluation indices list information from compressed row format ...
static KOKKOS_INLINE_FUNCTION int getNN(const int m, const int dimension=3, const ReconstructionSpace r_space=ReconstructionSpace::ScalarTaylorPolynomial)
Returns number of neighbors needed for unisolvency for a given basis order and dimension.
Kokkos::View< double * > _P
P*sqrt(w) matrix for all problems.
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 setPolynomialOrder(const int poly_order)
Sets basis order to be used when reconstructing any function.
NeighborLists< Kokkos::View< int * > > neighbor_lists_type
GMLS(ReconstructionSpace reconstruction_space, const SamplingFunctional polynomial_sampling_strategy, const SamplingFunctional data_sampling_strategy, const int poly_order, const int dimensions, const DenseSolverType dense_solver_type, const ProblemType problem_type, const ConstraintType constraint_type, const int manifold_curvature_poly_order)
Maximal constructor, not intended for users.
SamplingFunctional getDataSamplingFunctional() const
Get the data sampling functional specified at instantiation (often the same as the polynomial samplin...
decltype(_T) * getTangentDirections()
Get a view (device) of all tangent direction bundles.
void resetCoefficientData()
const GMLSBasisData extractBasisData() const
Get GMLS basis data.
void setTargetSites(view_type target_coordinates)
Sets target coordinate information. Rows of this 2D-array should correspond to rows of the neighbor l...
static DenseSolverType parseSolverType(const std::string &dense_solver_type)
Parses a string to determine solver type.
Kokkos::View< double **, layout_right > coordinates_type
point_connections_type _pc
connections between points and neighbors
SamplingFunctional getPolynomialSamplingFunctional() const
Get the polynomial sampling functional specified at instantiation.
void setSourceSites(view_type source_coordinates)
Sets source coordinate information.
Kokkos::View< double *****, layout_right > _prestencil_weights
generated weights for nontraditional samples required to transform data into expected sampling functi...
Kokkos::View< double * >::HostMirror _host_ref_N
reference outward normal vectors information (host)
std::enable_if< view_type::rank==1 &&std::is_same< neighbor_lists_type::internal_view_type, view_type >::value==1, void >::type setAuxiliaryEvaluationIndicesLists(view_type additional_evaluation_indices, view_type number_of_neighbors_list)
(OPTIONAL) Sets the additional target evaluation indices list information from compressed row format ...
decltype(_additional_pc) * getAdditionalPointConnections()
(OPTIONAL) Get a view (device) of all additional evaluation point connection info
ProblemType getProblemType() const
Get problem type.
ConstraintType getConstraintType() const
Get constraint type.
void setCurvatureWeightingType(const std::string &wt)
Type for weighting kernel for curvature.
Kokkos::View< double * > _w
contains weights for all problems
host_managed_local_index_type getPolynomialCoefficientsDomainRangeSize() const
Returns (size of the basis used in instance's polynomial reconstruction) x (data input dimension)
void setWindowSizes(view_type epsilons)
Sets window sizes, also called the support of the kernel.
int _curvature_weighting_n
second parameter to be used for weighting kernel for curvature
bool _entire_batch_computed_at_once
whether entire calculation was computed at once the alternative is that it was broken up over many sm...
static KOKKOS_INLINE_FUNCTION double Wab(const double r, const double h, const WeightingFunctionType &weighting_type, const int p, const int n)
Evaluates the weighting kernel.
int _local_dimensions
dimension of the problem, set at class instantiation only. For manifolds, generally _global_dimension...
GMLS(ReconstructionSpace reconstruction_space, const SamplingFunctional polynomial_sampling_strategy, const SamplingFunctional data_sampling_strategy, const int poly_order, const int dimensions=3, const std::string dense_solver_type=std::string("QR"), const std::string problem_type=std::string("STANDARD"), const std::string constraint_type=std::string("NO_CONSTRAINT"), const int manifold_curvature_poly_order=2)
Maximal constructor, but with string arguments.
Kokkos::View< double * > _manifold_curvature_gradient
_dimension-1 gradient values for curvature for all problems
Kokkos::View< double * >::HostMirror _host_T
tangent vectors information (host)
bool verifyAdditionalPointConnections(bool assert_valid=false)
Verify whether _additional_pc is valid.
bool _store_PTWP_inv_PTW
whether polynomial coefficients were requested to be stored (in a state not yet applied to data)
int _dimension_of_quadrature_points
dimension of quadrature rule
void setTargetExtraData(view_type extra_data)
(OPTIONAL) Sets extra data to be used by target operations in certain instances.
double getTangentBundle(const int target_index, const int direction, const int component) const
Get component of tangent or normal directions for manifold problems.
void setTargetSites(coordinates_type target_coordinates)
Sets target coordinate information. Rows of this 2D-array should correspond to rows of the neighbor l...
Kokkos::View< double * > _ref_N
Rank 2 tensor for high order approximation of tangent vectors for all problems.
int getNumberOfQuadraturePoints() const
Number of quadrature points.
void setAuxiliaryEvaluationCoordinates(coordinates_type evaluation_coordinates)
(OPTIONAL) Sets additional points for evaluation of target operation on polynomial reconstruction.
Quadrature _qm
manages and calculates quadrature
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...
KOKKOS_INLINE_FUNCTION int getNumberOfQuadraturePoints() const
KOKKOS_INLINE_FUNCTION int getSize(const int degree, const int dimension)
Returns size of basis.
KOKKOS_INLINE_FUNCTION int getSize(const int degree, const int dimension)
Returns size of basis.
KOKKOS_INLINE_FUNCTION int getSize(const int degree, const int dimension)
Returns size of basis.
ProblemType
Problem type, that optionally can handle manifolds.
@ STANDARD
Standard GMLS problem type.
@ MANIFOLD
Solve GMLS problem on a manifold (will use QR or SVD to solve the resultant GMLS problem dependent on...
KOKKOS_INLINE_FUNCTION void getRHSDims(DenseSolverType dense_solver_type, ConstraintType constraint_type, ReconstructionSpace reconstruction_space, const int dimension, const int M, const int N, int &RHS_row, int &RHS_col)
KOKKOS_INLINE_FUNCTION int getActualReconstructionSpaceRank(const int &index)
Number of actual components in the ReconstructionSpace.
ConstraintType
Constraint type.
@ NEUMANN_GRAD_SCALAR
Neumann Gradient Scalar Type.
@ NO_CONSTRAINT
No constraint.
constexpr SamplingFunctional PointSample
Available sampling functionals.
TargetOperation
Available target functionals.
@ GradientOfScalarPointEvaluation
Point evaluation of the gradient of a scalar.
@ DifferentEachNeighbor
Each target applies a different transform for each neighbor.
@ DifferentEachTarget
Each target applies a different data transform, but the same to each neighbor.
DenseSolverType
Dense solver type.
@ LU
LU factorization performed on P^T*W*P matrix.
@ QR
QR+Pivoting factorization performed on P*sqrt(w) matrix.
WeightingFunctionType
Available weighting kernel function types.
NeighborLists< view_type_1d > Convert2DToCompressedRowNeighborLists(view_type_2d neighbor_lists)
Converts 2D neighbor lists to compressed row neighbor lists.
KOKKOS_INLINE_FUNCTION void getPDims(DenseSolverType dense_solver_type, ConstraintType constraint_type, ReconstructionSpace reconstruction_space, const int dimension, const int M, const int N, int &out_row, int &out_col)
constexpr SamplingFunctional ManifoldVectorPointSample
Point evaluations of the entire vector source function (but on a manifold, so it includes a transform...
constexpr SamplingFunctional VectorPointSample
Point evaluations of the entire vector source function.
ReconstructionSpace
Space in which to reconstruct polynomial.
@ DivergenceFreeVectorTaylorPolynomial
Divergence-free vector polynomial basis.
@ BernsteinPolynomial
Bernstein polynomial basis.
@ ScalarTaylorPolynomial
Scalar polynomial basis centered at the target site and scaled by sum of basis powers e....
@ VectorOfScalarClonesTaylorPolynomial
Scalar basis reused as many times as there are components in the vector resulting in a much cheaper p...
NeighborLists assists in accessing entries of compressed row neighborhood lists.
KOKKOS_INLINE_FUNCTION int getMaxNumNeighbors() const
Get the maximum number of neighbors of all targets' neighborhoods (host/device)
KOKKOS_INLINE_FUNCTION int getNumberOfTargets() const
Get number of total targets having neighborhoods (host/device).
Combines NeighborLists with the PointClouds from which it was derived Assumed that memory_space is th...
void setSourceCoordinates(view_type_2 source_coordinates)
Update only source coordinates.
void setNeighborLists(nla_type nla)
Update only target coordinates.
void setTargetCoordinates(view_type_1 target_coordinates)
Update only target coordinates.
device_mirror_source_view_type _source_coordinates
device_mirror_target_view_type _target_coordinates
int transform_type
Describes the SamplingFunction relationship to targets, neighbors.
All vairables and functionality related to the layout and storage of GMLS solutions (alpha values)
void addTargets(TargetOperation lro)
Adds a target to the vector of target functional to be applied to the reconstruction.
void clearTargets()
Empties the vector of target functionals to apply to the reconstruction.
int _max_evaluation_sites_per_target
maximum number of evaluation sites for each target (includes target site)
bool _contains_valid_alphas
whether internal alpha values are valid (set externally on a solve)
void copyAlphas(SolutionSet< other_memory_space > &other)
Copies alphas between two instances of SolutionSet Copying of alphas is intentionally omitted in copy...