Compadre 1.5.5
Loading...
Searching...
No Matches
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"
6
7namespace 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*/
20namespace 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\
762ports 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\
815ports 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\
937ports 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\
1042ports 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\
1147ports 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\
1203ports 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\
1243ports 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
Kokkos::View< double *, Kokkos::MemoryTraits< Kokkos::Unmanaged > > scratch_vector_type
#define compadre_kernel_assert_release(condition)
compadre_kernel_assert_release is similar to compadre_assert_release, but is a call on the device,...
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 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...
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_...
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...
KOKKOS_INLINE_FUNCTION int getSize(const int degree, const int dimension)
Returns size of basis.
KOKKOS_INLINE_FUNCTION int getSize(const int degree, const int dimension)
Returns size of basis.