9 #ifndef _COMPADRE_GMLS_BERNSTEIN_BASIS_HPP_
10 #define _COMPADRE_GMLS_BERNSTEIN_BASIS_HPP_
23 namespace BernsteinPolynomialBasis {
29 KOKKOS_INLINE_FUNCTION
30 int getSize(
const int degree,
const int dimension) {
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};
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));
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];
90 }
else if (dimension==2) {
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));
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];
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));
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];
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};
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));
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) {
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];
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];
194 }
else if (partial_direction==1) {
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];
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];
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];
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];
209 *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * new_contribution;
214 }
else if (dimension==2) {
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));
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) {
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];
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];
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];
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];
250 *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * new_contribution;
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));
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;
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];
271 if (max_degree-n>0) {
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];
275 *(delta+i) = weight_of_original_value * *(delta+i) + weight_of_new_value * new_contribution;
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), [&] () {
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...