if (i < max_q_indices)
return tensor_polys.template compute_derivative<order>(i, p);
- const unsigned int comp = i - tensor_polys.n();
+ [[maybe_unused]] const unsigned int comp = i - tensor_polys.n();
- Tensor<order, dim> derivative;
- switch (order)
+ if constexpr (order == 1)
{
- case 1:
+ Tensor<1, dim> derivative;
+ for (unsigned int d = 0; d < dim; ++d)
{
- Tensor<1, dim> &derivative_1 =
- *reinterpret_cast<Tensor<1, dim> *>(&derivative);
-
- for (unsigned int d = 0; d < dim; ++d)
- {
- derivative_1[d] = 1.;
- // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p)
- for (unsigned j = 0; j < dim; ++j)
- derivative_1[d] *=
- (d == j ? 4 * (1 - 2 * p[j]) : 4 * p[j] * (1 - p[j]));
- // and multiply with (2*x_i-1)^{r-1}
- for (unsigned int i = 0; i < q_degree - 1; ++i)
- derivative_1[d] *= 2 * p[comp] - 1;
- }
-
- if (q_degree >= 2)
- {
- // add \prod_{i=1}^d 4*(x_i(1-x_i))(p)
- double value = 1.;
- for (unsigned int j = 0; j < dim; ++j)
- value *= 4 * p[j] * (1 - p[j]);
- // and multiply with grad(2*x_i-1)^{r-1}
- double tmp = value * 2 * (q_degree - 1);
- for (unsigned int i = 0; i < q_degree - 2; ++i)
- tmp *= 2 * p[comp] - 1;
- derivative_1[comp] += tmp;
- }
-
- return derivative;
+ derivative[d] = 1.;
+ // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p)
+ for (unsigned j = 0; j < dim; ++j)
+ derivative[d] *=
+ (d == j ? 4 * (1 - 2 * p[j]) : 4 * p[j] * (1 - p[j]));
+ // and multiply with (2*x_i-1)^{r-1}
+ for (unsigned int i = 0; i < q_degree - 1; ++i)
+ derivative[d] *= 2 * p[comp] - 1;
}
- case 2:
+
+ if (q_degree >= 2)
{
- Tensor<2, dim> &derivative_2 =
- *reinterpret_cast<Tensor<2, dim> *>(&derivative);
+ // add \prod_{i=1}^d 4*(x_i(1-x_i))(p)
+ double value = 1.;
+ for (unsigned int j = 0; j < dim; ++j)
+ value *= 4 * p[j] * (1 - p[j]);
+ // and multiply with grad(2*x_i-1)^{r-1}
+ double tmp = value * 2 * (q_degree - 1);
+ for (unsigned int i = 0; i < q_degree - 2; ++i)
+ tmp *= 2 * p[comp] - 1;
+ derivative[comp] += tmp;
+ }
- double v[dim + 1][3];
- {
- for (unsigned int c = 0; c < dim; ++c)
- {
- v[c][0] = 4 * p[c] * (1 - p[c]);
- v[c][1] = 4 * (1 - 2 * p[c]);
- v[c][2] = -8;
- }
+ return derivative;
+ }
+ else if constexpr (order == 2)
+ {
+ Tensor<2, dim> derivative;
- double tmp = 1.;
- for (unsigned int i = 0; i < q_degree - 1; ++i)
- tmp *= 2 * p[comp] - 1;
- v[dim][0] = tmp;
+ double v[dim + 1][3];
+ {
+ for (unsigned int c = 0; c < dim; ++c)
+ {
+ v[c][0] = 4 * p[c] * (1 - p[c]);
+ v[c][1] = 4 * (1 - 2 * p[c]);
+ v[c][2] = -8;
+ }
- if (q_degree >= 2)
- {
- double tmp = 2 * (q_degree - 1);
- for (unsigned int i = 0; i < q_degree - 2; ++i)
- tmp *= 2 * p[comp] - 1;
- v[dim][1] = tmp;
- }
- else
- v[dim][1] = 0.;
+ double tmp = 1.;
+ for (unsigned int i = 0; i < q_degree - 1; ++i)
+ tmp *= 2 * p[comp] - 1;
+ v[dim][0] = tmp;
- if (q_degree >= 3)
- {
- double tmp = 4 * (q_degree - 2) * (q_degree - 1);
- for (unsigned int i = 0; i < q_degree - 3; ++i)
- tmp *= 2 * p[comp] - 1;
- v[dim][2] = tmp;
- }
- else
- v[dim][2] = 0.;
+ if (q_degree >= 2)
+ {
+ double tmp = 2 * (q_degree - 1);
+ for (unsigned int i = 0; i < q_degree - 2; ++i)
+ tmp *= 2 * p[comp] - 1;
+ v[dim][1] = tmp;
}
+ else
+ v[dim][1] = 0.;
- // calculate (\partial_j \partial_k \psi) * monomial
- Tensor<2, dim> grad_grad_1;
- for (unsigned int d1 = 0; d1 < dim; ++d1)
- for (unsigned int d2 = 0; d2 < dim; ++d2)
+ if (q_degree >= 3)
+ {
+ double tmp = 4 * (q_degree - 2) * (q_degree - 1);
+ for (unsigned int i = 0; i < q_degree - 3; ++i)
+ tmp *= 2 * p[comp] - 1;
+ v[dim][2] = tmp;
+ }
+ else
+ v[dim][2] = 0.;
+ }
+
+ // calculate (\partial_j \partial_k \psi) * monomial
+ Tensor<2, dim> grad_grad_1;
+ for (unsigned int d1 = 0; d1 < dim; ++d1)
+ for (unsigned int d2 = 0; d2 < dim; ++d2)
+ {
+ grad_grad_1[d1][d2] = v[dim][0];
+ for (unsigned int x = 0; x < dim; ++x)
{
- grad_grad_1[d1][d2] = v[dim][0];
- for (unsigned int x = 0; x < dim; ++x)
+ unsigned int derivative = 0;
+ if (d1 == x || d2 == x)
{
- unsigned int derivative = 0;
- if (d1 == x || d2 == x)
- {
- if (d1 == d2)
- derivative = 2;
- else
- derivative = 1;
- }
- grad_grad_1[d1][d2] *= v[x][derivative];
+ if (d1 == d2)
+ derivative = 2;
+ else
+ derivative = 1;
}
+ grad_grad_1[d1][d2] *= v[x][derivative];
}
+ }
- // calculate (\partial_j \psi) *(\partial_k monomial)
- // and (\partial_k \psi) *(\partial_j monomial)
- Tensor<2, dim> grad_grad_2;
- Tensor<2, dim> grad_grad_3;
- for (unsigned int d = 0; d < dim; ++d)
+ // calculate (\partial_j \psi) *(\partial_k monomial)
+ // and (\partial_k \psi) *(\partial_j monomial)
+ Tensor<2, dim> grad_grad_2;
+ Tensor<2, dim> grad_grad_3;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ grad_grad_2[d][comp] = v[dim][1];
+ grad_grad_3[comp][d] = v[dim][1];
+ for (unsigned int x = 0; x < dim; ++x)
{
- grad_grad_2[d][comp] = v[dim][1];
- grad_grad_3[comp][d] = v[dim][1];
- for (unsigned int x = 0; x < dim; ++x)
- {
- grad_grad_2[d][comp] *= v[x][d == x];
- grad_grad_3[comp][d] *= v[x][d == x];
- }
+ grad_grad_2[d][comp] *= v[x][d == x];
+ grad_grad_3[comp][d] *= v[x][d == x];
}
+ }
- // calculate \psi *(\partial j \partial_k monomial) and sum
- double psi_value = 1.;
- for (unsigned int x = 0; x < dim; ++x)
- psi_value *= v[x][0];
+ // calculate \psi *(\partial j \partial_k monomial) and sum
+ double psi_value = 1.;
+ for (unsigned int x = 0; x < dim; ++x)
+ psi_value *= v[x][0];
- for (unsigned int d1 = 0; d1 < dim; ++d1)
- for (unsigned int d2 = 0; d2 < dim; ++d2)
- derivative_2[d1][d2] =
- grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2];
- derivative_2[comp][comp] += psi_value * v[dim][2];
+ for (unsigned int d1 = 0; d1 < dim; ++d1)
+ for (unsigned int d2 = 0; d2 < dim; ++d2)
+ derivative[d1][d2] =
+ grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2];
+ derivative[comp][comp] += psi_value * v[dim][2];
- return derivative;
- }
- default:
- {
- DEAL_II_NOT_IMPLEMENTED();
- return derivative;
- }
+ return derivative;
+ }
+ else
+ {
+ DEAL_II_NOT_IMPLEMENTED();
+ return {};
}
}