Compadre  1.3.3
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"
5 #include "Compadre_Typedefs.hpp"
6 
7 #include "Compadre_Misc.hpp"
8 #include "Compadre_Operators.hpp"
11 #include "Compadre_Quadrature.hpp"
15 
16 namespace Compadre {
17 
18 //! Generalized Moving Least Squares (GMLS)
19 /*!
20 * This class sets up a batch of GMLS problems from a given set of neighbor lists, target sites, and source sites.
21 * GMLS requires a target functional, reconstruction space, and sampling functional to be specified.
22 * For a given choice of reconstruction space and sampling functional, multiple targets can be generated with very little
23 * additional computation, which is why this class allows for multiple target functionals to be specified.
24 */
25 class GMLS {
26 protected:
27 
28  // random numbe generator pool
30 
31  // matrices that may be needed for matrix factorization on the device
32  // supports batched matrix factorization dispatch
33 
34  //! contains weights for all problems
35  Kokkos::View<double*> _w;
36 
37  //! P*sqrt(w) matrix for all problems
38  Kokkos::View<double*> _P;
39 
40  //! sqrt(w)*Identity matrix for all problems, later holds polynomial coefficients for all problems
41  Kokkos::View<double*> _RHS;
42 
43  //! Rank 3 tensor for high order approximation of tangent vectors for all problems. First rank is
44  //! for the target index, the second is for the local direction to the manifolds 0..(_dimensions-1)
45  //! are tangent, _dimensions is the normal, and the third is for the spatial dimension (_dimensions)
46  Kokkos::View<double*> _T;
47 
48  //! Rank 2 tensor for high order approximation of tangent vectors for all problems. First rank is
49  //! for the target index, the second is for the spatial dimension (_dimensions)
50  Kokkos::View<double*> _ref_N;
51 
52  //! tangent vectors information (host)
53  Kokkos::View<double*>::HostMirror _host_T;
54 
55  //! reference outward normal vectors information (host)
56  Kokkos::View<double*>::HostMirror _host_ref_N;
57 
58  //! metric tensor inverse for all problems
59  Kokkos::View<double*> _manifold_metric_tensor_inverse;
60 
61  //! curvature polynomial coefficients for all problems
62  Kokkos::View<double*> _manifold_curvature_coefficients;
63 
64  //! _dimension-1 gradient values for curvature for all problems
65  Kokkos::View<double*> _manifold_curvature_gradient;
66 
67  //! Extra data available to basis functions (optional)
68  Kokkos::View<double**, layout_right> _source_extra_data;
69 
70  //! Extra data available to target operations (optional)
71  Kokkos::View<double**, layout_right> _target_extra_data;
72 
73  //! Accessor to get neighbor list data, offset data, and number of neighbors per target
75 
76  //! convenient copy on host of number of neighbors
77  Kokkos::View<int*, host_memory_space> _host_number_of_neighbors_list;
78 
79  //! all coordinates for the source for which _neighbor_lists refers (device)
80  Kokkos::View<double**, layout_right> _source_coordinates;
81 
82  //! coordinates for target sites for reconstruction (device)
83  Kokkos::View<double**, layout_right> _target_coordinates;
84 
85  //! h supports determined through neighbor search (device)
86  Kokkos::View<double*> _epsilons;
87 
88  //! h supports determined through neighbor search (host)
89  Kokkos::View<double*>::HostMirror _host_epsilons;
90 
91  //! generated alpha coefficients (device)
92  Kokkos::View<double*, layout_right> _alphas;
93 
94  //! generated alpha coefficients (host)
95  Kokkos::View<const double*, layout_right>::HostMirror _host_alphas;
96 
97  //! generated weights for nontraditional samples required to transform data into expected sampling
98  //! functional form (device).
99  Kokkos::View<double*****, layout_right> _prestencil_weights;
100 
101  //! generated weights for nontraditional samples required to transform data into expected sampling
102  //! functional form (host)
103  Kokkos::View<const double*****, layout_right>::HostMirror _host_prestencil_weights;
104 
105  //! (OPTIONAL) user provided additional coordinates for target operation evaluation (device)
106  Kokkos::View<double**, layout_right> _additional_evaluation_coordinates;
107 
108  //! (OPTIONAL) contains indices of entries in the _additional_evaluation_coordinates view (device)
109  Kokkos::View<int**, layout_right> _additional_evaluation_indices;
110 
111  //! (OPTIONAL) contains indices of entries in the _additional_evaluation_coordinates view (host)
112  Kokkos::View<int**, layout_right>::HostMirror _host_additional_evaluation_indices;
113 
114  //! (OPTIONAL) contains the # of additional coordinate indices for each target
116 
117 
118  //! order of basis for polynomial reconstruction
120 
121  //! order of basis for curvature reconstruction
123 
124  //! dimension of basis for polynomial reconstruction
125  int _NP;
126 
127  //! spatial dimension of the points, set at class instantiation only
129 
130  //! dimension of the problem, set at class instantiation only. For manifolds, generally _global_dimensions-1
132 
133  //! dimension of the problem, set at class instantiation only
135 
136  //! reconstruction space for GMLS problems, set at GMLS class instantiation
138 
139  //! actual rank of reconstruction basis
141 
142  //! solver type for GMLS problem - can be QR, SVD or LU
144 
145  //! problem type for GMLS problem, can also be set to STANDARD for normal or MANIFOLD for manifold problems
147 
148  //! constraint type for GMLS problem
150 
151  //! polynomial sampling functional used to construct P matrix, set at GMLS class instantiation
153 
154  //! generally the same as _polynomial_sampling_functional, but can differ if specified at
155  //! GMLS class instantiation
157 
158  //! vector containing target functionals to be applied for curvature
159  Kokkos::View<TargetOperation*> _curvature_support_operations;
160 
161  //! vector containing target functionals to be applied for reconstruction problem (device)
162  Kokkos::View<TargetOperation*> _operations;
163 
164  //! vector containing target functionals to be applied for reconstruction problem (host)
165  Kokkos::View<TargetOperation*>::HostMirror _host_operations;
166 
167 
168 
169  //! weighting kernel type for GMLS
171 
172  //! weighting kernel type for curvature problem
174 
175  //! power to be used for weighting kernel
177 
178  //! power to be used for weighting kernel for curvature
180 
181  //! dimension of the reconstructed function
182  //! e.g. reconstruction of vector on a 2D manifold in 3D would have _basis_multiplier of 2
184 
185  //! actual dimension of the sampling functional
186  //! e.g. reconstruction of vector on a 2D manifold in 3D would have _basis_multiplier of 2
187  //! e.g. in 3D, a scalar will be 1, a vector will be 3, and a vector of reused scalars will be 1
189 
190  //! effective dimension of the data sampling functional
191  //! e.g. in 3D, a scalar will be 1, a vector will be 3, and a vector of reused scalars will be 3
193 
194  //! whether or not operator to be inverted for GMLS problem has a nontrivial nullspace (requiring SVD)
196 
197  //! whether or not the orthonormal tangent directions were provided by the user. If they are not,
198  //! then for the case of calculations on manifolds, a GMLS approximation of the tangent space will
199  //! be made and stored for use.
201 
202  //! whether or not the reference outward normal directions were provided by the user.
204 
205  //! whether or not to use reference outward normal directions to orient the surface in a manifold problem.
207 
208  //! whether entire calculation was computed at once
209  //! the alternative is that it was broken up over many smaller groups, in which case
210  //! this is false, and so the _RHS matrix can not be stored or requested
212 
213  //! whether polynomial coefficients were requested to be stored (in a state not yet applied to data)
215 
216  //! initial index for current batch
218 
219  //! maximum number of neighbors over all target sites
221 
222  //! maximum number of evaluation sites for each target (includes target site)
224 
225  //! vector of user requested target operations
226  std::vector<TargetOperation> _lro;
227 
228  //! vector containing a mapping from a target functionals enum value to the its place in the list
229  //! of target functionals to be applied
230  std::vector<int> _lro_lookup;
231 
232  //! index for where this operation begins the for _alpha coefficients (device)
233  Kokkos::View<int*> _lro_total_offsets;
234 
235  //! index for where this operation begins the for _alpha coefficients (host)
236  Kokkos::View<int*>::HostMirror _host_lro_total_offsets;
237 
238  //! dimensions ^ rank of tensor of output for each target functional (device)
239  Kokkos::View<int*> _lro_output_tile_size;
240 
241  //! dimensions ^ rank of tensor of output for each target functional (host)
242  Kokkos::View<int*>::HostMirror _host_lro_output_tile_size;
243 
244  //! dimensions ^ rank of tensor of output for each sampling functional (device)
245  Kokkos::View<int*> _lro_input_tile_size;
246 
247  //! dimensions ^ rank of tensor of output for each sampling functional (host)
248  Kokkos::View<int*>::HostMirror _host_lro_input_tile_size;
249 
250  //! tensor rank of target functional (device)
251  Kokkos::View<int*> _lro_output_tensor_rank;
252 
253  //! tensor rank of target functional (host)
254  Kokkos::View<int*>::HostMirror _host_lro_output_tensor_rank;
255 
256  //! tensor rank of sampling functional (device)
257  Kokkos::View<int*> _lro_input_tensor_rank;
258 
259  //! tensor rank of sampling functional (host)
260  Kokkos::View<int*>::HostMirror _host_lro_input_tensor_rank;
261 
262  //! used for sizing P_target_row and the _alphas view
264 
265  //! additional alpha coefficients due to constraints
267 
268  //! determines scratch level spaces and is used to call kernels
270 
271  //! order of exact polynomial integration for quadrature rule
273 
274  //! dimension of quadrature rule
276 
277  //! quadrature rule type
278  std::string _quadrature_type;
279 
280  //! manages and calculates quadrature
282 
283 
284 /** @name Private Modifiers
285  * Private function because information lives on the device
286  */
287 ///@{
288 
289  /*! \brief Evaluates the polynomial basis under a particular sampling function. Generally used to fill a row of P.
290  \param delta [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large as the _basis_multipler*the dimension of the polynomial basis.
291  \param thread_workspace [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large as the _poly_order*the spatial dimension of the polynomial basis.
292  \param target_index [in] - target number
293  \param neighbor_index [in] - index of neighbor for this target with respect to local numbering [0,...,number of neighbors for target]
294  \param alpha [in] - double to determine convex combination of target and neighbor site at which to evaluate polynomials. (1-alpha)*neighbor + alpha*target
295  \param dimension [in] - spatial dimension of basis to evaluate. e.g. dimension two basis of order one is 1, x, y, whereas for dimension 3 it is 1, x, y, z
296  \param poly_order [in] - polynomial basis degree
297  \param specific_order_only [in] - boolean for only evaluating one degree of polynomial when true
298  \param V [in] - orthonormal basis matrix size _dimensions * _dimensions whose first _dimensions-1 columns are an approximation of the tangent plane
299  \param reconstruction_space [in] - space of polynomial that a sampling functional is to evaluate
300  \param sampling_strategy [in] - sampling functional specification
301  \param additional_evaluation_local_index [in] - local index for evaluation sites
302  */
303  KOKKOS_INLINE_FUNCTION
304  void calcPij(const member_type& teamMember, double* delta, double* thread_workspace, const int target_index, int neighbor_index, const double alpha, const int dimension, const int poly_order, bool specific_order_only = false, const scratch_matrix_right_type* V = NULL, const ReconstructionSpace reconstruction_space = ReconstructionSpace::ScalarTaylorPolynomial, const SamplingFunctional sampling_strategy = PointSample, const int additional_evaluation_local_index = 0) const;
305 
306  /*! \brief Evaluates the gradient of a polynomial basis under the Dirac Delta (pointwise) sampling function.
307  \param delta [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large is the _basis_multipler*the dimension of the polynomial basis.
308  \param thread_workspace [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large as the _poly_order*the spatial dimension of the polynomial basis.
309  \param target_index [in] - target number
310  \param neighbor_index [in] - index of neighbor for this target with respect to local numbering [0,...,number of neighbors for target]
311  \param alpha [in] - double to determine convex combination of target and neighbor site at which to evaluate polynomials. (1-alpha)*neighbor + alpha*target
312  \param partial_direction [in] - direction that partial is taken with respect to, e.g. 0 is x direction, 1 is y direction
313  \param dimension [in] - spatial dimension of basis to evaluate. e.g. dimension two basis of order one is 1, x, y, whereas for dimension 3 it is 1, x, y, z
314  \param poly_order [in] - polynomial basis degree
315  \param specific_order_only [in] - boolean for only evaluating one degree of polynomial when true
316  \param V [in] - orthonormal basis matrix size _dimensions * _dimensions whose first _dimensions-1 columns are an approximation of the tangent plane
317  \param reconstruction_space [in] - space of polynomial that a sampling functional is to evaluate
318  \param sampling_strategy [in] - sampling functional specification
319  \param additional_evaluation_local_index [in] - local index for evaluation sites
320  */
321  KOKKOS_INLINE_FUNCTION
322  void calcGradientPij(const member_type& teamMember, double* delta, double* thread_workspace, const int target_index, int neighbor_index, const double alpha, const int partial_direction, const int dimension, const int poly_order, bool specific_order_only, const scratch_matrix_right_type* V, const ReconstructionSpace reconstruction_space, const SamplingFunctional sampling_strategy, const int additional_evaluation_local_index = 0) const;
323 
324  /*! \brief Evaluates the Hessian of a polynomial basis under the Dirac Delta (pointwise) sampling function.
325  \param delta [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large is the _basis_multipler*the dimension of the polynomial basis.
326  \param thread_workspace [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large as the _poly_order*the spatial dimension of the polynomial basis.
327  \param target_index [in] - target number
328  \param neighbor_index [in] - index of neighbor for this target with respect to local numbering [0,...,number of neighbors for target]
329  \param alpha [in] - double to determine convex combination of target and neighbor site at which to evaluate polynomials. (1-alpha)*neighbor + alpha*target
330  \param partial_direction_1 [in] - first direction that partial is taken with respect to, e.g. 0 is x direction, 1 is y direction
331  \param partial_direction_2 [in] - second direction that partial is taken with respect to, e.g. 0 is x direction, 1 is y direction
332  \param dimension [in] - spatial dimension of basis to evaluate. e.g. dimension two basis of order one is 1, x, y, whereas for dimension 3 it is 1, x, y, z
333  \param poly_order [in] - polynomial basis degree
334  \param specific_order_only [in] - boolean for only evaluating one degree of polynomial when true
335  \param V [in] - orthonormal basis matrix size _dimensions * _dimensions whose first _dimensions-1 columns are an approximation of the tangent plane
336  \param reconstruction_space [in] - space of polynomial that a sampling functional is to evaluate
337  \param sampling_strategy [in] - sampling functional specification
338  \param additional_evaluation_local_index [in] - local index for evaluation sites
339  */
340  KOKKOS_INLINE_FUNCTION
341  void calcHessianPij(const member_type& teamMember, double* delta, double* thread_workspace, const int target_index, int neighbor_index, const double alpha, const int partial_direction_1, const int partial_direction_2, const int dimension, const int poly_order, bool specific_order_only, const scratch_matrix_right_type* V, const ReconstructionSpace reconstruction_space, const SamplingFunctional sampling_strategy, const int additional_evaluation_local_index = 0) const;
342 
343  /*! \brief Fills the _P matrix with either P or P*sqrt(w)
344  \param teamMember [in] - Kokkos::TeamPolicy member type (created by parallel_for)
345  \param delta [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large is the _basis_multipler*the dimension of the polynomial basis.
346  \param thread_workspace [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large as the _poly_order*the spatial dimension of the polynomial basis.
347  \param P [out] - 2D Kokkos View which will contain evaluation of sampling functional on polynomial basis for each neighbor the target has (stored column major)
348  \param w [out] - 1D Kokkos View which will contain weighting kernel values for the target with each neighbor if weight_p = true
349  \param dimension [in] - spatial dimension of basis to evaluate. e.g. dimension two basis of order one is 1, x, y, whereas for dimension 3 it is 1, x, y, z
350  \param polynomial_order [in] - polynomial basis degree
351  \param weight_p [in] - boolean whether to fill w with kernel weights
352  \param V [in] - orthonormal basis matrix size _dimensions * _dimensions whose first _dimensions-1 columns are an approximation of the tangent plane
353  \param reconstruction_space [in] - space of polynomial that a sampling functional is to evaluate
354  \param sampling_strategy [in] - sampling functional specification
355  */
356  KOKKOS_INLINE_FUNCTION
357  void createWeightsAndP(const member_type& teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P, scratch_vector_type w, const int dimension, int polynomial_order, bool weight_p = false, scratch_matrix_right_type* V = NULL, const ReconstructionSpace reconstruction_space = ReconstructionSpace::ScalarTaylorPolynomial, const SamplingFunctional sampling_strategy = PointSample) const;
358 
359  /*! \brief Fills the _P matrix with P*sqrt(w) for use in solving for curvature
360 
361  Uses _curvature_poly_order as the polynomial order of the basis
362 
363  \param teamMember [in] - Kokkos::TeamPolicy member type (created by parallel_for)
364  \param delta [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large is the _basis_multipler*the dimension of the polynomial basis.
365  \param thread_workspace [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large as the _poly_order*the spatial dimension of the polynomial basis.
366  \param P [out] - 2D Kokkos View which will contain evaluation of sampling functional on polynomial basis for each neighbor the target has (stored column major)
367  \param w [out] - 1D Kokkos View which will contain weighting kernel values for the target with each neighbor if weight_p = true
368  \param dimension [in] - spatial dimension of basis to evaluate. e.g. dimension two basis of order one is 1, x, y, whereas for dimension 3 it is 1, x, y, z
369  \param only_specific_order [in] - boolean for only evaluating one degree of polynomial when true
370  \param V [in] - orthonormal basis matrix size _dimensions * _dimensions whose first _dimensions-1 columns are an approximation of the tangent plane
371  */
372  KOKKOS_INLINE_FUNCTION
373  void createWeightsAndPForCurvature(const member_type& teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P, scratch_vector_type w, const int dimension, bool only_specific_order, scratch_matrix_right_type* V = NULL) const;
374 
375  /*! \brief Evaluates a polynomial basis with a target functional applied to each member of the basis
376  \param teamMember [in] - Kokkos::TeamPolicy member type (created by parallel_for)
377  \param delta [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large is the _basis_multipler*the dimension of the polynomial basis.
378  \param thread_workspace [in/out] - scratch space that is allocated so that each team has its own copy. Must be at least as large is the _poly_order*_global_dimensions.
379  \param P_target_row [out] - 1D Kokkos View where the evaluation of the polynomial basis is stored
380  */
381  KOKKOS_INLINE_FUNCTION
383 
384  /*! \brief Evaluates a polynomial basis for the curvature with a gradient target functional applied
385 
386  _operations is used by this function which is set through a modifier function
387 
388  \param teamMember [in] - Kokkos::TeamPolicy member type (created by parallel_for)
389  \param delta [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large is the _basis_multipler*the dimension of the polynomial basis.
390  \param thread_workspace [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large as the _curvature_poly_order*the spatial dimension of the polynomial basis.
391  \param P_target_row [out] - 1D Kokkos View where the evaluation of the polynomial basis is stored
392  \param V [in] - orthonormal basis matrix size _dimensions * _dimensions whose first _dimensions-1 columns are an approximation of the tangent plane
393  */
394  KOKKOS_INLINE_FUNCTION
395  void computeCurvatureFunctionals(const member_type& teamMember, scratch_vector_type t1, scratch_vector_type t2, scratch_matrix_right_type P_target_row, const scratch_matrix_right_type* V, const local_index_type local_neighbor_index = -1) const;
396 
397  /*! \brief Evaluates a polynomial basis with a target functional applied, using information from the manifold curvature
398 
399  _operations is used by this function which is set through a modifier function
400 
401  \param teamMember [in] - Kokkos::TeamPolicy member type (created by parallel_for)
402  \param delta [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large is the _basis_multipler*the dimension of the polynomial basis.
403  \param thread_workspace [in/out] - scratch space that is allocated so that each thread has its own copy. Must be at least as large as the _curvature_poly_order*the spatial dimension of the polynomial basis.
404  \param P_target_row [out] - 1D Kokkos View where the evaluation of the polynomial basis is stored
405  \param V [in] - orthonormal basis matrix size _dimensions * _dimensions whose first _dimensions-1 columns are an approximation of the tangent plane
406  \param G_inv [in] - (_dimensions-1)*(_dimensions-1) Kokkos View containing inverse of metric tensor
407  \param curvature_coefficients [in] - polynomial coefficients for curvature
408  \param curvature_gradients [in] - approximation of gradient of curvature, Kokkos View of size (_dimensions-1)
409  */
410  KOKKOS_INLINE_FUNCTION
412 
413  //! Helper function for applying the evaluations from a target functional to the polynomial coefficients
414  KOKKOS_INLINE_FUNCTION
416 
417 ///@}
418 
419 
420 
421 /** @name Private Accessors
422  * Private function because information lives on the device
423  */
424 ///@{
425 
426  //! Returns number of neighbors for a particular target
427  KOKKOS_INLINE_FUNCTION
428  int getNNeighbors(const int target_index) const {
429  return _neighbor_lists.getNumberOfNeighborsDevice(target_index);
430  }
431 
432  //! Mapping from [0,number of neighbors for a target] to the row that contains the source coordinates for
433  //! that neighbor
434  KOKKOS_INLINE_FUNCTION
435  int getNeighborIndex(const int target_index, const int neighbor_list_num) const {
436  return _neighbor_lists.getNeighborDevice(target_index, neighbor_list_num);
437  }
438 
439  //! Returns the maximum neighbor lists size over all target sites
440  KOKKOS_INLINE_FUNCTION
441  int getMaxNNeighbors() const {
443  }
444 
445  //! Returns the maximum number of evaluation sites over all target sites (target sites are included in total)
446  KOKKOS_INLINE_FUNCTION
449  }
450 
451  //! (OPTIONAL)
452  //! Returns number of additional evaluation sites for a particular target
453  KOKKOS_INLINE_FUNCTION
454  int getNEvaluationSitesPerTarget(const int target_index) const {
455  return _number_of_additional_evaluation_indices(target_index)+1;
456  }
457 
458  //! (OPTIONAL)
459  //! Mapping from [0,number of additional evaluation sites for a target] to the row that contains the coordinates for
460  //! that evaluation
461  KOKKOS_INLINE_FUNCTION
462  int getAdditionalEvaluationIndex(const int target_index, const int additional_list_num) const {
463  compadre_kernel_assert_debug((additional_list_num >= 1)
464  && "additional_list_num must be greater than or equal to 1, unlike neighbor lists which begin indexing at 0.");
465  return _additional_evaluation_indices(target_index, additional_list_num);
466  }
467 
468  //! Returns one component of the target coordinate for a particular target. Whether global or local coordinates
469  //! depends upon V being specified
470  KOKKOS_INLINE_FUNCTION
471  double getTargetCoordinate(const int target_index, const int dim, const scratch_matrix_right_type* V = NULL) const {
472  compadre_kernel_assert_debug((_target_coordinates.extent(0) >= (size_t)target_index) && "Target index is out of range for _target_coordinates.");
473  if (V==NULL) {
474  return _target_coordinates(target_index, dim);
475  } else {
476  XYZ target_coord = XYZ( _target_coordinates(target_index, 0),
477  _target_coordinates(target_index, 1),
478  _target_coordinates(target_index, 2));
479  return this->convertGlobalToLocalCoordinate(target_coord, dim, V);
480  }
481  }
482 
483  //! (OPTIONAL)
484  //! Returns one component of the additional evaluation coordinates. Whether global or local coordinates
485  //! depends upon V being specified
486  KOKKOS_INLINE_FUNCTION
487  double getTargetAuxiliaryCoordinate(const int target_index, const int additional_list_num, const int dim, const scratch_matrix_right_type* V = NULL) const {
488  auto additional_evaluation_index = getAdditionalEvaluationIndex(target_index, additional_list_num);
489  compadre_kernel_assert_debug((_additional_evaluation_coordinates.extent(0) >= (size_t)additional_evaluation_index) && "Additional evaluation index is out of range for _additional_evaluation_coordinates.");
490  if (V==NULL) {
491  return _additional_evaluation_coordinates(additional_evaluation_index, dim);
492  } else {
493  XYZ additional_target_coord = XYZ( _additional_evaluation_coordinates(additional_evaluation_index, 0),
494  _additional_evaluation_coordinates(additional_evaluation_index, 1),
495  _additional_evaluation_coordinates(additional_evaluation_index, 2));
496  return this->convertGlobalToLocalCoordinate(additional_target_coord, dim, V);
497  }
498  }
499 
500  //! Returns one component of the neighbor coordinate for a particular target. Whether global or local coordinates
501  //! depends upon V being specified
502  KOKKOS_INLINE_FUNCTION
503  double getNeighborCoordinate(const int target_index, const int neighbor_list_num, const int dim, const scratch_matrix_right_type* V = NULL) const {
504  compadre_kernel_assert_debug((_source_coordinates.extent(0) >= (size_t)(this->getNeighborIndex(target_index, neighbor_list_num))) && "Source index is out of range for _source_coordinates.");
505  if (V==NULL) {
506  return _source_coordinates(this->getNeighborIndex(target_index, neighbor_list_num), dim);
507  } else {
508  XYZ neighbor_coord = XYZ(_source_coordinates(this->getNeighborIndex(target_index, neighbor_list_num), 0), _source_coordinates(this->getNeighborIndex(target_index, neighbor_list_num), 1), _source_coordinates(this->getNeighborIndex(target_index, neighbor_list_num), 2));
509  return this->convertGlobalToLocalCoordinate(neighbor_coord, dim, V);
510  }
511  }
512 
513  //! Returns the relative coordinate as a vector between the target site and the neighbor site.
514  //! Whether global or local coordinates depends upon V being specified
515  KOKKOS_INLINE_FUNCTION
516  XYZ getRelativeCoord(const int target_index, const int neighbor_list_num, const int dimension, const scratch_matrix_right_type* V = NULL) const {
517  XYZ coordinate_delta;
518 
519  coordinate_delta.x = this->getNeighborCoordinate(target_index, neighbor_list_num, 0, V) - this->getTargetCoordinate(target_index, 0, V);
520  if (dimension>1) coordinate_delta.y = this->getNeighborCoordinate(target_index, neighbor_list_num, 1, V) - this->getTargetCoordinate(target_index, 1, V);
521  if (dimension>2) coordinate_delta.z = this->getNeighborCoordinate(target_index, neighbor_list_num, 2, V) - this->getTargetCoordinate(target_index, 2, V);
522 
523  return coordinate_delta;
524  }
525 
526  //! Returns a component of the local coordinate after transformation from global to local under the orthonormal basis V.
527  KOKKOS_INLINE_FUNCTION
528  double convertGlobalToLocalCoordinate(const XYZ global_coord, const int dim, const scratch_matrix_right_type* V) const {
529  // only written for 2d manifold in 3d space
530  double val = 0;
531  val += global_coord.x * (*V)(dim, 0);
532  val += global_coord.y * (*V)(dim, 1); // can't be called from dimension 1 problem
533  if (_dimensions>2) val += global_coord.z * (*V)(dim, 2);
534  return val;
535  }
536 
537  //! Returns a component of the global coordinate after transformation from local to global under the orthonormal basis V^T.
538  KOKKOS_INLINE_FUNCTION
539  double convertLocalToGlobalCoordinate(const XYZ local_coord, const int dim, const scratch_matrix_right_type* V) const {
540  // only written for 2d manifold in 3d space
541  double val = 0.0;
542  if (dim == 0 && _dimensions==2) { // 2D problem with 1D manifold
543  val = local_coord.x * (*V)(0, dim);
544  } else if (dim == 0) { // 3D problem with 2D manifold
545  val = local_coord.x * ((*V)(0, dim) + (*V)(1, dim));
546  } else if (dim == 1) { // 3D problem with 2D manifold
547  val = local_coord.y * ((*V)(0, dim) + (*V)(1, dim));
548  }
549  return val;
550  }
551 
552  //! Handles offset from operation input/output + extra evaluation sites
553  int getTargetOffsetIndexHost(const int lro_num, const int input_component, const int output_component, const int additional_evaluation_local_index = 0) const {
554  return ( _total_alpha_values*additional_evaluation_local_index
555  + _host_lro_total_offsets[lro_num]
556  + input_component*_host_lro_output_tile_size[lro_num]
557  + output_component );
558  }
559 
560  //! Handles offset from operation input/output + extra evaluation sites
561  KOKKOS_INLINE_FUNCTION
562  int getTargetOffsetIndexDevice(const int lro_num, const int input_component, const int output_component, const int additional_evaluation_local_index = 0) const {
563  return ( _total_alpha_values*additional_evaluation_local_index
564  + _lro_total_offsets[lro_num]
565  + input_component*_lro_output_tile_size[lro_num]
566  + output_component );
567  }
568 
569 ///@}
570 
571 /** @name Private Utility
572  *
573  */
574 ///@{
575 
576  //! Parses a string to determine solver type
577  static DenseSolverType parseSolverType(const std::string& dense_solver_type) {
578  std::string solver_type_to_lower = dense_solver_type;
579  transform(solver_type_to_lower.begin(), solver_type_to_lower.end(), solver_type_to_lower.begin(), ::tolower);
580  if (solver_type_to_lower == "lu") {
581  return DenseSolverType::LU;
582  } else {
583  return DenseSolverType::QR;
584  }
585  }
586 
587  //! Parses a string to determine problem type
588  static ProblemType parseProblemType(const std::string& problem_type) {
589  std::string problem_type_to_lower = problem_type;
590  transform(problem_type_to_lower.begin(), problem_type_to_lower.end(), problem_type_to_lower.begin(), ::tolower);
591  if (problem_type_to_lower == "standard") {
592  return ProblemType::STANDARD;
593  } else if (problem_type_to_lower == "manifold") {
594  return ProblemType::MANIFOLD;
595  } else {
596  return ProblemType::STANDARD;
597  }
598  }
599 
600  //! Parses a string to determine constraint type
601  static ConstraintType parseConstraintType(const std::string& constraint_type) {
602  std::string constraint_type_to_lower = constraint_type;
603  transform(constraint_type_to_lower.begin(), constraint_type_to_lower.end(), constraint_type_to_lower.begin(), ::tolower);
604  if (constraint_type_to_lower == "none") {
606  } else if (constraint_type_to_lower == "neumann_grad_scalar") {
608  } else {
610  }
611  }
612 
613 ///@}
614 
615 public:
616 
617 /** @name Instantiation / Destruction
618  *
619  */
620 ///@{
621 
622  //! Minimal constructor providing no data (neighbor lists, source sites, target sites)
623  GMLS(ReconstructionSpace reconstruction_space,
624  const SamplingFunctional polynomial_sampling_strategy,
625  const SamplingFunctional data_sampling_strategy,
626  const int poly_order,
627  const int dimensions = 3,
628  const std::string dense_solver_type = std::string("QR"),
629  const std::string problem_type = std::string("STANDARD"),
630  const std::string constraint_type = std::string("NO_CONSTRAINT"),
631  const int manifold_curvature_poly_order = 2) :
632  _poly_order(poly_order),
633  _curvature_poly_order(manifold_curvature_poly_order),
634  _dimensions(dimensions),
635  _reconstruction_space(reconstruction_space),
636  _dense_solver_type(parseSolverType(dense_solver_type)),
637  _problem_type(parseProblemType(problem_type)),
638  _constraint_type(parseConstraintType(constraint_type)),
640  && (polynomial_sampling_strategy == VectorPointSample)) ? ManifoldVectorPointSample : polynomial_sampling_strategy),
642  && (data_sampling_strategy == VectorPointSample)) ? ManifoldVectorPointSample : data_sampling_strategy)
643  {
644 
645  // seed random number generator pool
647 
648  _NP = this->getNP(_poly_order, dimensions, _reconstruction_space);
649  Kokkos::fence();
650 
651  // register curvature operations for manifold problems
653  _curvature_support_operations = Kokkos::View<TargetOperation*>
654  ("operations needed for manifold gradient reconstruction", 1);
655  auto curvature_support_operations_mirror =
656  Kokkos::create_mirror_view(_curvature_support_operations);
657  curvature_support_operations_mirror(0) =
659  Kokkos::deep_copy(_curvature_support_operations, curvature_support_operations_mirror);
660  }
661 
662  _lro_lookup = std::vector<int>(TargetOperation::COUNT,-1); // hard coded against number of operations defined
663  _lro = std::vector<TargetOperation>();
664 
665  // various initializations
667 
670  _weighting_power = 2;
672 
674 
675  _basis_multiplier = 1;
677 
678  _nontrivial_nullspace = false;
683  _store_PTWP_inv_PTW = false;
684 
686 
687  _max_num_neighbors = 0;
689 
690  _global_dimensions = dimensions;
692  _local_dimensions = dimensions-1;
693  } else {
694  _local_dimensions = dimensions;
695  }
696 
699  }
700 
701  //! Constructor for the case when the data sampling functional does not match the polynomial
702  //! sampling functional. Only case anticipated is staggered Laplacian.
703  GMLS(const int poly_order,
704  const int dimensions = 3,
705  const std::string dense_solver_type = std::string("QR"),
706  const std::string problem_type = std::string("STANDARD"),
707  const std::string constraint_type = std::string("NO_CONSTRAINT"),
708  const int manifold_curvature_poly_order = 2)
709  : GMLS(ReconstructionSpace::VectorOfScalarClonesTaylorPolynomial, VectorPointSample, VectorPointSample, poly_order, dimensions, dense_solver_type, problem_type, constraint_type, manifold_curvature_poly_order) {}
710 
711  //! Constructor for the case when nonstandard sampling functionals or reconstruction spaces
712  //! are to be used. Reconstruction space and sampling strategy can only be set at instantiation.
713  GMLS(ReconstructionSpace reconstruction_space,
714  SamplingFunctional dual_sampling_strategy,
715  const int poly_order,
716  const int dimensions = 3,
717  const std::string dense_solver_type = std::string("QR"),
718  const std::string problem_type = std::string("STANDARD"),
719  const std::string constraint_type = std::string("NO_CONSTRAINT"),
720  const int manifold_curvature_poly_order = 2)
721  : GMLS(reconstruction_space, dual_sampling_strategy, dual_sampling_strategy, poly_order, dimensions, dense_solver_type, problem_type, constraint_type, manifold_curvature_poly_order) {}
722 
723  //! Destructor
724  ~GMLS(){
725  };
726 
727 ///@}
728 
729 /** @name Functors
730  * Member functions that perform operations on the entire batch
731  */
732 ///@{
733 
734  //! Tag for functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity
736 
737  //! Tag for functor to evaluate targets, apply target evaluation to polynomial coefficients to
738  //! store in _alphas
740 
741  //! Tag for functor to create a coarse tangent approximation from a given neighborhood of points
743 
744  //! Tag for functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity for curvature
746 
747  //! Tag for functor to evaluate curvature targets and construct accurate tangent direction approximation for manifolds
749 
750  //! Tag for functor to determine if tangent directions need reordered, and to reorder them if needed
752 
753  //! Tag for functor to evaluate curvature targets and apply to coefficients of curvature reconstruction
755 
756  //! Tag for functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity
758 
759  //! Tag for functor to evaluate targets, apply target evaluation to polynomial coefficients to store in _alphas
761 
762  //! Tag for functor to calculate prestencil weights to apply to data to transform into a format expected by a GMLS stencil
764 
765 
766 
767  //! Functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity
768  KOKKOS_INLINE_FUNCTION
769  void operator() (const AssembleStandardPsqrtW&, const member_type& teamMember) const;
770 
771  //! Functor to evaluate targets, apply target evaluation to polynomial coefficients to store in _alphas
772  KOKKOS_INLINE_FUNCTION
773  void operator() (const ApplyStandardTargets&, const member_type& teamMember) const;
774 
775  //! Functor to create a coarse tangent approximation from a given neighborhood of points
776  KOKKOS_INLINE_FUNCTION
777  void operator() (const ComputeCoarseTangentPlane&, const member_type& teamMember) const;
778 
779  //! Functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity for curvature
780  KOKKOS_INLINE_FUNCTION
781  void operator() (const AssembleCurvaturePsqrtW&, const member_type& teamMember) const;
782 
783  //! Functor to evaluate curvature targets and construct accurate tangent direction approximation for manifolds
784  KOKKOS_INLINE_FUNCTION
785  void operator() (const GetAccurateTangentDirections&, const member_type& teamMember) const;
786 
787  //! Functor to determine if tangent directions need reordered, and to reorder them if needed
788  //! We require that the normal is consistent with a right hand rule on the tangent vectors
789  KOKKOS_INLINE_FUNCTION
790  void operator() (const FixTangentDirectionOrdering&, const member_type& teamMember) const;
791 
792  //! Functor to evaluate curvature targets and apply to coefficients of curvature reconstruction
793  KOKKOS_INLINE_FUNCTION
794  void operator() (const ApplyCurvatureTargets&, const member_type& teamMember) const;
795 
796  //! Functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity
797  KOKKOS_INLINE_FUNCTION
798  void operator() (const AssembleManifoldPsqrtW&, const member_type& teamMember) const;
799 
800  //! Functor to evaluate targets, apply target evaluation to polynomial coefficients to store in _alphas
801  KOKKOS_INLINE_FUNCTION
802  void operator() (const ApplyManifoldTargets&, const member_type& teamMember) const;
803 
804  //! Functor to calculate prestencil weights to apply to data to transform into a format expected by a GMLS stencil
805  KOKKOS_INLINE_FUNCTION
806  void operator() (const ComputePrestencilWeights&, const member_type& teamMember) const;
807 
808 ///@}
809 
810 /** @name Public Utility
811  *
812  */
813 ///@{
814 
815  //! Returns size of the basis for a given polynomial order and dimension
816  //! General to dimension 1..3 and polynomial order m
817  //! The divfree options will return the divergence-free basis if true
818  KOKKOS_INLINE_FUNCTION
819  static int getNP(const int m, const int dimension = 3, const ReconstructionSpace r_space = ReconstructionSpace::ScalarTaylorPolynomial) {
821  return ScalarTaylorPolynomialBasis::getSize(m, dimension);
822  } else {
823  return DivergenceFreePolynomialBasis::getSize(m, dimension);
824  }
825  }
826 
827  //! Returns number of neighbors needed for unisolvency for a given basis order and dimension
828  KOKKOS_INLINE_FUNCTION
829  static int getNN(const int m, const int dimension = 3, const ReconstructionSpace r_space = ReconstructionSpace::ScalarTaylorPolynomial) {
830  // may need div-free argument in the future
831  const int np = getNP(m, dimension, r_space);
832  int nn = np;
833  switch (dimension) {
834  case 3:
835  nn = np * (1.7 + m*0.1);
836  break;
837  case 2:
838  nn = np * (1.4 + m*0.03);
839  break;
840  case 1:
841  nn = np * 1.1;
842  }
843  return nn;
844  }
845 
846  /*! \brief Evaluates the weighting kernel
847  \param r [in] - Euclidean distance of relative vector. Euclidean distance of (target - neighbor) in some basis.
848  \param h [in] - window size. Kernel is guaranteed to take on a value of zero if it exceeds h.
849  \param weighting_type [in] - weighting type to be evaluated as the kernel. e,g. power, Gaussian, etc..
850  \param power [in] - power parameter to be given to the kernel.
851  */
852  KOKKOS_INLINE_FUNCTION
853  static double Wab(const double r, const double h, const WeightingFunctionType& weighting_type, const int power) {
854  if (weighting_type == WeightingFunctionType::Power) {
855  return std::pow(1.0-std::abs(r/h), power) * double(1.0-std::abs(r/h)>0.0);
856  } else if (weighting_type == WeightingFunctionType::Gaussian) {
857  // 2.5066282746310002416124 = sqrt(2*pi)
858  double h_over_3 = h/3.0;
859  return 1./( h_over_3 * 2.5066282746310002416124 ) * std::exp(-.5*r*r/(h_over_3*h_over_3));
860  } else if (weighting_type == WeightingFunctionType::CubicSpline) {
861  double x = r/h;
862  return ((1-x)+x*(1-x)*(1-2*x)) * double(r<=h);
863  } else { // no weighting power specified
864  compadre_kernel_assert_release(false && "Invalid WeightingFunctionType selected.");
865  return 0;
866  }
867  }
868 
869  //! Returns Euclidean norm of a vector
870  KOKKOS_INLINE_FUNCTION
871  static double EuclideanVectorLength(const XYZ& delta_vector, const int dimension) {
872  double inside_val = delta_vector.x*delta_vector.x;
873  switch (dimension) {
874  case 3:
875  inside_val += delta_vector.z*delta_vector.z;
876  // no break is intentional
877  case 2:
878  inside_val += delta_vector.y*delta_vector.y;
879  // no break is intentional
880  default:
881  break;
882  }
883  return std::sqrt(inside_val);
884  }
885 
886 
887  //! Helper function for finding alpha coefficients
888  static int getTargetOutputIndex(const int operation_num, const int output_component_axis_1, const int output_component_axis_2, const int dimensions) {
889  const int axis_1_size = (TargetOutputTensorRank[operation_num] > 1) ? dimensions : 1;
890  return axis_1_size*output_component_axis_1 + output_component_axis_2; // 0 for scalar, 0 for vector;
891  }
892 
893  //! Helper function for finding alpha coefficients
894  static int getSamplingOutputIndex(const SamplingFunctional sf, const int output_component_axis_1, const int output_component_axis_2) {
895  const int axis_1_size = (sf.output_rank > 1) ? sf.output_rank : 1;
896  return axis_1_size*output_component_axis_1 + output_component_axis_2; // 0 for scalar, 0 for vector;
897  }
898 
899  //! Output rank for sampling operation
901  return sro.output_rank;
902  }
903 
904  //! Input rank for sampling operation
906  return sro.input_rank;
907  }
908 
909 
910 ///@}
911 
912 /** @name Accessors
913  * Retrieve member variables through public member functions
914  */
915 ///@{
916 
917 
918  //! Returns (size of the basis used in instance's polynomial reconstruction) x (data input dimension)
920  host_managed_local_index_type sizes("sizes", 2);
921  sizes(0) = _basis_multiplier*_NP;
923  return sizes;
924  }
925 
926  //! Returns size of the basis used in instance's polynomial reconstruction
928  auto sizes = this->getPolynomialCoefficientsDomainRangeSize();
929  return sizes(0);
930  }
931 
932  //! Returns 2D array size in memory on which coefficients are stored
934  auto M_by_N = this->getPolynomialCoefficientsDomainRangeSize();
935  compadre_assert_release(_entire_batch_computed_at_once
936  && "Entire batch not computed at once, so getFullPolynomialCoefficientsBasis() can not be called.");
938  && "generateAlphas() called with keep_coefficients set to false.");
939  host_managed_local_index_type sizes("sizes", 2);
941  getRHSDims(_dense_solver_type, _constraint_type, _reconstruction_space, _dimensions, M_by_N[1], M_by_N[0], sizes(0), sizes(1));
942  } else {
943  getPDims(_dense_solver_type, _constraint_type, _reconstruction_space, _dimensions, M_by_N[1], M_by_N[0], sizes(1), sizes(0));
944  }
945  return sizes;
946  }
947 
948  //! Dimension of the GMLS problem, set only at class instantiation
949  int getDimensions() const { return _dimensions; }
950 
951  //! Dimension of the GMLS problem's point data (spatial description of points in ambient space), set only at class instantiation
952  int getGlobalDimensions() const { return _global_dimensions; }
953 
954  //! Local dimension of the GMLS problem (less than global dimension if on a manifold), set only at class instantiation
955  int getLocalDimensions() const { return _local_dimensions; }
956 
957  //! Get dense solver type
959 
960  //! Get problem type
962 
963  //! Get constraint type
965 
966  //! Type for weighting kernel for GMLS problem
968 
969  //! Type for weighting kernel for curvature
971 
972  //! Power for weighting kernel for GMLS problem
973  int getWeightingPower() const { return _weighting_power; }
974 
975  //! Power for weighting kernel for curvature
977 
978  //! Number of quadrature points
980 
981  //! Order of quadrature points
983 
984  //! Dimensions of quadrature points
986 
987  //! Type of quadrature points
988  std::string getQuadratureType() const { return _quadrature_type; }
989 
990  //! Get neighbor list accessor
992 
993  //! Get a view (device) of all tangent direction bundles.
994  decltype(_T) getTangentDirections() const { return _T; }
995 
996  //! Get a view (device) of all reference outward normal directions.
997  decltype(_ref_N) getReferenceNormalDirections() const { return _ref_N; }
998 
999  //! Get component of tangent or normal directions for manifold problems
1000  double getTangentBundle(const int target_index, const int direction, const int component) const {
1001  // Component index 0.._dimensions-2 will return tangent direction
1002  // Component index _dimensions-1 will return the normal direction
1003  scratch_matrix_right_type::HostMirror
1004  T(_host_T.data() + target_index*_dimensions*_dimensions, _dimensions, _dimensions);
1005  return T(direction, component);
1006  }
1007 
1008  //! Get component of tangent or normal directions for manifold problems
1009  double getReferenceNormalDirection(const int target_index, const int component) const {
1011  "getRefenceNormalDirection called, but reference outwrad normal directions were never provided.");
1012  scratch_vector_type::HostMirror
1013  ref_N(_host_ref_N.data() + target_index*_dimensions, _dimensions);
1014  return ref_N(component);
1015  }
1016 
1017  //! Get the local index (internal) to GMLS for a particular TargetOperation
1018  //! Every TargetOperation has a global index which can be readily found in Compadre::TargetOperation
1019  //! but this function returns the index used inside of the GMLS class
1021  return _lro_lookup[(int)lro];
1022  }
1023 
1024  //! Get a view (device) of all rank 2 preprocessing tensors
1025  //! This is a rank 5 tensor that is able to provide data transformation
1026  //! into a form that GMLS is able to operate on. The ranks are as follows:
1027  //!
1028  //! 1 - Either size 2 if it operates on the target site and neighbor site (staggered schemes)
1029  //! or 1 if it operates only on the neighbor sites (almost every scheme)
1030  //!
1031  //! 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
1032  //!
1033  //! 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)
1034  //!
1035  //! 4 - Data transform resulting in rank 1 data for the GMLS operator will have size _local_dimensions, otherwise 1
1036  //!
1037  //! 5 - Data transform taking in rank 1 data will have size _global_dimensions, otherwise 1
1039  return _prestencil_weights;
1040  }
1041 
1042  //! Retrieves the offset for an operator based on input and output component, generic to row
1043  //! (but still multiplied by the number of neighbors for each row and then needs a neighbor number added
1044  //! to this returned value to be meaningful)
1045  int getAlphaColumnOffset(TargetOperation lro, const int output_component_axis_1,
1046  const int output_component_axis_2, const int input_component_axis_1,
1047  const int input_component_axis_2, const int additional_evaluation_local_index = 0) const {
1048 
1049  const int lro_number = _lro_lookup[(int)lro];
1050  compadre_assert_debug((lro_number >= 0) && "getAlphaColumnOffset called for a TargetOperation that was not registered.");
1051 
1052  // the target functional input indexing is sized based on the output rank of the sampling
1053  // functional used, which can not be inferred unless a specification of target functional,
1054  // reconstruction space, and sampling functional are all known (as was the case at the
1055  // construction of this class)
1056  const int input_index = getSamplingOutputIndex(_data_sampling_functional, input_component_axis_1, input_component_axis_2);
1057  const int output_index = getTargetOutputIndex((int)lro, output_component_axis_1, output_component_axis_2, _dimensions);
1058 
1059  return getTargetOffsetIndexHost(lro_number, input_index, output_index, additional_evaluation_local_index);
1060  }
1061 
1062  //! Get a view (device) of all alphas
1063  decltype(_alphas) getAlphas() const { return _alphas; }
1064 
1065  //! Get a view (device) of all polynomial coefficients basis
1068  && "Entire batch not computed at once, so getFullPolynomialCoefficientsBasis() can not be called.");
1070  && "generateAlphas() called with keep_coefficients set to false.");
1072  return _RHS;
1073  } else {
1074  return _P;
1075  }
1076  }
1077 
1078  //! Get the polynomial sampling functional specified at instantiation
1080 
1081  //! Get the data sampling functional specified at instantiation (often the same as the polynomial sampling functional)
1083 
1084  //! Get the reconstruction space specified at instantiation
1086 
1087  //! Helper function for getting alphas for scalar reconstruction from scalar data
1088  double getAlpha0TensorTo0Tensor(TargetOperation lro, const int target_index, const int neighbor_index, const int additional_evaluation_site = 0) const {
1089  // e.g. Dirac Delta target of a scalar field
1090  return getAlpha(lro, target_index, 0, 0, neighbor_index, 0, 0, additional_evaluation_site);
1091  }
1092 
1093  //! Helper function for getting alphas for vector reconstruction from scalar data
1094  double getAlpha0TensorTo1Tensor(TargetOperation lro, const int target_index, const int output_component, const int neighbor_index, const int additional_evaluation_site = 0) const {
1095  // e.g. gradient of a scalar field
1096  return getAlpha(lro, target_index, output_component, 0, neighbor_index, 0, 0, additional_evaluation_site);
1097  }
1098 
1099  //! Helper function for getting alphas for matrix reconstruction from scalar data
1100  double getAlpha0TensorTo2Tensor(TargetOperation lro, const int target_index, const int output_component_axis_1, const int output_component_axis_2, const int neighbor_index, const int additional_evaluation_site = 0) const {
1101  return getAlpha(lro, target_index, output_component_axis_1, output_component_axis_2, neighbor_index, 0, 0, additional_evaluation_site);
1102  }
1103 
1104  //! Helper function for getting alphas for scalar reconstruction from vector data
1105  double getAlpha1TensorTo0Tensor(TargetOperation lro, const int target_index, const int neighbor_index, const int input_component, const int additional_evaluation_site = 0) const {
1106  // e.g. divergence of a vector field
1107  return getAlpha(lro, target_index, 0, 0, neighbor_index, input_component, 0, additional_evaluation_site);
1108  }
1109 
1110  //! Helper function for getting alphas for vector reconstruction from vector data
1111  double getAlpha1TensorTo1Tensor(TargetOperation lro, const int target_index, const int output_component, const int neighbor_index, const int input_component, const int additional_evaluation_site = 0) const {
1112  // e.g. curl of a vector field
1113  return getAlpha(lro, target_index, output_component, 0, neighbor_index, input_component, 0, additional_evaluation_site);
1114  }
1115 
1116  //! Helper function for getting alphas for matrix reconstruction from vector data
1117  double getAlpha1TensorTo2Tensor(TargetOperation lro, const int target_index, const int output_component_axis_1, const int output_component_axis_2, const int neighbor_index, const int input_component, const int additional_evaluation_site = 0) const {
1118  // e.g. gradient of a vector field
1119  return getAlpha(lro, target_index, output_component_axis_1, output_component_axis_2, neighbor_index, input_component, 0, additional_evaluation_site);
1120  }
1121 
1122  //! Helper function for getting alphas for scalar reconstruction from matrix data
1123  double getAlpha2TensorTo0Tensor(TargetOperation lro, const int target_index, const int neighbor_index, const int input_component_axis_1, const int input_component_axis_2, const int additional_evaluation_site = 0) const {
1124  return getAlpha(lro, target_index, 0, 0, neighbor_index, input_component_axis_1, input_component_axis_2, additional_evaluation_site);
1125  }
1126 
1127  //! Helper function for getting alphas for vector reconstruction from matrix data
1128  double getAlpha2TensorTo1Tensor(TargetOperation lro, const int target_index, const int output_component, const int neighbor_index, const int input_component_axis_1, const int input_component_axis_2, const int additional_evaluation_site = 0) const {
1129  return getAlpha(lro, target_index, output_component, 0, neighbor_index, input_component_axis_1, input_component_axis_2, additional_evaluation_site);
1130  }
1131 
1132  //! Helper function for getting alphas for matrix reconstruction from matrix data
1133  double getAlpha2TensorTo2Tensor(TargetOperation lro, const int target_index, const int output_component_axis_1, const int output_component_axis_2, const int neighbor_index, const int input_component_axis_1, const int input_component_axis_2, const int additional_evaluation_site = 0) const {
1134  return getAlpha(lro, target_index, output_component_axis_1, output_component_axis_2, neighbor_index, input_component_axis_1, input_component_axis_2, additional_evaluation_site);
1135  }
1136 
1137  //! Gives index into alphas given two axes, which when incremented by the neighbor number transforms access into
1138  //! alphas from a rank 1 view into a rank 3 view.
1139  KOKKOS_INLINE_FUNCTION
1140  global_index_type getAlphaIndexDevice(const int target_index, const int alpha_column_offset) const {
1141 
1142  global_index_type total_neighbors_before_target = _neighbor_lists.getRowOffsetDevice(target_index);
1143  int total_added_alphas_before_target = target_index*_added_alpha_size;
1144 
1145  int alphas_per_tile_per_target = _neighbor_lists.getNumberOfNeighborsDevice(target_index) + _added_alpha_size;
1146 
1147  return (total_neighbors_before_target+TO_GLOBAL(total_added_alphas_before_target))
1149  + TO_GLOBAL(alpha_column_offset*alphas_per_tile_per_target);
1150 
1151  }
1152 
1153  //! Gives index into alphas given two axes, which when incremented by the neighbor number transforms access into
1154  //! alphas from a rank 1 view into a rank 3 view.
1155  global_index_type getAlphaIndexHost(const int target_index, const int alpha_column_offset) const {
1156 
1157  global_index_type total_neighbors_before_target = _neighbor_lists.getRowOffsetHost(target_index);
1158  int total_added_alphas_before_target = target_index*_added_alpha_size;
1159 
1160  int alphas_per_tile_per_target = _neighbor_lists.getNumberOfNeighborsHost(target_index) + _added_alpha_size;
1161 
1162  return (total_neighbors_before_target+TO_GLOBAL(total_added_alphas_before_target))
1164  + TO_GLOBAL(alpha_column_offset*alphas_per_tile_per_target);
1165 
1166  }
1167 
1168  //! Underlying function all interface helper functions call to retrieve alpha values
1169  double getAlpha(TargetOperation lro, const int target_index, const int output_component_axis_1, const int output_component_axis_2, const int neighbor_index, const int input_component_axis_1, const int input_component_axis_2, const int additional_evaluation_site = 0) const {
1170  // lro - the operator from TargetOperations
1171  // target_index - the # for the target site where information is required
1172  // neighbor_index - the # for the neighbor of the target
1173  //
1174  // This code support up to rank 2 tensors for inputs and outputs
1175  //
1176  // scalar reconstruction from scalar data: rank 0 to rank 0
1177  // provides 1 piece of information for each neighbor
1178  // scalar reconstruction from vector data (e.g. divergence): rank 1 to rank 0
1179  // provides 'd' pieces of information for each neighbor
1180  // vector reconstruction from scalar data (e.g. gradient): rank 0 to rank 1
1181  // provides 'd' piece of information for each neighbor
1182  // vector reconstruction from vector data (e.g. curl): rank 1 to rank 1
1183  // provides 'd'x'd' pieces of information for each neighbor
1184  //
1185  // This function would more reasonably be called from one of the getAlphaNTensorFromNTensor
1186  // which is much easier to understand with respect to indexing and only requesting indices
1187  // that are relavent to the operator in question.
1188  //
1189 
1190  const int alpha_column_offset = this->getAlphaColumnOffset( lro, output_component_axis_1,
1191  output_component_axis_2, input_component_axis_1, input_component_axis_2, additional_evaluation_site);
1192 
1193  auto alphas_index = this->getAlphaIndexHost(target_index, alpha_column_offset);
1194  return _host_alphas(alphas_index + neighbor_index);
1195  }
1196 
1197  //! Returns a stencil to transform data from its existing state into the input expected
1198  //! for some sampling functionals.
1199  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 {
1200  // for certain sampling strategies, linear combinations of the neighbor and target value are needed
1201  // for the traditional PointSample, this value is 1 for the neighbor and 0 for the target
1202  if (sro == PointSample ) {
1203  if (for_target) return 0; else return 1;
1204  }
1205 
1206  // these check conditions on the sampling operator and change indexing on target and neighbors
1207  // in order to reuse information, such as if the same data transformation is used, regardless
1208  // of target site or neighbor site
1209  const int target_index_in_weights =
1212  target_index : 0;
1213  const int neighbor_index_in_weights =
1215  neighbor_index : 0;
1216 
1217  return _host_prestencil_weights((int)for_target, target_index_in_weights, neighbor_index_in_weights,
1218  output_component, input_component);
1219  }
1220 
1221  //! Dimensions ^ output rank for target operation
1222  int getOutputDimensionOfOperation(TargetOperation lro, bool ambient = false) const {
1223  int return_val;
1224  if (ambient) return_val = std::pow(_global_dimensions, TargetOutputTensorRank[(int)lro]);
1225  else return_val = std::pow(_local_dimensions, TargetOutputTensorRank[(int)lro]);
1226  return return_val;
1227  }
1228 
1229  //! Dimensions ^ input rank for target operation (always in local chart if on a manifold, never ambient space)
1231  return this->_host_lro_input_tile_size[_lro_lookup[(int)lro]];
1232  // this is the same return values as the OutputDimensionOfSampling for the GMLS class's SamplingFunctional
1233  }
1234 
1235  //! Dimensions ^ output rank for sampling operation
1236  //! (always in local chart if on a manifold, never ambient space)
1238  return std::pow(_local_dimensions, sro.output_rank);
1239  }
1240 
1241  //! Dimensions ^ output rank for sampling operation
1242  //! (always in ambient space, never local chart on a manifold)
1244  return std::pow(_global_dimensions, sro.input_rank);
1245  }
1246 
1247  //! Calculate basis_multiplier
1249  // calculate the dimension of the basis
1250  // (a vector space on a manifold requires two components, for example)
1251  return std::pow(_local_dimensions, ActualReconstructionSpaceRank[(int)rs]);
1252  }
1253 
1254  //! Calculate sampling_multiplier
1256  // this would normally be SamplingOutputTensorRank[_data_sampling_functional], but we also want to handle the
1257  // case of reconstructions where a scalar basis is reused as a vector, and this handles everything
1258  // this handles scalars, vectors, and scalars that are reused as vectors
1259  int bm = this->calculateBasisMultiplier(rs);
1260  int sm = this->getOutputDimensionOfSampling(sro);
1262  // storage and computational efficiency by reusing solution to scalar problem for
1263  // a vector problem (in 3d, 27x cheaper computation, 9x cheaper storage)
1264  sm = std::min(bm,sm);
1265  }
1266  return sm;
1267  }
1268 
1269 
1270 ///@}
1271 
1272 
1273 /** @name Modifiers
1274  * Changed member variables through public member functions
1275  */
1276 ///@{
1277 
1279  if (_RHS.extent(0) > 0)
1280  _RHS = Kokkos::View<double*>("RHS",0);
1281  }
1282 
1283  //! Sets basic problem data (neighbor lists, source coordinates, and target coordinates)
1284  template<typename view_type_1, typename view_type_2, typename view_type_3, typename view_type_4>
1286  view_type_1 neighbor_lists,
1287  view_type_2 source_coordinates,
1288  view_type_3 target_coordinates,
1289  view_type_4 epsilons) {
1290  this->setNeighborLists<view_type_1>(neighbor_lists);
1291  this->setSourceSites<view_type_2>(source_coordinates);
1292  this->setTargetSites<view_type_3>(target_coordinates);
1293  this->setWindowSizes<view_type_4>(epsilons);
1294  }
1295 
1296  //! Sets basic problem data (neighbor lists data, number of neighbors list, source coordinates, and target coordinates)
1297  template<typename view_type_1, typename view_type_2, typename view_type_3, typename view_type_4>
1299  view_type_1 cr_neighbor_lists,
1300  view_type_1 number_of_neighbors_list,
1301  view_type_2 source_coordinates,
1302  view_type_3 target_coordinates,
1303  view_type_4 epsilons) {
1304  this->setNeighborLists<view_type_1>(cr_neighbor_lists, number_of_neighbors_list);
1305  this->setSourceSites<view_type_2>(source_coordinates);
1306  this->setTargetSites<view_type_3>(target_coordinates);
1307  this->setWindowSizes<view_type_4>(epsilons);
1308  }
1309 
1310  //! (OPTIONAL) Sets additional evaluation sites for each target site
1311  template<typename view_type_1, typename view_type_2>
1313  view_type_1 additional_evaluation_indices,
1314  view_type_2 additional_evaluation_coordinates) {
1315  this->setAuxiliaryEvaluationIndicesLists<view_type_1>(additional_evaluation_indices);
1316  this->setAuxiliaryEvaluationCoordinates<view_type_2>(additional_evaluation_coordinates);
1317  }
1318 
1319  //! Sets neighbor list information from compressed row neighborhood lists data (if same view_type).
1320  template <typename view_type>
1321  typename std::enable_if<view_type::rank==1&&std::is_same<decltype(_neighbor_lists)::internal_view_type,view_type>::value==1, void>::type
1322  setNeighborLists(view_type neighbor_lists, view_type number_of_neighbors_list) {
1323 
1324  _neighbor_lists = NeighborLists<view_type>(neighbor_lists, number_of_neighbors_list);
1327  Kokkos::parallel_for("copy neighbor list sizes", Kokkos::RangePolicy<host_execution_space>(0, _host_number_of_neighbors_list.extent(0)), KOKKOS_LAMBDA(const int i) {
1329  });
1330  Kokkos::fence();
1331  this->resetCoefficientData();
1332 
1333  }
1334 
1335  //! Sets neighbor list information from compressed row neighborhood lists data (if different view_type).
1336  template <typename view_type>
1337  typename std::enable_if<view_type::rank==1&&std::is_same<decltype(_neighbor_lists)::internal_view_type,view_type>::value==0, void>::type
1338  setNeighborLists(view_type neighbor_lists, view_type number_of_neighbors_list) {
1339 
1340  typedef decltype(_neighbor_lists)::internal_view_type gmls_view_type;
1341  gmls_view_type d_neighbor_lists("compressed row neighbor lists data", neighbor_lists.extent(0));
1342  gmls_view_type d_number_of_neighbors_list("number of neighbors list", number_of_neighbors_list.extent(0));
1343  Kokkos::deep_copy(d_neighbor_lists, neighbor_lists);
1344  Kokkos::deep_copy(d_number_of_neighbors_list, number_of_neighbors_list);
1345  Kokkos::fence();
1346  _neighbor_lists = NeighborLists<gmls_view_type>(d_neighbor_lists, d_number_of_neighbors_list);
1349  Kokkos::parallel_for("copy neighbor list sizes", Kokkos::RangePolicy<host_execution_space>(0, _host_number_of_neighbors_list.extent(0)), KOKKOS_LAMBDA(const int i) {
1351  });
1352  Kokkos::fence();
1353  this->resetCoefficientData();
1354 
1355  }
1356 
1357  //! Sets neighbor list information. Should be # targets x maximum number of neighbors for any target + 1.
1358  //! first entry in ever row should be the number of neighbors for the corresponding target.
1359  template <typename view_type>
1360  typename std::enable_if<view_type::rank==2, void>::type setNeighborLists(view_type neighbor_lists) {
1361 
1362  _neighbor_lists = Convert2DToCompressedRowNeighborLists<decltype(neighbor_lists), Kokkos::View<int*> >(neighbor_lists);
1365  Kokkos::parallel_for("copy neighbor list sizes", Kokkos::RangePolicy<host_execution_space>(0, _host_number_of_neighbors_list.extent(0)), KOKKOS_LAMBDA(const int i) {
1367  });
1368  Kokkos::fence();
1369  this->resetCoefficientData();
1370 
1371  }
1372 
1373  //! Sets source coordinate information. Rows of this 2D-array should correspond to neighbor IDs contained in the entries
1374  //! of the neighbor lists 2D array.
1375  template<typename view_type>
1376  void setSourceSites(view_type source_coordinates) {
1377 
1378  // allocate memory on device
1379  _source_coordinates = decltype(_source_coordinates)("device neighbor coordinates",
1380  source_coordinates.extent(0), source_coordinates.extent(1));
1381 
1382  typedef typename view_type::memory_space input_array_memory_space;
1383  if (std::is_same<input_array_memory_space, device_memory_space>::value) {
1384  // check if on the device, then copy directly
1385  // if it is, then it doesn't match the internal layout we use
1386  // then copy to the host mirror
1387  // switches potential layout mismatches
1388  Kokkos::deep_copy(_source_coordinates, source_coordinates);
1389  } else {
1390  // if is on the host, copy to the host mirror
1391  // then copy to the device
1392  // switches potential layout mismatches
1393  auto host_source_coordinates = Kokkos::create_mirror_view(_source_coordinates);
1394  Kokkos::deep_copy(host_source_coordinates, source_coordinates);
1395  // switches memory spaces
1396  Kokkos::deep_copy(_source_coordinates, host_source_coordinates);
1397  }
1398  this->resetCoefficientData();
1399  }
1400 
1401  //! Sets source coordinate information. Rows of this 2D-array should correspond to neighbor IDs contained in the entries
1402  //! of the neighbor lists 2D array.
1403  template<typename view_type>
1404  void setSourceSites(decltype(_source_coordinates) source_coordinates) {
1405  // allocate memory on device
1406  _source_coordinates = source_coordinates;
1407  this->resetCoefficientData();
1408  }
1409 
1410  //! Sets target coordinate information. Rows of this 2D-array should correspond to rows of the neighbor lists.
1411  template<typename view_type>
1412  void setTargetSites(view_type target_coordinates) {
1413  // allocate memory on device
1414  _target_coordinates = decltype(_target_coordinates)("device target coordinates",
1415  target_coordinates.extent(0), target_coordinates.extent(1));
1416 
1417  typedef typename view_type::memory_space input_array_memory_space;
1418  if (std::is_same<input_array_memory_space, device_memory_space>::value) {
1419  // check if on the device, then copy directly
1420  // if it is, then it doesn't match the internal layout we use
1421  // then copy to the host mirror
1422  // switches potential layout mismatches
1423  Kokkos::deep_copy(_target_coordinates, target_coordinates);
1424  } else {
1425  // if is on the host, copy to the host mirror
1426  // then copy to the device
1427  // switches potential layout mismatches
1428  auto host_target_coordinates = Kokkos::create_mirror_view(_target_coordinates);
1429  Kokkos::deep_copy(host_target_coordinates, target_coordinates);
1430  // switches memory spaces
1431  Kokkos::deep_copy(_target_coordinates, host_target_coordinates);
1432  }
1434  = decltype(_number_of_additional_evaluation_indices)("number of additional evaluation indices", target_coordinates.extent(0));
1435  Kokkos::deep_copy(_number_of_additional_evaluation_indices, 0);
1436  this->resetCoefficientData();
1437  }
1438 
1439  //! Sets target coordinate information. Rows of this 2D-array should correspond to rows of the neighbor lists.
1440  template<typename view_type>
1441  void setTargetSites(decltype(_target_coordinates) target_coordinates) {
1442  // allocate memory on device
1443  _target_coordinates = target_coordinates;
1445  = decltype(_number_of_additional_evaluation_indices)("number of additional evaluation indices", target_coordinates.extent(0));
1446  Kokkos::deep_copy(_number_of_additional_evaluation_indices, 0);
1447  this->resetCoefficientData();
1448  }
1449 
1450  //! Sets window sizes, also called the support of the kernel
1451  template<typename view_type>
1452  void setWindowSizes(view_type epsilons) {
1453 
1454  // allocate memory on device
1455  _epsilons = decltype(_epsilons)("device epsilons", epsilons.extent(0));
1456 
1457  _host_epsilons = Kokkos::create_mirror_view(_epsilons);
1458  Kokkos::deep_copy(_host_epsilons, epsilons);
1459  // copy data from host to device
1460  Kokkos::deep_copy(_epsilons, _host_epsilons);
1461  this->resetCoefficientData();
1462  }
1463 
1464  //! Sets window sizes, also called the support of the kernel (device)
1465  template<typename view_type>
1466  void setWindowSizes(decltype(_epsilons) epsilons) {
1467  // allocate memory on device
1468  _epsilons = epsilons;
1469  this->resetCoefficientData();
1470  }
1471 
1472  //! (OPTIONAL)
1473  //! Sets orthonormal tangent directions for reconstruction on a manifold. The first rank of this 2D array
1474  //! corresponds to the target indices, i.e., rows of the neighbor lists 2D array. The second rank is the
1475  //! ordinal of the tangent direction (spatial dimensions-1 are tangent, last one is normal), and the third
1476  //! rank is indices into the spatial dimension.
1477  template<typename view_type>
1478  void setTangentBundle(view_type tangent_directions) {
1479  // accept input from user as a rank 3 tensor
1480  // but convert data to a rank 2 tensor with the last rank of dimension = _dimensions x _dimensions
1481  // this allows for nonstrided views on the device later
1482 
1483  // allocate memory on device
1484  _T = decltype(_T)("device tangent directions", _target_coordinates.extent(0)*_dimensions*_dimensions);
1485 
1486  compadre_assert_release( (std::is_same<decltype(_T)::memory_space, typename view_type::memory_space>::value) &&
1487  "Memory space does not match between _T and tangent_directions");
1488 
1489  auto this_dimensions = _dimensions;
1490  auto this_T = _T;
1491  // rearrange data on device from data given on host
1492  Kokkos::parallel_for("copy tangent vectors", Kokkos::RangePolicy<device_execution_space>(0, _target_coordinates.extent(0)), KOKKOS_LAMBDA(const int i) {
1493  scratch_matrix_right_type T(this_T.data() + i*this_dimensions*this_dimensions, this_dimensions, this_dimensions);
1494  for (int j=0; j<this_dimensions; ++j) {
1495  for (int k=0; k<this_dimensions; ++k) {
1496  T(j,k) = tangent_directions(i, j, k);
1497  }
1498  }
1499  });
1501 
1502  // copy data from device back to host in rearranged format
1503  _host_T = Kokkos::create_mirror_view(_T);
1504  Kokkos::deep_copy(_host_T, _T);
1505  this->resetCoefficientData();
1506  }
1507 
1508  //! (OPTIONAL)
1509  //! Sets outward normal direction. For manifolds this may be used for orienting surface. It is also accessible
1510  //! for sampling operators that require a normal direction.
1511  template<typename view_type>
1512  void setReferenceOutwardNormalDirection(view_type outward_normal_directions, bool use_to_orient_surface = true) {
1513  // accept input from user as a rank 2 tensor
1514 
1515  // allocate memory on device
1516  _ref_N = decltype(_ref_N)("device normal directions", _target_coordinates.extent(0)*_dimensions);
1517  // to assist LAMBDA capture
1518  auto this_ref_N = this->_ref_N;
1519  auto this_dimensions = this->_dimensions;
1520 
1521  // rearrange data on device from data given on host
1522  Kokkos::parallel_for("copy normal vectors", Kokkos::RangePolicy<device_execution_space>(0, _target_coordinates.extent(0)), KOKKOS_LAMBDA(const int i) {
1523  for (int j=0; j<this_dimensions; ++j) {
1524  this_ref_N(i*this_dimensions + j) = outward_normal_directions(i, j);
1525  }
1526  });
1527  Kokkos::fence();
1530 
1531  // copy data from device back to host in rearranged format
1532  _host_ref_N = Kokkos::create_mirror_view(_ref_N);
1533  Kokkos::deep_copy(_host_ref_N, _ref_N);
1534  this->resetCoefficientData();
1535  }
1536 
1537  //! (OPTIONAL)
1538  //! Sets extra data to be used by sampling functionals in certain instances.
1539  template<typename view_type>
1540  void setSourceExtraData(view_type extra_data) {
1541 
1542  // allocate memory on device
1543  _source_extra_data = decltype(_source_extra_data)("device source extra data", extra_data.extent(0), extra_data.extent(1));
1544 
1545  auto host_extra_data = Kokkos::create_mirror_view(_source_extra_data);
1546  Kokkos::deep_copy(host_extra_data, extra_data);
1547  // copy data from host to device
1548  Kokkos::deep_copy(_source_extra_data, host_extra_data);
1549  this->resetCoefficientData();
1550  }
1551 
1552  //! (OPTIONAL)
1553  //! Sets extra data to be used by sampling functionals in certain instances. (device)
1554  template<typename view_type>
1555  void setSourceExtraData(decltype(_source_extra_data) extra_data) {
1556  // allocate memory on device
1557  _source_extra_data = extra_data;
1558  this->resetCoefficientData();
1559  }
1560 
1561  //! (OPTIONAL)
1562  //! Sets extra data to be used by target operations in certain instances.
1563  template<typename view_type>
1564  void setTargetExtraData(view_type extra_data) {
1565 
1566  // allocate memory on device
1567  _target_extra_data = decltype(_target_extra_data)("device target extra data", extra_data.extent(0), extra_data.extent(1));
1568 
1569  auto host_extra_data = Kokkos::create_mirror_view(_target_extra_data);
1570  Kokkos::deep_copy(host_extra_data, extra_data);
1571  // copy data from host to device
1572  Kokkos::deep_copy(_target_extra_data, host_extra_data);
1573  this->resetCoefficientData();
1574  }
1575 
1576  //! (OPTIONAL)
1577  //! Sets extra data to be used by target operations in certain instances. (device)
1578  template<typename view_type>
1579  void setTargetExtraData(decltype(_target_extra_data) extra_data) {
1580  // allocate memory on device
1581  _target_extra_data = extra_data;
1582  this->resetCoefficientData();
1583  }
1584 
1585  //! (OPTIONAL)
1586  //! Sets additional points for evaluation of target operation on polynomial reconstruction.
1587  //! If this is never called, then the target sites are the only locations where the target
1588  //! operations will be evaluated and applied to polynomial reconstructions.
1589  template <typename view_type>
1590  void setAuxiliaryEvaluationCoordinates(view_type evaluation_coordinates) {
1591  // allocate memory on device
1592  _additional_evaluation_coordinates = decltype(_additional_evaluation_coordinates)("device additional evaluation coordinates",
1593  evaluation_coordinates.extent(0), evaluation_coordinates.extent(1));
1594 
1595  typedef typename view_type::memory_space input_array_memory_space;
1596  if (std::is_same<input_array_memory_space, device_memory_space>::value) {
1597  // check if on the device, then copy directly
1598  // if it is, then it doesn't match the internal layout we use
1599  // then copy to the host mirror
1600  // switches potential layout mismatches
1601  Kokkos::deep_copy(_additional_evaluation_coordinates, evaluation_coordinates);
1602  } else {
1603  // if is on the host, copy to the host mirror
1604  // then copy to the device
1605  // switches potential layout mismatches
1606  auto host_additional_evaluation_coordinates = Kokkos::create_mirror_view(_additional_evaluation_coordinates);
1607  Kokkos::deep_copy(host_additional_evaluation_coordinates, evaluation_coordinates);
1608  // switches memory spaces
1609  Kokkos::deep_copy(_additional_evaluation_coordinates, host_additional_evaluation_coordinates);
1610  }
1611  this->resetCoefficientData();
1612  }
1613 
1614  //! (OPTIONAL)
1615  //! Sets additional points for evaluation of target operation on polynomial reconstruction.
1616  //! If this is never called, then the target sites are the only locations where the target
1617  //! operations will be evaluated and applied to polynomial reconstructions. (device)
1618  template <typename view_type>
1620  _additional_evaluation_coordinates = evaluation_coordinates;
1621  this->resetCoefficientData();
1622  }
1623 
1624  //! (OPTIONAL)
1625  //! Sets the additional target evaluation coordinate indices list information. Should be # targets x maximum number of indices
1626  //! evaluation indices for any target + 1. first entry in every row should be the number of indices for the corresponding target.
1627  template <typename view_type>
1628  void setAuxiliaryEvaluationIndicesLists(view_type indices_lists) {
1629  // allocate memory on device
1630  _additional_evaluation_indices = decltype(_additional_evaluation_indices)("device additional evaluation indices",
1631  indices_lists.extent(0), indices_lists.extent(1));
1632 
1634 
1635  typedef typename view_type::memory_space input_array_memory_space;
1636  if (std::is_same<input_array_memory_space, device_memory_space>::value) {
1637  // check if on the device, then copy directly
1638  // if it is, then it doesn't match the internal layout we use
1639  // then copy to the host mirror
1640  // switches potential layout mismatches
1641  Kokkos::deep_copy(_additional_evaluation_indices, indices_lists);
1643  } else {
1644  // if is on the host, copy to the host mirror
1645  // then copy to the device
1646  // switches potential layout mismatches
1647  Kokkos::deep_copy(_host_additional_evaluation_indices, indices_lists);
1648  // copy data from host to device
1650  }
1651 
1653  auto number_of_additional_evaluation_indices = _number_of_additional_evaluation_indices;
1654  auto additional_evaluation_indices = _additional_evaluation_indices;
1655  Kokkos::parallel_reduce("additional evaluation indices",
1656  Kokkos::RangePolicy<device_execution_space>(0, _additional_evaluation_indices.extent(0)),
1657  KOKKOS_LAMBDA(const int i, int& t_max_evaluation_sites_per_target) {
1658  number_of_additional_evaluation_indices(i) = additional_evaluation_indices(i,0);
1659  t_max_evaluation_sites_per_target = (t_max_evaluation_sites_per_target > number_of_additional_evaluation_indices(i)+1)
1660  ? t_max_evaluation_sites_per_target : number_of_additional_evaluation_indices(i)+1;
1661  }, Kokkos::Max<int>(_max_evaluation_sites_per_target));
1662  Kokkos::fence();
1663  this->resetCoefficientData();
1664  }
1665 
1666  //! (OPTIONAL)
1667  //! Sets the additional target evaluation coordinate indices list information. Should be # targets x maximum number of indices
1668  //! evaluation indices for any target + 1. first entry in every row should be the number of indices for the corresponding target.
1669  template <typename view_type>
1671  // allocate memory on device
1672  _additional_evaluation_indices = indices_lists;
1673 
1675  // copy data from host to device
1677 
1679  auto number_of_additional_evaluation_indices = _number_of_additional_evaluation_indices;
1680  auto additional_evaluation_indices = _additional_evaluation_indices;
1681  Kokkos::parallel_reduce("additional evaluation indices",
1682  Kokkos::RangePolicy<device_execution_space>(0, _additional_evaluation_indices.extent(0)),
1683  KOKKOS_LAMBDA(const int i, int& t_max_evaluation_sites_per_target) {
1684  number_of_additional_evaluation_indices(i) = additional_evaluation_indices(i,0);
1685  t_max_evaluation_sites_per_target = (t_max_evaluation_sites_per_target > number_of_additional_evaluation_indices(i)+1)
1686  ? t_max_evaluation_sites_per_target : number_of_additional_evaluation_indices(i)+1;
1687  }, Kokkos::Max<int>(_max_evaluation_sites_per_target));
1688  Kokkos::fence();
1689  this->resetCoefficientData();
1690  }
1691 
1692 
1693  //! Type for weighting kernel for GMLS problem
1694  void setWeightingType( const std::string &wt) {
1695  std::string wt_to_lower = wt;
1696  transform(wt_to_lower.begin(), wt_to_lower.end(), wt_to_lower.begin(), ::tolower);
1697  if (wt_to_lower == "power") {
1699  } else if (wt_to_lower == "gaussian") {
1701  } else if (wt_to_lower == "cubicspline") {
1703  } else {
1704  // Power is default
1706  }
1707  this->resetCoefficientData();
1708  }
1709 
1710  //! Type for weighting kernel for GMLS problem
1712  _weighting_type = wt;
1713  this->resetCoefficientData();
1714  }
1715 
1716  //! Type for weighting kernel for curvature
1717  void setCurvatureWeightingType( const std::string &wt) {
1718  std::string wt_to_lower = wt;
1719  transform(wt_to_lower.begin(), wt_to_lower.end(), wt_to_lower.begin(), ::tolower);
1720  if (wt_to_lower == "power") {
1722  } else if (wt_to_lower == "gaussian") {
1724  } else if (wt_to_lower == "cubicspline") {
1726  } else {
1727  // Power is default
1729  }
1730  this->resetCoefficientData();
1731  }
1732 
1733  //! Type for weighting kernel for curvature
1736  this->resetCoefficientData();
1737  }
1738 
1739  //! Sets basis order to be used when reoncstructing any function
1740  void setPolynomialOrder(const int poly_order) {
1741  _poly_order = poly_order;
1742  _NP = this->getNP(_poly_order, _dimensions, _reconstruction_space);
1743  this->resetCoefficientData();
1744  }
1745 
1746  //! Sets basis order to be used when reoncstructing curvature
1747  void setCurvaturePolynomialOrder(const int manifold_poly_order) {
1748  _curvature_poly_order = manifold_poly_order;
1749  this->resetCoefficientData();
1750  }
1751 
1752  //! Power for weighting kernel for GMLS problem
1753  void setWeightingPower(int wp) {
1754  _weighting_power = wp;
1755  this->resetCoefficientData();
1756  }
1757 
1758  //! Power for weighting kernel for curvature
1761  this->resetCoefficientData();
1762  }
1763 
1764  //! Number quadrature points to use
1765  void setOrderOfQuadraturePoints(int order) {
1767  this->resetCoefficientData();
1768  }
1769 
1770  //! Dimensions of quadrature points to use
1773  this->resetCoefficientData();
1774  }
1775 
1776  //! Type of quadrature points
1777  void setQuadratureType(std::string quadrature_type) {
1778  _quadrature_type = quadrature_type;
1779  this->resetCoefficientData();
1780  }
1781 
1782  //! Adds a target to the vector of target functional to be applied to the reconstruction
1784  std::vector<TargetOperation> temporary_lro_vector(1, lro);
1785  this->addTargets(temporary_lro_vector);
1786  this->resetCoefficientData();
1787  }
1788 
1789  //! Adds a vector of target functionals to the vector of target functionals already to be applied to the reconstruction
1790  void addTargets(std::vector<TargetOperation> lro) {
1791  // if called multiple times with different dimensions, only the last
1792  // dimension called with is used for all
1793 
1794  // loop over requested targets
1795  for (size_t i=0; i<lro.size(); ++i) {
1796 
1797  bool operation_found = false;
1798  // loop over existing targets registered
1799  for (size_t j=0; j<_lro.size(); ++j) {
1800 
1801  // if found
1802  if (_lro[j]==lro[i]) {
1803 
1804  operation_found = true;
1805 
1806  // the operation should now point to where the operation is stored
1807  _lro_lookup[(int)lro[i]] = j;
1808 
1809  break;
1810 
1811  }
1812  }
1813 
1814  if (!operation_found) {
1815  _lro_lookup[(int)lro[i]] = _lro.size();
1816  _lro.push_back(lro[i]);
1817  }
1818  }
1819 
1820  _lro_total_offsets = Kokkos::View<int*>("total offsets for alphas", _lro.size());
1821  _lro_output_tile_size = Kokkos::View<int*>("output tile size for each operation", _lro.size());
1822  _lro_input_tile_size = Kokkos::View<int*>("output tile size for each operation", _lro.size());
1823  _lro_output_tensor_rank = Kokkos::View<int*>("output tensor rank", _lro.size());
1824  _lro_input_tensor_rank = Kokkos::View<int*>("input tensor rank", _lro.size());
1825 
1826  _host_lro_total_offsets = create_mirror_view(_lro_total_offsets);
1827  _host_lro_output_tile_size = create_mirror_view(_lro_output_tile_size);
1828  _host_lro_input_tile_size = create_mirror_view(_lro_input_tile_size);
1831 
1832  int total_offset = 0; // need total offset
1833  int output_offset = 0;
1834  int input_offset = 0;
1835  for (size_t i=0; i<_lro.size(); ++i) {
1836  _host_lro_total_offsets(i) = total_offset;
1837 
1838  // allows for a tile of the product of dimension^input_tensor_rank * dimension^output_tensor_rank * the number of neighbors
1839  int output_tile_size = getOutputDimensionOfOperation(_lro[i], false /* uses _local_dimension */);
1840 
1841  // the target functional input indexing is sized based on the output rank of the sampling
1842  // functional used
1844  _host_lro_output_tile_size(i) = output_tile_size;
1845  _host_lro_input_tile_size(i) = input_tile_size;
1846 
1847  total_offset += input_tile_size * output_tile_size;
1848  output_offset += output_tile_size;
1849  input_offset += input_tile_size;
1850 
1851  // the target functional output rank is based on the output rank of the sampling
1852  // functional used
1855  }
1856 
1857  _total_alpha_values = total_offset;
1858 
1860  // if on a manifold, the total alphas values must be large enough to hold the gradient
1861  // of the geometry reconstruction
1863  _total_alpha_values : std::pow(_local_dimensions, 1);
1864  }
1865 
1866  Kokkos::deep_copy(_lro_total_offsets, _host_lro_total_offsets);
1868  Kokkos::deep_copy(_lro_input_tile_size, _host_lro_input_tile_size);
1871  this->resetCoefficientData();
1872  }
1873 
1874  //! Empties the vector of target functionals to apply to the reconstruction
1875  void clearTargets() {
1876  _lro.clear();
1877  for (int i=0; i<TargetOperation::COUNT; ++i) {
1878  _lro_lookup[i] = -1;
1879  }
1880  this->resetCoefficientData();
1881  }
1882 
1883  /*! \brief Generates polynomial coefficients by setting up and solving least squares problems
1884  //! Sets up the batch of GMLS problems to be solved for. Provides alpha values
1885  //! that can later be contracted against data or degrees of freedom to form a
1886  //! global linear system.
1887  //! \param number_of_batches [in] - how many batches to break up the total workload into (for storage)
1888  //! \param keep_coefficients [in] - whether to store (P^T W P)^-1 * P^T * W
1889  */
1890  void generatePolynomialCoefficients(const int number_of_batches = 1, const bool keep_coefficients = false);
1891 
1892  /*! \brief Meant to calculate target operations and apply the evaluations to the previously
1893  //! constructed polynomial coefficients. But now that is inside of generatePolynomialCoefficients because
1894  //! it must be to handle number_of_batches>1. Effectively, this just calls generatePolynomialCoefficients.
1895  //! \param number_of_batches [in] - how many batches to break up the total workload into (for storage)
1896  //! \param keep_coefficients [in] - whether to store (P^T W P)^-1 * P^T * W
1897  */
1898  void generateAlphas(const int number_of_batches = 1, const bool keep_coefficients = false);
1899 
1900 ///@}
1901 
1902 
1903 
1904 
1905 }; // GMLS Class
1906 } // Compadre
1907 
1908 #endif
1909 
1910 
Kokkos::View< double * > _manifold_metric_tensor_inverse
metric tensor inverse for all problems
KOKKOS_INLINE_FUNCTION double convertLocalToGlobalCoordinate(const XYZ local_coord, const int dim, const scratch_matrix_right_type *V) const
Returns a component of the global coordinate after transformation from local to global under the orth...
Divergence-free vector polynomial basis.
Kokkos::View< int * > _lro_input_tile_size
dimensions ^ rank of tensor of output for each sampling functional (device)
Kokkos::View< double *, layout_right > _alphas
generated alpha coefficients (device)
KOKKOS_INLINE_FUNCTION int getNumberOfNeighborsDevice(int target_index) const
Get number of neighbors for a given target (device)
decltype(_RHS) getFullPolynomialCoefficientsBasis() const
Get a view (device) of all polynomial coefficients basis.
Standard GMLS problem type.
std::size_t global_index_type
Kokkos::View< const double *****, layout_right >::HostMirror _host_prestencil_weights
generated weights for nontraditional samples required to transform data into expected sampling functi...
Kokkos::View< TargetOperation * > _operations
vector containing target functionals to be applied for reconstruction problem (device) ...
bool _entire_batch_computed_at_once
whether entire calculation was computed at once the alternative is that it was broken up over many sm...
KOKKOS_INLINE_FUNCTION double getNeighborCoordinate(const int target_index, const int neighbor_list_num, const int dim, const scratch_matrix_right_type *V=NULL) const
Returns one component of the neighbor coordinate for a particular target.
int _max_evaluation_sites_per_target
maximum number of evaluation sites for each target (includes target site)
ProblemType getProblemType()
Get problem type.
Kokkos::View< double * > _manifold_curvature_coefficients
curvature polynomial coefficients for all problems
int _curvature_weighting_power
power to be used for weighting kernel for curvature
KOKKOS_INLINE_FUNCTION void createWeightsAndP(const member_type &teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P, scratch_vector_type w, const int dimension, int polynomial_order, bool weight_p=false, scratch_matrix_right_type *V=NULL, const ReconstructionSpace reconstruction_space=ReconstructionSpace::ScalarTaylorPolynomial, const SamplingFunctional sampling_strategy=PointSample) const
Fills the _P matrix with either P or P*sqrt(w)
int output_rank
Rank of sampling functional output for each SamplingFunctional.
void generatePolynomialCoefficients(const int number_of_batches=1, const bool keep_coefficients=false)
Generates polynomial coefficients by setting up and solving least squares problems ! Sets up the batc...
static KOKKOS_INLINE_FUNCTION double EuclideanVectorLength(const XYZ &delta_vector, const int dimension)
Returns Euclidean norm of a vector.
Kokkos::View< int **, layout_right > _additional_evaluation_indices
(OPTIONAL) contains indices of entries in the _additional_evaluation_coordinates view (device) ...
Tag for functor to evaluate curvature targets and apply to coefficients of curvature reconstruction...
Kokkos::View< double * > _ref_N
Rank 2 tensor for high order approximation of tangent vectors for all problems.
KOKKOS_INLINE_FUNCTION void calcGradientPij(const member_type &teamMember, double *delta, double *thread_workspace, const int target_index, int neighbor_index, const double alpha, const int partial_direction, const int dimension, const int poly_order, bool specific_order_only, const scratch_matrix_right_type *V, const ReconstructionSpace reconstruction_space, const SamplingFunctional sampling_strategy, const int additional_evaluation_local_index=0) const
Evaluates the gradient of a polynomial basis under the Dirac Delta (pointwise) sampling function...
KOKKOS_INLINE_FUNCTION int getSize(const int degree, const int dimension)
Returns size of basis.
#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...
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...
bool _reference_outward_normal_direction_provided
whether or not the reference outward normal directions were provided by the user. ...
KOKKOS_INLINE_FUNCTION int getMaxNumNeighbors() const
Get the maximum number of neighbors of all targets&#39; neighborhoods (host/device)
team_policy::member_type member_type
int getOrderOfQuadraturePoints() const
Order of quadrature points.
Neumann Gradient Scalar Type.
double getAlpha2TensorTo0Tensor(TargetOperation lro, const int target_index, const int neighbor_index, const int input_component_axis_1, const int input_component_axis_2, const int additional_evaluation_site=0) const
Helper function for getting alphas for scalar reconstruction from matrix data.
int _NP
dimension of basis for polynomial reconstruction
int _total_alpha_values
used for sizing P_target_row and the _alphas view
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)
Kokkos::View< int *, host_memory_space > _host_number_of_neighbors_list
convenient copy on host of number of neighbors
void setQuadratureType(std::string quadrature_type)
Type of quadrature points.
static int getTargetOutputIndex(const int operation_num, const int output_component_axis_1, const int output_component_axis_2, const int dimensions)
Helper function for finding alpha coefficients.
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 ...
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
int getWeightingPower() const
Power for weighting kernel for GMLS problem.
Kokkos::View< double **, layout_right > _target_coordinates
coordinates for target sites for reconstruction (device)
Each target applies a different data transform, but the same to each neighbor.
KOKKOS_INLINE_FUNCTION int getAdditionalEvaluationIndex(const int target_index, const int additional_list_num) const
(OPTIONAL) Mapping from [0,number of additional evaluation sites for a target] to the row that contai...
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.
std::vector< int > _lro_lookup
vector containing a mapping from a target functionals enum value to the its place in the list of targ...
int getInputDimensionOfSampling(SamplingFunctional sro) const
Dimensions ^ output rank for sampling operation (always in ambient space, never local chart on a mani...
int getNumberOfQuadraturePoints() const
Number of quadrature points.
Each target applies a different transform for each neighbor.
NeighborLists< Kokkos::View< int * > > _neighbor_lists
Accessor to get neighbor list data, offset data, and number of neighbors per target.
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...
Quadrature _qm
manages and calculates quadrature
Kokkos::View< double * >::HostMirror _host_T
tangent vectors information (host)
decltype(_ref_N) getReferenceNormalDirections() const
Get a view (device) of all reference outward normal directions.
Kokkos::View< int * >::HostMirror _host_lro_input_tensor_rank
tensor rank of sampling functional (host)
double getAlpha0TensorTo0Tensor(TargetOperation lro, const int target_index, const int neighbor_index, const int additional_evaluation_site=0) const
Helper function for getting alphas for scalar reconstruction from scalar data.
KOKKOS_INLINE_FUNCTION double getTargetCoordinate(const int target_index, const int dim, const scratch_matrix_right_type *V=NULL) const
Returns one component of the target coordinate for a particular target.
void clearTargets()
Empties the vector of target functionals to apply to the reconstruction.
std::vector< TargetOperation > _lro
vector of user requested target operations
void setTargetExtraData(view_type extra_data)
(OPTIONAL) Sets extra data to be used by target operations in certain instances.
int getTargetOperationLocalIndex(TargetOperation lro) const
Get the local index (internal) to GMLS for a particular TargetOperation Every TargetOperation has a g...
void setReferenceOutwardNormalDirection(view_type outward_normal_directions, bool use_to_orient_surface=true)
(OPTIONAL) Sets outward normal direction.
void setAuxiliaryEvaluationIndicesLists(decltype(_additional_evaluation_indices) indices_lists)
(OPTIONAL) Sets the additional target evaluation coordinate indices list information.
void setDimensionOfQuadraturePoints(int dim)
Dimensions of quadrature points to use.
Scalar polynomial basis centered at the target site and scaled by sum of basis powers e...
double getAlpha1TensorTo0Tensor(TargetOperation lro, const int target_index, const int neighbor_index, const int input_component, const int additional_evaluation_site=0) const
Helper function for getting alphas for scalar reconstruction from vector data.
Solve GMLS problem on a manifold (will use QR or SVD to solve the resultant GMLS problem dependent on...
KOKKOS_INLINE_FUNCTION void calcPij(const member_type &teamMember, double *delta, double *thread_workspace, const int target_index, int neighbor_index, const double alpha, const int dimension, const int poly_order, bool specific_order_only=false, const scratch_matrix_right_type *V=NULL, const ReconstructionSpace reconstruction_space=ReconstructionSpace::ScalarTaylorPolynomial, const SamplingFunctional sampling_strategy=PointSample, const int additional_evaluation_local_index=0) const
Evaluates the polynomial basis under a particular sampling function. Generally used to fill a row of ...
int _max_num_neighbors
maximum number of neighbors over all target sites
Kokkos::View< int * > _lro_input_tensor_rank
tensor rank of sampling functional (device)
Should be the total count of all available target functionals.
KOKKOS_INLINE_FUNCTION double convertGlobalToLocalCoordinate(const XYZ global_coord, const int dim, const scratch_matrix_right_type *V) const
Returns a component of the local coordinate after transformation from global to local under the ortho...
SamplingFunctional getPolynomialSamplingFunctional() const
Get the polynomial sampling functional specified at instantiation.
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 * >::HostMirror _host_lro_output_tensor_rank
tensor rank of target functional (host)
WeightingFunctionType _curvature_weighting_type
weighting kernel type for curvature problem
WeightingFunctionType _weighting_type
weighting kernel type for GMLS
Kokkos::View< int * > _lro_output_tile_size
dimensions ^ rank of tensor of output for each target functional (device)
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.
const SamplingFunctional _data_sampling_functional
generally the same as _polynomial_sampling_functional, but can differ if specified at GMLS class inst...
constexpr int TargetOutputTensorRank[]
Rank of target functional output for each TargetOperation Rank of target functional input for each Ta...
Tag for functor to evaluate targets, apply target evaluation to polynomial coefficients to store in _...
double getAlpha0TensorTo2Tensor(TargetOperation lro, const int target_index, const int output_component_axis_1, const int output_component_axis_2, const int neighbor_index, const int additional_evaluation_site=0) const
Helper function for getting alphas for matrix reconstruction from scalar data.
static ConstraintType parseConstraintType(const std::string &constraint_type)
Parses a string to determine constraint type.
KOKKOS_INLINE_FUNCTION int getNNeighbors(const int target_index) const
Returns number of neighbors for a particular target.
int calculateBasisMultiplier(const ReconstructionSpace rs) const
Calculate basis_multiplier.
DenseSolverType _dense_solver_type
solver type for GMLS problem - can be QR, SVD or LU
static ProblemType parseProblemType(const std::string &problem_type)
Parses a string to determine problem type.
KOKKOS_INLINE_FUNCTION int getNEvaluationSitesPerTarget(const int target_index) const
(OPTIONAL) Returns number of additional evaluation sites for a particular target
QR+Pivoting factorization performed on P*sqrt(w) matrix.
void addTargets(TargetOperation lro)
Adds a target to the vector of target functional to be applied to the reconstruction.
Tag for functor to create a coarse tangent approximation from a given neighborhood of points...
decltype(_alphas) getAlphas() const
Get a view (device) of all alphas.
pool_type _random_number_pool
Kokkos::View< int * > _lro_total_offsets
index for where this operation begins the for _alpha coefficients (device)
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.
ProblemType
Problem type, that optionally can handle manifolds.
ConstraintType getConstraintType()
Get constraint type.
Kokkos::View< int *, host_execution_space > host_managed_local_index_type
void setAuxiliaryEvaluationCoordinates(decltype(_additional_evaluation_coordinates) evaluation_coordinates)
(OPTIONAL) Sets additional points for evaluation of target operation on polynomial reconstruction...
void setAuxiliaryEvaluationIndicesLists(view_type indices_lists)
(OPTIONAL) Sets the additional target evaluation coordinate indices list information.
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 ...
#define TO_GLOBAL(variable)
KOKKOS_INLINE_FUNCTION int getSize(const int degree, const int dimension)
Returns size of basis.
void setCurvatureWeightingType(const std::string &wt)
Type for weighting kernel for curvature.
std::enable_if< view_type::rank==1 &&std::is_same< decltype(_neighbor_lists)::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 _initial_index_for_batch
initial index for current batch
int input_rank
Rank of sampling functional input for each SamplingFunctional.
KOKKOS_INLINE_FUNCTION int getTargetOffsetIndexDevice(const int lro_num, const int input_component, const int output_component, const int additional_evaluation_local_index=0) const
Handles offset from operation input/output + extra evaluation sites.
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...
SamplingFunctional getDataSamplingFunctional() const
Get the data sampling functional specified at instantiation (often the same as the polynomial samplin...
void setSourceSites(decltype(_source_coordinates) source_coordinates)
Sets source coordinate information.
Kokkos::View< double **, layout_right > _source_extra_data
Extra data available to basis functions (optional)
KOKKOS_INLINE_FUNCTION int getMaxEvaluationSitesPerTarget() const
Returns the maximum number of evaluation sites over all target sites (target sites are included in to...
void setTargetSites(decltype(_target_coordinates) target_coordinates)
Sets target coordinate information. Rows of this 2D-array should correspond to rows of the neighbor l...
Scalar basis reused as many times as there are components in the vector resulting in a much cheaper p...
Kokkos::View< int **, layout_right >::HostMirror _host_additional_evaluation_indices
(OPTIONAL) contains indices of entries in the _additional_evaluation_coordinates view (host) ...
KOKKOS_INLINE_FUNCTION int getMaxNNeighbors() const
Returns the maximum neighbor lists size over all target sites.
void setAdditionalEvaluationSitesData(view_type_1 additional_evaluation_indices, view_type_2 additional_evaluation_coordinates)
(OPTIONAL) Sets additional evaluation sites for each target site
std::string getQuadratureType() const
Type of quadrature points.
ReconstructionSpace
Space in which to reconstruct polynomial.
const SamplingFunctional _polynomial_sampling_functional
polynomial sampling functional used to construct P matrix, set at GMLS class instantiation ...
Tag for functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity.
std::enable_if< view_type::rank==1 &&std::is_same< decltype(_neighbor_lists)::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)...
DenseSolverType getDenseSolverType()
Get dense solver 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...
decltype(_T) getTangentDirections() const
Get a view (device) of all tangent direction bundles.
Tag for functor to calculate prestencil weights to apply to data to transform into a format expected ...
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)
Kokkos::View< double * >::HostMirror _host_epsilons
h supports determined through neighbor search (host)
host_managed_local_index_type getPolynomialCoefficientsMemorySize() const
Returns 2D array size in memory on which coefficients are stored.
double getAlpha0TensorTo1Tensor(TargetOperation lro, const int target_index, const int output_component, const int neighbor_index, const int additional_evaluation_site=0) const
Helper function for getting alphas for vector reconstruction from scalar data.
Kokkos::View< double **, layout_right, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_matrix_right_type
Kokkos::View< double *, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_vector_type
double getAlpha2TensorTo1Tensor(TargetOperation lro, const int target_index, const int output_component, const int neighbor_index, const int input_component_axis_1, const int input_component_axis_2, const int additional_evaluation_site=0) const
Helper function for getting alphas for vector reconstruction from matrix data.
KOKKOS_INLINE_FUNCTION int getNumberOfQuadraturePoints() const
int getAlphaColumnOffset(TargetOperation lro, const int output_component_axis_1, const int output_component_axis_2, const int input_component_axis_1, const int input_component_axis_2, const int additional_evaluation_local_index=0) const
Retrieves the offset for an operator based on input and output component, generic to row (but still m...
void setWeightingType(const WeightingFunctionType wt)
Type for weighting kernel for GMLS problem.
double getAlpha1TensorTo2Tensor(TargetOperation lro, const int target_index, const int output_component_axis_1, const int output_component_axis_2, const int neighbor_index, const int input_component, const int additional_evaluation_site=0) const
Helper function for getting alphas for matrix reconstruction from vector data.
WeightingFunctionType
Available weighting kernel function types.
int getOutputDimensionOfSampling(SamplingFunctional sro) const
Dimensions ^ output rank for sampling operation (always in local chart if on a manifold, never ambient space)
NeighborLists assists in accessing entries of compressed row neighborhood lists.
Kokkos::View< double **, layout_right > _source_coordinates
all coordinates for the source for which _neighbor_lists refers (device)
Kokkos::View< int * > _lro_output_tensor_rank
tensor rank of target functional (device)
void setCurvaturePolynomialOrder(const int manifold_poly_order)
Sets basis order to be used when reoncstructing curvature.
double getAlpha(TargetOperation lro, const int target_index, const int output_component_axis_1, const int output_component_axis_2, const int neighbor_index, const int input_component_axis_1, const int input_component_axis_2, const int additional_evaluation_site=0) const
Underlying function all interface helper functions call to retrieve alpha values. ...
void generateAlphas(const int number_of_batches=1, const bool keep_coefficients=false)
Meant to calculate target operations and apply the evaluations to the previously ! constructed polyno...
Kokkos::View< double **, layout_right > _target_extra_data
Extra data available to target operations (optional)
int _dimension_of_quadrature_points
dimension of quadrature rule
global_index_type getAlphaIndexHost(const int target_index, const int alpha_column_offset) const
Gives index into alphas given two axes, which when incremented by the neighbor number transforms acce...
double getReferenceNormalDirection(const int target_index, const int component) const
Get component of tangent or normal directions for manifold problems.
Kokkos::View< double * >::HostMirror _host_ref_N
reference outward normal vectors information (host)
KOKKOS_INLINE_FUNCTION double getTargetAuxiliaryCoordinate(const int target_index, const int additional_list_num, const int dim, const scratch_matrix_right_type *V=NULL) const
(OPTIONAL) Returns one component of the additional evaluation coordinates.
int getManifoldWeightingPower() const
Power for weighting kernel for curvature.
Point evaluation of the gradient of a scalar.
KOKKOS_INLINE_FUNCTION int getNeighborDevice(int target_index, int neighbor_num) const
Offers N(i,j) indexing where N(i,j) is the index of the jth neighbor of i (device) ...
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)
Minimal constructor providing no data (neighbor lists, source sites, target sites) ...
int getOutputDimensionOfOperation(TargetOperation lro, bool ambient=false) const
Dimensions ^ output rank for target operation.
double getAlpha1TensorTo1Tensor(TargetOperation lro, const int target_index, const int output_component, const int neighbor_index, const int input_component, const int additional_evaluation_site=0) const
Helper function for getting alphas for vector reconstruction from vector data.
Kokkos::View< int * >::HostMirror _host_lro_input_tile_size
dimensions ^ rank of tensor of output for each sampling functional (host)
ParallelManager _pm
determines scratch level spaces and is used to call kernels
host_managed_local_index_type getPolynomialCoefficientsDomainRangeSize() const
Returns (size of the basis used in instance&#39;s polynomial reconstruction) x (data input dimension) ...
void setWindowSizes(decltype(_epsilons) epsilons)
Sets window sizes, also called the support of the kernel (device)
int getDimensions() const
Dimension of the GMLS problem, set only at class instantiation.
int getPolynomialCoefficientsSize() const
Returns size of the basis used in instance&#39;s polynomial reconstruction.
Kokkos::View< TargetOperation * > _curvature_support_operations
vector containing target functionals to be applied for curvature
static int getInputRankOfSampling(SamplingFunctional sro)
Input rank for sampling operation.
Generalized Moving Least Squares (GMLS)
~GMLS()
Destructor.
KOKKOS_INLINE_FUNCTION XYZ getRelativeCoord(const int target_index, const int neighbor_list_num, const int dimension, const scratch_matrix_right_type *V=NULL) const
Returns the relative coordinate as a vector between the target site and the neighbor site...
KOKKOS_INLINE_FUNCTION global_index_type getRowOffsetDevice(int target_index) const
Get offset into compressed row neighbor lists (device)
void setTargetExtraData(decltype(_target_extra_data) extra_data)
(OPTIONAL) Sets extra data to be used by target operations in certain instances.
int _added_alpha_size
additional alpha coefficients due to constraints
ProblemType _problem_type
problem type for GMLS problem, can also be set to STANDARD for normal or MANIFOLD for manifold proble...
int getInputDimensionOfOperation(TargetOperation lro) const
Dimensions ^ input rank for target operation (always in local chart if on a manifold, never ambient space)
int _basis_multiplier
dimension of the reconstructed function e.g.
WeightingFunctionType getWeightingType() const
Type for weighting kernel for GMLS problem.
int _global_dimensions
spatial dimension of the points, set at class instantiation only
Kokkos::View< double **, layout_right > _additional_evaluation_coordinates
(OPTIONAL) user provided additional coordinates for target operation evaluation (device) ...
DenseSolverType
Dense solver type.
Kokkos::View< double * > _T
Rank 3 tensor for high order approximation of tangent vectors for all problems.
double getTangentBundle(const int target_index, const int direction, const int component) const
Get component of tangent or normal directions for manifold problems.
void resetCoefficientData()
KOKKOS_INLINE_FUNCTION void computeTargetFunctionals(const member_type &teamMember, scratch_vector_type t1, scratch_vector_type t2, scratch_matrix_right_type P_target_row) const
Evaluates a polynomial basis with a target functional applied to each member of the basis...
int getDimensionOfQuadraturePoints() const
Dimensions of quadrature points.
ReconstructionSpace getReconstructionSpace() const
Get the reconstruction space specified at instantiation.
void setSourceSites(view_type source_coordinates)
Sets source coordinate information.
LU factorization performed on P^T*W*P matrix.
bool _store_PTWP_inv_PTW
whether polynomial coefficients were requested to be stored (in a state not yet applied to data) ...
static int getOutputRankOfSampling(SamplingFunctional sro)
Output rank for sampling operation.
int transform_type
Describes the SamplingFunction relationship to targets, neighbors.
Kokkos::View< int * > _number_of_additional_evaluation_indices
(OPTIONAL) contains the # of additional coordinate indices for each target
void setPolynomialOrder(const int poly_order)
Sets basis order to be used when reoncstructing any function.
constexpr int ActualReconstructionSpaceRank[]
Number of actual components in the ReconstructionSpace.
void setTargetSites(view_type target_coordinates)
Sets target coordinate information. Rows of this 2D-array should correspond to rows of the neighbor l...
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< const double *, layout_right >::HostMirror _host_alphas
generated alpha coefficients (host)
int getTargetOffsetIndexHost(const int lro_num, const int input_component, const int output_component, const int additional_evaluation_local_index=0) const
Handles offset from operation input/output + extra evaluation sites.
Kokkos::Random_XorShift64_Pool pool_type
void setCurvatureWeightingPower(int wp)
Power for weighting kernel for curvature.
int _weighting_power
power to be used for weighting kernel
Kokkos::View< double * > _P
P*sqrt(w) matrix for all problems.
Tag for functor to evaluate targets, apply target evaluation to polynomial coefficients to store in _...
Tag for functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity.
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
KOKKOS_INLINE_FUNCTION global_index_type getAlphaIndexDevice(const int target_index, const int alpha_column_offset) const
Gives index into alphas given two axes, which when incremented by the neighbor number transforms acce...
decltype(_neighbor_lists) * getNeighborLists()
Get neighbor list accessor.
std::string _quadrature_type
quadrature rule type
Kokkos::View< TargetOperation * >::HostMirror _host_operations
vector containing target functionals to be applied for reconstruction problem (host) ...
static KOKKOS_INLINE_FUNCTION double Wab(const double r, const double h, const WeightingFunctionType &weighting_type, const int power)
Evaluates the weighting kernel.
bool _nontrivial_nullspace
whether or not operator to be inverted for GMLS problem has a nontrivial nullspace (requiring SVD) ...
double getAlpha2TensorTo2Tensor(TargetOperation lro, const int target_index, const int output_component_axis_1, const int output_component_axis_2, const int neighbor_index, const int input_component_axis_1, const int input_component_axis_2, const int additional_evaluation_site=0) const
Helper function for getting alphas for matrix reconstruction from matrix data.
KOKKOS_INLINE_FUNCTION void calcHessianPij(const member_type &teamMember, double *delta, double *thread_workspace, const int target_index, int neighbor_index, const double alpha, const int partial_direction_1, const int partial_direction_2, const int dimension, const int poly_order, bool specific_order_only, const scratch_matrix_right_type *V, const ReconstructionSpace reconstruction_space, const SamplingFunctional sampling_strategy, const int additional_evaluation_local_index=0) const
Evaluates the Hessian of a polynomial basis under the Dirac Delta (pointwise) sampling function...
KOKKOS_INLINE_FUNCTION void applyTargetsToCoefficients(const member_type &teamMember, scratch_vector_type t1, scratch_vector_type t2, scratch_matrix_right_type Q, scratch_vector_type w, scratch_matrix_right_type P_target_row, const int target_NP) const
Helper function for applying the evaluations from a target functional to the polynomial coefficients...
Tag for functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity for curva...
#define compadre_assert_debug(condition)
compadre_assert_debug is used for assertions that are checked in loops, as these significantly impact...
Kokkos::View< int * >::HostMirror _host_lro_total_offsets
index for where this operation begins the for _alpha coefficients (host)
KOKKOS_INLINE_FUNCTION void operator()(const AssembleStandardPsqrtW &, const member_type &teamMember) const
Functor to assemble the P*sqrt(weights) matrix and construct sqrt(weights)*Identity.
int getNumberOfNeighborsHost(int target_index) const
Get number of neighbors for a given target (host)
void setWeightingPower(int wp)
Power for weighting kernel for GMLS problem.
void setCurvatureWeightingType(const WeightingFunctionType wt)
Type for weighting kernel for curvature.
int getGlobalDimensions() const
Dimension of the GMLS problem&#39;s point data (spatial description of points in ambient space)...
int calculateSamplingMultiplier(const ReconstructionSpace rs, const SamplingFunctional sro) const
Calculate sampling_multiplier.
int local_index_type
ConstraintType _constraint_type
constraint type for GMLS problem
KOKKOS_INLINE_FUNCTION int getNeighborIndex(const int target_index, const int neighbor_list_num) const
Mapping from [0,number of neighbors for a target] to the row that contains the source coordinates for...
KOKKOS_INLINE_FUNCTION void computeCurvatureFunctionals(const member_type &teamMember, scratch_vector_type t1, scratch_vector_type t2, scratch_matrix_right_type P_target_row, const scratch_matrix_right_type *V, const local_index_type local_neighbor_index=-1) const
Evaluates a polynomial basis for the curvature with a gradient target functional applied.
int _local_dimensions
dimension of the problem, set at class instantiation only. For manifolds, generally _global_dimension...
Kokkos::View< int * >::HostMirror _host_lro_output_tile_size
dimensions ^ rank of tensor of output for each target functional (host)
global_index_type getRowOffsetHost(int target_index) const
Get offset into compressed row neighbor lists (host)
int _order_of_quadrature_points
order of exact polynomial integration for quadrature rule
#define compadre_assert_release(condition)
compadre_assert_release is used for assertions that should always be checked, but generally are not e...
Tag for functor to determine if tangent directions need reordered, and to reorder them if needed...
Tag for functor to evaluate curvature targets and construct accurate tangent direction approximation ...
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 void computeTargetFunctionalsOnManifold(const member_type &teamMember, scratch_vector_type t1, scratch_vector_type t2, scratch_matrix_right_type P_target_row, scratch_matrix_right_type V, scratch_matrix_right_type G_inv, scratch_vector_type curvature_coefficients, scratch_vector_type curvature_gradients) const
Evaluates a polynomial basis with a target functional applied, using information from the manifold cu...
WeightingFunctionType getManifoldWeightingType() const
Type for weighting kernel for curvature.
constexpr SamplingFunctional VectorPointSample
Point evaluations of the entire vector source function.
void setOrderOfQuadraturePoints(int order)
Number quadrature points to use.
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
TargetOperation
Available target functionals.
void setSourceExtraData(decltype(_source_extra_data) extra_data)
(OPTIONAL) Sets extra data to be used by sampling functionals in certain instances.
KOKKOS_INLINE_FUNCTION int getNumberOfTargets() const
Get number of total targets having neighborhoods (host/device).
void setTangentBundle(view_type tangent_directions)
(OPTIONAL) Sets orthonormal tangent directions for reconstruction on a manifold.
#define compadre_kernel_assert_debug(condition)
ConstraintType
Constraint type.
static int getSamplingOutputIndex(const SamplingFunctional sf, const int output_component_axis_1, const int output_component_axis_2)
Helper function for finding alpha coefficients.
KOKKOS_INLINE_FUNCTION void createWeightsAndPForCurvature(const member_type &teamMember, scratch_vector_type delta, scratch_vector_type thread_workspace, scratch_matrix_right_type P, scratch_vector_type w, const int dimension, bool only_specific_order, scratch_matrix_right_type *V=NULL) const
Fills the _P matrix with P*sqrt(w) for use in solving for curvature.