Compadre  1.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Compadre_ScalarTaylorPolynomial.hpp
Go to the documentation of this file.
1 #ifndef _COMPADRE_SCALAR_TAYLORPOLYNOMIAL_BASIS_HPP_
2 #define _COMPADRE_SCALAR_TAYLORPOLYNOMIAL_BASIS_HPP_
3 
4 #include "Compadre_GMLS.hpp"
5 
6 namespace Compadre {
7 /*! \brief Definition of scalar Taylor polynomial basis
8 
9  For order P, the sum of the basis is defined by:
10  \li in 1D) \f[\sum_{n=0}^{n=P} \frac{(x/h)^n}{n!}\f]
11  \li in 2D) \f[\sum_{n=0}^{n=P} \sum_{k=0}^{k=n} \frac{(x/h)^{n-k}*(y/h)^k}{(n-k)!k!}\f]
12  \li in 3D) \f[\sum_{p=0}^{p=P} \sum_{k_1+k_2+k_3=n} \frac{(x/h)^{k_1}*(y/h)^{k_2}*(z/h)^{k_3}}{k_1!k_2!k_3!}\f]
13 */
14 namespace ScalarTaylorPolynomialBasis {
15 
16  /*! \brief Returns size of basis
17  \param degree [in] - highest degree of polynomial
18  \param dimension [in] - spatial dimension to evaluate
19  */
20  KOKKOS_INLINE_FUNCTION
21  int getSize(const int degree, const int dimension) {
22  if (dimension == 3) return (degree+1)*(degree+2)*(degree+3)/6;
23  else if (dimension == 2) return (degree+1)*(degree+2)/2;
24  else return degree+1;
25  }
26 
27  /*! \brief Evaluates the scalar Taylor polynomial basis
28  * delta[j] = weight_of_original_value * delta[j] + weight_of_new_value * (calculation of this function)
29  \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.
30  \param 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.
31  \param dimension [in] - spatial dimension to evaluate
32  \param max_degree [in] - highest degree of polynomial
33  \param h [in] - epsilon/window size
34  \param x [in] - x coordinate (already shifted by target)
35  \param y [in] - y coordinate (already shifted by target)
36  \param z [in] - z coordinate (already shifted by target)
37  \param starting_order [in] - polynomial order to start with (default=0)
38  \param weight_of_original_value [in] - weighting to assign to original value in delta (default=0, replace, 1-increment current value)
39  \param weight_of_new_value [in] - weighting to assign to new contribution (default=1, add to/replace)
40  */
41  KOKKOS_INLINE_FUNCTION
42  void evaluate(const member_type& teamMember, double* delta, double* workspace, const int dimension, const int max_degree, const double h, const double x, const double y, const double z, const int starting_order = 0, const double weight_of_original_value = 0.0, const double weight_of_new_value = 1.0) {
43  Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
44  const double factorial[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
45  if (dimension==3) {
46  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
47  scratch_vector_type y_over_h_to_i(workspace+1*(max_degree+1), max_degree+1);
48  scratch_vector_type z_over_h_to_i(workspace+2*(max_degree+1), max_degree+1);
49  x_over_h_to_i[0] = 1;
50  y_over_h_to_i[0] = 1;
51  z_over_h_to_i[0] = 1;
52  for (int i=1; i<=max_degree; ++i) {
53  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
54  y_over_h_to_i[i] = y_over_h_to_i[i-1]*(y/h);
55  z_over_h_to_i[i] = z_over_h_to_i[i-1]*(z/h);
56  }
57  // (in 3D) \sum_{p=0}^{p=P} \sum_{k1+k2+k3=n} (x/h)^k1*(y/h)^k2*(z/h)^k3 / (k1!k2!k3!)
58  int alphax, alphay, alphaz;
59  double alphaf;
60  int i=0, s=0;
61  for (int n = starting_order; n <= max_degree; n++){
62  for (alphaz = 0; alphaz <= n; alphaz++){
63  s = n - alphaz;
64  for (alphay = 0; alphay <= s; alphay++){
65  alphax = s - alphay;
66  alphaf = factorial[alphax]*factorial[alphay]*factorial[alphaz];
67  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * x_over_h_to_i[alphax]*y_over_h_to_i[alphay]*z_over_h_to_i[alphaz]/alphaf;
68  i++;
69  }
70  }
71  }
72  } else if (dimension==2) {
73  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
74  scratch_vector_type y_over_h_to_i(workspace+1*(max_degree+1), max_degree+1);
75  x_over_h_to_i[0] = 1;
76  y_over_h_to_i[0] = 1;
77  for (int i=1; i<=max_degree; ++i) {
78  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
79  y_over_h_to_i[i] = y_over_h_to_i[i-1]*(y/h);
80  }
81  // (in 2D) \sum_{n=0}^{n=P} \sum_{k=0}^{k=n} (x/h)^(n-k)*(y/h)^k / ((n-k)!k!)
82  int alphax, alphay;
83  double alphaf;
84  int i = 0;
85  for (int n = starting_order; n <= max_degree; n++){
86  for (alphay = 0; alphay <= n; alphay++){
87  alphax = n - alphay;
88  alphaf = factorial[alphax]*factorial[alphay];
89  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * x_over_h_to_i[alphax]*y_over_h_to_i[alphay]/alphaf;
90  i++;
91  }
92  }
93  } else {
94  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
95  x_over_h_to_i[0] = 1;
96  for (int i=1; i<=max_degree; ++i) {
97  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
98  }
99  // (in 1D) \sum_{n=0}^{n=P} (x/h)^n / n!
100  for (int i=starting_order; i<=max_degree; ++i) {
101  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * x_over_h_to_i[i]/factorial[i];
102  }
103  }
104  });
105  }
106 
107  /*! \brief Evaluates the first partial derivatives of scalar Taylor polynomial basis
108  * delta[j] = weight_of_original_value * delta[j] + weight_of_new_value * (calculation of this function)
109  \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.
110  \param 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.
111  \param dimension [in] - spatial dimension to evaluate
112  \param max_degree [in] - highest degree of polynomial
113  \param partial_direction [in] - direction in which to take partial derivative
114  \param h [in] - epsilon/window size
115  \param x [in] - x coordinate (already shifted by target)
116  \param y [in] - y coordinate (already shifted by target)
117  \param z [in] - z coordinate (already shifted by target)
118  \param starting_order [in] - polynomial order to start with (default=0)
119  \param weight_of_original_value [in] - weighting to assign to original value in delta (default=0, replace, 1-increment current value)
120  \param weight_of_new_value [in] - weighting to assign to new contribution (default=1, add to/replace)
121  */
122  KOKKOS_INLINE_FUNCTION
123  void evaluatePartialDerivative(const member_type& teamMember, double* delta, double* workspace, const int dimension, const int max_degree, const int partial_direction, const double h, const double x, const double y, const double z, const int starting_order = 0, const double weight_of_original_value = 0.0, const double weight_of_new_value = 1.0) {
124  Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
125  const double factorial[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
126  if (dimension==3) {
127  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
128  scratch_vector_type y_over_h_to_i(workspace+1*(max_degree+1), max_degree+1);
129  scratch_vector_type z_over_h_to_i(workspace+2*(max_degree+1), max_degree+1);
130  x_over_h_to_i[0] = 1;
131  y_over_h_to_i[0] = 1;
132  z_over_h_to_i[0] = 1;
133  for (int i=1; i<=max_degree; ++i) {
134  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
135  y_over_h_to_i[i] = y_over_h_to_i[i-1]*(y/h);
136  z_over_h_to_i[i] = z_over_h_to_i[i-1]*(z/h);
137  }
138  // (in 3D) \sum_{p=0}^{p=P} \sum_{k1+k2+k3=n} (x/h)^k1*(y/h)^k2*(z/h)^k3 / (k1!k2!k3!)
139  int alphax, alphay, alphaz;
140  double alphaf;
141  int i=0, s=0;
142  for (int n = starting_order; n <= max_degree; n++){
143  for (alphaz = 0; alphaz <= n; alphaz++){
144  s = n - alphaz;
145  for (alphay = 0; alphay <= s; alphay++){
146  alphax = s - alphay;
147 
148  int var_pow[3] = {alphax, alphay, alphaz};
149  var_pow[partial_direction]--;
150 
151  if (var_pow[0]<0 || var_pow[1]<0 || var_pow[2]<0) {
152  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
153  } else {
154  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]]*factorial[var_pow[2]];
155  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 1./h * x_over_h_to_i[var_pow[0]]*y_over_h_to_i[var_pow[1]]*z_over_h_to_i[var_pow[2]]/alphaf;
156  }
157  i++;
158  }
159  }
160  }
161  } else if (dimension==2) {
162  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
163  scratch_vector_type y_over_h_to_i(workspace+1*(max_degree+1), max_degree+1);
164  x_over_h_to_i[0] = 1;
165  y_over_h_to_i[0] = 1;
166  for (int i=1; i<=max_degree; ++i) {
167  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
168  y_over_h_to_i[i] = y_over_h_to_i[i-1]*(y/h);
169  }
170  // (in 3D) \sum_{p=0}^{p=P} \sum_{k1+k2+k3=n} (x/h)^k1*(y/h)^k2*(z/h)^k3 / (k1!k2!k3!)
171  int alphax, alphay;
172  double alphaf;
173  int i=0;
174  for (int n = starting_order; n <= max_degree; n++){
175  for (alphay = 0; alphay <= n; alphay++){
176  alphax = n - alphay;
177 
178  int var_pow[2] = {alphax, alphay};
179  var_pow[partial_direction]--;
180 
181  if (var_pow[0]<0 || var_pow[1]<0) {
182  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
183  } else {
184  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]];
185  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 1./h * x_over_h_to_i[var_pow[0]]*y_over_h_to_i[var_pow[1]]/alphaf;
186  }
187  i++;
188  }
189  }
190  } else {
191  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
192  x_over_h_to_i[0] = 1;
193  for (int i=1; i<=max_degree; ++i) {
194  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
195  }
196  int alphax;
197  double alphaf;
198  int i=0;
199  for (int n = starting_order; n <= max_degree; n++){
200  alphax = n;
201 
202  int var_pow[1] = {alphax};
203  var_pow[partial_direction]--;
204 
205  if (var_pow[0]<0) {
206  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
207  } else {
208  alphaf = factorial[var_pow[0]];
209  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 1./h * x_over_h_to_i[var_pow[0]]/alphaf;
210  }
211  i++;
212  }
213  }
214  });
215  }
216 
217  /*! \brief Evaluates the second partial derivatives of scalar Taylor polynomial basis
218  * delta[j] = weight_of_original_value * delta[j] + weight_of_new_value * (calculation of this function)
219  \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.
220  \param 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.
221  \param dimension [in] - spatial dimension to evaluate
222  \param max_degree [in] - highest degree of polynomial
223  \param partial_direction_1 [in] - direction in which to take first partial derivative
224  \param partial_direction_2 [in] - direction in which to take second partial derivative
225  \param h [in] - epsilon/window size
226  \param x [in] - x coordinate (already shifted by target)
227  \param y [in] - y coordinate (already shifted by target)
228  \param z [in] - z coordinate (already shifted by target)
229  \param starting_order [in] - polynomial order to start with (default=0)
230  \param weight_of_original_value [in] - weighting to assign to original value in delta (default=0, replace, 1-increment current value)
231  \param weight_of_new_value [in] - weighting to assign to new contribution (default=1, add to/replace)
232  */
233  KOKKOS_INLINE_FUNCTION
234  void evaluateSecondPartialDerivative(const member_type& teamMember, double* delta, double* workspace, const int dimension, const int max_degree, const int partial_direction_1, const int partial_direction_2, const double h, const double x, const double y, const double z, const int starting_order = 0, const double weight_of_original_value = 0.0, const double weight_of_new_value = 1.0) {
235  Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
236  const double factorial[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
237  if (dimension==3) {
238  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
239  scratch_vector_type y_over_h_to_i(workspace+1*(max_degree+1), max_degree+1);
240  scratch_vector_type z_over_h_to_i(workspace+2*(max_degree+1), max_degree+1);
241  x_over_h_to_i[0] = 1;
242  y_over_h_to_i[0] = 1;
243  z_over_h_to_i[0] = 1;
244  for (int i=1; i<=max_degree; ++i) {
245  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
246  y_over_h_to_i[i] = y_over_h_to_i[i-1]*(y/h);
247  z_over_h_to_i[i] = z_over_h_to_i[i-1]*(z/h);
248  }
249  // (in 3D) \sum_{p=0}^{p=P} \sum_{k1+k2+k3=n} (x/h)^k1*(y/h)^k2*(z/h)^k3 / (k1!k2!k3!)
250  int alphax, alphay, alphaz;
251  double alphaf;
252  int i=0, s=0;
253  for (int n = starting_order; n <= max_degree; n++){
254  for (alphaz = 0; alphaz <= n; alphaz++){
255  s = n - alphaz;
256  for (alphay = 0; alphay <= s; alphay++){
257  alphax = s - alphay;
258 
259  int var_pow[3] = {alphax, alphay, alphaz};
260  var_pow[partial_direction_1]--;
261  var_pow[partial_direction_2]--;
262 
263  if (var_pow[0]<0 || var_pow[1]<0 || var_pow[2]<0) {
264  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
265  } else {
266  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]]*factorial[var_pow[2]];
267  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 1./h * x_over_h_to_i[var_pow[0]]*y_over_h_to_i[var_pow[1]]*z_over_h_to_i[var_pow[2]]/alphaf;
268  }
269  i++;
270  }
271  }
272  }
273  } else if (dimension==2) {
274  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
275  scratch_vector_type y_over_h_to_i(workspace+1*(max_degree+1), max_degree+1);
276  x_over_h_to_i[0] = 1;
277  y_over_h_to_i[0] = 1;
278  for (int i=1; i<=max_degree; ++i) {
279  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
280  y_over_h_to_i[i] = y_over_h_to_i[i-1]*(y/h);
281  }
282  // (in 3D) \sum_{p=0}^{p=P} \sum_{k1+k2+k3=n} (x/h)^k1*(y/h)^k2*(z/h)^k3 / (k1!k2!k3!)
283  int alphax, alphay;
284  double alphaf;
285  int i=0;
286  for (int n = starting_order; n <= max_degree; n++){
287  for (alphay = 0; alphay <= n; alphay++){
288  alphax = n - alphay;
289 
290  int var_pow[2] = {alphax, alphay};
291  var_pow[partial_direction_1]--;
292  var_pow[partial_direction_2]--;
293 
294  if (var_pow[0]<0 || var_pow[1]<0) {
295  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
296  } else {
297  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]];
298  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 1./h * x_over_h_to_i[var_pow[0]]*y_over_h_to_i[var_pow[1]]/alphaf;
299  }
300  i++;
301  }
302  }
303  } else {
304  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
305  x_over_h_to_i[0] = 1;
306  for (int i=1; i<=max_degree; ++i) {
307  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
308  }
309  int alphax;
310  double alphaf;
311  int i=0;
312  for (int n = starting_order; n <= max_degree; n++){
313  alphax = n;
314 
315  int var_pow[1] = {alphax};
316  var_pow[partial_direction_1]--;
317  var_pow[partial_direction_2]--;
318 
319  if (var_pow[0]<0) {
320  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
321  } else {
322  alphaf = factorial[var_pow[0]];
323  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 1./h * x_over_h_to_i[var_pow[0]]/alphaf;
324  }
325  i++;
326  }
327  }
328  });
329  }
330 
331 } // ScalarTaylorPolynomialBasis
332 
333 } // Compadre
334 
335 #endif
KOKKOS_INLINE_FUNCTION int getSize(const int degree, const int dimension)
Returns size of basis.
KOKKOS_INLINE_FUNCTION void evaluatePartialDerivative(const member_type &teamMember, double *delta, double *workspace, const int dimension, const int max_degree, const int partial_direction, const double h, const double x, const double y, const double z, const int starting_order=0, const double weight_of_original_value=0.0, const double weight_of_new_value=1.0)
Evaluates the first partial derivatives of scalar Taylor polynomial basis delta[j] = weight_of_origin...
team_policy::member_type member_type
KOKKOS_INLINE_FUNCTION void evaluate(const member_type &teamMember, double *delta, double *workspace, const int dimension, const int max_degree, const double h, const double x, const double y, const double z, const int starting_order=0, const double weight_of_original_value=0.0, const double weight_of_new_value=1.0)
Evaluates the scalar Taylor polynomial basis delta[j] = weight_of_original_value * delta[j] + weight_...
KOKKOS_INLINE_FUNCTION void evaluateSecondPartialDerivative(const member_type &teamMember, double *delta, double *workspace, const int dimension, const int max_degree, const int partial_direction_1, const int partial_direction_2, const double h, const double x, const double y, const double z, const int starting_order=0, const double weight_of_original_value=0.0, const double weight_of_new_value=1.0)
Evaluates the second partial derivatives of scalar Taylor polynomial basis delta[j] = weight_of_origi...
Kokkos::View< double *, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_vector_type