17 #ifndef dealii_matrix_free_tensor_product_kernels_h
18 #define dealii_matrix_free_tensor_product_kernels_h
116 bool transpose_matrix,
128 static_assert(n_rows > 0 || n_columns > 0,
129 "Specialization only for n_rows, n_columns > 0");
130 Assert(n_rows > 0 && n_columns > 0,
135 "This function should only use EvaluatorQuantity::value");
137 constexpr
int mm = transpose_matrix ? n_rows : n_columns,
138 nn = transpose_matrix ? n_columns : n_rows;
140 std::array<Number, mm> x;
141 for (
int i = 0; i < mm; ++i)
142 x[i] = in[stride_in * i];
143 for (
int col = 0; col < nn; ++col)
146 if (transpose_matrix ==
true)
148 res0 =
matrix[col] * x[0];
149 for (
int i = 1; i < mm; ++i)
150 res0 +=
matrix[i * n_columns + col] * x[i];
154 res0 =
matrix[col * n_columns] * x[0];
155 for (
int i = 1; i < mm; ++i)
156 res0 +=
matrix[col * n_columns + i] * x[i];
159 out[stride_out * col] += res0;
161 out[stride_out * col] = res0;
173 bool transpose_matrix,
184 const int stride_out)
186 const int mm = transpose_matrix ? n_rows : n_columns,
187 nn = transpose_matrix ? n_columns : n_rows;
189 Assert(n_rows > 0 && n_columns > 0,
191 Assert(n_rows > 0 && n_columns > 0,
197 "This function should only use EvaluatorQuantity::value");
204 if (transpose_matrix && n_rows == 2)
206 const Number2 *matrix_1 =
matrix + n_columns;
207 const Number x0 = in[0], x1 = in[stride_in];
208 for (
int col = 0; col < nn; ++col)
210 const Number result =
matrix[col] * x0 + matrix_1[col] * x1;
212 out[stride_out * col] += result;
214 out[stride_out * col] = result;
217 else if (transpose_matrix && n_rows == 3)
219 const Number2 *matrix_1 =
matrix + n_columns;
220 const Number2 *matrix_2 = matrix_1 + n_columns;
221 const Number x0 = in[0], x1 = in[stride_in], x2 = in[2 * stride_in];
222 for (
int col = 0; col < nn; ++col)
224 const Number result =
225 matrix[col] * x0 + matrix_1[col] * x1 + matrix_2[col] * x2;
227 out[stride_out * col] += result;
229 out[stride_out * col] = result;
234 std::array<Number, 129> x;
235 for (
int i = 0; i < mm; ++i)
236 x[i] = in[stride_in * i];
239 for (
int col = 0; col < nn; ++col)
241 if (transpose_matrix ==
true)
243 res0 =
matrix[col] * x[0];
244 for (
int i = 1; i < mm; ++i)
245 res0 +=
matrix[i * n_columns + col] * x[i];
249 res0 =
matrix[col * n_columns] * x[0];
250 for (
int i = 1; i < mm; ++i)
251 res0 +=
matrix[col * n_columns + i] * x[i];
254 out[stride_out * col] += res0;
256 out[stride_out * col] = res0;
275 bool transpose_matrix,
287 static_assert(n_rows > 0 || n_columns > 0,
288 "Specialization only for n_rows, n_columns > 0");
289 Assert(n_rows > 0 && n_columns > 0,
294 constexpr
int mm = transpose_matrix ? n_rows : n_columns,
295 nn = transpose_matrix ? n_columns : n_rows;
296 constexpr
int n_cols = nn / 2;
297 constexpr
int mid = mm / 2;
299 std::array<Number, mm> x;
300 for (
int i = 0; i < mm; ++i)
301 x[i] = in[stride_in * i];
324 for (
int col = 0; col < n_cols; ++col)
328 if (transpose_matrix ==
true)
331 val1 =
matrix[nn - 1 - col];
335 val0 =
matrix[col * n_columns];
336 val1 =
matrix[(col + 1) * n_columns - 1];
342 res0 += val1 * x[mm - 1];
343 res1 += val0 * x[mm - 1];
344 for (
int ind = 1; ind < mid; ++ind)
346 if (transpose_matrix ==
true)
348 val0 =
matrix[ind * n_columns + col];
349 val1 =
matrix[ind * n_columns + nn - 1 - col];
353 val0 =
matrix[col * n_columns + ind];
354 val1 =
matrix[(col + 1) * n_columns - 1 - ind];
356 res0 += val0 * x[ind];
357 res1 += val1 * x[ind];
358 res0 += val1 * x[mm - 1 - ind];
359 res1 += val0 * x[mm - 1 - ind];
363 res0 = res1 = Number();
364 if (transpose_matrix ==
true)
368 const Number tmp =
matrix[mid * n_columns + col] * x[mid];
375 if (mm % 2 == 1 && nn % 2 == 0)
377 const Number tmp =
matrix[col * n_columns + mid] * x[mid];
384 out[stride_out * col] += res0;
385 out[stride_out * (nn - 1 - col)] += res1;
389 out[stride_out * col] = res0;
390 out[stride_out * (nn - 1 - col)] = res1;
393 if (transpose_matrix ==
true && nn % 2 == 1 && mm % 2 == 1)
396 out[stride_out * n_cols] += x[mid];
398 out[stride_out * n_cols] = x[mid];
400 else if (transpose_matrix ==
true && nn % 2 == 1)
405 res0 =
matrix[n_cols] * (x[0] + x[mm - 1]);
406 for (
int ind = 1; ind < mid; ++ind)
408 const Number2 val0 =
matrix[ind * n_columns + n_cols];
409 res0 += val0 * (x[ind] + in[mm - 1 - ind]);
415 out[stride_out * n_cols] += res0;
417 out[stride_out * n_cols] = res0;
419 else if (transpose_matrix ==
false && nn % 2 == 1)
424 res0 =
matrix[n_cols * n_columns] * (x[0] + x[mm - 1]);
425 for (
int ind = 1; ind < mid; ++ind)
427 const Number2 val0 =
matrix[n_cols * n_columns + ind];
428 res0 += val0 * (x[ind] + x[mm - 1 - ind]);
437 out[stride_out * n_cols] += res0;
439 out[stride_out * n_cols] = res0;
460 for (
int col = 0; col < n_cols; ++col)
464 if (transpose_matrix ==
true)
467 val1 =
matrix[nn - 1 - col];
471 val0 =
matrix[col * n_columns];
472 val1 =
matrix[(nn - col - 1) * n_columns];
478 res0 -= val1 * x[mm - 1];
479 res1 -= val0 * x[mm - 1];
480 for (
int ind = 1; ind < mid; ++ind)
482 if (transpose_matrix ==
true)
484 val0 =
matrix[ind * n_columns + col];
485 val1 =
matrix[ind * n_columns + nn - 1 - col];
489 val0 =
matrix[col * n_columns + ind];
490 val1 =
matrix[(nn - col - 1) * n_columns + ind];
492 res0 += val0 * x[ind];
493 res1 += val1 * x[ind];
494 res0 -= val1 * x[mm - 1 - ind];
495 res1 -= val0 * x[mm - 1 - ind];
499 res0 = res1 = Number();
502 if (transpose_matrix ==
true)
503 val0 =
matrix[mid * n_columns + col];
505 val0 =
matrix[col * n_columns + mid];
506 const Number tmp = val0 * x[mid];
512 out[stride_out * col] += res0;
513 out[stride_out * (nn - 1 - col)] += res1;
517 out[stride_out * col] = res0;
518 out[stride_out * (nn - 1 - col)] = res1;
525 if (transpose_matrix ==
true)
528 val0 =
matrix[n_cols * n_columns];
529 res0 = val0 * (x[0] - x[mm - 1]);
530 for (
int ind = 1; ind < mid; ++ind)
532 if (transpose_matrix ==
true)
533 val0 =
matrix[ind * n_columns + n_cols];
535 val0 =
matrix[n_cols * n_columns + ind];
536 Number in1 = val0 * (x[ind] - x[mm - 1 - ind]);
540 out[stride_out * n_cols] += res0;
542 out[stride_out * n_cols] = res0;
549 for (
int col = 0; col < n_cols; ++col)
553 if (transpose_matrix ==
true)
556 val1 =
matrix[nn - 1 - col];
560 val0 =
matrix[col * n_columns];
561 val1 =
matrix[(col + 1) * n_columns - 1];
567 res0 += val1 * x[mm - 1];
568 res1 += val0 * x[mm - 1];
569 for (
int ind = 1; ind < mid; ++ind)
571 if (transpose_matrix ==
true)
573 val0 =
matrix[ind * n_columns + col];
574 val1 =
matrix[ind * n_columns + nn - 1 - col];
578 val0 =
matrix[col * n_columns + ind];
579 val1 =
matrix[(col + 1) * n_columns - 1 - ind];
581 res0 += val0 * x[ind];
582 res1 += val1 * x[ind];
583 res0 += val1 * x[mm - 1 - ind];
584 res1 += val0 * x[mm - 1 - ind];
588 res0 = res1 = Number();
591 if (transpose_matrix ==
true)
592 val0 =
matrix[mid * n_columns + col];
594 val0 =
matrix[col * n_columns + mid];
595 const Number tmp = val0 * x[mid];
601 out[stride_out * col] += res0;
602 out[stride_out * (nn - 1 - col)] += res1;
606 out[stride_out * col] = res0;
607 out[stride_out * (nn - 1 - col)] = res1;
614 if (transpose_matrix ==
true)
617 val0 =
matrix[n_cols * n_columns];
620 res0 = val0 * (x[0] + x[mm - 1]);
621 for (
int ind = 1; ind < mid; ++ind)
623 if (transpose_matrix ==
true)
624 val0 =
matrix[ind * n_columns + n_cols];
626 val0 =
matrix[n_cols * n_columns + ind];
627 Number in1 = val0 * (x[ind] + x[mm - 1 - ind]);
635 if (transpose_matrix ==
true)
636 val0 =
matrix[mid * n_columns + n_cols];
638 val0 =
matrix[n_cols * n_columns + mid];
639 res0 += val0 * x[mid];
642 out[stride_out * n_cols] += res0;
644 out[stride_out * n_cols] = res0;
672 int n_columns_static,
673 int stride_in_static,
674 int stride_out_static,
675 bool transpose_matrix,
686 int n_rows_runtime = 0,
687 int n_columns_runtime = 0,
688 int stride_in_runtime = 0,
689 int stride_out_runtime = 0)
691 static_assert(n_rows_static >= 0 && n_columns_static >= 0,
692 "Negative loop ranges are not allowed!");
694 const int n_rows = n_rows_static == 0 ? n_rows_runtime : n_rows_static;
695 const int n_columns =
696 n_rows_static == 0 ? n_columns_runtime : n_columns_static;
697 const int stride_in =
698 n_rows_static == 0 ? stride_in_runtime : stride_in_static;
699 const int stride_out =
700 n_rows_static == 0 ? stride_out_runtime : stride_out_static;
702 Assert(n_rows > 0 && n_columns > 0,
707 const int mm = transpose_matrix ? n_rows : n_columns,
708 nn = transpose_matrix ? n_columns : n_rows;
709 const int n_half = nn / 2;
710 const int m_half = mm / 2;
712 constexpr
int array_length =
713 (n_rows_static == 0) ?
716 (1 + (transpose_matrix ? n_rows_static : n_columns_static) / 2);
717 const int offset = (n_columns + 1) / 2;
721 std::array<Number, array_length> xp, xm;
722 for (
int i = 0; i < m_half; ++i)
726 xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
727 xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
731 xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
732 xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
735 Number xmid = in[stride_in * m_half];
736 for (
int col = 0; col < n_half; ++col)
741 if (transpose_matrix ==
true)
744 r1 =
matrix[(n_rows - 1) * offset + col] * xm[0];
748 r0 =
matrix[col * offset] * xp[0];
749 r1 =
matrix[(n_rows - 1 - col) * offset] * xm[0];
751 for (
int ind = 1; ind < m_half; ++ind)
753 if (transpose_matrix ==
true)
755 r0 +=
matrix[ind * offset + col] * xp[ind];
756 r1 +=
matrix[(n_rows - 1 - ind) * offset + col] * xm[ind];
760 r0 +=
matrix[col * offset + ind] * xp[ind];
761 r1 +=
matrix[(n_rows - 1 - col) * offset + ind] * xm[ind];
767 if (mm % 2 == 1 && transpose_matrix ==
true)
770 r1 +=
matrix[m_half * offset + col] * xmid;
772 r0 +=
matrix[m_half * offset + col] * xmid;
774 else if (mm % 2 == 1 &&
777 r0 +=
matrix[col * offset + m_half] * xmid;
781 out[stride_out * col] += r0 + r1;
783 transpose_matrix ==
false)
784 out[stride_out * (nn - 1 - col)] += r1 - r0;
786 out[stride_out * (nn - 1 - col)] += r0 - r1;
790 out[stride_out * col] = r0 + r1;
792 transpose_matrix ==
false)
793 out[stride_out * (nn - 1 - col)] = r1 - r0;
795 out[stride_out * (nn - 1 - col)] = r0 - r1;
799 nn % 2 == 1 && mm % 2 == 1 && mm > 3)
802 out[stride_out * n_half] +=
matrix[m_half * offset + n_half] * xmid;
804 out[stride_out * n_half] =
matrix[m_half * offset + n_half] * xmid;
806 else if (transpose_matrix ==
true && nn % 2 == 1)
811 r0 =
matrix[n_half] * xp[0];
812 for (
int ind = 1; ind < m_half; ++ind)
813 r0 +=
matrix[ind * offset + n_half] * xp[ind];
818 r0 +=
matrix[m_half * offset + n_half] * xmid;
821 out[stride_out * n_half] += r0;
823 out[stride_out * n_half] = r0;
825 else if (transpose_matrix ==
false && nn % 2 == 1)
832 r0 =
matrix[n_half * offset] * xm[0];
833 for (
int ind = 1; ind < m_half; ++ind)
834 r0 +=
matrix[n_half * offset + ind] * xm[ind];
838 r0 =
matrix[n_half * offset] * xp[0];
839 for (
int ind = 1; ind < m_half; ++ind)
840 r0 +=
matrix[n_half * offset + ind] * xp[ind];
847 r0 +=
matrix[n_half * offset + m_half] * xmid;
850 out[stride_out * n_half] += r0;
852 out[stride_out * n_half] = r0;
864 bool transpose_matrix,
885 matrix, in, out, n_rows, n_columns, stride_in, stride_out);
911 bool transpose_matrix,
920 static_assert(n_rows > 0 && n_columns > 0,
921 "Specialization requires n_rows, n_columns > 0");
923 constexpr
bool evaluate_antisymmetric =
926 constexpr
int mm = transpose_matrix ? n_rows : n_columns,
927 nn = transpose_matrix ? n_columns : n_rows;
928 constexpr
int n_half = nn / 2;
929 constexpr
int m_half = mm / 2;
931 if (transpose_matrix)
933 std::array<Number, mm> x;
934 for (
unsigned int i = 0; i < mm; ++i)
935 x[i] = in[stride_in * i];
936 for (
unsigned int col = 0; col < n_half; ++col)
942 r1 =
matrix[col + n_columns] * x[1];
943 for (
unsigned int ind = 1; ind < m_half; ++ind)
945 r0 +=
matrix[col + 2 * ind * n_columns] * x[2 * ind];
947 matrix[col + (2 * ind + 1) * n_columns] * x[2 * ind + 1];
953 r0 +=
matrix[col + (mm - 1) * n_columns] * x[mm - 1];
956 out[stride_out * col] += r0 + r1;
957 if (evaluate_antisymmetric)
958 out[stride_out * (nn - 1 - col)] += r1 - r0;
960 out[stride_out * (nn - 1 - col)] += r0 - r1;
964 out[stride_out * col] = r0 + r1;
965 if (evaluate_antisymmetric)
966 out[stride_out * (nn - 1 - col)] = r1 - r0;
968 out[stride_out * (nn - 1 - col)] = r0 - r1;
974 const unsigned int shift = evaluate_antisymmetric ? 1 : 0;
978 for (
unsigned int ind = 1; ind < m_half; ++ind)
979 r0 +=
matrix[n_half + (2 * ind +
shift) * n_columns] *
984 if (!evaluate_antisymmetric && mm % 2 == 1)
985 r0 +=
matrix[n_half + (mm - 1) * n_columns] * x[mm - 1];
987 out[stride_out * n_half] += r0;
989 out[stride_out * n_half] = r0;
994 std::array<Number, m_half + 1> xp, xm;
995 for (
int i = 0; i < m_half; ++i)
996 if (!evaluate_antisymmetric)
998 xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
999 xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1003 xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1004 xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1007 xp[m_half] = in[stride_in * m_half];
1008 for (
unsigned int col = 0; col < n_half; ++col)
1013 r0 =
matrix[2 * col * n_columns] * xp[0];
1014 r1 =
matrix[(2 * col + 1) * n_columns] * xm[0];
1015 for (
unsigned int ind = 1; ind < m_half; ++ind)
1017 r0 +=
matrix[2 * col * n_columns + ind] * xp[ind];
1018 r1 +=
matrix[(2 * col + 1) * n_columns + ind] * xm[ind];
1025 if (evaluate_antisymmetric)
1026 r1 +=
matrix[(2 * col + 1) * n_columns + m_half] * xp[m_half];
1028 r0 +=
matrix[2 * col * n_columns + m_half] * xp[m_half];
1032 out[stride_out * (2 * col)] += r0;
1033 out[stride_out * (2 * col + 1)] += r1;
1037 out[stride_out * (2 * col)] = r0;
1038 out[stride_out * (2 * col + 1)] = r1;
1046 r0 =
matrix[(nn - 1) * n_columns] * xp[0];
1047 for (
unsigned int ind = 1; ind < m_half; ++ind)
1048 r0 +=
matrix[(nn - 1) * n_columns + ind] * xp[ind];
1052 if (mm % 2 == 1 && !evaluate_antisymmetric)
1053 r0 +=
matrix[(nn - 1) * n_columns + m_half] * xp[m_half];
1055 out[stride_out * (nn - 1)] += r0;
1057 out[stride_out * (nn - 1)] = r0;
1091 typename Number2 = Number>
1115 const unsigned int = 0,
1116 const unsigned int = 0)
1125 n_rows * ((n_columns + 1) / 2));
1128 n_rows * ((n_columns + 1) / 2));
1131 n_rows * ((n_columns + 1) / 2));
1141 n_rows * n_columns));
1145 n_rows * n_columns));
1155 const unsigned int dummy1 = 0,
1156 const unsigned int dummy2 = 0)
1190 template <
int direction,
bool contract_over_rows,
bool add,
int str
ide = 1>
1192 values(
const Number in[], Number out[])
const
1195 apply<direction, contract_over_rows, add, false, value_type, stride>(
1204 template <
int direction,
bool contract_over_rows,
bool add,
int str
ide = 1>
1211 apply<direction, contract_over_rows, add, false, gradient_type, stride>(
1220 template <
int direction,
bool contract_over_rows,
bool add>
1229 apply<direction, contract_over_rows, add, false, hessian_type>(
1240 template <
int direction,
bool contract_over_rows,
bool add>
1245 apply<direction, contract_over_rows, add, true, EvaluatorQuantity::value>(
1256 template <
int direction,
bool contract_over_rows,
bool add>
1264 apply<direction, contract_over_rows, add, true, gradient_type>(
1275 template <
int direction,
bool contract_over_rows,
bool add>
1285 apply<direction, contract_over_rows, add, true, hessian_type>(
1325 template <
int direction,
1326 bool contract_over_rows,
1328 bool one_line =
false,
1350 template <
int direction,
1351 bool contract_over_rows,
1362 static_assert(one_line ==
false || direction == dim - 1,
1363 "Single-line evaluation only works for direction=dim-1.");
1364 Assert(shape_data !=
nullptr,
1366 "The given array shape_data must not be the null pointer!"));
1367 Assert(dim == direction + 1 || one_line ==
true || n_rows == n_columns ||
1369 ExcMessage(
"In-place operation only supported for "
1370 "n_rows==n_columns or single-line interpolation"));
1372 constexpr
int mm = contract_over_rows ? n_rows : n_columns,
1373 nn = contract_over_rows ? n_columns : n_rows;
1375 constexpr
int stride_operation =
Utilities::pow(n_columns, direction);
1376 constexpr
int n_blocks1 = one_line ? 1 : stride_operation;
1377 constexpr
int n_blocks2 =
1378 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1380 constexpr
int stride_in = !contract_over_rows ? stride : 1;
1381 constexpr
int stride_out = contract_over_rows ? stride : 1;
1382 for (
int i2 = 0; i2 < n_blocks2; ++i2)
1384 for (
int i1 = 0; i1 < n_blocks1; ++i1)
1390 stride_operation * stride_in,
1391 stride_operation * stride_out,
1393 add>(shape_data, in, out);
1395 if (one_line ==
false)
1401 if (one_line ==
false)
1403 in += stride_operation * (mm - 1) * stride_in;
1404 out += stride_operation * (nn - 1) * stride_out;
1453 const unsigned int n_rows = 0,
1454 const unsigned int n_columns = 0)
1459 , n_columns(n_columns)
1465 n_rows * ((n_columns + 1) / 2));
1468 n_rows * ((n_columns + 1) / 2));
1471 n_rows * ((n_columns + 1) / 2));
1481 n_rows * n_columns));
1485 n_rows * n_columns));
1495 const unsigned int n_rows = 0,
1496 const unsigned int n_columns = 0)
1501 , n_columns(n_columns)
1504 template <
int direction,
bool contract_over_rows,
bool add,
int str
ide = 1>
1509 apply<direction, contract_over_rows, add, false, value_type, stride>(
1513 template <
int direction,
bool contract_over_rows,
bool add,
int str
ide = 1>
1520 apply<direction, contract_over_rows, add, false, gradient_type, stride>(
1524 template <
int direction,
bool contract_over_rows,
bool add>
1531 apply<direction, contract_over_rows, add, false, hessian_type>(
1535 template <
int direction,
bool contract_over_rows,
bool add>
1540 apply<direction, contract_over_rows, add, true, EvaluatorQuantity::value>(
1544 template <
int direction,
bool contract_over_rows,
bool add>
1552 apply<direction, contract_over_rows, add, true, gradient_type>(
1556 template <
int direction,
bool contract_over_rows,
bool add>
1564 apply<direction, contract_over_rows, add, true, hessian_type>(
1568 template <
int direction,
1569 bool contract_over_rows,
1571 bool one_line =
false,
1592 template <
int direction,
1593 bool contract_over_rows,
1604 static_assert(one_line ==
false || direction == dim - 1,
1605 "Single-line evaluation only works for direction=dim-1.");
1606 Assert(shape_data !=
nullptr,
1608 "The given array shape_data must not be the null pointer!"));
1609 Assert(dim == direction + 1 || one_line ==
true || n_rows == n_columns ||
1611 ExcMessage(
"In-place operation only supported for "
1612 "n_rows==n_columns or single-line interpolation"));
1614 const int mm = contract_over_rows ? n_rows : n_columns,
1615 nn = contract_over_rows ? n_columns : n_rows;
1617 const int stride_operation =
1618 direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
1619 const int n_blocks1 = one_line ? 1 : stride_operation;
1620 const int n_blocks2 = direction >= dim - 1 ?
1625 constexpr
int stride_in = !contract_over_rows ? stride : 1;
1626 constexpr
int stride_out = contract_over_rows ? stride : 1;
1627 for (
int i2 = 0; i2 < n_blocks2; ++i2)
1629 for (
int i1 = 0; i1 < n_blocks1; ++i1)
1643 stride_operation * stride_in,
1644 stride_operation * stride_out);
1646 if (one_line ==
false)
1652 if (one_line ==
false)
1654 in += stride_operation * (mm - 1) * stride_in;
1655 out += stride_operation * (nn - 1) * stride_out;
1665 bool contract_over_rows,
1666 bool symmetric_evaluate =
true>
1669 template <
int direction,
1671 typename Number =
double,
1672 typename Number2 =
double>
1677 const bool add_into_result =
false,
1678 const int subface_index_1d = 0)
1683 constexpr
int n_rows = fe_degree + 1;
1684 constexpr
int n_columns = n_q_points_1d;
1685 constexpr
int mm = contract_over_rows ? n_rows : n_columns;
1686 constexpr
int nn = contract_over_rows ? n_columns : n_rows;
1687 const Number2 *shape_data =
1688 symmetric_evaluate ?
1692 Assert(contract_over_rows ==
false || !add_into_result,
1693 ExcMessage(
"Cannot add into result if contract_over_rows = true"));
1696 constexpr
int n_blocks2 =
Utilities::pow(fe_degree, dim - direction - 1);
1697 constexpr
int stride_in = contract_over_rows ? 1 : stride;
1698 constexpr
int stride_out = contract_over_rows ? stride : 1;
1702 for (
int i2 = 0; i2 < n_blocks2; ++i2)
1704 for (
int i1 = 0; i1 < n_blocks1; ++i1)
1706 if (contract_over_rows ==
false && add_into_result)
1711 n_blocks1 * stride_in,
1712 n_blocks1 * stride_out,
1714 true>(shape_data, in, out);
1720 n_blocks1 * stride_in,
1721 n_blocks1 * stride_out,
1723 false>(shape_data, in, out);
1728 in += n_blocks1 * (mm - 1) * stride_in;
1729 out += n_blocks1 * (nn - 1) * stride_out;
1733 template <
int direction,
1734 int normal_direction,
1736 typename Number =
double,
1737 typename Number2 =
double>
1742 const int subface_index_1d = 0)
1747 static_assert(direction != normal_direction,
1748 "Cannot interpolate tangentially in normal direction");
1750 constexpr
int n_rows =
std::max(fe_degree, 0);
1751 constexpr
int n_columns = n_q_points_1d;
1752 const Number2 *shape_data =
1753 symmetric_evaluate ?
1758 constexpr
int n_blocks1 =
1759 (direction > normal_direction) ?
1764 constexpr
int n_blocks2 =
1765 (direction > normal_direction) ?
1767 ((direction + 1 < dim) ?
1768 (
Utilities::pow(fe_degree, dim - 2 - direction) * n_q_points_1d) :
1777 if (contract_over_rows)
1779 in += (n_blocks2 - 1) * n_blocks1 * n_rows + n_blocks1 - 1;
1781 stride * ((n_blocks2 - 1) * n_blocks1 * n_columns + n_blocks1 - 1);
1782 for (
int i2 = 0; i2 < n_blocks2; ++i2)
1784 for (
int i1 = 0; i1 < n_blocks1; ++i1)
1793 false>(shape_data, in, out);
1798 in -= n_blocks1 * (n_rows - 1);
1799 out -= n_blocks1 * (n_columns - 1) * stride;
1804 for (
int i2 = 0; i2 < n_blocks2; ++i2)
1806 for (
int i1 = 0; i1 < n_blocks1; ++i1)
1815 false>(shape_data, in, out);
1820 in += n_blocks1 * (n_columns - 1) * stride;
1821 out += n_blocks1 * (n_rows - 1);
1875 template <
int n_rows_template,
1876 int stride_template,
1877 bool contract_onto_face,
1882 inline std::enable_if_t<contract_onto_face, void>
1884 const std::array<int, 2> &
n_blocks,
1885 const std::array<int, 2> &steps,
1886 const Number *input,
1888 const int n_rows_runtime = 0,
1889 const int stride_runtime = 1)
1891 const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
1892 const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
1896 for (
int i2 = 0; i2 <
n_blocks[1]; ++i2)
1898 for (
int i1 = 0; i1 <
n_blocks[0]; ++i1)
1900 Number res0 = shape_values[0] * input[0];
1902 if (max_derivative > 0)
1903 res1 = shape_values[n_rows] * input[0];
1904 if (max_derivative > 1)
1905 res2 = shape_values[2 * n_rows] * input[0];
1906 for (
int ind = 1; ind < n_rows; ++ind)
1908 res0 += shape_values[ind] * input[stride * ind];
1909 if (max_derivative > 0)
1910 res1 += shape_values[ind + n_rows] * input[stride * ind];
1911 if (max_derivative > 1)
1912 res2 += shape_values[ind + 2 * n_rows] * input[stride * ind];
1917 if (max_derivative > 0)
1918 output1[i1] += res1;
1919 if (max_derivative > 1)
1920 output2[i2] += res2;
1925 if (max_derivative > 0)
1927 if (max_derivative > 1)
1933 if (max_derivative > 0)
1935 if (max_derivative > 1)
1952 const unsigned int n_q_points_1d)
1954 return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
1955 (n_q_points_1d <= 3 * fe_degree / 2 + 1);
1965 template <
int n_rows_template,
1966 int stride_template,
1967 bool contract_onto_face,
1972 inline std::enable_if_t<!contract_onto_face, void>
1974 const std::array<int, 2> &
n_blocks,
1975 const std::array<int, 2> &steps,
1976 const Number *input,
1978 const int n_rows_runtime = 0,
1979 const int stride_runtime = 1)
1981 const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
1982 const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
1986 for (
int i2 = 0; i2 <
n_blocks[1]; ++i2)
1988 for (
int i1 = 0; i1 <
n_blocks[0]; ++i1)
1990 const Number in = input[i1];
1992 if (max_derivative > 0)
1994 if (max_derivative > 1)
1996 for (
int col = 0; col < n_rows; ++col)
1999 add ? (output[col * stride] + shape_values[col] * in) :
2000 (shape_values[col] * in);
2001 if (max_derivative > 0)
2002 result += shape_values[col + n_rows] * in1;
2003 if (max_derivative > 1)
2004 result += shape_values[col + 2 * n_rows] * in2;
2006 output[col * stride] = result;
2011 if (max_derivative > 0)
2013 if (max_derivative > 1)
2027 template <
typename Number,
typename Number2>
2033 template <
int dim,
typename Number,
typename Number2>
2045 template <
int dim,
typename Number>
2051 const unsigned int derivative = 1)
2053 const int n_shapes = poly.size();
2056 std::array<Number, dim>
point;
2057 for (
unsigned int d = 0;
d < dim; ++
d)
2059 for (
int i = 0; i < n_shapes; ++i)
2060 poly[i].values_of_array(
point, derivative, shapes[i].data());
2068 template <
typename Number>
2088 bool do_renumber =
true>
2093 std::array<typename ProductTypeNoPoint<Number, Number2>::type,
2096 const std::vector<unsigned int> &renumber,
2097 const ::ndarray<Number2, 2, dim> *shapes,
2098 const int n_shapes_runtime,
2101 static_assert(0 <= dim && dim <= 3,
"Only dim=0,1,2,3 implemented");
2102 static_assert(1 <= n_values && n_values <= 2,
2103 "Only n_values=1,2 implemented");
2105 const int n_shapes = length > 0 ? length : n_shapes_runtime;
2110 const Number *values_2 =
2113 Utilities::fixed_power<dim>(n_shapes_runtime)) :
2116 std::array<Number3, 2 + n_values> result = {};
2117 for (
int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
2120 std::array<Number3, 1 + n_values> inner_result = {};
2124 if (do_renumber && !renumber.empty())
2125 for (
int i0 = 0; i0 < n_shapes; ++i0, ++i)
2128 inner_result[0] += shapes[i0][1][0] *
values[renumber[i]];
2130 inner_result[1] += shapes[i0][0][0] *
values[renumber[i]];
2132 inner_result[2] += shapes[i0][0][0] * values_2[renumber[i]];
2135 for (
int i0 = 0; i0 < n_shapes; ++i0, ++i)
2138 inner_result[0] += shapes[i0][1][0] *
values[i];
2140 inner_result[1] += shapes[i0][0][0] *
values[i];
2142 inner_result[2] += shapes[i0][0][0] * values_2[i];
2149 result[0] += inner_result[0] * shapes[i1][0][1];
2150 result[1] += inner_result[1] * shapes[i1][1][1];
2152 result[2] += inner_result[1] * shapes[i1][0][1];
2154 result[3] += inner_result[2] * shapes[i1][0][1];
2159 result[0] = inner_result[0];
2161 result[1] = inner_result[1];
2163 result[2] = inner_result[2];
2179 bool do_renumber =
true>
2180 inline std::array<typename ProductTypeNoPoint<Number, Number2>::type,
2183 const ::ndarray<Number2, 2, dim> *shapes,
2186 const std::vector<unsigned int> &renumber = {})
2188 static_assert(0 <= dim && dim <= 3,
"Only dim=0,1,2,3 implemented");
2189 static_assert(1 <= n_values && n_values <= 2,
2190 "Only n_values=1,2 implemented");
2194 std::array<Number3, dim + n_values> result = {};
2208 for (
int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
2210 std::array<Number3, 2 + n_values> inner_result;
2215 do_interpolate_xy<dim, 2, Number2, Number, n_values, do_renumber>(
2216 values, renumber, shapes, n_shapes, i);
2217 else if (n_shapes == 3)
2219 do_interpolate_xy<dim, 3, Number2, Number, n_values, do_renumber>(
2220 values, renumber, shapes, n_shapes, i);
2221 else if (n_shapes == 4)
2223 do_interpolate_xy<dim, 4, Number2, Number, n_values, do_renumber>(
2224 values, renumber, shapes, n_shapes, i);
2225 else if (n_shapes == 5)
2227 do_interpolate_xy<dim, 5, Number2, Number, n_values, do_renumber>(
2228 values, renumber, shapes, n_shapes, i);
2229 else if (n_shapes == 6)
2231 do_interpolate_xy<dim, 6, Number2, Number, n_values, do_renumber>(
2232 values, renumber, shapes, n_shapes, i);
2236 values, renumber, shapes, n_shapes, i);
2241 result[0] += inner_result[0] * shapes[i2][0][2];
2242 result[1] += inner_result[1] * shapes[i2][0][2];
2243 result[2] += inner_result[2] * shapes[i2][1][2];
2245 result[3] += inner_result[2] * shapes[i2][0][2];
2247 result[4] += inner_result[3] * shapes[i2][0][2];
2252 result[0] = inner_result[0];
2253 result[1] = inner_result[1];
2255 result[2] = inner_result[2];
2257 result[3] = inner_result[3];
2262 result[0] = inner_result[0];
2264 result[1] = inner_result[1];
2266 result[2] = inner_result[2];
2284 inline std::array<typename ProductTypeNoPoint<Number, Number2>::type,
2290 static_assert(0 <= dim && dim <= 3,
"Only dim=0,1,2,3 implemented");
2291 static_assert(1 <= n_values && n_values <= 2,
2292 "Only n_values=1,2 implemented");
2297 n_values == 1 || stride == 1,
2298 "Either n_values or stride has to be one for correct data access!");
2303 std::array<Number3, dim + n_values> result;
2316 result[1] = Number3(
values[0]) + p[0] * result[0];
2322 const Number3 val10 = Number3(
values[stride] -
values[0]);
2323 const Number3 val32 = Number3(
values[3 * stride] -
values[2 * stride]);
2324 const Number3 tmp0 = Number3(
values[0]) + p[0] * val10;
2325 const Number3 tmp1 = Number3(
values[2 * stride]) + p[0] * val32;
2328 result[0] = val10 + p[1] * (val32 - val10);
2329 result[1] = tmp1 - tmp0;
2332 result[2] = tmp0 + p[1] * result[1];
2336 const Number3 tmp0_2 =
2338 const Number3 tmp1_2 =
2340 result[3] = tmp0_2 + p[1] * (tmp1_2 - tmp0_2);
2345 const Number3 val10 = Number3(
values[stride] -
values[0]);
2346 const Number3 val32 = Number3(
values[3 * stride] -
values[2 * stride]);
2347 const Number3 tmp0 = Number3(
values[0]) + p[0] * val10;
2348 const Number3 tmp1 = Number3(
values[2 * stride]) + p[0] * val32;
2349 const Number3 tmp10 = tmp1 - tmp0;
2350 const Number3 tmpy0 = tmp0 + p[1] * tmp10;
2352 const Number3 val54 = Number3(
values[5 * stride] -
values[4 * stride]);
2353 const Number3 val76 = Number3(
values[7 * stride] -
values[6 * stride]);
2354 const Number3 tmp2 = Number3(
values[4 * stride]) + p[0] * val54;
2355 const Number3 tmp3 = Number3(
values[6 * stride]) + p[0] * val76;
2356 const Number3 tmp32 = tmp3 - tmp2;
2357 const Number3 tmpy1 = tmp2 + p[1] * tmp32;
2360 result[2] = tmpy1 - tmpy0;
2361 result[1] = tmp10 + p[2] * (tmp32 - tmp10);
2362 const Number3 tmpz0 = val10 + p[1] * (val32 - val10);
2363 result[0] = tmpz0 + p[2] * (val54 + p[1] * (val76 - val54) - tmpz0);
2366 result[3] = tmpy0 + p[2] * result[2];
2409 template <
int dim,
typename Number,
typename Number2>
2415 const std::vector<Number> &
values,
2417 const bool d_linear =
false,
2418 const std::vector<unsigned int> &renumber = {})
2422 std::array<Number3, dim + 1> result;
2431 std::array<::ndarray<Number2, 2, dim>, 200> shapes;
2436 shapes.data(), poly.size(),
values.data(), renumber);
2438 return std::make_pair(result[dim],
2449 bool do_renumber =
true>
2456 const std::vector<unsigned int> &renumber,
2457 const ::ndarray<Number2, 2, dim> *shapes,
2458 const int n_shapes_runtime,
2461 const int n_shapes = length > 0 ? length : n_shapes_runtime;
2463 Number3 result = {};
2464 for (
int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
2471 if (do_renumber && !renumber.empty())
2472 for (
int i0 = 0; i0 < n_shapes; ++i0, ++i)
2475 for (
int i0 = 0; i0 < n_shapes; ++i0, ++i)
2479 result +=
value * shapes[i1][0][1];
2488 template <
int dim,
typename Number,
typename Number2,
bool do_renumber = true>
2491 const ::ndarray<Number2, 2, dim> *shapes,
2494 const std::vector<unsigned int> &renumber = {})
2496 static_assert(dim >= 0 && dim <= 3,
"Only dim=0,1,2,3 implemented");
2508 Number3 result = {};
2509 for (
int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
2511 Number3 inner_result;
2516 do_interpolate_xy_value<dim, 2, Number2, Number, do_renumber>(
2517 values, renumber, shapes, n_shapes, i);
2518 else if (n_shapes == 3)
2520 do_interpolate_xy_value<dim, 3, Number2, Number, do_renumber>(
2521 values, renumber, shapes, n_shapes, i);
2522 else if (n_shapes == 4)
2524 do_interpolate_xy_value<dim, 4, Number2, Number, do_renumber>(
2525 values, renumber, shapes, n_shapes, i);
2526 else if (n_shapes == 5)
2528 do_interpolate_xy_value<dim, 5, Number2, Number, do_renumber>(
2529 values, renumber, shapes, n_shapes, i);
2530 else if (n_shapes == 6)
2532 do_interpolate_xy_value<dim, 6, Number2, Number, do_renumber>(
2533 values, renumber, shapes, n_shapes, i);
2537 values, renumber, shapes, n_shapes, i);
2541 result += inner_result * shapes[i2][0][2];
2545 result = inner_result;
2554 template <
int dim,
typename Number,
typename Number2,
int str
ide = 1>
2559 static_assert(dim >= 0 && dim <= 3,
"Only dim=0,1,2,3 implemented");
2574 const Number3 val10 = Number3(
values[stride] -
values[0]);
2575 const Number3 val32 = Number3(
values[3 * stride] -
values[2 * stride]);
2576 const Number3 tmp0 = Number3(
values[0]) + p[0] * val10;
2577 const Number3 tmp1 = Number3(
values[2 * stride]) + p[0] * val32;
2578 return tmp0 + p[1] * (tmp1 - tmp0);
2582 const Number3 val10 = Number3(
values[stride] -
values[0]);
2583 const Number3 val32 = Number3(
values[3 * stride] -
values[2 * stride]);
2584 const Number3 tmp0 = Number3(
values[0]) + p[0] * val10;
2585 const Number3 tmp1 = Number3(
values[2 * stride]) + p[0] * val32;
2586 const Number3 tmpy0 = tmp0 + p[1] * (tmp1 - tmp0);
2588 const Number3 val54 = Number3(
values[5 * stride] -
values[4 * stride]);
2589 const Number3 val76 = Number3(
values[7 * stride] -
values[6 * stride]);
2590 const Number3 tmp2 = Number3(
values[4 * stride]) + p[0] * val54;
2591 const Number3 tmp3 = Number3(
values[6 * stride]) + p[0] * val76;
2592 const Number3 tmpy1 = tmp2 + p[1] * (tmp3 - tmp2);
2594 return tmpy0 + p[2] * (tmpy1 - tmpy0);
2603 template <
int dim,
typename Number,
typename Number2>
2607 const std::vector<Number> &
values,
2609 const bool d_linear =
false,
2610 const std::vector<unsigned int> &renumber = {})
2620 std::array<::ndarray<Number2, 2, dim>, 200> shapes;
2621 const int n_shapes = poly.size();
2622 std::array<Number2, dim>
point;
2623 for (
unsigned int d = 0;
d < dim; ++
d)
2625 for (
int i = 0; i < n_shapes; ++i)
2626 poly[i].values_of_array(
point, 0, shapes[i].data());
2627 result = evaluate_tensor_product_value_shapes<dim, Number, Number2>(
2628 shapes.data(), n_shapes,
values.data(), renumber);
2639 template <
int derivative_order,
typename Number,
typename Number2>
2643 const std::vector<Number> &
values,
2645 const std::vector<unsigned int> &renumber = {})
2649 const int n_shapes = poly.size();
2651 Assert(renumber.empty() || renumber.size() ==
values.size(),
2654 std::array<Number2, derivative_order + 1> shapes;
2656 if (renumber.empty())
2657 for (
int i = 0; i < n_shapes; ++i)
2659 poly[i].value(p[0], derivative_order, shapes.data());
2660 result[0] += shapes[derivative_order] *
values[i];
2663 for (
int i = 0; i < n_shapes; ++i)
2665 poly[i].value(p[0], derivative_order, shapes.data());
2666 result[0] += shapes[derivative_order] *
values[renumber[i]];
2677 template <
int derivative_order,
typename Number,
typename Number2>
2679 derivative_order + 1,
2683 const std::vector<Number> &
values,
2685 const std::vector<unsigned int> &renumber = {})
2688 constexpr
int dim = 2;
2690 const int n_shapes = poly.size();
2692 Assert(renumber.empty() || renumber.size() ==
values.size(),
2698 std::array<Number2, dim>
point;
2699 for (
unsigned int d = 0;
d < dim; ++
d)
2701 for (
int i = 0; i < n_shapes; ++i)
2702 poly[i].values_of_array(
point, derivative_order, &shapes[i][0]);
2705 for (
int i1 = 0, i = 0; i1 < n_shapes; ++i1)
2708 if (renumber.empty())
2709 for (
int i0 = 0; i0 < n_shapes; ++i0, ++i)
2710 for (
unsigned int d = 0;
d <= derivative_order; ++
d)
2711 result_x[
d] += shapes[i0][
d][0] *
values[i];
2713 for (
int i0 = 0; i0 < n_shapes; ++i0, ++i)
2714 for (
unsigned int d = 0;
d <= derivative_order; ++
d)
2715 result_x[
d] += shapes[i0][
d][0] *
values[renumber[i]];
2717 for (
unsigned int d = 0;
d <= derivative_order; ++
d)
2718 result[
d] += shapes[i1][
d][1] * result_x[derivative_order -
d];
2729 template <
int derivative_order,
typename Number,
typename Number2>
2731 ((derivative_order + 1) * (derivative_order + 2)) / 2,
2735 const std::vector<Number> &
values,
2737 const std::vector<unsigned int> &renumber = {})
2740 constexpr
int dim = 3;
2741 constexpr
int n_derivatives =
2742 ((derivative_order + 1) * (derivative_order + 2)) / 2;
2744 const int n_shapes = poly.size();
2746 Assert(renumber.empty() || renumber.size() ==
values.size(),
2752 std::array<Number2, dim>
point;
2753 for (
unsigned int d = 0;
d < dim; ++
d)
2755 for (
int i = 0; i < n_shapes; ++i)
2756 poly[i].values_of_array(
point, derivative_order, &shapes[i][0]);
2759 for (
int i2 = 0, i = 0; i2 < n_shapes; ++i2)
2762 for (
int i1 = 0; i1 < n_shapes; ++i1)
2766 if (renumber.empty())
2767 for (
int i0 = 0; i0 < n_shapes; ++i0, ++i)
2768 for (
unsigned int d = 0;
d <= derivative_order; ++
d)
2769 result_x[
d] += shapes[i0][
d][0] *
values[i];
2771 for (
int i0 = 0; i0 < n_shapes; ++i0, ++i)
2772 for (
unsigned int d = 0;
d <= derivative_order; ++
d)
2773 result_x[
d] += shapes[i0][
d][0] *
values[renumber[i]];
2779 for (
unsigned int d = 0, c = 0;
d <= derivative_order; ++
d)
2780 for (
unsigned int e =
d;
e <= derivative_order; ++
e, ++c)
2782 shapes[i1][
e -
d][1] * result_x[derivative_order -
e];
2786 for (
unsigned int d = 0, c = 0;
d <= derivative_order; ++
d)
2787 for (
unsigned int e =
d;
e <= derivative_order; ++
e, ++c)
2788 result[c] += shapes[i2][
d][2] * result_xy[c];
2795 template <
int dim,
typename Number,
typename Number2>
2799 const std::vector<Number> &
values,
2801 const std::vector<unsigned int> &renumber = {})
2803 static_assert(dim >= 1 && dim <= 3,
"Only dim=1,2,3 implemented");
2805 const auto hessian =
2806 evaluate_tensor_product_higher_derivatives<2>(poly,
values, p, renumber);
2811 result[0][0] = hessian[0];
2816 for (
unsigned int d = 0, c = 0;
d < 2; ++
d)
2817 for (
unsigned int e =
d;
e < 2; ++
e, ++c)
2818 result[
d][
e] = hessian[c];
2821 for (
unsigned int d = 0;
d < 2; ++
d)
2822 result[
d][2] = hessian[3 +
d];
2823 result[2][2] = hessian[5];
2848 const ::ndarray<Number, 2, dim> *shapes,
2849 const std::array<Number2, 2 + n_values> &test_grads_value,
2850 const int n_shapes_runtime,
2853 static_assert(0 <= dim && dim <= 3,
"Only dim=0,1,2,3 implemented");
2854 static_assert(1 <= n_values && n_values <= 2,
2855 "Only n_values=1,2 implemented");
2861 constexpr
unsigned int array_size = length > 0 ? length : 1;
2862 std::array<Number, array_size> shape_values_x;
2863 std::array<Number, array_size> shape_derivs_x;
2864 for (
unsigned int j = 0; j < array_size; ++j)
2866 shape_values_x[j] = shapes[j][0][0];
2867 shape_derivs_x[j] = shapes[j][1][0];
2869 for (
int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1)
2871 const Number2 test_value_y =
2872 dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] +
2873 test_grads_value[1] * shapes[i1][1][1]) :
2874 test_grads_value[2];
2875 const Number2 test_grad_xy =
2876 dim > 1 ? test_grads_value[0] * shapes[i1][0][1] :
2877 test_grads_value[0];
2878 Number2 test_value_y_2;
2880 test_value_y_2 = dim > 1 ?
2881 test_grads_value[3] * shapes[i1][0][1] :
2882 test_grads_value[3];
2884 Number2 *values_ptr =
values + i + i1 * length;
2885 Number2 *values_ptr_2 =
2886 n_values > 1 ? values_ptr +
Utilities::pow(length, dim) :
nullptr;
2887 for (
int i0 = 0; i0 < length; ++i0)
2890 values_ptr[i0] += shape_values_x[i0] * test_value_y;
2892 values_ptr[i0] = shape_values_x[i0] * test_value_y;
2893 values_ptr[i0] += shape_derivs_x[i0] * test_grad_xy;
2897 values_ptr_2[i0] += shape_values_x[i0] * test_value_y_2;
2899 values_ptr_2[i0] = shape_values_x[i0] * test_value_y_2;
2903 i += (dim > 1 ? length * length : length);
2907 for (
int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1)
2909 const Number2 test_value_y =
2910 dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] +
2911 test_grads_value[1] * shapes[i1][1][1]) :
2912 test_grads_value[2];
2913 const Number2 test_grad_xy =
2914 dim > 1 ? test_grads_value[0] * shapes[i1][0][1] :
2915 test_grads_value[0];
2916 Number2 test_value_y_2;
2918 test_value_y_2 = dim > 1 ?
2919 test_grads_value[3] * shapes[i1][0][1] :
2920 test_grads_value[3];
2922 Number2 *values_ptr =
values + i + i1 * n_shapes_runtime;
2923 Number2 *values_ptr_2 =
2925 values_ptr + Utilities::fixed_power<dim>(n_shapes_runtime) :
2927 for (
int i0 = 0; i0 < n_shapes_runtime; ++i0)
2930 values_ptr[i0] += shapes[i0][0][0] * test_value_y;
2932 values_ptr[i0] = shapes[i0][0][0] * test_value_y;
2933 values_ptr[i0] += shapes[i0][1][0] * test_grad_xy;
2937 values_ptr_2[i0] += shapes[i0][0][0] * test_value_y_2;
2939 values_ptr_2[i0] = shapes[i0][0][0] * test_value_y_2;
2943 i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime);
2960 const ::ndarray<Number, 2, dim> *shapes,
2962 const Number2 *value,
2966 static_assert(0 <= dim && dim <= 3,
"Only dim=0,1,2,3 implemented");
2967 static_assert(1 <= n_values && n_values <= 2,
2968 "Only n_values=1,2 implemented");
2990 std::array<Number2, 2 + n_values> test_grads_value;
2991 for (
int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
2994 test_grads_value[0] =
2997 test_grads_value[1] = dim > 2 ?
gradient[1] * shapes[i2][0][2] :
2998 (dim > 1 ?
gradient[1] : Number2());
3000 test_grads_value[2] =
3002 (
value[0] * shapes[i2][0][2] +
gradient[2] * shapes[i2][1][2]) :
3006 test_grads_value[3] =
3007 dim > 2 ?
value[1] * shapes[i2][0][2] :
value[1];
3011 do_apply_test_functions_xy<dim, 2, Number2, Number, add, n_values>(
3012 values, shapes, test_grads_value, n_shapes, i);
3013 else if (n_shapes == 3)
3014 do_apply_test_functions_xy<dim, 3, Number2, Number, add, n_values>(
3015 values, shapes, test_grads_value, n_shapes, i);
3016 else if (n_shapes == 4)
3017 do_apply_test_functions_xy<dim, 4, Number2, Number, add, n_values>(
3018 values, shapes, test_grads_value, n_shapes, i);
3019 else if (n_shapes == 5)
3020 do_apply_test_functions_xy<dim, 5, Number2, Number, add, n_values>(
3021 values, shapes, test_grads_value, n_shapes, i);
3022 else if (n_shapes == 6)
3023 do_apply_test_functions_xy<dim, 6, Number2, Number, add, n_values>(
3024 values, shapes, test_grads_value, n_shapes, i);
3027 values, shapes, test_grads_value, n_shapes, i);
3044 const Number2 *value,
3049 static_assert(0 <= dim && dim <= 3,
"Only dim=0,1,2,3 implemented");
3050 static_assert(1 <= n_values && n_values <= 2,
3051 "Only n_values=1,2 implemented");
3084 const Number2 product =
value[1] * p[0];
3099 const Number2 test_value_y1 =
value[0] * p[1] +
gradient[1];
3100 const Number2 test_value_y0 =
value[0] - test_value_y1;
3101 const Number2 test_grad_xy1 =
gradient[0] * p[1];
3102 const Number2 test_grad_xy0 =
gradient[0] - test_grad_xy1;
3103 const Number2 value0 = p[0] * test_value_y0 + test_grad_xy0;
3104 const Number2 value1 = p[0] * test_value_y1 + test_grad_xy1;
3108 values[0] += test_value_y0 - value0;
3110 values[2] += test_value_y1 - value1;
3115 values[0] = test_value_y0 - value0;
3117 values[2] = test_value_y1 - value1;
3123 const Number2 test_value_y1_2 =
value[1] * p[1];
3124 const Number2 test_value_y0_2 =
value[1] - test_value_y1_2;
3125 const Number2 value0_2 = p[0] * test_value_y0_2;
3126 const Number2 value1_2 = p[0] * test_value_y1_2;
3130 values[4] += test_value_y0_2 - value0_2;
3132 values[6] += test_value_y1_2 - value1_2;
3137 values[4] = test_value_y0_2 - value0_2;
3139 values[6] = test_value_y1_2 - value1_2;
3148 const Number2 test_value_z1 =
value[0] * p[2] +
gradient[2];
3149 const Number2 test_value_z0 =
value[0] - test_value_z1;
3150 const Number2 test_grad_x1 =
gradient[0] * p[2];
3151 const Number2 test_grad_x0 =
gradient[0] - test_grad_x1;
3152 const Number2 test_grad_y1 =
gradient[1] * p[2];
3153 const Number2 test_grad_y0 =
gradient[1] - test_grad_y1;
3155 const Number2 test_value_y01 = test_value_z0 * p[1] + test_grad_y0;
3156 const Number2 test_value_y00 = test_value_z0 - test_value_y01;
3157 const Number2 test_grad_xy01 = test_grad_x0 * p[1];
3158 const Number2 test_grad_xy00 = test_grad_x0 - test_grad_xy01;
3159 const Number2 test_value_y11 = test_value_z1 * p[1] + test_grad_y1;
3160 const Number2 test_value_y10 = test_value_z1 - test_value_y11;
3161 const Number2 test_grad_xy11 = test_grad_x1 * p[1];
3162 const Number2 test_grad_xy10 = test_grad_x1 - test_grad_xy11;
3164 const Number2 value00 = p[0] * test_value_y00 + test_grad_xy00;
3165 const Number2 value01 = p[0] * test_value_y01 + test_grad_xy01;
3166 const Number2 value10 = p[0] * test_value_y10 + test_grad_xy10;
3167 const Number2 value11 = p[0] * test_value_y11 + test_grad_xy11;
3171 values[0] += test_value_y00 - value00;
3173 values[2] += test_value_y01 - value01;
3175 values[4] += test_value_y10 - value10;
3177 values[6] += test_value_y11 - value11;
3182 values[0] = test_value_y00 - value00;
3184 values[2] = test_value_y01 - value01;
3186 values[4] = test_value_y10 - value10;
3188 values[6] = test_value_y11 - value11;
3201 template <
bool is_linear,
3208 const ::ndarray<Number, 2, dim> *shapes,
3209 const unsigned int n_shapes,
3210 const Number2 *value,
3257 template <
int dim,
int length,
typename Number2,
typename Number,
bool add>
3265 const ::ndarray<Number, 2, dim> *shapes,
3266 const Number2 &test_value,
3267 const int n_shapes_runtime,
3272 constexpr
unsigned int array_size = length > 0 ? length : 1;
3273 std::array<Number, array_size> shape_values_x;
3274 for (
unsigned int i1 = 0; i1 < array_size; ++i1)
3275 shape_values_x[i1] = shapes[i1][0][0];
3276 for (
unsigned int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1)
3278 const Number2 test_value_y =
3279 dim > 1 ? test_value * shapes[i1][0][1] : test_value;
3281 Number2 *values_ptr =
values + i + i1 * length;
3282 for (
unsigned int i0 = 0; i0 < length; ++i0)
3285 values_ptr[i0] += shape_values_x[i0] * test_value_y;
3287 values_ptr[i0] = shape_values_x[i0] * test_value_y;
3290 i += (dim > 1 ? length * length : length);
3294 for (
int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1)
3296 const Number2 test_value_y =
3297 dim > 1 ? test_value * shapes[i1][0][1] : test_value;
3299 Number2 *values_ptr =
values + i + i1 * n_shapes_runtime;
3300 for (
int i0 = 0; i0 < n_shapes_runtime; ++i0)
3303 values_ptr[i0] += shapes[i0][0][0] * test_value_y;
3305 values_ptr[i0] = shapes[i0][0][0] * test_value_y;
3308 i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime);
3317 template <
int dim,
typename Number,
typename Number2,
bool add>
3320 const ::ndarray<Number, 2, dim> *shapes,
3322 const Number2 &value,
3325 static_assert(dim >= 0 && dim <= 3,
"Only dim=0,1,2,3 implemented");
3340 for (
int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
3343 test_value = dim > 2 ?
value * shapes[i2][0][2] :
value;
3348 do_apply_test_functions_xy_value<dim, 2, Number2, Number, add>(
3349 values, shapes, test_value, n_shapes, i);
3350 else if (n_shapes == 3)
3351 do_apply_test_functions_xy_value<dim, 3, Number2, Number, add>(
3352 values, shapes, test_value, n_shapes, i);
3353 else if (n_shapes == 4)
3354 do_apply_test_functions_xy_value<dim, 4, Number2, Number, add>(
3355 values, shapes, test_value, n_shapes, i);
3356 else if (n_shapes == 5)
3357 do_apply_test_functions_xy_value<dim, 5, Number2, Number, add>(
3358 values, shapes, test_value, n_shapes, i);
3359 else if (n_shapes == 6)
3360 do_apply_test_functions_xy_value<dim, 6, Number2, Number, add>(
3361 values, shapes, test_value, n_shapes, i);
3364 values, shapes, test_value, n_shapes, i);
3374 template <
int dim,
typename Number,
typename Number2,
bool add>
3380 static_assert(dim >= 0 && dim <= 3,
"Only dim=0,1,2,3 implemented");
3391 const auto x0 = 1. - p[0], x1 = p[0];
3406 const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1];
3408 const auto test_value_y0 =
value * y0;
3409 const auto test_value_y1 =
value * y1;
3413 values[0] += x0 * test_value_y0;
3414 values[1] += x1 * test_value_y0;
3415 values[2] += x0 * test_value_y1;
3416 values[3] += x1 * test_value_y1;
3420 values[0] = x0 * test_value_y0;
3421 values[1] = x1 * test_value_y0;
3422 values[2] = x0 * test_value_y1;
3423 values[3] = x1 * test_value_y1;
3428 const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
3429 z0 = 1. - p[2], z1 = p[2];
3431 const auto test_value_z0 =
value * z0;
3432 const auto test_value_z1 =
value * z1;
3434 const auto test_value_y00 = test_value_z0 * y0;
3435 const auto test_value_y01 = test_value_z0 * y1;
3436 const auto test_value_y10 = test_value_z1 * y0;
3437 const auto test_value_y11 = test_value_z1 * y1;
3441 values[0] += x0 * test_value_y00;
3442 values[1] += x1 * test_value_y00;
3443 values[2] += x0 * test_value_y01;
3444 values[3] += x1 * test_value_y01;
3445 values[4] += x0 * test_value_y10;
3446 values[5] += x1 * test_value_y10;
3447 values[6] += x0 * test_value_y11;
3448 values[7] += x1 * test_value_y11;
3452 values[0] = x0 * test_value_y00;
3453 values[1] = x1 * test_value_y00;
3454 values[2] = x0 * test_value_y01;
3455 values[3] = x1 * test_value_y01;
3456 values[4] = x0 * test_value_y10;
3457 values[5] = x1 * test_value_y10;
3458 values[6] = x0 * test_value_y11;
3459 values[7] = x1 * test_value_y11;
3471 template <
bool is_linear,
int dim,
typename Number,
typename Number2>
3474 const unsigned int n_shapes,
3475 const Number2 &value,
3520 template <
int dim,
int n_po
ints_1d_
template,
typename Number>
3523 const unsigned int n_components,
3524 const int n_points_1d_non_template,
3527 const int n_points_1d = n_points_1d_template != -1 ?
3528 n_points_1d_template :
3529 n_points_1d_non_template;
3534 unsigned int compressed_index[100];
3535 compressed_index[0] = 0;
3536 for (
int i = 1; i < n_points_1d - 1; ++i)
3537 compressed_index[i] = 1;
3538 compressed_index[n_points_1d - 1] = 2;
3540 for (
unsigned int c = 0; c < n_components; ++c)
3541 for (
int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3542 for (
int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
3544 const unsigned int shift =
3545 9 * compressed_index[k] + 3 * compressed_index[j];
3546 data[0] *= weights[
shift];
3549 const Number weight = weights[
shift + 1];
3550 for (
int i = 1; i < n_points_1d - 1; ++i)
3552 data[n_points_1d - 1] *= weights[
shift + 2];
3553 data += n_points_1d;
3558 template <
int dim,
int n_po
ints_1d_
template,
typename Number>
3561 const unsigned int n_components,
3562 const int n_points_1d_non_template,
3565 const int n_points_1d = n_points_1d_template != -1 ?
3566 n_points_1d_template :
3567 n_points_1d_non_template;
3569 Assert((n_points_1d % 2) == 1,
3570 ExcMessage(
"The function can only with add number of points"));
3574 const unsigned int n_inside_1d = n_points_1d / 2;
3576 unsigned int compressed_index[100];
3579 for (
int i = 0; i < n_inside_1d; ++i)
3580 compressed_index[c++] = 0;
3581 compressed_index[c++] = 1;
3582 for (
int i = 0; i < n_inside_1d; ++i)
3583 compressed_index[c++] = 2;
3585 for (
unsigned int c = 0; c < n_components; ++c)
3586 for (
int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3587 for (
int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
3589 const unsigned int shift =
3590 9 * compressed_index[k] + 3 * compressed_index[j];
3593 const Number weight1 = weights[
shift];
3594 for (
int i = 0; i < n_inside_1d; ++i)
3595 data[c++] *= weight1;
3596 data[c++] *= weights[
shift + 1];
3597 const Number weight2 = weights[
shift + 2];
3598 for (
int i = 0; i < n_inside_1d; ++i)
3599 data[c++] *= weight2;
3600 data += n_points_1d;
3605 template <
int dim,
int n_po
ints_1d_
template,
typename Number>
3608 const unsigned int n_components,
3609 const int n_points_1d_non_template,
3612 const int n_points_1d = n_points_1d_template != -1 ?
3613 n_points_1d_template :
3614 n_points_1d_non_template;
3619 unsigned int compressed_index[100];
3620 compressed_index[0] = 0;
3621 for (
int i = 1; i < n_points_1d - 1; ++i)
3622 compressed_index[i] = 1;
3623 compressed_index[n_points_1d - 1] = 2;
3630 const auto check_and_set = [](Number &weight,
const Number &data) {
3631 if (weight == Number(-1.0) || weight == data)
3640 for (
unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
3641 weights[c] = Number(-1.0);
3643 for (
unsigned int c = 0; c < n_components; ++c)
3644 for (
int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3645 for (
int j = 0; j < (dim > 1 ? n_points_1d : 1);
3646 ++j, data += n_points_1d)
3648 const unsigned int shift =
3649 9 * compressed_index[k] + 3 * compressed_index[j];
3651 if (!check_and_set(weights[
shift], data[0]))
3654 for (
int i = 1; i < n_points_1d - 1; ++i)
3655 if (!check_and_set(weights[
shift + 1], data[i]))
3658 if (!check_and_set(weights[
shift + 2], data[n_points_1d - 1]))
3666 template <
int dim,
int n_po
ints_1d_
template,
typename Number>
3670 const unsigned int n_components,
3671 const int n_points_1d_non_template,
3674 const int n_points_1d = n_points_1d_template != -1 ?
3675 n_points_1d_template :
3676 n_points_1d_non_template;
3678 Assert((n_points_1d % 2) == 1,
3679 ExcMessage(
"The function can only with add number of points"));
3683 const unsigned int n_inside_1d = n_points_1d / 2;
3685 unsigned int compressed_index[100];
3688 for (
int i = 0; i < n_inside_1d; ++i)
3689 compressed_index[c++] = 0;
3690 compressed_index[c++] = 1;
3691 for (
int i = 0; i < n_inside_1d; ++i)
3692 compressed_index[c++] = 2;
3699 const auto check_and_set = [](Number &weight,
const Number &data) {
3700 if (weight == Number(-1.0) || weight == data)
3709 for (
unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
3710 weights[c] = Number(-1.0);
3712 for (
unsigned int comp = 0; comp < n_components; ++comp)
3713 for (
int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3714 for (
int j = 0; j < (dim > 1 ? n_points_1d : 1);
3715 ++j, data += n_points_1d)
3717 const unsigned int shift =
3718 9 * compressed_index[k] + 3 * compressed_index[j];
3722 for (
int i = 0; i < n_inside_1d; ++i)
3723 if (!check_and_set(weights[
shift], data[c++]))
3726 if (!check_and_set(weights[
shift + 1], data[c++]))
3729 for (
int i = 0; i < n_inside_1d; ++i)
3730 if (!check_and_set(weights[
shift + 2], data[c++]))
#define DEAL_II_ALWAYS_INLINE
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcMessage(std::string arg1)
@ matrix
Contents is actually a matrix.
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
constexpr T pow(const T base, const int iexp)
std::enable_if_t<(variant==evaluate_general), void > apply_matrix_vector_product(const Number2 *matrix, const Number *in, Number *out)
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, 2+n_values > do_interpolate_xy(const Number *values, const std::vector< unsigned int > &renumber, const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes_runtime, int &i)
void integrate_add_tensor_product_value_linear(const Number2 &value, Number2 *values, const Point< dim, Number > &p)
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_shapes(const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes, const Number *values, const std::vector< unsigned int > &renumber={})
void weight_fe_q_dofs_by_entity(const Number *weights, const unsigned int n_components, const int n_points_1d_non_template, Number *data)
ProductTypeNoPoint< Number, Number2 >::type do_interpolate_xy_value(const Number *values, const std::vector< unsigned int > &renumber, const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes_runtime, int &i)
SymmetricTensor< 2, dim, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_hessian(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const std::vector< unsigned int > &renumber={})
void integrate_tensor_product_value_and_gradient(const ::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 *value, const Tensor< 1, dim, Number2 > &gradient, Number2 *values, const Point< dim, Number > &p, const bool do_add)
void do_apply_test_functions_xy_value(Number2 *values, const ::ndarray< Number, 2, dim > *shapes, const Number2 &test_value, const int n_shapes_runtime, int &i)
void integrate_add_tensor_product_value_and_gradient_linear(const Number2 *value, const Tensor< 1, dim, Number2 > &gradient, Number2 *values, const Point< dim, Number > &p)
void integrate_add_tensor_product_value_shapes(const ::ndarray< Number, 2, dim > *shapes, const int n_shapes, const Number2 &value, Number2 *values)
Tensor< 1, 1, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_higher_derivatives(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< 1, Number2 > &p, const std::vector< unsigned int > &renumber={})
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, dim+n_values > evaluate_tensor_product_value_and_gradient_linear(const Number *values, const Point< dim, Number2 > &p)
std::enable_if_t< contract_onto_face, void > interpolate_to_face(const Number2 *shape_values, const std::array< int, 2 > &n_blocks, const std::array< int, 2 > &steps, const Number *input, Number *DEAL_II_RESTRICT output, const int n_rows_runtime=0, const int stride_runtime=1)
std::pair< typename ProductTypeNoPoint< Number, Number2 >::type, Tensor< 1, dim, typename ProductTypeNoPoint< Number, Number2 >::type > > evaluate_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
void compute_values_of_array(::ndarray< Number, 2, dim > *shapes, const std::vector< Polynomials::Polynomial< double >> &poly, const Point< dim, Number > &p, const unsigned int derivative=1)
void integrate_add_tensor_product_value_and_gradient_shapes(const ::ndarray< Number, 2, dim > *shapes, const int n_shapes, const Number2 *value, const Tensor< 1, dim, Number2 > &gradient, Number2 *values)
void weight_fe_q_dofs_by_entity_shifted(const Number *weights, const unsigned int n_components, const int n_points_1d_non_template, Number *data)
void do_apply_test_functions_xy(Number2 *values, const ::ndarray< Number, 2, dim > *shapes, const std::array< Number2, 2+n_values > &test_grads_value, const int n_shapes_runtime, int &i)
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, dim+n_values > evaluate_tensor_product_value_and_gradient_shapes(const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes, const Number *values, const std::vector< unsigned int > &renumber={})
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_linear(const Number *values, const Point< dim, Number2 > &p)
bool compute_weights_fe_q_dofs_by_entity(const Number *data, const unsigned int n_components, const int n_points_1d_non_template, Number *weights)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
bool compute_weights_fe_q_dofs_by_entity_shifted(const Number *data, const unsigned int n_components, const int n_points_1d_non_template, Number *weights)
void integrate_tensor_product_value(const ::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 &value, Number2 *values, const Point< dim, Number > &p, const bool do_add)
@ evaluate_symmetric_hierarchical
static const unsigned int invalid_unsigned_int
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type
static void normal(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const Number *in, Number *out, const bool add_into_result=false, const int subface_index_1d=0)
static void tangential(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const Number *in, Number *out, const int subface_index_1d=0)
const Number2 * shape_hessians
const Number2 * shape_gradients
const unsigned int n_rows
void hessians_one_line(const Number in[], Number out[]) const
void gradients(const Number *in, Number *out) const
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int n_rows=0, const unsigned int n_columns=0)
const Number2 * shape_values
void values(const Number *in, Number *out) const
void hessians(const Number *in, Number *out) const
void gradients_one_line(const Number in[], Number out[]) const
void values_one_line(const Number in[], Number out[]) const
EvaluatorTensorProduct(const Number2 *shape_values, const Number2 *shape_gradients, const Number2 *shape_hessians, const unsigned int n_rows=0, const unsigned int n_columns=0)
const unsigned int n_columns
void values(const Number in[], Number out[]) const
EvaluatorTensorProduct(const Number2 *shape_values, const Number2 *shape_gradients, const Number2 *shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
void values_one_line(const Number in[], Number out[]) const
void gradients(const Number in[], Number out[]) const
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int=0, const unsigned int=0)
const Number2 * shape_gradients
static constexpr unsigned int n_rows_of_product
void hessians(const Number in[], Number out[]) const
const Number2 * shape_values
const Number2 * shape_hessians
static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out)
static constexpr unsigned int n_columns_of_product
void gradients_one_line(const Number in[], Number out[]) const
void hessians_one_line(const Number in[], Number out[]) const
AlignedVector< Number > shape_values_eo
unsigned int n_q_points_1d
std::array< AlignedVector< Number >, 2 > values_within_subface
typename ProductType< Tensor< 1, dim, Number >, Number2 >::type type
typename ProductType< Number, Number2 >::type type