Compadre 1.5.5
Loading...
Searching...
No Matches
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
6namespace 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*/
14namespace 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::View< double *, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_vector_type
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 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...
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...