Compadre  1.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Compadre_DivergenceFreePolynomial.hpp
Go to the documentation of this file.
1 #ifndef _COMPADRE_GMLS_DIVFREE_BASIS_HPP_
2 #define _COMPADRE_GMLS_DIVFREE_BASIS_HPP_
3 
4 #include "Compadre_GMLS.hpp"
5 #include "Compadre_ScalarTaylorPolynomial.hpp"
6 
7 namespace Compadre {
8 /*! \brief Definition of the divergence-free polynomial basis
9  *
10  Let \f$NB_d=\f$ the number of basis components in the scalar Taylor polynomial basis, for the same dimension \f$d\f$ and up to the same degree of polynomial.
11  Let \f$NB_{(d-1)/v}=\f$ the number of basis components in the scalar Taylor polynomial basis, for one dimension lower \f$(d-1)\f$ that does not use same contain the direction v.
12  \li in 1D) No definition.
13  \li in 2D) The first \f$NB_d\f$ elements are \f[\begin{bmatrix}b_i\\ -\int \frac{\partial}{\partial x}(b_i)\; dy\end{bmatrix},\f] where \f$i\f$ is the basis element number for the divergence-free polynomial basis.\newline
14  The next \f$NB_{(d-1)/y}\f$ elements are \f[\begin{bmatrix}0\\b_j\end{bmatrix},\f] where \f$j\f$ is the basis element number for the divergence-free polynomial basis of one dimension lower using only the variable \f$x\f$.
15  \li in 3D) The first \f$NB_d\f$ elements are \f[\begin{bmatrix}b_i\\ 0 \\-\int \frac{\partial}{\partial x}(b_i)\; dz\end{bmatrix},\f] where \f$i\f$ is the basis element number for the divergence-free polynomial basis.\newline
16  The next \f$NB_d\f$ elements are \f[\begin{bmatrix}0\\ b_i \\-\int \frac{\partial}{\partial y}(b_i)\; dz\end{bmatrix},\f] where \f$i\f$ is the basis element number for the divergence-free polynomial basis.\newline
17  The next \f$NB_{(d-1)/z}\f$ elements are \f[\begin{bmatrix}0\\0\\b_j\end{bmatrix},\f] where \f$j\f$ is the basis element number for the divergence-free polynomial basis of one dimension lower using only the variables \f$x\f$ and \f$y\f$.
18 
19 */
20 namespace DivergenceFreePolynomialBasis {
21 
22  /*! \brief Returns size of basis
23  \param degree [in] - highest degree of polynomial
24  \param dimension [in] - spatial dimension to evaluate
25  */
26  KOKKOS_INLINE_FUNCTION
27  int getSize(const int degree, const int dimension) {
28  return (dimension-1)*ScalarTaylorPolynomialBasis::getSize(degree, dimension)
29  + ScalarTaylorPolynomialBasis::getSize(degree, dimension-1);
30  }
31 
32  /*! \brief Evaluates the divergence-free polynomial basis
33  * delta[j] = weight_of_original_value * delta[j] + weight_of_new_value * (calculation of this function)
34  \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.
35  \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.
36  \param dimension [in] - spatial dimension to evaluate
37  \param max_degree [in] - highest degree of polynomial
38  \param h [in] - epsilon/window size
39  \param x [in] - x coordinate (already shifted by target)
40  \param y [in] - y coordinate (already shifted by target)
41  \param z [in] - z coordinate (already shifted by target)
42  \param starting_order [in] - polynomial order to start with (default=0)
43  \param weight_of_original_value [in] - weighting to assign to original value in delta (default=0, replace, 1-increment current value)
44  \param weight_of_new_value [in] - weighting to assign to new contribution (default=1, add to/replace)
45  */
46  KOKKOS_INLINE_FUNCTION
47  void evaluate(const member_type& teamMember, double* delta, double* workspace, const int dimension, const int max_degree, const int component, 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) {
48  Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
49  const double factorial[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
50  if (dimension==3) {
51  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
52  scratch_vector_type y_over_h_to_i(workspace+1*(max_degree+1), max_degree+1);
53  scratch_vector_type z_over_h_to_i(workspace+2*(max_degree+1), max_degree+1);
54  x_over_h_to_i[0] = 1;
55  y_over_h_to_i[0] = 1;
56  z_over_h_to_i[0] = 1;
57  for (int i=1; i<=max_degree; ++i) {
58  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
59  y_over_h_to_i[i] = y_over_h_to_i[i-1]*(y/h);
60  z_over_h_to_i[i] = z_over_h_to_i[i-1]*(z/h);
61  }
62  int i=0;
63  for (int d=0; d<dimension; ++d) {
64  if ((d+1)==dimension) {
65  if (component==d) {
66  // use 2D scalar basis definition
67  // (in 2D) \sum_{n=0}^{n=P} \sum_{k=0}^{k=n} (x/h)^(n-k)*(y/h)^k / ((n-k)!k!)
68  int alphax, alphay;
69  double alphaf;
70  for (int n = starting_order; n <= max_degree; n++){
71  for (alphay = 0; alphay <= n; alphay++){
72  alphax = n - alphay;
73  alphaf = factorial[alphax]*factorial[alphay];
74  *(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;
75  i++;
76  }
77  }
78  } else {
79  for (int j=0; j<ScalarTaylorPolynomialBasis::getSize(max_degree, dimension-1); ++j) {
80  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0;
81  i++;
82  }
83  }
84  } else {
85  if (component==d) {
86  // use 3D scalar basis definition
87  // (in 3D) \sum_{p=0}^{p=P} \sum_{k1+k2+k3=n} (x/h)^k1*(y/h)^k2*(z/h)^k3 / (k1!k2!k3!)
88  int alphax, alphay, alphaz;
89  double alphaf;
90  int s=0;
91  for (int n = starting_order; n <= max_degree; n++){
92  for (alphaz = 0; alphaz <= n; alphaz++){
93  s = n - alphaz;
94  for (alphay = 0; alphay <= s; alphay++){
95  alphax = s - alphay;
96  alphaf = factorial[alphax]*factorial[alphay]*factorial[alphaz];
97  *(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;
98  i++;
99  }
100  }
101  }
102  } else if (component==(d+1)%3) {
103  for (int j=0; j<ScalarTaylorPolynomialBasis::getSize(max_degree, dimension); ++j) {
104  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0;
105  i++;
106  }
107  } else {
108  // (in 3D)
109  int alphax, alphay, alphaz;
110  double alphaf;
111  int s=0;
112  for (int n = starting_order; n <= max_degree; n++){
113  for (alphaz = 0; alphaz <= n; alphaz++){
114  s = n - alphaz;
115  for (alphay = 0; alphay <= s; alphay++){
116  alphax = s - alphay;
117 
118  int var_pow[3] = {(d == 0) ? alphax-1 : alphax, (d == 1) ? alphay-1 : alphay, (d == 2) ? alphaz-1 : alphaz};
119 
120  if (var_pow[0]<0 || var_pow[1]<0 || var_pow[2]<0) {
121  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
122  } else {
123  var_pow[component]++;
124  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]]*factorial[var_pow[2]];
125  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * -1 * 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;
126  }
127  i++;
128  }
129  }
130  }
131  }
132  }
133  }
134  } else if (dimension==2) {
135  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
136  scratch_vector_type y_over_h_to_i(workspace+1*(max_degree+1), max_degree+1);
137  x_over_h_to_i[0] = 1;
138  y_over_h_to_i[0] = 1;
139  for (int i=1; i<=max_degree; ++i) {
140  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
141  y_over_h_to_i[i] = y_over_h_to_i[i-1]*(y/h);
142  }
143  int i=0;
144  for (int d=0; d<dimension; ++d) {
145  if ((d+1)==dimension) {
146  if (component==d) {
147  // use 1D scalar basis definition
148  // (in 1D) \sum_{n=0}^{n=P} (x/h)^n / n!
149  for (int j=starting_order; j<=max_degree; ++j) {
150  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * x_over_h_to_i[j]/factorial[j];
151  i++;
152  }
153  } else {
154  for (int j=0; j<ScalarTaylorPolynomialBasis::getSize(max_degree, dimension-1); ++j) {
155  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0;
156  i++;
157  }
158  }
159  } else {
160  if (component==d) {
161  // use 2D scalar basis definition
162  // (in 2D) \sum_{n=0}^{n=P} \sum_{k=0}^{k=n} (x/h)^(n-k)*(y/h)^k / ((n-k)!k!)
163  int alphax, alphay;
164  double alphaf;
165  for (int n = starting_order; n <= max_degree; n++){
166  for (alphay = 0; alphay <= n; alphay++){
167  alphax = n - alphay;
168  alphaf = factorial[alphax]*factorial[alphay];
169  *(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;
170  i++;
171  }
172  }
173  } else {
174  // (in 2D)
175  int alphax, alphay;
176  double alphaf;
177  for (int n = starting_order; n <= max_degree; n++){
178  for (alphay = 0; alphay <= n; alphay++){
179  alphax = n - alphay;
180 
181  int var_pow[2] = {(d == 0) ? alphax-1 : alphax, (d == 1) ? alphay-1 : alphay};
182 
183  if (var_pow[0]<0 || var_pow[1]<0) {
184  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
185  } else {
186  var_pow[component]++;
187  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]];
188  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * -1 * x_over_h_to_i[var_pow[0]]*y_over_h_to_i[var_pow[1]]/alphaf;
189  }
190  i++;
191  }
192  }
193  }
194  }
195  }
196 
197  } else {
198  compadre_kernel_assert_release((false) && "Divergence-free basis only defined or dimensions 2 and 3.");
199  }
200  });
201  }
202 
203  /*! \brief Evaluates the first partial derivatives of the divergence-free polynomial basis
204  * delta[j] = weight_of_original_value * delta[j] + weight_of_new_value * (calculation of this function)
205  \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.
206  \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.
207  \param dimension [in] - spatial dimension to evaluate
208  \param max_degree [in] - highest degree of polynomial
209  \param partial_direction [in] - direction in which to take partial derivative
210  \param h [in] - epsilon/window size
211  \param x [in] - x coordinate (already shifted by target)
212  \param y [in] - y coordinate (already shifted by target)
213  \param z [in] - z coordinate (already shifted by target)
214  \param starting_order [in] - polynomial order to start with (default=0)
215  \param weight_of_original_value [in] - weighting to assign to original value in delta (default=0, replace, 1-increment current value)
216  \param weight_of_new_value [in] - weighting to assign to new contribution (default=1, add to/replace)
217  */
218  KOKKOS_INLINE_FUNCTION
219  void evaluatePartialDerivative(const member_type& teamMember, double* delta, double* workspace, const int dimension, const int max_degree, const int component, 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) {
220  Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
221  const double factorial[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
222  if (dimension==3) {
223  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
224  scratch_vector_type y_over_h_to_i(workspace+1*(max_degree+1), max_degree+1);
225  scratch_vector_type z_over_h_to_i(workspace+2*(max_degree+1), max_degree+1);
226  x_over_h_to_i[0] = 1;
227  y_over_h_to_i[0] = 1;
228  z_over_h_to_i[0] = 1;
229  for (int i=1; i<=max_degree; ++i) {
230  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
231  y_over_h_to_i[i] = y_over_h_to_i[i-1]*(y/h);
232  z_over_h_to_i[i] = z_over_h_to_i[i-1]*(z/h);
233  }
234  int i=0;
235  for (int d=0; d<dimension; ++d) {
236  if ((d+1)==dimension) {
237  if (component==d) {
238  // use 2D partial derivative of scalar basis definition
239  // (in 2D) \sum_{n=0}^{n=P} \sum_{k=0}^{k=n} (x/h)^(n-k)*(y/h)^k / ((n-k)!k!)
240  int alphax, alphay;
241  double alphaf;
242  for (int n = starting_order; n <= max_degree; n++){
243  for (alphay = 0; alphay <= n; alphay++){
244  alphax = n - alphay;
245 
246  int var_pow[3] = {alphax, alphay, 0};
247  var_pow[partial_direction]--;
248 
249  if (var_pow[0]<0 || var_pow[1]<0 || var_pow[2]<0) {
250  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
251  } else {
252  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]];
253  *(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;
254  }
255  i++;
256  }
257  }
258  } else {
259  for (int j=0; j<ScalarTaylorPolynomialBasis::getSize(max_degree, dimension-1); ++j) {
260  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0;
261  i++;
262  }
263  }
264  } else {
265  if (component==d) {
266  // use 3D partial derivative of scalar basis definition
267  // (in 3D) \sum_{p=0}^{p=P} \sum_{k1+k2+k3=n} (x/h)^k1*(y/h)^k2*(z/h)^k3 / (k1!k2!k3!)
268  int alphax, alphay, alphaz;
269  double alphaf;
270  int s=0;
271  for (int n = starting_order; n <= max_degree; n++){
272  for (alphaz = 0; alphaz <= n; alphaz++){
273  s = n - alphaz;
274  for (alphay = 0; alphay <= s; alphay++){
275  alphax = s - alphay;
276 
277  int var_pow[3] = {alphax, alphay, alphaz};
278  var_pow[partial_direction]--;
279 
280  if (var_pow[0]<0 || var_pow[1]<0 || var_pow[2]<0) {
281  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
282  } else {
283  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]]*factorial[var_pow[2]];
284  *(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;
285  }
286  i++;
287  }
288  }
289  }
290  } else if (component==(d+1)%3) {
291  for (int j=0; j<ScalarTaylorPolynomialBasis::getSize(max_degree, dimension); ++j) {
292  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0;
293  i++;
294  }
295  } else {
296  // (in 3D)
297  int alphax, alphay, alphaz;
298  double alphaf;
299  int s=0;
300  for (int n = starting_order; n <= max_degree; n++){
301  for (alphaz = 0; alphaz <= n; alphaz++){
302  s = n - alphaz;
303  for (alphay = 0; alphay <= s; alphay++){
304  alphax = s - alphay;
305 
306  int var_pow[3] = {alphax, alphay, alphaz};
307  var_pow[d]--;
308 
309  if (var_pow[0]<0 || var_pow[1]<0 || var_pow[2]<0) {
310  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
311  } else {
312  var_pow[component]++;
313  var_pow[partial_direction]--;
314  if (var_pow[0]<0 || var_pow[1]<0 || var_pow[2]<0) {
315  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
316  } else {
317  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]]*factorial[var_pow[2]];
318  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * -1.0/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;
319  }
320  }
321  i++;
322  }
323  }
324  }
325  }
326  }
327  }
328  } else if (dimension==2) {
329  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
330  scratch_vector_type y_over_h_to_i(workspace+1*(max_degree+1), max_degree+1);
331  x_over_h_to_i[0] = 1;
332  y_over_h_to_i[0] = 1;
333  for (int i=1; i<=max_degree; ++i) {
334  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
335  y_over_h_to_i[i] = y_over_h_to_i[i-1]*(y/h);
336  }
337  int i=0;
338  for (int d=0; d<dimension; ++d) {
339  if ((d+1)==dimension) {
340  if (component==d) {
341  // use 1D partial derivative of scalar basis definition
342  int alphax;
343  double alphaf;
344  for (int n = starting_order; n <= max_degree; n++){
345  alphax = n;
346 
347  int var_pow[2] = {alphax, 0};
348  var_pow[partial_direction]--;
349 
350  if (var_pow[0]<0 || var_pow[1]<0) {
351  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
352  } else {
353  alphaf = factorial[var_pow[0]];
354  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 1./h * x_over_h_to_i[var_pow[0]]/alphaf;
355  }
356  i++;
357  }
358  } else {
359  for (int j=0; j<ScalarTaylorPolynomialBasis::getSize(max_degree, dimension-1); ++j) {
360  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0;
361  i++;
362  }
363  }
364  } else {
365  if (component==d) {
366  // use 2D partial derivative of scalar basis definition
367  int alphax, alphay;
368  double alphaf;
369  for (int n = starting_order; n <= max_degree; n++){
370  for (alphay = 0; alphay <= n; alphay++){
371  alphax = n - alphay;
372 
373  int var_pow[2] = {alphax, alphay};
374  var_pow[partial_direction]--;
375 
376  if (var_pow[0]<0 || var_pow[1]<0) {
377  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
378  } else {
379  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]];
380  *(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;
381  }
382  i++;
383  }
384  }
385  } else {
386  // (in 2D)
387  int alphax, alphay;
388  double alphaf;
389  for (int n = starting_order; n <= max_degree; n++){
390  for (alphay = 0; alphay <= n; alphay++){
391  alphax = n - alphay;
392 
393  int var_pow[2] = {alphax, alphay};
394  var_pow[d]--;
395 
396  if (var_pow[0]<0 || var_pow[1]<0) {
397  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
398  } else {
399  var_pow[component]++;
400  var_pow[partial_direction]--;
401  if (var_pow[0]<0 || var_pow[1]<0) {
402  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
403  } else {
404  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]];
405  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * -1.0/h * x_over_h_to_i[var_pow[0]]*y_over_h_to_i[var_pow[1]]/alphaf;
406  }
407  }
408  i++;
409  }
410  }
411  }
412  }
413  }
414 
415  } else {
416  compadre_kernel_assert_release((false) && "Divergence-free basis only defined or dimensions 2 and 3.");
417  }
418  });
419  }
420 
421  /*! \brief Evaluates the second partial derivatives of the divergence-free polynomial basis
422  * delta[j] = weight_of_original_value * delta[j] + weight_of_new_value * (calculation of this function)
423  \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.
424  \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.
425  \param dimension [in] - spatial dimension to evaluate
426  \param max_degree [in] - highest degree of polynomial
427  \param partial_direction_1 [in] - direction in which to take first partial derivative
428  \param partial_direction_2 [in] - direction in which to take second partial derivative
429  \param h [in] - epsilon/window size
430  \param x [in] - x coordinate (already shifted by target)
431  \param y [in] - y coordinate (already shifted by target)
432  \param z [in] - z coordinate (already shifted by target)
433  \param starting_order [in] - polynomial order to start with (default=0)
434  \param weight_of_original_value [in] - weighting to assign to original value in delta (default=0, replace, 1-increment current value)
435  \param weight_of_new_value [in] - weighting to assign to new contribution (default=1, add to/replace)
436  */
437  KOKKOS_INLINE_FUNCTION
438  void evaluateSecondPartialDerivative(const member_type& teamMember, double* delta, double* workspace, const int dimension, const int max_degree, const int component, 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) {
439  Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
440  const double factorial[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
441  if (dimension==3) {
442  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
443  scratch_vector_type y_over_h_to_i(workspace+1*(max_degree+1), max_degree+1);
444  scratch_vector_type z_over_h_to_i(workspace+2*(max_degree+1), max_degree+1);
445  x_over_h_to_i[0] = 1;
446  y_over_h_to_i[0] = 1;
447  z_over_h_to_i[0] = 1;
448  for (int i=1; i<=max_degree; ++i) {
449  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
450  y_over_h_to_i[i] = y_over_h_to_i[i-1]*(y/h);
451  z_over_h_to_i[i] = z_over_h_to_i[i-1]*(z/h);
452  }
453  int i=0;
454  for (int d=0; d<dimension; ++d) {
455  if ((d+1)==dimension) {
456  if (component==d) {
457  // use 2D partial derivative of scalar basis definition
458  // (in 2D) \sum_{n=0}^{n=P} \sum_{k=0}^{k=n} (x/h)^(n-k)*(y/h)^k / ((n-k)!k!)
459  int alphax, alphay;
460  double alphaf;
461  for (int n = starting_order; n <= max_degree; n++){
462  for (alphay = 0; alphay <= n; alphay++){
463  alphax = n - alphay;
464 
465  int var_pow[3] = {alphax, alphay, 0};
466  var_pow[partial_direction_1]--;
467  var_pow[partial_direction_2]--;
468 
469  if (var_pow[0]<0 || var_pow[1]<0 || var_pow[2]<0) {
470  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
471  } else {
472  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]];
473  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 1./h/h * x_over_h_to_i[var_pow[0]]*y_over_h_to_i[var_pow[1]]/alphaf;
474  }
475  i++;
476  }
477  }
478  } else {
479  for (int j=0; j<ScalarTaylorPolynomialBasis::getSize(max_degree, dimension-1); ++j) {
480  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0;
481  i++;
482  }
483  }
484  } else {
485  if (component==d) {
486  // use 3D partial derivative of scalar basis definition
487  // (in 3D) \sum_{p=0}^{p=P} \sum_{k1+k2+k3=n} (x/h)^k1*(y/h)^k2*(z/h)^k3 / (k1!k2!k3!)
488  int alphax, alphay, alphaz;
489  double alphaf;
490  int s=0;
491  for (int n = starting_order; n <= max_degree; n++){
492  for (alphaz = 0; alphaz <= n; alphaz++){
493  s = n - alphaz;
494  for (alphay = 0; alphay <= s; alphay++){
495  alphax = s - alphay;
496 
497  int var_pow[3] = {alphax, alphay, alphaz};
498  var_pow[partial_direction_1]--;
499  var_pow[partial_direction_2]--;
500 
501  if (var_pow[0]<0 || var_pow[1]<0 || var_pow[2]<0) {
502  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
503  } else {
504  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]]*factorial[var_pow[2]];
505  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 1./h/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;
506  }
507  i++;
508  }
509  }
510  }
511  } else if (component==(d+1)%3) {
512  for (int j=0; j<ScalarTaylorPolynomialBasis::getSize(max_degree, dimension); ++j) {
513  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0;
514  i++;
515  }
516  } else {
517  // (in 3D)
518  int alphax, alphay, alphaz;
519  double alphaf;
520  int s=0;
521  for (int n = starting_order; n <= max_degree; n++){
522  for (alphaz = 0; alphaz <= n; alphaz++){
523  s = n - alphaz;
524  for (alphay = 0; alphay <= s; alphay++){
525  alphax = s - alphay;
526 
527  int var_pow[3] = {alphax, alphay, alphaz};
528  var_pow[d]--;
529 
530  if (var_pow[0]<0 || var_pow[1]<0 || var_pow[2]<0) {
531  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
532  } else {
533  var_pow[component]++;
534  var_pow[partial_direction_1]--;
535  var_pow[partial_direction_2]--;
536  if (var_pow[0]<0 || var_pow[1]<0 || var_pow[2]<0) {
537  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
538  } else {
539  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]]*factorial[var_pow[2]];
540  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * -1.0/h/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;
541  }
542  }
543  i++;
544  }
545  }
546  }
547  }
548  }
549  }
550  } else if (dimension==2) {
551  scratch_vector_type x_over_h_to_i(workspace, max_degree+1);
552  scratch_vector_type y_over_h_to_i(workspace+1*(max_degree+1), max_degree+1);
553  x_over_h_to_i[0] = 1;
554  y_over_h_to_i[0] = 1;
555  for (int i=1; i<=max_degree; ++i) {
556  x_over_h_to_i[i] = x_over_h_to_i[i-1]*(x/h);
557  y_over_h_to_i[i] = y_over_h_to_i[i-1]*(y/h);
558  }
559  int i=0;
560  for (int d=0; d<dimension; ++d) {
561  if ((d+1)==dimension) {
562  if (component==d) {
563  // use 1D partial derivative of scalar basis definition
564  int alphax;
565  double alphaf;
566  for (int n = starting_order; n <= max_degree; n++){
567  alphax = n;
568 
569  int var_pow[2] = {alphax, 0};
570  var_pow[partial_direction_1]--;
571  var_pow[partial_direction_2]--;
572 
573  if (var_pow[0]<0 || var_pow[1]<0) {
574  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
575  } else {
576  alphaf = factorial[var_pow[0]];
577  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 1./h/h * x_over_h_to_i[var_pow[0]]/alphaf;
578  }
579  i++;
580  }
581  } else {
582  for (int j=0; j<ScalarTaylorPolynomialBasis::getSize(max_degree, dimension-1); ++j) {
583  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0;
584  i++;
585  }
586  }
587  } else {
588  if (component==d) {
589  // use 2D partial derivative of scalar basis definition
590  int alphax, alphay;
591  double alphaf;
592  for (int n = starting_order; n <= max_degree; n++){
593  for (alphay = 0; alphay <= n; alphay++){
594  alphax = n - alphay;
595 
596  int var_pow[2] = {alphax, alphay};
597  var_pow[partial_direction_1]--;
598  var_pow[partial_direction_2]--;
599 
600  if (var_pow[0]<0 || var_pow[1]<0) {
601  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
602  } else {
603  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]];
604  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 1./h/h * x_over_h_to_i[var_pow[0]]*y_over_h_to_i[var_pow[1]]/alphaf;
605  }
606  i++;
607  }
608  }
609  } else {
610  // (in 2D)
611  int alphax, alphay;
612  double alphaf;
613  for (int n = starting_order; n <= max_degree; n++){
614  for (alphay = 0; alphay <= n; alphay++){
615  alphax = n - alphay;
616 
617  int var_pow[2] = {alphax, alphay};
618  var_pow[d]--;
619 
620  if (var_pow[0]<0 || var_pow[1]<0) {
621  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
622  } else {
623  var_pow[component]++;
624  var_pow[partial_direction_1]--;
625  var_pow[partial_direction_2]--;
626  if (var_pow[0]<0 || var_pow[1]<0) {
627  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * 0.0;
628  } else {
629  alphaf = factorial[var_pow[0]]*factorial[var_pow[1]];
630  *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * -1.0/h/h * x_over_h_to_i[var_pow[0]]*y_over_h_to_i[var_pow[1]]/alphaf;
631  }
632  }
633  i++;
634  }
635  }
636  }
637  }
638  }
639 
640  } else {
641  compadre_kernel_assert_release((false) && "Divergence-free basis only defined or dimensions 2 and 3.");
642  }
643  });
644  }
645 
646  /*! \brief Evaluates the hard-coded divergence-free polynomial basis up to order 4 for 3D
647  *
648  Polynomial degree is inferred from np
649  \warning This function is deprecated.
650  \param np [in] - basis component number
651  \param sx [in] - x coordinate (already shifted by target and scaled)
652  \param sy [in] - y coordinate (already shifted by target and scaled)
653  \param sz [in] - z coordinate (already shifted by target and scaled)
654  */
655  KOKKOS_INLINE_FUNCTION
656  XYZ evaluate(int np, const double sx, const double sy, const double sz) {
657  // define one-third
658  const double th = 1./3.;
659 
660  switch (np) {
661  // P0
662  case 0: return XYZ(1.0,0.0,0.0);
663  case 1: return XYZ(0.0,1.0,0.0);
664  case 2: return XYZ(0.0,0.0,1.0);
665 
666  // P1
667  case 3 : return XYZ( sy, 0, 0);
668  case 4 : return XYZ( sz, 0, 0);
669  case 5 : return XYZ( 0, sx, 0);
670  case 6 : return XYZ( 0, sz, 0);
671  case 7 : return XYZ( 0, 0, sx);
672  case 8 : return XYZ( 0, 0, sy);
673  case 9 : return XYZ( sx,-sy, 0);
674  case 10: return XYZ( sx, 0,-sz);
675 
676  // P2
677  case 11: return XYZ( sy*sy, 0, 0);
678  case 12: return XYZ( sz*sz, 0, 0);
679  case 13: return XYZ( sy*sz, 0, 0);
680  case 14: return XYZ( 0, sx*sx, 0);
681  case 15: return XYZ( 0, sz*sz, 0);
682  case 16: return XYZ( 0, sx*sz, 0);
683  case 17: return XYZ( 0, 0, sx*sx);
684  case 18: return XYZ( 0, 0, sy*sy);
685  case 19: return XYZ( 0, 0, sx*sy);
686 
687  case 20: return XYZ( sx*sx, -2.0*sx*sy, 0);
688  case 21: return XYZ( sx*sy, 0,-1.0*sy*sz);
689  case 22: return XYZ( sx*sz, -1.0*sy*sz, 0);
690  case 23: return XYZ(-2.0*sx*sy, sy*sy, 0);
691  case 24: return XYZ( 0, sx*sy, -1.0*sx*sz);
692  case 25: return XYZ(-2.0*sx*sz, 0, sz*sz);
693 
694  // P3
695  case 26: return XYZ(sy*sy*sy , 0, 0);
696  case 27: return XYZ(sy*sy*sz , 0, 0);
697  case 28: return XYZ(sy*sz*sz , 0, 0);
698  case 29: return XYZ(sz*sz*sz , 0, 0);
699  case 30: return XYZ( 0, sx*sx*sx , 0);
700  case 31: return XYZ( 0, sx*sx*sz , 0);
701  case 32: return XYZ( 0, sx*sz*sz , 0);
702  case 33: return XYZ( 0, sz*sz*sz , 0);
703  case 34: return XYZ( 0, 0, sx*sx*sx );
704  case 35: return XYZ( 0, 0, sx*sx*sy );
705  case 36: return XYZ( 0, 0, sx*sy*sy );
706  case 37: return XYZ( 0, 0, sy*sy*sy );
707 
708  case 38: return XYZ( sx*sx*sx ,-3.0*sx*sx*sy, 0);
709  case 39: return XYZ( sx*sx*sx , 0,-3.0*sx*sx*sz);
710  case 40: return XYZ( sx*sx*sy ,-1.0*sx*sy*sy, 0);
711  case 41: return XYZ( sx*sx*sy , 0,-2.0*sx*sy*sz);
712  case 42: return XYZ( sx*sx*sz ,-2.0*sx*sy*sz, 0);
713  case 43: return XYZ( sx*sx*sz , 0,-1.0*sx*sz*sz);
714  case 44: return XYZ( sx*sy*sy ,- th*sy*sy*sy, 0);
715  case 45: return XYZ( sx*sy*sy , 0,-1.0*sy*sy*sz);
716  case 46: return XYZ( sx*sy*sz ,-0.5*sy*sy*sz, 0);
717  case 47: return XYZ( sx*sy*sz , 0,-0.5*sy*sz*sz);
718  case 48: return XYZ( sx*sz*sz ,-1.0*sy*sz*sz, 0);
719  case 49: return XYZ( sx*sz*sz , 0,- th*sz*sz*sz);
720 
721  // P4
722  case 50: return XYZ(sy*sy*sy*sy , 0, 0);
723  case 51: return XYZ(sy*sy*sy*sz , 0, 0);
724  case 52: return XYZ(sy*sy*sz*sz , 0, 0);
725  case 53: return XYZ(sy*sz*sz*sz , 0, 0);
726  case 54: return XYZ(sz*sz*sz*sz , 0, 0);
727  case 55: return XYZ( 0, sx*sx*sx*sx , 0);
728  case 56: return XYZ( 0, sx*sx*sx*sz , 0);
729  case 57: return XYZ( 0, sx*sx*sz*sz , 0);
730  case 58: return XYZ( 0, sx*sz*sz*sz , 0);
731  case 59: return XYZ( 0, sz*sz*sz*sz , 0);
732  case 60: return XYZ( 0, 0, sx*sx*sx*sx );
733  case 61: return XYZ( 0, 0, sx*sx*sx*sy );
734  case 62: return XYZ( 0, 0, sx*sx*sy*sy );
735  case 63: return XYZ( 0, 0, sx*sy*sy*sy );
736  case 64: return XYZ( 0, 0, sy*sy*sy*sy );
737 
738  case 65: return XYZ(sx*sx*sx*sx ,-4.0*sx*sx*sx*sy , 0);
739  case 66: return XYZ(sx*sx*sx*sx , 0, -4.0*sx*sx*sx*sz );
740  case 67: return XYZ(sx*sx*sx*sy ,-1.5*sx*sx*sy*sy , 0);
741  case 68: return XYZ(sx*sx*sx*sy , 0, -3.0*sx*sx*sy*sz );
742  case 69: return XYZ(sx*sx*sx*sz ,-3.0*sx*sx*sy*sz , 0);
743  case 70: return XYZ(sx*sx*sx*sz , 0, -1.5*sx*sx*sz*sz );
744 
745  case 71: return XYZ(sx*sx*sy*sy , -2.0*th*sx*sy*sy*sy, 0);
746  case 72: return XYZ(sx*sx*sy*sy , 0, -2.0*sx*sy*sy*sz );
747  case 73: return XYZ(sx*sx*sy*sz , -sx*sy*sy*sz , 0);
748  case 74: return XYZ(sx*sx*sy*sz , 0, -sx*sy*sz*sz );
749  case 75: return XYZ(sx*sx*sz*sz , -2.0*sx*sy*sz*sz , 0);
750  case 76: return XYZ(sx*sx*sz*sz , 0, -2.0*th*sx*sz*sz*sz);
751  case 77: return XYZ(sx*sy*sy*sy , -0.25*sy*sy*sy*sy , 0);
752  case 78: return XYZ(sx*sy*sy*sy , 0, -sy*sy*sy*sz );
753  case 79: return XYZ(sx*sy*sy*sz , -th*sy*sy*sy*sz , 0);
754  case 80: return XYZ(sx*sy*sy*sz , 0, -0.5*sy*sy*sz*sz );
755  case 81: return XYZ(sx*sy*sz*sz , -0.5*sy*sy*sz*sz , 0);
756  case 82: return XYZ(sx*sy*sz*sz , 0, -th*sy*sz*sz*sz );
757  case 83: return XYZ(sx*sz*sz*sz , -sy*sz*sz*sz , 0);
758  case 84: return XYZ(sx*sz*sz*sz , 0,-0.25*sz*sz*sz*sz );
759 
760  // default
761  default: compadre_kernel_assert_release((false) && "Divergence-free basis only sup\
762 ports up to 4th-order polynomials for now.");
763  return XYZ(0,0,0); // avoid warning about no return
764  }
765  }
766 
767  /*! \brief Evaluates the hard-coded divergence-free polynomial basis up to order 4 for 2D
768  *
769  Polynomial degree is inferred from np
770  \warning This function is deprecated.
771  \param np [in] - basis component number
772  \param sx [in] - x coordinate (already shifted by target and scaled)
773  \param sy [in] - y coordinate (already shifted by target and scaled)
774  */
775  KOKKOS_INLINE_FUNCTION
776  XYZ evaluate(int np, const double sx, const double sy) {
777  // define one-third
778 
779  switch (np) {
780  // P0
781  case 0: return XYZ(1.0,0.0,0.0);
782  case 1: return XYZ(0.0,1.0,0.0);
783 
784  // P1
785  case 2: return XYZ(sy, 0, 0);
786  case 3: return XYZ( 0, sx, 0);
787  case 4: return XYZ(sx,-sy, 0);
788 
789  // P2
790  case 5: return XYZ( sy*sy, 0, 0);
791  case 6: return XYZ( 0, sx*sx, 0);
792 
793  case 7: return XYZ( sx*sx, -2.0*sx*sy, 0);
794  case 8: return XYZ( -2.0*sx*sy, sy*sy, 0);
795 
796  // P3
797  case 9 : return XYZ( sy*sy*sy, 0, 0);
798  case 10: return XYZ( 0, sx*sx*sx, 0);
799 
800  case 11: return XYZ( sx*sx*sx, -3.0*sx*sx*sy, 0);
801  case 12: return XYZ( sx*sx*sy, -1.0*sx*sy*sy, 0);
802  case 13: return XYZ(-3.0*sx*sy*sy, sy*sy*sy, 0);
803 
804  // P4
805  case 14: return XYZ( sy*sy*sy*sy, 0, 0);
806  case 15: return XYZ( 0, sx*sx*sx*sx, 0);
807 
808  case 16: return XYZ( sx*sx*sx*sx, -4.0*sx*sx*sx*sy, 0);
809  case 17: return XYZ( 2.0*sx*sx*sx*sy, -3.0*sx*sx*sy*sy, 0);
810  case 18: return XYZ(-3.0*sx*sx*sy*sy, 2.0*sx*sy*sy*sy, 0);
811  case 19: return XYZ(-4.0*sx*sy*sy*sy, sy*sy*sy*sy, 0);
812 
813  // default
814  default: compadre_kernel_assert_release((false) && "Divergence-free basis only sup\
815 ports up to 4th-order polynomials for now.");
816  return XYZ(0,0,0); // avoid warning about no return
817  }
818  }
819 
820  /*! \brief Evaluates the hard-coded first partial derivative of the divergence-free polynomial basis up to order 4 for 3D
821  *
822  Polynomial degree is inferred from np
823  \warning This function is deprecated.
824  \param np [in] - basis component number
825  \param partial_direction [in] - partial derivative direction
826  \param sx [in] - x coordinate (already shifted by target and scaled)
827  \param sy [in] - y coordinate (already shifted by target and scaled)
828  \param sz [in] - z coordinate (already shifted by target and scaled)
829  */
830  KOKKOS_INLINE_FUNCTION
831  XYZ evaluatePartialDerivative(int np, const int partial_direction, const double h, const double sx, const double sy, const double sz) {
832  // define one-third
833  const double th = 1./3.;
834 
835  if (partial_direction==0) {
836  switch (np) {
837  // P0
838  case 0: return XYZ(0.0,0.0,0.0);
839  case 1: return XYZ(0.0,0.0,0.0);
840  case 2: return XYZ(0.0,0.0,0.0);
841 
842  // P1
843  case 3 : return XYZ( 0, 0, 0);
844  case 4 : return XYZ( 0, 0, 0);
845  case 5 : return XYZ( 0, 1.0/h, 0);
846  case 6 : return XYZ( 0, 0, 0);
847  case 7 : return XYZ( 0, 0, 1.0/h);
848  case 8 : return XYZ( 0, 0, 0);
849  case 9 : return XYZ( 1.0/h, 0, 0);
850  case 10: return XYZ( 1.0/h, 0, 0);
851 
852  // P2
853  case 11: return XYZ( 0, 0, 0);
854  case 12: return XYZ( 0, 0, 0);
855  case 13: return XYZ( 0, 0, 0);
856  case 14: return XYZ( 0, 2*sx/h, 0);
857  case 15: return XYZ( 0, 0, 0);
858  case 16: return XYZ( 0, sz/h, 0);
859  case 17: return XYZ( 0, 0, 2*sx/h);
860  case 18: return XYZ( 0, 0, 0);
861  case 19: return XYZ( 0, 0, sy/h);
862 
863  case 20: return XYZ( 2*sx/h, -2.0*sy/h, 0);
864  case 21: return XYZ( sy/h, 0, 0);
865  case 22: return XYZ( sz/h, 0, 0);
866  case 23: return XYZ( -2.0*sy/h, 0, 0);
867  case 24: return XYZ( 0, sy/h, -1.0*sz/h);
868  case 25: return XYZ( -2.0*sz/h, 0, 0);
869 
870  // P3
871  case 26: return XYZ( 0, 0, 0);
872  case 27: return XYZ( 0, 0, 0);
873  case 28: return XYZ( 0, 0, 0);
874  case 29: return XYZ( 0, 0, 0);
875  case 30: return XYZ( 0, 3*sx*sx/h, 0);
876  case 31: return XYZ( 0, 2*sx*sz/h, 0);
877  case 32: return XYZ( 0, sz*sz/h, 0);
878  case 33: return XYZ( 0, 0, 0);
879  case 34: return XYZ( 0, 0, 3*sx*sx/h);
880  case 35: return XYZ( 0, 0, 2*sx*sy/h);
881  case 36: return XYZ( 0, 0, sy*sy/h);
882  case 37: return XYZ( 0, 0, 0);
883 
884  case 38: return XYZ( 3*sx*sx/h, -6.0*sx*sy/h, 0);
885  case 39: return XYZ( 3*sx*sx/h, 0, -6.0*sx*sz/h);
886  case 40: return XYZ( 2*sx*sy/h, -1.0*sy*sy/h, 0);
887  case 41: return XYZ( 2*sx*sy/h, 0, -2.0*sy*sz/h);
888  case 42: return XYZ( 2*sx*sz/h, -2.0*sy*sz/h, 0);
889  case 43: return XYZ( 2*sx*sz/h, 0, -1.0*sz*sz/h);
890  case 44: return XYZ( sy*sy/h, 0, 0);
891  case 45: return XYZ( sy*sy/h, 0, 0);
892  case 46: return XYZ( sy*sz/h, 0, 0);
893  case 47: return XYZ( sy*sz/h, 0, 0);
894  case 48: return XYZ( sz*sz/h, 0, 0);
895  case 49: return XYZ( sz*sz/h, 0, 0);
896 
897  // P4
898  case 50: return XYZ( 0, 0, 0);
899  case 51: return XYZ( 0, 0, 0);
900  case 52: return XYZ( 0, 0, 0);
901  case 53: return XYZ( 0, 0, 0);
902  case 54: return XYZ( 0, 0, 0);
903  case 55: return XYZ( 0, 4*sx*sx*sx/h, 0);
904  case 56: return XYZ( 0, 3*sx*sx*sz/h, 0);
905  case 57: return XYZ( 0, 2*sx*sz*sz/h, 0);
906  case 58: return XYZ( 0, sz*sz*sz/h, 0);
907  case 59: return XYZ( 0, 0, 0);
908  case 60: return XYZ( 0, 0, 4*sx*sx*sx/h);
909  case 61: return XYZ( 0, 0, 3*sx*sx*sy/h);
910  case 62: return XYZ( 0, 0, 2*sx*sy*sy/h);
911  case 63: return XYZ( 0, 0, sy*sy*sy/h);
912  case 64: return XYZ( 0, 0, 0);
913 
914  case 65: return XYZ( 4*sx*sx*sx/h ,-12.0*sx*sx*sy/h , 0);
915  case 66: return XYZ( 4*sx*sx*sx/h , 0, -12.0*sx*sx*sz/h );
916  case 67: return XYZ( 3*sx*sx*sy/h ,-3.0*sx*sy*sy/h , 0);
917  case 68: return XYZ( 3*sx*sx*sy/h , 0, -6.0*sx*sy*sz/h );
918  case 69: return XYZ( 3*sx*sx*sz/h ,-6.0*sx*sy*sz/h , 0);
919  case 70: return XYZ( 3*sx*sx*sz/h , 0, -3.0*sx*sz*sz/h );
920  case 71: return XYZ( 2*sx*sy*sy/h , -2.0*th*sy*sy*sy/h, 0);
921  case 72: return XYZ( 2*sx*sy*sy/h , 0, -2.0*sy*sy*sz/h );
922  case 73: return XYZ( 2*sx*sy*sz/h , -sy*sy*sz/h , 0);
923  case 74: return XYZ( 2*sx*sy*sz/h , 0, -sy*sz*sz/h );
924  case 75: return XYZ( 2*sx*sz*sz/h , -2.0*sy*sz*sz/h , 0);
925  case 76: return XYZ( 2*sx*sz*sz/h , 0, -2.0*th*sz*sz*sz/h);
926  case 77: return XYZ( sy*sy*sy/h , 0, 0);
927  case 78: return XYZ( sy*sy*sy/h , 0, 0);
928  case 79: return XYZ( sy*sy*sz/h , 0, 0);
929  case 80: return XYZ( sy*sy*sz/h , 0, 0);
930  case 81: return XYZ( sy*sz*sz/h , 0, 0);
931  case 82: return XYZ( sy*sz*sz/h , 0, 0);
932  case 83: return XYZ( sz*sz*sz/h , 0, 0);
933  case 84: return XYZ( sz*sz*sz/h , 0, 0);
934 
935  // default
936  default: compadre_kernel_assert_release((false) && "Divergence-free basis only sup\
937 ports up to 4th-order polynomials for now.");
938  return XYZ(0,0,0); // avoid warning about no return
939  }
940  } else if (partial_direction==1) {
941  switch (np) {
942  // P0
943  case 0: return XYZ(0.0,0.0,0.0);
944  case 1: return XYZ(0.0,0.0,0.0);
945  case 2: return XYZ(0.0,0.0,0.0);
946 
947  // P1
948  case 3 : return XYZ( 1.0/h, 0, 0);
949  case 4 : return XYZ( 0, 0, 0);
950  case 5 : return XYZ( 0, 0, 0);
951  case 6 : return XYZ( 0, 0, 0);
952  case 7 : return XYZ( 0, 0, 0);
953  case 8 : return XYZ( 0, 0, 1.0/h);
954  case 9 : return XYZ( 0, -1.0/h, 0);
955  case 10: return XYZ( 0, 0, 0);
956 
957  // P2
958  case 11: return XYZ( 2*sy/h, 0, 0);
959  case 12: return XYZ( 0, 0, 0);
960  case 13: return XYZ( sz/h, 0, 0);
961  case 14: return XYZ( 0, 0, 0);
962  case 15: return XYZ( 0, 0, 0);
963  case 16: return XYZ( 0, 0, 0);
964  case 17: return XYZ( 0, 0, 0);
965  case 18: return XYZ( 0, 0, 2*sy/h);
966  case 19: return XYZ( 0, 0, sx/h);
967 
968  case 20: return XYZ( 0, -2.0*sx/h, 0);
969  case 21: return XYZ( sx/h, 0, -1.0*sz/h);
970  case 22: return XYZ( 0, -1.0*sz/h, 0);
971  case 23: return XYZ( -2.0*sx/h, 2*sy/h, 0);
972  case 24: return XYZ( 0, sx/h, 0);
973  case 25: return XYZ( 0, 0, 0);
974 
975  // P3
976  case 26: return XYZ( 3*sy*sy/h, 0, 0);
977  case 27: return XYZ( 2*sy*sz/h, 0, 0);
978  case 28: return XYZ( sz*sz/h, 0, 0);
979  case 29: return XYZ( 0, 0, 0);
980  case 30: return XYZ( 0, 0, 0);
981  case 31: return XYZ( 0, 0, 0);
982  case 32: return XYZ( 0, 0, 0);
983  case 33: return XYZ( 0, 0, 0);
984  case 34: return XYZ( 0, 0, 0);
985  case 35: return XYZ( 0, 0, sx*sx/h);
986  case 36: return XYZ( 0, 0, 2*sx*sy/h);
987  case 37: return XYZ( 0, 0, 3*sy*sy/h);
988 
989  case 38: return XYZ( 0, -3.0*sx*sx/h, 0);
990  case 39: return XYZ( 0, 0, 0);
991  case 40: return XYZ( sx*sx/h, -2.0*sx*sy/h, 0);
992  case 41: return XYZ( sx*sx/h, 0, -2.0*sx*sz/h);
993  case 42: return XYZ( 0, -2.0*sx*sz/h, 0);
994  case 43: return XYZ( 0, 0, 0);
995  case 44: return XYZ( 2*sx*sy/h,-3*th*sy*sy/h, 0);
996  case 45: return XYZ( 2*sx*sy/h, 0, -2.0*sy*sz/h);
997  case 46: return XYZ( sx*sz/h, -1.0*sy*sz/h, 0);
998  case 47: return XYZ( sx*sz/h, 0, -0.5*sz*sz/h);
999  case 48: return XYZ( 0, -1.0*sz*sz/h, 0);
1000  case 49: return XYZ( 0, 0, 0);
1001 
1002  // P4
1003  case 50: return XYZ( 4*sy*sy*sy/h, 0, 0);
1004  case 51: return XYZ( 3*sy*sy*sz/h, 0, 0);
1005  case 52: return XYZ( 2*sy*sz*sz/h, 0, 0);
1006  case 53: return XYZ( sz*sz*sz/h, 0, 0);
1007  case 54: return XYZ( 0, 0, 0);
1008  case 55: return XYZ( 0, 0, 0);
1009  case 56: return XYZ( 0, 0, 0);
1010  case 57: return XYZ( 0, 0, 0);
1011  case 58: return XYZ( 0, 0, 0);
1012  case 59: return XYZ( 0, 0, 0);
1013  case 60: return XYZ( 0, 0, 0);
1014  case 61: return XYZ( 0, 0, sx*sx*sx/h);
1015  case 62: return XYZ( 0, 0, 2*sx*sx*sy/h);
1016  case 63: return XYZ( 0, 0, 3*sx*sy*sy/h);
1017  case 64: return XYZ( 0, 0, 4*sy*sy*sy/h);
1018 
1019  case 65: return XYZ( 0, -4.0*sx*sx*sx/h, 0);
1020  case 66: return XYZ( 0, 0, 0);
1021  case 67: return XYZ( sx*sx*sx/h, -3.0*sx*sx*sy/h, 0);
1022  case 68: return XYZ( sx*sx*sx/h, 0, -3.0*sx*sx*sz/h);
1023  case 69: return XYZ( 0, -3.0*sx*sx*sz/h, 0);
1024  case 70: return XYZ( 0, 0, 0);
1025  case 71: return XYZ( 2*sx*sx*sy/h, -6.0*th*sx*sy*sy/h, 0);
1026  case 72: return XYZ( 2*sx*sx*sy/h, 0, -4.0*sx*sy*sz/h);
1027  case 73: return XYZ( sx*sx*sz/h, -2*sx*sy*sz, 0);
1028  case 74: return XYZ( sx*sx*sz/h, 0, -sx*sz*sz/h);
1029  case 75: return XYZ( 0, -2.0*sx*sz*sz/h, 0);
1030  case 76: return XYZ( 0, 0, 0);
1031  case 77: return XYZ( 3*sx*sy*sy/h, -1.0*sy*sy*sy/h, 0);
1032  case 78: return XYZ( 3*sx*sy*sy/h, 0, -3*sy*sy*sz/h);
1033  case 79: return XYZ( 2*sx*sy*sz/h, -3*th*sy*sy*sz/h, 0);
1034  case 80: return XYZ( 2*sx*sy*sz/h, 0, -1.0*sy*sz*sz/h);
1035  case 81: return XYZ( sx*sz*sz/h, -1.0*sy*sz*sz/h, 0);
1036  case 82: return XYZ( sx*sz*sz/h, 0, -th*sz*sz*sz/h);
1037  case 83: return XYZ( 0, -sz*sz*sz/h, 0);
1038  case 84: return XYZ( 0, 0, 0);
1039 
1040  // default
1041  default: compadre_kernel_assert_release((false) && "Divergence-free basis only sup\
1042 ports up to 4th-order polynomials for now.");
1043  return XYZ(0,0,0); // avoid warning about no return
1044  }
1045  } else {
1046  switch (np) {
1047  // P0
1048  case 0: return XYZ(0.0,0.0,0.0);
1049  case 1: return XYZ(0.0,0.0,0.0);
1050  case 2: return XYZ(0.0,0.0,0.0);
1051 
1052  // P1
1053  case 3 : return XYZ( 0, 0, 0);
1054  case 4 : return XYZ( 1/h, 0, 0);
1055  case 5 : return XYZ( 0, 0, 0);
1056  case 6 : return XYZ( 0, 1/h, 0);
1057  case 7 : return XYZ( 0, 0, 0);
1058  case 8 : return XYZ( 0, 0, 0);
1059  case 9 : return XYZ( 0, 0, 0);
1060  case 10: return XYZ( 0, 0,-1/h);
1061 
1062  // P2
1063  case 11: return XYZ( 0, 0, 0);
1064  case 12: return XYZ( 2*sz/h, 0, 0);
1065  case 13: return XYZ( sy/h, 0, 0);
1066  case 14: return XYZ( 0, 0, 0);
1067  case 15: return XYZ( 0, 2*sz/h, 0);
1068  case 16: return XYZ( 0, sx/h, 0);
1069  case 17: return XYZ( 0, 0, 0);
1070  case 18: return XYZ( 0, 0, 0);
1071  case 19: return XYZ( 0, 0, 0);
1072 
1073  case 20: return XYZ( 0, 0, 0);
1074  case 21: return XYZ( 0, 0, -1.0*sy/h);
1075  case 22: return XYZ( sx/h, -1.0*sy/h, 0);
1076  case 23: return XYZ( 0, 0, 0);
1077  case 24: return XYZ( 0, 0, -1.0*sx/h);
1078  case 25: return XYZ( -2.0*sx/h, 0, 2*sz/h);
1079 
1080  // P3
1081  case 26: return XYZ( 0, 0, 0);
1082  case 27: return XYZ( sy*sy/h, 0, 0);
1083  case 28: return XYZ( 2*sy*sz/h, 0, 0);
1084  case 29: return XYZ( 3*sz*sz/h, 0, 0);
1085  case 30: return XYZ( 0, 0, 0);
1086  case 31: return XYZ( 0, sx*sx/h, 0);
1087  case 32: return XYZ( 0, 2*sx*sz/h, 0);
1088  case 33: return XYZ( 0, 3*sz*sz/h, 0);
1089  case 34: return XYZ( 0, 0, 0);
1090  case 35: return XYZ( 0, 0, 0);
1091  case 36: return XYZ( 0, 0, 0);
1092  case 37: return XYZ( 0, 0, 0);
1093 
1094  case 38: return XYZ( 0, 0, 0);
1095  case 39: return XYZ( 0, 0, -3.0*sx*sx/h);
1096  case 40: return XYZ( 0, 0, 0);
1097  case 41: return XYZ( 0, 0, -2.0*sx*sy/h);
1098  case 42: return XYZ( sx*sx/h, -2.0*sx*sy/h, 0);
1099  case 43: return XYZ( sx*sx/h, 0, -2.0*sx*sz/h);
1100  case 44: return XYZ( 0, 0, 0);
1101  case 45: return XYZ( 0, 0, -1.0*sy*sy/h);
1102  case 46: return XYZ( sx*sy/h, -0.5*sy*sy/h, 0);
1103  case 47: return XYZ( sx*sy/h, 0, -1.0*sy*sz/h);
1104  case 48: return XYZ( 2*sx*sz/h, -2.0*sy*sz/h, 0);
1105  case 49: return XYZ( 2*sx*sz/h, 0,-3*th*sz*sz/h);
1106 
1107  // P4
1108  case 50: return XYZ( 0, 0, 0);
1109  case 51: return XYZ( sy*sy*sy/h, 0, 0);
1110  case 52: return XYZ( 2*sy*sy*sz/h, 0, 0);
1111  case 53: return XYZ( 3*sy*sz*sz/h, 0, 0);
1112  case 54: return XYZ( 4*sz*sz*sz/h, 0, 0);
1113  case 55: return XYZ( 0, 0, 0);
1114  case 56: return XYZ( 0, sx*sx*sx/h, 0);
1115  case 57: return XYZ( 0, 2*sx*sx*sz/h, 0);
1116  case 58: return XYZ( 0, 3*sx*sz*sz/h, 0);
1117  case 59: return XYZ( 0, 4*sz*sz*sz/h, 0);
1118  case 60: return XYZ( 0, 0, 0);
1119  case 61: return XYZ( 0, 0, 0);
1120  case 62: return XYZ( 0, 0, 0);
1121  case 63: return XYZ( 0, 0, 0);
1122  case 64: return XYZ( 0, 0, 0);
1123 
1124  case 65: return XYZ( 0, 0, 0);
1125  case 66: return XYZ( 0, 0, -4.0*sx*sx*sx/h);
1126  case 67: return XYZ( 0, 0, 0);
1127  case 68: return XYZ( 0, 0, -3.0*sx*sx*sy/h);
1128  case 69: return XYZ( sx*sx*sx/h, -3.0*sx*sx*sy/h, 0);
1129  case 70: return XYZ( sx*sx*sx/h, 0, -3.0*sx*sx*sz/h);
1130  case 71: return XYZ( 0, 0, 0);
1131  case 72: return XYZ( 0, 0, -2.0*sx*sy*sy/h);
1132  case 73: return XYZ( sx*sx*sy/h, -sx*sy*sy/h, 0);
1133  case 74: return XYZ( sx*sx*sy/h, 0, -2*sx*sy*sz/h);
1134  case 75: return XYZ( 2*sx*sx*sz/h, -4.0*sx*sy*sz/h, 0);
1135  case 76: return XYZ( 2*sx*sx*sz/h, 0, -6.0*th*sx*sz*sz/h);
1136  case 77: return XYZ( 0, 0, 0);
1137  case 78: return XYZ( 0, 0, -sy*sy*sy/h);
1138  case 79: return XYZ( sx*sy*sy/h, -th*sy*sy*sy/h, 0);
1139  case 80: return XYZ( sx*sy*sy/h, 0, -1.0*sy*sy*sz/h);
1140  case 81: return XYZ( 2*sx*sy*sz/h, -1.0*sy*sy*sz/h, 0);
1141  case 82: return XYZ( 2*sx*sy*sz/h, 0, -3*th*sy*sz*sz/h);
1142  case 83: return XYZ( 3*sx*sz*sz/h, -3*sy*sz*sz/h, 0);
1143  case 84: return XYZ( 3*sx*sz*sz/h, 0, -1.0*sz*sz*sz/h);
1144 
1145  // default
1146  default: compadre_kernel_assert_release((false) && "Divergence-free basis only sup\
1147 ports up to 4th-order polynomials for now.");
1148  return XYZ(0,0,0); // avoid warning about no return
1149  }
1150  }
1151  }
1152 
1153  /*! \brief Evaluates the hard-coded first partial derivative of the divergence-free polynomial basis up to order 4 for 2D
1154  *
1155  Polynomial degree is inferred from np
1156  \warning This function is deprecated.
1157  \param np [in] - basis component number
1158  \param partial_direction [in] - partial derivative direction
1159  \param sx [in] - x coordinate (already shifted by target and scaled)
1160  \param sy [in] - y coordinate (already shifted by target and scaled)
1161  */
1162  KOKKOS_INLINE_FUNCTION
1163  XYZ evaluatePartialDerivative(int np, const int partial_direction, const double h, const double sx, const double sy) {
1164  // define one-third
1165 
1166  if (partial_direction==0) {
1167  switch (np) {
1168  // P0
1169  case 0: return XYZ(0.0,0.0,0.0);
1170  case 1: return XYZ(0.0,0.0,0.0);
1171 
1172  // P1
1173  case 2: return XYZ( 0, 0, 0);
1174  case 3: return XYZ( 0, 1.0/h, 0);
1175  case 4: return XYZ( 1.0/h, 0, 0);
1176 
1177  // P2
1178  case 5: return XYZ( 0, 0, 0);
1179  case 6: return XYZ( 0, 2*sx/h, 0);
1180 
1181  case 7: return XYZ( 2*sx/h, -2.0*sy/h, 0);
1182  case 8: return XYZ( -2.0*sy/h, 0, 0);
1183 
1184  // P3
1185  case 9 : return XYZ( 0, 0, 0);
1186  case 10: return XYZ( 0, 3*sx*sx/h, 0);
1187 
1188  case 11: return XYZ( 3*sx*sx/h, -6.0*sx*sy/h, 0);
1189  case 12: return XYZ( 2.0*sx*sy/h, -1.0*sy*sy/h, 0);
1190  case 13: return XYZ( -3.0*sy*sy/h, 0, 0);
1191 
1192  // P4
1193  case 14: return XYZ( 0, 0, 0);
1194  case 15: return XYZ( 0, 4*sx*sx*sx/h, 0);
1195 
1196  case 16: return XYZ( 4*sx*sx*sx/h, -12.0*sx*sx*sy/h, 0);
1197  case 17: return XYZ( 6.0*sx*sx*sy/h, -6.0*sx*sy*sy/h, 0);
1198  case 18: return XYZ( -6.0*sx*sy*sy/h, 2.0*sy*sy*sy/h, 0);
1199  case 19: return XYZ( -4.0*sy*sy*sy/h, 0, 0);
1200 
1201  // default
1202  default: compadre_kernel_assert_release((false) && "Divergence-free basis only sup\
1203 ports up to 4th-order polynomials for now.");
1204  return XYZ(0,0,0); // avoid warning about no return
1205  }
1206  } else {
1207  switch (np) {
1208  // P0
1209  case 0: return XYZ(0.0,0.0,0.0);
1210  case 1: return XYZ(0.0,0.0,0.0);
1211 
1212  // P1
1213  case 2: return XYZ( 1.0/h, 0, 0);
1214  case 3: return XYZ( 0, 0, 0);
1215  case 4: return XYZ( 0, -1.0/h, 0);
1216 
1217  // P2
1218  case 5: return XYZ( 2*sy/h, 0, 0);
1219  case 6: return XYZ( 0, 0, 0);
1220 
1221  case 7: return XYZ( 0, -2.0*sx/h, 0);
1222  case 8: return XYZ( -2.0*sx/h, 2*sy/h, 0);
1223 
1224  // P3
1225  case 9 : return XYZ( 3*sy*sy/h, 0, 0);
1226  case 10: return XYZ( 0, 0, 0);
1227 
1228  case 11: return XYZ( 0, -3.0*sx*sx/h, 0);
1229  case 12: return XYZ( sx*sx/h, -2.0*sx*sy/h, 0);
1230  case 13: return XYZ( -6.0*sx*sy/h, 3*sy*sy/h, 0);
1231 
1232  // P4
1233  case 14: return XYZ( 4*sy*sy*sy/h, 0, 0);
1234  case 15: return XYZ( 0, 0, 0);
1235 
1236  case 16: return XYZ( 0, -4.0*sx*sx*sx/h, 0);
1237  case 17: return XYZ( 2.0*sx*sx*sx/h, -6.0*sx*sx*sy/h, 0);
1238  case 18: return XYZ( -6.0*sx*sx*sy/h, 6.0*sx*sy*sy/h, 0);
1239  case 19: return XYZ( -12.0*sx*sy*sy/h, 4*sy*sy*sy/h, 0);
1240 
1241  // default
1242  default: compadre_kernel_assert_release((false) && "Divergence-free basis only sup\
1243 ports up to 4th-order polynomials for now.");
1244  return XYZ(0,0,0); // avoid warning about no return
1245  }
1246  }
1247  }
1248 } // DivergenceFreePolynomialBasis
1249 } // Compadre
1250 
1251 #endif
Compadre::ScalarTaylorPolynomialBasis::getSize
KOKKOS_INLINE_FUNCTION int getSize(const int degree, const int dimension)
Returns size of basis.
Definition: Compadre_ScalarTaylorPolynomial.hpp:21
compadre_kernel_assert_release
#define compadre_kernel_assert_release(condition)
compadre_kernel_assert_release is similar to compadre_assert_release, but is a call on the device...
Definition: Compadre_Typedefs.hpp:129
Compadre::DivergenceFreePolynomialBasis::evaluate
KOKKOS_INLINE_FUNCTION void evaluate(const member_type &teamMember, double *delta, double *workspace, const int dimension, const int max_degree, const int component, 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 divergence-free polynomial basis delta[j] = weight_of_original_value * delta[j] + weigh...
Definition: Compadre_DivergenceFreePolynomial.hpp:47
member_type
team_policy::member_type member_type
Definition: Compadre_Typedefs.hpp:48
Compadre::DivergenceFreePolynomialBasis::evaluateSecondPartialDerivative
KOKKOS_INLINE_FUNCTION void evaluateSecondPartialDerivative(const member_type &teamMember, double *delta, double *workspace, const int dimension, const int max_degree, const int component, 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 the divergence-free polynomial basis delta[j] = weight_of...
Definition: Compadre_DivergenceFreePolynomial.hpp:438
Compadre_GMLS.hpp
Compadre::DivergenceFreePolynomialBasis::getSize
KOKKOS_INLINE_FUNCTION int getSize(const int degree, const int dimension)
Returns size of basis.
Definition: Compadre_DivergenceFreePolynomial.hpp:27
Compadre::DivergenceFreePolynomialBasis::evaluatePartialDerivative
KOKKOS_INLINE_FUNCTION void evaluatePartialDerivative(const member_type &teamMember, double *delta, double *workspace, const int dimension, const int max_degree, const int component, 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 the divergence-free polynomial basis delta[j] = weight_of_...
Definition: Compadre_DivergenceFreePolynomial.hpp:219
scratch_vector_type
Kokkos::View< double *, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_vector_type
Definition: Compadre_Typedefs.hpp:63
Compadre_ScalarTaylorPolynomial.hpp
d
import subprocess import os import re import math import sys import argparse d d d(20 *num_target_sites *pow(4, grid_num))
Generated on Fri Nov 22 2024 09:26:40 for Compadre by   doxygen 1.8.5