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