18 constexpr int edge_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
19 {0, 3}, {1, 3}, {2, 3}};
21 std::fill(alpha, &alpha[(
D + 1) *
NBEDGE_H1(
N)], 0);
22 for (
int n =
N - 1;
n != 0; --
n) {
26 alpha[edge_nodes[Side][0]] =
n;
27 alpha[edge_nodes[Side][1]] =
N -
n;
36 constexpr int tri_nodes[4][3] = {{0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
39 for (
int n0 = 1; n0 !=
N - 1; ++n0) {
40 for (
int n1 = 1; n1 !=
N - n0; ++n1) {
49 alpha[tri_nodes[Side][0]] = n0;
50 alpha[tri_nodes[Side][1]] = n1;
51 alpha[tri_nodes[Side][2]] =
N - n0 - n1;
113 return generateIndicesEdgeOnSimplex<2, 0>(
N, alpha);
115 return generateIndicesEdgeOnSimplex<2, 1>(
N, alpha);
117 return generateIndicesEdgeOnSimplex<2, 2>(
N, alpha);
120 "Wrong side number, on triangle are edges from 0 to 2");
141 CHKERR generateIndicesEdgeOnSimplex<3, 0>(
N[0], alpha[0]);
142 CHKERR generateIndicesEdgeOnSimplex<3, 1>(
N[1], alpha[1]);
143 CHKERR generateIndicesEdgeOnSimplex<3, 2>(
N[2], alpha[2]);
144 CHKERR generateIndicesEdgeOnSimplex<3, 3>(
N[3], alpha[3]);
145 CHKERR generateIndicesEdgeOnSimplex<3, 4>(
N[4], alpha[4]);
146 CHKERR generateIndicesEdgeOnSimplex<3, 5>(
N[5], alpha[5]);
155 return generateIndicesEdgeOnSimplex<3, 0>(
N, alpha);
157 return generateIndicesEdgeOnSimplex<3, 1>(
N, alpha);
159 return generateIndicesEdgeOnSimplex<3, 2>(
N, alpha);
161 return generateIndicesEdgeOnSimplex<3, 3>(
N, alpha);
163 return generateIndicesEdgeOnSimplex<3, 4>(
N, alpha);
165 return generateIndicesEdgeOnSimplex<3, 5>(
N, alpha);
168 "Wrong side number, on tetrahedron are edges from 0 to 5 "
169 "in generateIndicesEdgeTet");
176 CHKERR generateIndicesTriOnSimplex<3, 0>(
N[0], alpha[0]);
177 CHKERR generateIndicesTriOnSimplex<3, 1>(
N[1], alpha[1]);
178 CHKERR generateIndicesTriOnSimplex<3, 2>(
N[2], alpha[2]);
179 CHKERR generateIndicesTriOnSimplex<3, 3>(
N[3], alpha[3]);
184 const int N,
int *alpha) {
187 return generateIndicesTriOnSimplex<3, 0>(
N, alpha);
189 return generateIndicesTriOnSimplex<3, 1>(
N, alpha);
191 return generateIndicesTriOnSimplex<3, 2>(
N, alpha);
193 return generateIndicesTriOnSimplex<3, 3>(
N, alpha);
196 "Wrong side number, on tetrahedron are from from 0 to 3 "
197 "in generateIndicesTriTet");
207 const int N,
const int n_alpha,
const int *alpha,
const int *diff,
208 const int n_alpha_diff,
const int *alpha_diff,
double *
c) {
211 std::fill(
c, &
c[n_alpha * n_alpha_diff], 0);
212 int abs_diff = diff[0];
213 for (
int d = 1; d !=
D + 1; ++d)
216 const double b = boost::math::factorial<double>(
N) /
217 boost::math::factorial<double>(
N - abs_diff);
219 for (
int i = 0;
i != n_alpha; ++
i) {
220 for (
int j = 0;
j != n_alpha_diff; ++
j) {
222 for (; d !=
D + 1; ++d)
223 if (alpha_diff[(
D + 1) *
j + d] + diff[d] != alpha[(
D + 1) *
i + d])
226 c[
i * n_alpha_diff +
j] = b;
236 const int N,
const int n_alpha,
const int *alpha,
const int *diff,
237 const int n_alpha_diff,
const int *alpha_diff,
double *
c) {
238 return genrateDerivativeIndices<1>(
N, n_alpha, alpha, diff, n_alpha_diff,
243 const int N,
const int n_alpha,
const int *alpha,
const int *diff,
244 const int n_alpha_diff,
const int *alpha_diff,
double *
c) {
245 return genrateDerivativeIndices<2>(
N, n_alpha, alpha, diff, n_alpha_diff,
250 const int N,
const int n_alpha,
const int *alpha,
const int *diff,
251 const int n_alpha_diff,
const int *alpha_diff,
double *
c) {
252 return genrateDerivativeIndices<3>(
N, n_alpha, alpha, diff, n_alpha_diff,
292BernsteinBezier::baseFunctions(
const int N,
const int gdim,
const int n_alpha,
293 const int *alpha,
const double *
lambda,
294 const double *grad_lambda,
double *base,
298 const int *
const alpha0 = alpha;
299 const double fN = boost::math::factorial<double>(
N);
300 constexpr int MAX_ALPHA = 12;
301 int max_alpha = *std::max_element(alpha, &alpha[(
D + 1) * n_alpha]);
302 if(max_alpha > MAX_ALPHA)
304 "Is assumed maximal order not to be bigger than %d", MAX_ALPHA);
305 std::array<
double, (
D + 1) * (MAX_ALPHA + 1)> pow_alpha;
306 std::array<
double, (MAX_ALPHA + 1)> factorial_alpha;
307 std::array<double, D + 1> terms, diff_terms;
308 const double *
const grad_lambda0 = grad_lambda;
310 factorial_alpha[0] = 1;
312 factorial_alpha[1] = 1;
314 for (
int a = 2;
a <= max_alpha; ++
a)
315 factorial_alpha[
a] = factorial_alpha[
a - 1] *
a;
317 for (
int g = 0;
g != gdim; ++
g) {
320 const size_t shift = (MAX_ALPHA + 1) *
n;
321 double *pow_alpha_ptr = &pow_alpha[shift];
324 if (max_alpha >= 1) {
330 for (
int a = 2;
a <= max_alpha; ++
a) {
331 const double p = (*pow_alpha_ptr) * (*lambda);
337 for (
int n0 = 0; n0 != n_alpha; ++n0) {
339 grad_lambda = grad_lambda0;
341 double f = factorial_alpha[(*alpha)];
342 double *terms_ptr = terms.data();
343 double *diff_terms_ptr = diff_terms.data();
344 *terms_ptr = pow_alpha[(*alpha)];
347 *diff_terms_ptr = (*alpha) * pow_alpha[(*alpha) - 1];
356 for (
int n1 = 1; n1 <
D + 1;
357 ++n1, ++alpha, ++terms_ptr, ++diff_terms_ptr) {
358 f *= factorial_alpha[(*alpha)];
359 const size_t shift = (MAX_ALPHA + 1) * n1;
360 *terms_ptr = pow_alpha[shift + (*alpha)];
363 *diff_terms_ptr = (*alpha) * pow_alpha[shift + (*alpha) - 1];
370 const double b = fN / f;
375 double *terms_ptr = terms.data();
376 double *diff_terms_ptr = diff_terms.data();
377 double z = *diff_terms_ptr;
380 for (
int n2 = 1; n2 !=
D + 1; ++n2, ++terms_ptr)
383 for (
int d = 0; d !=
D; ++d, ++grad_lambda)
384 grad_base[d] = z * (*grad_lambda);
386 for (
int n1 = 1; n1 <
D + 1; ++n1) {
390 for (terms_ptr = terms.data(); n2 != n1; ++n2, ++terms_ptr)
396 for (; n2 <
D + 1; ++n2, ++terms_ptr)
399 for (
int d = 0; d !=
D; ++d, ++grad_lambda)
400 grad_base[d] += z * (*grad_lambda);
404 for (
int d = 0; d !=
D; ++d)