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)