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