Compadre  1.3.3
Compadre_Misc.hpp
Go to the documentation of this file.
1 #ifndef _COMPADRE_MISC_HPP_
2 #define _COMPADRE_MISC_HPP_
3 
4 #include "Compadre_Operators.hpp"
5 
6 namespace Compadre {
7 
8 struct XYZ {
9 
10  double x;
11  double y;
12  double z;
13 
14  KOKKOS_INLINE_FUNCTION
15  XYZ(double _x = 0.0, double _y = 0.0, double _z = 0.0) : x(_x), y(_y), z(_z) {}
16 
17  KOKKOS_INLINE_FUNCTION
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 
22  KOKKOS_INLINE_FUNCTION
23  XYZ operator += (const XYZ& other)
24  { x += other.x;
25  y += other.y;
26  z += other.z;
27  return *this; }
28  KOKKOS_INLINE_FUNCTION
30  { x += other.x;
31  y += other.y;
32  z += other.z;
33  return *this; }
34  KOKKOS_INLINE_FUNCTION
35  XYZ operator -= (const XYZ& other)
36  { x -= other.x;
37  y -= other.y;
38  z -= other.z;
39  return *this; }
40  KOKKOS_INLINE_FUNCTION
42  { x -= other.x;
43  y -= other.y;
44  z -= other.z;
45  return *this; }
46  KOKKOS_INLINE_FUNCTION
47  XYZ operator *= (const double& scaling)
48  { x *= scaling;
49  y *= scaling;
50  z *= scaling;
51  return *this; }
52  KOKKOS_INLINE_FUNCTION
53  scalar_type& operator [](const int i) {
54  switch (i) {
55  case 0:
56  return x;
57  case 1:
58  return y;
59  default:
60  return z;
61  }
62  }
63  KOKKOS_INLINE_FUNCTION
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  }
74  KOKKOS_INLINE_FUNCTION
75  XYZ operator *(double scalar) {
76  XYZ result;
77  result.x = scalar*x;
78  result.y = scalar*y;
79  result.z = scalar*z;
80  return result;
81  }
82 }; // XYZ
83 
84 
85 KOKKOS_INLINE_FUNCTION
86 XYZ operator + ( const XYZ& vecA, const XYZ& vecB ) {
87  return XYZ( vecA.x + vecB.x, vecA.y + vecB.y, vecA.z + vecB.z); }
88 
89  KOKKOS_INLINE_FUNCTION
90 XYZ operator - ( const XYZ& vecA, const XYZ& vecB ) {
91  return XYZ( vecA.x - vecB.x, vecA.y - vecB.y, vecA.z - vecB.z); }
92 
93 KOKKOS_INLINE_FUNCTION
94 XYZ operator * ( const XYZ& vecA, const XYZ& vecB ) {
95  return XYZ( vecA.x * vecB.x, vecA.y * vecB.y, vecA.z * vecB.z); }
96 
97 KOKKOS_INLINE_FUNCTION
98 XYZ operator + ( const XYZ& vecA, const scalar_type& constant ) {
99  return XYZ( vecA.x + constant, vecA.y + constant, vecA.z + constant); }
100 
101 KOKKOS_INLINE_FUNCTION
102 XYZ operator + ( const scalar_type& constant, const XYZ& vecA ) {
103  return XYZ( vecA.x + constant, vecA.y + constant, vecA.z + constant); }
104 
105 KOKKOS_INLINE_FUNCTION
106 XYZ operator - ( const XYZ& vecA, const scalar_type& constant ) {
107  return XYZ( vecA.x - constant, vecA.y - constant, vecA.z - constant); }
108 
109 KOKKOS_INLINE_FUNCTION
110 XYZ operator - ( const scalar_type& constant, const XYZ& vecA ) {
111  return XYZ( constant - vecA.x, constant - vecA.y, constant - vecA.z); }
112 
113 KOKKOS_INLINE_FUNCTION
114 XYZ operator * ( const XYZ& vecA, const scalar_type& constant ) {
115  return XYZ( vecA.x * constant, vecA.y * constant, vecA.z * constant); }
116 
117 KOKKOS_INLINE_FUNCTION
118 XYZ operator * (const scalar_type& constant, const XYZ& vecA) {
119  return XYZ( vecA.x * constant, vecA.y * constant, vecA.z * constant); }
120 
121 KOKKOS_INLINE_FUNCTION
122 XYZ operator / ( const XYZ& vecA, const scalar_type& constant ) {
123  return XYZ( vecA.x / constant, vecA.y / constant, vecA.z / constant); }
124 
125 inline std::ostream& operator << ( std::ostream& os, const XYZ& vec ) {
126  os << "(" << vec.x << ", " << vec.y << ", " << vec.z << ")" ; return os; }
127 
128 KOKKOS_INLINE_FUNCTION
129 int getAdditionalAlphaSizeFromConstraint(DenseSolverType dense_solver_type, ConstraintType constraint_type) {
130  // Return the additional constraint size
131  if (constraint_type == NEUMANN_GRAD_SCALAR) {
132  return 1;
133  } else {
134  return 0;
135  }
136 }
137 
138 KOKKOS_INLINE_FUNCTION
139 int getAdditionalCoeffSizeFromConstraintAndSpace(DenseSolverType dense_solver_type, ConstraintType constraint_type, ReconstructionSpace reconstruction_space, const int dimension) {
140  // Return the additional constraint size
141  int added_alpha_size = getAdditionalAlphaSizeFromConstraint(dense_solver_type, constraint_type);
142  if (reconstruction_space == VectorTaylorPolynomial) {
143  return dimension*added_alpha_size;
144  } else {
145  return added_alpha_size;
146  }
147 }
148 
149 KOKKOS_INLINE_FUNCTION
150 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) {
151  // Return the appropriate size for _RHS. Since in LU, the system solves P^T*P against P^T*W.
152  // We store P^T*P in the RHS space, which means RHS can be much smaller compared to the
153  // case for QR/SVD where the system solves PsqrtW against sqrtW*Identity
154 
155  int added_coeff_size = getAdditionalCoeffSizeFromConstraintAndSpace(dense_solver_type, constraint_type, reconstruction_space, dimension);
156 
157  if (constraint_type == NEUMANN_GRAD_SCALAR) {
158  RHS_row = RHS_col = N + added_coeff_size;
159  } else {
160  if (dense_solver_type != LU) {
161  RHS_row = N;
162  RHS_col = M;
163  } else {
164  RHS_row = RHS_col = N + added_coeff_size;
165  }
166  }
167 }
168 
169 KOKKOS_INLINE_FUNCTION
170 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) {
171  // Return the appropriate size for _P.
172  // In the case of solving with LU and additional constraint is used, _P needs
173  // to be resized to include additional row(s) based on the type of constraint.
174 
175  int added_coeff_size = getAdditionalCoeffSizeFromConstraintAndSpace(dense_solver_type, constraint_type, reconstruction_space, dimension);
176  int added_alpha_size = getAdditionalAlphaSizeFromConstraint(dense_solver_type, constraint_type);
177 
178  if (constraint_type == NEUMANN_GRAD_SCALAR) {
179  out_row = M + added_alpha_size;
180  out_col = N + added_coeff_size;
181  } else {
182  if (dense_solver_type == LU) {
183  out_row = M + added_alpha_size;
184  out_col = N + added_coeff_size;
185  } else {
186  out_row = M;
187  out_col = N;
188  }
189  }
190 }
191 
192 } // Compadre
193 
194 #endif
KOKKOS_INLINE_FUNCTION int getAdditionalAlphaSizeFromConstraint(DenseSolverType dense_solver_type, ConstraintType constraint_type)
double scalar_type
KOKKOS_INLINE_FUNCTION scalar_type & operator[](const int i)
KOKKOS_INLINE_FUNCTION int getAdditionalCoeffSizeFromConstraintAndSpace(DenseSolverType dense_solver_type, ConstraintType constraint_type, ReconstructionSpace reconstruction_space, const int dimension)
Neumann Gradient Scalar Type.
KOKKOS_INLINE_FUNCTION XYZ operator-(const XYZ &vecA, const XYZ &vecB)
KOKKOS_INLINE_FUNCTION XYZ operator+(const XYZ &vecA, const XYZ &vecB)
KOKKOS_INLINE_FUNCTION XYZ operator+=(const XYZ &other)
KOKKOS_INLINE_FUNCTION XYZ(double _x=0.0, double _y=0.0, double _z=0.0)
KOKKOS_INLINE_FUNCTION XYZ operator*(double scalar)
KOKKOS_INLINE_FUNCTION XYZ operator*(const XYZ &vecA, const XYZ &vecB)
XYZ(const std::vector< scalar_type > &vec)
KOKKOS_INLINE_FUNCTION XYZ(const scalar_type *arry)
ReconstructionSpace
Space in which to reconstruct polynomial.
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 XYZ operator-=(const XYZ &other)
std::ostream & operator<<(std::ostream &os, const XYZ &vec)
DenseSolverType
Dense solver type.
LU factorization performed on P^T*W*P matrix.
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 double &scaling)
KOKKOS_INLINE_FUNCTION XYZ operator/(const XYZ &vecA, const scalar_type &constant)
Vector polynomial basis having # of components _dimensions, or (_dimensions-1) in the case of manifol...
ConstraintType
Constraint type.