1 #ifndef _COMPADRE_GMLS_BERNSTEIN_BASIS_HPP_
2 #define _COMPADRE_GMLS_BERNSTEIN_BASIS_HPP_
15 namespace BernsteinPolynomialBasis {
21 KOKKOS_INLINE_FUNCTION
22 int getSize(
const int degree,
const int dimension) {
40 KOKKOS_INLINE_FUNCTION
41 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) {
42 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
43 const double factorial[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
54 one_minus_x_over_h_to_i[0] = 1;
55 one_minus_y_over_h_to_i[0] = 1;
56 one_minus_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]*(0.5*x/h+0.5);
59 y_over_h_to_i[i] = y_over_h_to_i[i-1]*(0.5*y/h+0.5);
60 z_over_h_to_i[i] = z_over_h_to_i[i-1]*(0.5*z/h+0.5);
61 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));
62 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));
63 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));
71 for (
int n1 = starting_order; n1 <= max_degree; n1++){
72 int degree_choose_n1 = factorial[max_degree]/(factorial[n1]*factorial[max_degree-n1]);
73 for (
int n2 = starting_order; n2 <= max_degree; n2++){
74 int degree_choose_n2 = factorial[max_degree]/(factorial[n2]*factorial[max_degree-n2]);
75 for (
int n3 = starting_order; n3 <= max_degree; n3++){
76 int degree_choose_n3 = factorial[max_degree]/(factorial[n3]*factorial[max_degree-n3]);
77 *(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];
82 }
else if (dimension==2) {
89 one_minus_x_over_h_to_i[0] = 1;
90 one_minus_y_over_h_to_i[0] = 1;
91 for (
int i=1; i<=max_degree; ++i) {
92 x_over_h_to_i[i] = x_over_h_to_i[i-1]*(0.5*x/h+0.5);
93 y_over_h_to_i[i] = y_over_h_to_i[i-1]*(0.5*y/h+0.5);
94 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));
95 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));
102 for (
int n1 = starting_order; n1 <= max_degree; n1++){
103 int degree_choose_n1 = factorial[max_degree]/(factorial[n1]*factorial[max_degree-n1]);
104 for (
int n2 = starting_order; n2 <= max_degree; n2++){
105 int degree_choose_n2 = factorial[max_degree]/(factorial[n2]*factorial[max_degree-n2]);
106 *(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];
113 x_over_h_to_i[0] = 1;
114 one_minus_x_over_h_to_i[0] = 1;
115 for (
int i=1; i<=max_degree; ++i) {
116 x_over_h_to_i[i] = x_over_h_to_i[i-1]*(0.5*x/h+0.5);
117 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));
122 for (
int n=starting_order; n<=max_degree; ++n) {
123 int degree_choose_n = factorial[max_degree]/(factorial[n]*factorial[max_degree-n]);
124 *(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];
146 KOKKOS_INLINE_FUNCTION
147 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) {
148 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
149 const double factorial[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200};
157 x_over_h_to_i[0] = 1;
158 y_over_h_to_i[0] = 1;
159 z_over_h_to_i[0] = 1;
160 one_minus_x_over_h_to_i[0] = 1;
161 one_minus_y_over_h_to_i[0] = 1;
162 one_minus_z_over_h_to_i[0] = 1;
163 for (
int i=1; i<=max_degree; ++i) {
164 x_over_h_to_i[i] = x_over_h_to_i[i-1]*(0.5*x/h+0.5);
165 y_over_h_to_i[i] = y_over_h_to_i[i-1]*(0.5*y/h+0.5);
166 z_over_h_to_i[i] = z_over_h_to_i[i-1]*(0.5*z/h+0.5);
167 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));
168 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));
169 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));
172 for (
int n1 = starting_order; n1 <= max_degree; n1++){
173 int degree_choose_n1 = factorial[max_degree]/(factorial[n1]*factorial[max_degree-n1]);
174 for (
int n2 = starting_order; n2 <= max_degree; n2++){
175 int degree_choose_n2 = factorial[max_degree]/(factorial[n2]*factorial[max_degree-n2]);
176 for (
int n3 = starting_order; n3 <= max_degree; n3++){
177 int degree_choose_n3 = factorial[max_degree]/(factorial[n3]*factorial[max_degree-n3]);
178 double new_contribution = 0.0;
179 if (partial_direction==0) {
181 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];
183 if (max_degree-n1>0) {
184 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];
186 }
else if (partial_direction==1) {
188 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];
190 if (max_degree-n2>0) {
191 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];
195 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];
197 if (max_degree-n3>0) {
198 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];
201 *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * new_contribution;
206 }
else if (dimension==2) {
211 x_over_h_to_i[0] = 1;
212 y_over_h_to_i[0] = 1;
213 one_minus_x_over_h_to_i[0] = 1;
214 one_minus_y_over_h_to_i[0] = 1;
215 for (
int i=1; i<=max_degree; ++i) {
216 x_over_h_to_i[i] = x_over_h_to_i[i-1]*(0.5*x/h+0.5);
217 y_over_h_to_i[i] = y_over_h_to_i[i-1]*(0.5*y/h+0.5);
218 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));
219 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));
222 for (
int n1 = starting_order; n1 <= max_degree; n1++){
223 int degree_choose_n1 = factorial[max_degree]/(factorial[n1]*factorial[max_degree-n1]);
224 for (
int n2 = starting_order; n2 <= max_degree; n2++){
225 int degree_choose_n2 = factorial[max_degree]/(factorial[n2]*factorial[max_degree-n2]);
226 double new_contribution = 0.0;
227 if (partial_direction==0) {
229 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];
231 if (max_degree-n1>0) {
232 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];
236 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];
238 if (max_degree-n2>0) {
239 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];
242 *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * new_contribution;
249 x_over_h_to_i[0] = 1;
250 one_minus_x_over_h_to_i[0] = 1;
251 for (
int i=1; i<=max_degree; ++i) {
252 x_over_h_to_i[i] = x_over_h_to_i[i-1]*(0.5*x/h+0.5);
253 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));
256 for (
int n=starting_order; n<=max_degree; ++n) {
257 int degree_choose_n = factorial[max_degree]/(factorial[n]*factorial[max_degree-n]);
258 double new_contribution = 0.0;
261 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];
263 if (max_degree-n>0) {
265 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];
267 *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * new_contribution;
290 KOKKOS_INLINE_FUNCTION
291 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) {
292 Kokkos::single(Kokkos::PerThread(teamMember), [&] () {
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 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 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.
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...
KOKKOS_INLINE_FUNCTION int getSize(const int degree, const int dimension)
Returns size of basis.
KOKKOS_INLINE_FUNCTION int pown(int n, unsigned p)
n^p (n^0 returns 1, regardless of n)