Compadre 1.5.5
Loading...
Searching...
No Matches
Compadre_Misc.hpp
Go to the documentation of this file.
1#ifndef _COMPADRE_MISC_HPP_
2#define _COMPADRE_MISC_HPP_
3
5
6namespace Compadre {
7
8struct XYZ {
9
10 double x;
11 double y;
12 double z;
13
15 XYZ(double _x = 0.0, double _y = 0.0, double _z = 0.0) : x(_x), y(_y), z(_z) {}
16
18 XYZ(const scalar_type* arry) : x(arry[0]), y(arry[1]), z(arry[2]) {};
19
20 XYZ(const std::vector<scalar_type>& vec) : x(vec[0]), y(vec[1]), z(vec[2]) {};
21
24 { x += other.x;
25 y += other.y;
26 z += other.z;
27 return *this; }
30 { x += other.x;
31 y += other.y;
32 z += other.z;
33 return *this; }
36 { x -= other.x;
37 y -= other.y;
38 z -= other.z;
39 return *this; }
42 { x -= other.x;
43 y -= other.y;
44 z -= other.z;
45 return *this; }
47 XYZ operator *= (const double& scaling)
48 { x *= scaling;
49 y *= scaling;
50 z *= scaling;
51 return *this; }
54 switch (i) {
55 case 0:
56 return x;
57 case 1:
58 return y;
59 default:
60 return z;
61 }
62 }
64 scalar_type operator [](const int i) const {
65 switch (i) {
66 case 0:
67 return x;
68 case 1:
69 return y;
70 default:
71 return z;
72 }
73 }
76 XYZ result;
77 result.x = scalar*x;
78 result.y = scalar*y;
79 result.z = scalar*z;
80 return result;
81 }
82
84 size_t extent(int comp = 0) {
85 return 3;
86 }
87
89 int extent_int(int comp = 0) {
90 return 3;
91 }
92
93}; // XYZ
94
95
96KOKKOS_INLINE_FUNCTION
97XYZ operator + ( const XYZ& vecA, const XYZ& vecB ) {
98 return XYZ( vecA.x + vecB.x, vecA.y + vecB.y, vecA.z + vecB.z); }
99
100 KOKKOS_INLINE_FUNCTION
101XYZ operator - ( const XYZ& vecA, const XYZ& vecB ) {
102 return XYZ( vecA.x - vecB.x, vecA.y - vecB.y, vecA.z - vecB.z); }
103
104KOKKOS_INLINE_FUNCTION
105XYZ operator * ( const XYZ& vecA, const XYZ& vecB ) {
106 return XYZ( vecA.x * vecB.x, vecA.y * vecB.y, vecA.z * vecB.z); }
107
108KOKKOS_INLINE_FUNCTION
109XYZ operator + ( const XYZ& vecA, const scalar_type& constant ) {
110 return XYZ( vecA.x + constant, vecA.y + constant, vecA.z + constant); }
111
112KOKKOS_INLINE_FUNCTION
113XYZ operator + ( const scalar_type& constant, const XYZ& vecA ) {
114 return XYZ( vecA.x + constant, vecA.y + constant, vecA.z + constant); }
115
116KOKKOS_INLINE_FUNCTION
117XYZ operator - ( const XYZ& vecA, const scalar_type& constant ) {
118 return XYZ( vecA.x - constant, vecA.y - constant, vecA.z - constant); }
119
120KOKKOS_INLINE_FUNCTION
121XYZ operator - ( const scalar_type& constant, const XYZ& vecA ) {
122 return XYZ( constant - vecA.x, constant - vecA.y, constant - vecA.z); }
123
124KOKKOS_INLINE_FUNCTION
125XYZ operator * ( const XYZ& vecA, const scalar_type& constant ) {
126 return XYZ( vecA.x * constant, vecA.y * constant, vecA.z * constant); }
127
128KOKKOS_INLINE_FUNCTION
129XYZ operator * (const scalar_type& constant, const XYZ& vecA) {
130 return XYZ( vecA.x * constant, vecA.y * constant, vecA.z * constant); }
131
132KOKKOS_INLINE_FUNCTION
133XYZ operator / ( const XYZ& vecA, const scalar_type& constant ) {
134 return XYZ( vecA.x / constant, vecA.y / constant, vecA.z / constant); }
135
136inline std::ostream& operator << ( std::ostream& os, const XYZ& vec ) {
137 os << "(" << vec.x << ", " << vec.y << ", " << vec.z << ")" ; return os; }
138
139//! n^p (n^0 returns 1, regardless of n)
140KOKKOS_INLINE_FUNCTION
141int pown(int n, unsigned p) {
142 // O(p) implementation
143 int y = 1;
144 for (unsigned i=0; i<p; i++) {
145 y *= n;
146 }
147 return y;
148 // O(log(p)) implementation
149 //int result = 1;
150 //while (p) {
151 // if (p & 0x1) {
152 // result *= n;
153 // }
154 // n *= n;
155 // p >>= 1;
156 //}
157 //return result;
158}
159
160KOKKOS_INLINE_FUNCTION
162 // Return the additional constraint size
163 if (constraint_type == NEUMANN_GRAD_SCALAR) {
164 return 1;
165 } else {
166 return 0;
167 }
168}
169
170KOKKOS_INLINE_FUNCTION
171int getAdditionalCoeffSizeFromConstraintAndSpace(DenseSolverType dense_solver_type, ConstraintType constraint_type, ReconstructionSpace reconstruction_space, const int dimension) {
172 // Return the additional constraint size
173 int added_alpha_size = getAdditionalAlphaSizeFromConstraint(dense_solver_type, constraint_type);
174 if (reconstruction_space == VectorTaylorPolynomial) {
175 return dimension*added_alpha_size;
176 } else {
177 return added_alpha_size;
178 }
179}
180
181KOKKOS_INLINE_FUNCTION
182void 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) {
183 // Return the appropriate size for _RHS. Since in LU, the system solves P^T*P against P^T*W.
184 // We store P^T*P in the RHS space, which means RHS can be much smaller compared to the
185 // case for QR/SVD where the system solves PsqrtW against sqrtW*Identity
186
187 int added_coeff_size = getAdditionalCoeffSizeFromConstraintAndSpace(dense_solver_type, constraint_type, reconstruction_space, dimension);
188
189 if (constraint_type == NEUMANN_GRAD_SCALAR) {
190 RHS_row = RHS_col = N + added_coeff_size;
191 } else {
192 if (dense_solver_type != LU) {
193 RHS_row = N;
194 RHS_col = M;
195 } else {
196 RHS_row = RHS_col = N + added_coeff_size;
197 }
198 }
199}
200
201KOKKOS_INLINE_FUNCTION
202void 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) {
203 // Return the appropriate size for _P.
204 // In the case of solving with LU and additional constraint is used, _P needs
205 // to be resized to include additional row(s) based on the type of constraint.
206
207 int added_coeff_size = getAdditionalCoeffSizeFromConstraintAndSpace(dense_solver_type, constraint_type, reconstruction_space, dimension);
208 int added_alpha_size = getAdditionalAlphaSizeFromConstraint(dense_solver_type, constraint_type);
209
210 if (constraint_type == NEUMANN_GRAD_SCALAR) {
211 out_row = M + added_alpha_size;
212 out_col = N + added_coeff_size;
213 } else {
214 if (dense_solver_type == LU) {
215 out_row = M + added_alpha_size;
216 out_col = N + added_coeff_size;
217 } else {
218 out_row = M;
219 out_col = N;
220 }
221 }
222}
223
224//! Helper function for finding alpha coefficients
225KOKKOS_INLINE_FUNCTION
226int getTargetOutputIndex(const int operation_num, const int output_component_axis_1, const int output_component_axis_2, const int dimensions) {
227 const int axis_1_size = (getTargetOutputTensorRank(operation_num) > 1) ? dimensions : 1;
228 return axis_1_size*output_component_axis_1 + output_component_axis_2; // 0 for scalar, 0 for vector;
229}
230
231//! Helper function for finding alpha coefficients
232KOKKOS_INLINE_FUNCTION
233int getSamplingOutputIndex(const SamplingFunctional sf, const int output_component_axis_1, const int output_component_axis_2) {
234 const int axis_1_size = (sf.output_rank > 1) ? sf.output_rank : 1;
235 return axis_1_size*output_component_axis_1 + output_component_axis_2; // 0 for scalar, 0 for vector;
236}
237
238//! Input rank for sampling operation
239KOKKOS_INLINE_FUNCTION
243
244//! Dimensions ^ output rank for sampling operation
245//! (always in local chart if on a manifold, never ambient space)
246KOKKOS_INLINE_FUNCTION
247int getOutputDimensionOfSampling(SamplingFunctional sro, const int local_dimensions) {
248 return pown(local_dimensions, sro.output_rank);
249}
250
251//! Dimensions ^ output rank for sampling operation
252//! (always in ambient space, never local chart on a manifold)
253KOKKOS_INLINE_FUNCTION
254int getInputDimensionOfSampling(SamplingFunctional sro, const int global_dimensions) {
255 return pown(global_dimensions, sro.input_rank);
256}
257
258//! Calculate basis_multiplier
259KOKKOS_INLINE_FUNCTION
260int calculateBasisMultiplier(const ReconstructionSpace rs, const int local_dimensions) {
261 // calculate the dimension of the basis
262 // (a vector space on a manifold requires two components, for example)
263 return pown(local_dimensions, getActualReconstructionSpaceRank((int)rs));
264}
265
266//! Calculate sampling_multiplier
267KOKKOS_INLINE_FUNCTION
268int calculateSamplingMultiplier(const ReconstructionSpace rs, const SamplingFunctional sro, const int local_dimensions) {
269 // this would normally be SamplingOutputTensorRank[_data_sampling_functional], but we also want to handle the
270 // case of reconstructions where a scalar basis is reused as a vector, and this handles everything
271 // this handles scalars, vectors, and scalars that are reused as vectors
272 int bm = calculateBasisMultiplier(rs, local_dimensions);
273 int sm = getOutputDimensionOfSampling(sro, local_dimensions);
275 // storage and computational efficiency by reusing solution to scalar problem for
276 // a vector problem (in 3d, 27x cheaper computation, 9x cheaper storage)
277 sm =(bm < sm) ? bm : sm;
278 }
279 return sm;
280}
281
282//! Dimensions ^ output rank for target operation
283KOKKOS_INLINE_FUNCTION
284int getOutputDimensionOfOperation(TargetOperation lro, const int local_dimensions) {
285 return pown(local_dimensions, getTargetOutputTensorRank(lro));
286}
287
288//! Dimensions ^ input rank for target operation (always in local chart if on a manifold, never ambient space)
289KOKKOS_INLINE_FUNCTION
290int getInputDimensionOfOperation(TargetOperation lro, SamplingFunctional sro, const int local_dimensions) {
291 // this is the same return values as the OutputDimensionOfSampling for the GMLS class's SamplingFunctional
292 return getOutputDimensionOfSampling(sro, local_dimensions);
293}
294
295} // Compadre
296
297#endif
double scalar_type
KOKKOS_INLINE_FUNCTION int pown(int n, unsigned p)
n^p (n^0 returns 1, regardless of n)
KOKKOS_INLINE_FUNCTION int getTargetOutputTensorRank(const int &index)
Rank of target functional output for each TargetOperation Rank of target functional input for each Ta...
KOKKOS_INLINE_FUNCTION int getAdditionalCoeffSizeFromConstraintAndSpace(DenseSolverType dense_solver_type, ConstraintType constraint_type, ReconstructionSpace reconstruction_space, const int dimension)
KOKKOS_INLINE_FUNCTION void getRHSDims(DenseSolverType dense_solver_type, ConstraintType constraint_type, ReconstructionSpace reconstruction_space, const int dimension, const int M, const int N, int &RHS_row, int &RHS_col)
KOKKOS_INLINE_FUNCTION XYZ operator-(const XYZ &vecA, const XYZ &vecB)
KOKKOS_INLINE_FUNCTION int getActualReconstructionSpaceRank(const int &index)
Number of actual components in the ReconstructionSpace.
KOKKOS_INLINE_FUNCTION XYZ operator*(const XYZ &vecA, const XYZ &vecB)
ConstraintType
Constraint type.
@ NEUMANN_GRAD_SCALAR
Neumann Gradient Scalar Type.
std::ostream & operator<<(std::ostream &os, const XYZ &vec)
KOKKOS_INLINE_FUNCTION int calculateBasisMultiplier(const ReconstructionSpace rs, const int local_dimensions)
Calculate basis_multiplier.
TargetOperation
Available target functionals.
KOKKOS_INLINE_FUNCTION int calculateSamplingMultiplier(const ReconstructionSpace rs, const SamplingFunctional sro, const int local_dimensions)
Calculate sampling_multiplier.
KOKKOS_INLINE_FUNCTION XYZ operator+(const XYZ &vecA, const XYZ &vecB)
KOKKOS_INLINE_FUNCTION int getAdditionalAlphaSizeFromConstraint(DenseSolverType dense_solver_type, ConstraintType constraint_type)
DenseSolverType
Dense solver type.
@ LU
LU factorization performed on P^T*W*P matrix.
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_INLINE_FUNCTION int getOutputDimensionOfOperation(TargetOperation lro, const int local_dimensions)
Dimensions ^ output rank for target operation.
KOKKOS_INLINE_FUNCTION int getInputDimensionOfOperation(TargetOperation lro, SamplingFunctional sro, const int local_dimensions)
Dimensions ^ input rank for target operation (always in local chart if on a manifold,...
KOKKOS_INLINE_FUNCTION XYZ operator/(const XYZ &vecA, const scalar_type &constant)
KOKKOS_INLINE_FUNCTION int getInputDimensionOfSampling(SamplingFunctional sro, const int global_dimensions)
Dimensions ^ output rank for sampling operation (always in ambient space, never local chart on a mani...
KOKKOS_INLINE_FUNCTION 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.
ReconstructionSpace
Space in which to reconstruct polynomial.
@ VectorTaylorPolynomial
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...
@ VectorOfScalarClonesTaylorPolynomial
Scalar basis reused as many times as there are components in the vector resulting in a much cheaper p...
KOKKOS_INLINE_FUNCTION int getOutputDimensionOfSampling(SamplingFunctional sro, const int local_dimensions)
Dimensions ^ output rank for sampling operation (always in local chart if on a manifold,...
KOKKOS_INLINE_FUNCTION int getInputRankOfSampling(SamplingFunctional sro)
Input rank for sampling operation.
KOKKOS_INLINE_FUNCTION int getSamplingOutputIndex(const SamplingFunctional sf, const int output_component_axis_1, const int output_component_axis_2)
Helper function for finding alpha coefficients.
int output_rank
Rank of sampling functional output for each SamplingFunctional.
int input_rank
Rank of sampling functional input for each SamplingFunctional.
KOKKOS_INLINE_FUNCTION scalar_type & operator[](const int i)
KOKKOS_INLINE_FUNCTION XYZ(double _x=0.0, double _y=0.0, double _z=0.0)
KOKKOS_INLINE_FUNCTION XYZ operator+=(const XYZ &other)
KOKKOS_INLINE_FUNCTION XYZ operator*=(const double &scaling)
KOKKOS_INLINE_FUNCTION XYZ(const scalar_type *arry)
XYZ(const std::vector< scalar_type > &vec)
KOKKOS_INLINE_FUNCTION XYZ operator-=(const XYZ &other)
KOKKOS_INLINE_FUNCTION XYZ operator*(double scalar)
KOKKOS_INLINE_FUNCTION int extent_int(int comp=0)
KOKKOS_INLINE_FUNCTION size_t extent(int comp=0)