17 #ifndef dealii_matrix_free_tensor_product_kernels_h
18 #define dealii_matrix_free_tensor_product_kernels_h
122 typename Number2 = Number>
154 typename Number2 = Number>
189 static constexpr
unsigned int n_rows_of_product =
191 static constexpr
unsigned int n_columns_of_product =
199 : shape_values(nullptr)
200 , shape_gradients(nullptr)
201 , shape_hessians(nullptr)
210 const unsigned int dummy1 = 0,
211 const unsigned int dummy2 = 0)
212 : shape_values(shape_values.
begin())
213 , shape_gradients(shape_gradients.
begin())
214 , shape_hessians(shape_hessians.
begin())
221 shape_values.
size() == n_rows * n_columns ||
222 shape_values.
size() == 3 * n_rows,
225 shape_gradients.
size() == n_rows * n_columns,
228 shape_hessians.
size() == n_rows * n_columns,
234 template <
int direction,
bool contract_over_rows,
bool add>
236 values(
const Number in[], Number out[])
const
238 apply<direction, contract_over_rows, add>(shape_values, in, out);
241 template <
int direction,
bool contract_over_rows,
bool add>
245 apply<direction, contract_over_rows, add>(shape_gradients, in, out);
248 template <
int direction,
bool contract_over_rows,
bool add>
252 apply<direction, contract_over_rows, add>(shape_hessians, in, out);
255 template <
int direction,
bool contract_over_rows,
bool add>
260 apply<direction, contract_over_rows, add, true>(shape_values, in, out);
263 template <
int direction,
bool contract_over_rows,
bool add>
268 apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
271 template <
int direction,
bool contract_over_rows,
bool add>
276 apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
303 template <
int direction,
304 bool contract_over_rows,
306 bool one_line =
false>
341 template <
int face_direction,
342 bool contract_onto_face,
362 template <
int direction,
bool contract_over_rows,
bool add,
bool one_line>
374 static_assert(one_line ==
false || direction == dim - 1,
375 "Single-line evaluation only works for direction=dim-1.");
376 Assert(shape_data !=
nullptr,
378 "The given array shape_data must not be the null pointer!"));
379 Assert(dim == direction + 1 || one_line ==
true || n_rows == n_columns ||
381 ExcMessage(
"In-place operation only supported for "
382 "n_rows==n_columns or single-line interpolation"));
384 constexpr
int mm = contract_over_rows ? n_rows : n_columns,
385 nn = contract_over_rows ? n_columns : n_rows;
388 constexpr
int n_blocks1 = one_line ? 1 : stride;
389 constexpr
int n_blocks2 =
390 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
392 for (
int i2 = 0; i2 < n_blocks2; ++i2)
394 for (
int i1 = 0; i1 < n_blocks1; ++i1)
397 for (
int i = 0; i < mm; ++i)
398 x[i] = in[stride * i];
399 for (
int col = 0; col < nn; ++col)
402 if (contract_over_rows ==
true)
403 val0 = shape_data[col];
405 val0 = shape_data[col * n_columns];
406 Number res0 = val0 * x[0];
407 for (
int i = 1; i < mm; ++i)
409 if (contract_over_rows ==
true)
410 val0 = shape_data[i * n_columns + col];
412 val0 = shape_data[col * n_columns + i];
416 out[stride * col] += res0;
418 out[stride * col] = res0;
421 if (one_line ==
false)
427 if (one_line ==
false)
429 in += stride * (mm - 1);
430 out += stride * (nn - 1);
442 template <
int face_direction,
443 bool contract_onto_face,
457 static_assert(max_derivative >= 0 && max_derivative < 3,
458 "Only derivative orders 0-2 implemented");
459 Assert(shape_values !=
nullptr,
461 "The given array shape_values must not be the null pointer."));
463 constexpr
int n_blocks1 = (dim > 1 ? n_rows : 1);
464 constexpr
int n_blocks2 = (dim > 2 ? n_rows : 1);
471 for (
int i2 = 0; i2 < n_blocks2; ++i2)
473 for (
int i1 = 0; i1 < n_blocks1; ++i1)
475 if (contract_onto_face ==
true)
477 Number res0 = shape_values[0] * in[0];
479 if (max_derivative > 0)
480 res1 = shape_values[n_rows] * in[0];
481 if (max_derivative > 1)
482 res2 = shape_values[2 * n_rows] * in[0];
483 for (
int ind = 1; ind < n_rows; ++ind)
485 res0 += shape_values[ind] * in[in_stride * ind];
486 if (max_derivative > 0)
487 res1 += shape_values[ind + n_rows] * in[in_stride * ind];
488 if (max_derivative > 1)
490 shape_values[ind + 2 * n_rows] * in[in_stride * ind];
495 if (max_derivative > 0)
496 out[out_stride] += res1;
497 if (max_derivative > 1)
498 out[2 * out_stride] += res2;
503 if (max_derivative > 0)
504 out[out_stride] = res1;
505 if (max_derivative > 1)
506 out[2 * out_stride] = res2;
511 for (
int col = 0; col < n_rows; ++col)
514 out[col * in_stride] += shape_values[col] * in[0];
516 out[col * in_stride] = shape_values[col] * in[0];
517 if (max_derivative > 0)
518 out[col * in_stride] +=
519 shape_values[col + n_rows] * in[out_stride];
520 if (max_derivative > 1)
521 out[col * in_stride] +=
522 shape_values[col + 2 * n_rows] * in[2 * out_stride];
529 switch (face_direction)
532 in += contract_onto_face ? n_rows : 1;
533 out += contract_onto_face ? 1 : n_rows;
543 if (contract_onto_face)
559 if (face_direction == 1 && dim == 3)
561 if (contract_onto_face)
563 in += n_rows * (n_rows - 1);
564 out -= n_rows * n_rows - 1;
568 out += n_rows * (n_rows - 1);
569 in -= n_rows * n_rows - 1;
590 template <
int dim,
typename Number,
typename Number2>
593 static constexpr
unsigned int n_rows_of_product =
595 static constexpr
unsigned int n_columns_of_product =
603 : shape_values(nullptr)
604 , shape_gradients(nullptr)
605 , shape_hessians(nullptr)
616 const unsigned int n_rows,
617 const unsigned int n_columns)
618 : shape_values(shape_values.
begin())
619 , shape_gradients(shape_gradients.
begin())
620 , shape_hessians(shape_hessians.
begin())
622 , n_columns(n_columns)
629 shape_values.
size() == n_rows * n_columns ||
630 shape_values.
size() == n_rows * 3,
633 shape_gradients.
size() == n_rows * n_columns,
636 shape_hessians.
size() == n_rows * n_columns,
644 const Number2 * shape_gradients,
645 const Number2 * shape_hessians,
646 const unsigned int n_rows,
647 const unsigned int n_columns)
648 : shape_values(shape_values)
649 , shape_gradients(shape_gradients)
650 , shape_hessians(shape_hessians)
652 , n_columns(n_columns)
655 template <
int direction,
bool contract_over_rows,
bool add>
657 values(
const Number *in, Number *out)
const
659 apply<direction, contract_over_rows, add>(shape_values, in, out);
662 template <
int direction,
bool contract_over_rows,
bool add>
666 apply<direction, contract_over_rows, add>(shape_gradients, in, out);
669 template <
int direction,
bool contract_over_rows,
bool add>
673 apply<direction, contract_over_rows, add>(shape_hessians, in, out);
676 template <
int direction,
bool contract_over_rows,
bool add>
681 apply<direction, contract_over_rows, add, true>(shape_values, in, out);
684 template <
int direction,
bool contract_over_rows,
bool add>
689 apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
692 template <
int direction,
bool contract_over_rows,
bool add>
697 apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
700 template <
int direction,
701 bool contract_over_rows,
703 bool one_line =
false>
709 template <
int face_direction,
710 bool contract_onto_face,
726 template <
int dim,
typename Number,
typename Number2>
727 template <
int direction,
bool contract_over_rows,
bool add,
bool one_line>
734 static_assert(one_line ==
false || direction == dim - 1,
735 "Single-line evaluation only works for direction=dim-1.");
736 Assert(shape_data !=
nullptr,
738 "The given array shape_data must not be the null pointer!"));
739 Assert(dim == direction + 1 || one_line ==
true || n_rows == n_columns ||
741 ExcMessage(
"In-place operation only supported for "
742 "n_rows==n_columns or single-line interpolation"));
744 const int mm = contract_over_rows ? n_rows : n_columns,
745 nn = contract_over_rows ? n_columns : n_rows;
748 direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
749 const int n_blocks1 = one_line ? 1 : stride;
750 const int n_blocks2 = direction >= dim - 1 ?
760 if (contract_over_rows && n_rows == 2)
762 const Number2 *shape_data_1 = shape_data + n_columns;
763 for (
int i2 = 0; i2 < n_blocks2; ++i2)
765 for (
int i1 = 0; i1 < n_blocks1; ++i1)
767 const Number x0 = in[0], x1 = in[stride];
768 for (
int col = 0; col < nn; ++col)
770 const Number result =
771 shape_data[col] * x0 + shape_data_1[col] * x1;
773 out[stride * col] += result;
775 out[stride * col] = result;
778 if (one_line ==
false)
784 if (one_line ==
false)
786 in += stride * (mm - 1);
787 out += stride * (nn - 1);
792 else if (contract_over_rows && n_rows == 3)
794 const Number2 *shape_data_1 = shape_data + n_columns;
795 const Number2 *shape_data_2 = shape_data + 2 * n_columns;
796 for (
int i2 = 0; i2 < n_blocks2; ++i2)
798 for (
int i1 = 0; i1 < n_blocks1; ++i1)
800 const Number x0 = in[0], x1 = in[stride], x2 = in[2 * stride];
801 for (
int col = 0; col < nn; ++col)
803 const Number result = shape_data[col] * x0 +
804 shape_data_1[col] * x1 +
805 shape_data_2[col] * x2;
807 out[stride * col] += result;
809 out[stride * col] = result;
812 if (one_line ==
false)
818 if (one_line ==
false)
820 in += stride * (mm - 1);
821 out += stride * (nn - 1);
827 for (
int i2 = 0; i2 < n_blocks2; ++i2)
829 for (
int i1 = 0; i1 < n_blocks1; ++i1)
832 for (
int i = 0; i < mm; ++i)
833 x[i] = in[stride * i];
834 for (
int col = 0; col < nn; ++col)
837 if (contract_over_rows ==
true)
838 val0 = shape_data[col];
840 val0 = shape_data[col * n_columns];
841 Number res0 = val0 * x[0];
842 for (
int i = 1; i < mm; ++i)
844 if (contract_over_rows ==
true)
845 val0 = shape_data[i * n_columns + col];
847 val0 = shape_data[col * n_columns + i];
851 out[stride * col] += res0;
853 out[stride * col] = res0;
856 if (one_line ==
false)
862 if (one_line ==
false)
864 in += stride * (mm - 1);
865 out += stride * (nn - 1);
872 template <
int dim,
typename Number,
typename Number2>
873 template <
int face_direction,
874 bool contract_onto_face,
882 Assert(shape_values !=
nullptr,
884 "The given array shape_data must not be the null pointer!"));
885 static_assert(dim > 0 && dim < 4,
"Only dim=1,2,3 supported");
886 const int n_blocks1 = dim > 1 ? n_rows : 1;
887 const int n_blocks2 = dim > 2 ? n_rows : 1;
890 const int in_stride =
891 face_direction > 0 ? Utilities::fixed_power<face_direction>(n_rows) : 1;
892 const int out_stride =
895 for (
int i2 = 0; i2 < n_blocks2; ++i2)
897 for (
int i1 = 0; i1 < n_blocks1; ++i1)
899 if (contract_onto_face ==
true)
901 Number res0 = shape_values[0] * in[0];
903 if (max_derivative > 0)
904 res1 = shape_values[n_rows] * in[0];
905 if (max_derivative > 1)
906 res2 = shape_values[2 * n_rows] * in[0];
907 for (
unsigned int ind = 1; ind < n_rows; ++ind)
909 res0 += shape_values[ind] * in[in_stride * ind];
910 if (max_derivative > 0)
911 res1 += shape_values[ind + n_rows] * in[in_stride * ind];
912 if (max_derivative > 1)
914 shape_values[ind + 2 * n_rows] * in[in_stride * ind];
919 if (max_derivative > 0)
920 out[out_stride] += res1;
921 if (max_derivative > 1)
922 out[2 * out_stride] += res2;
927 if (max_derivative > 0)
928 out[out_stride] = res1;
929 if (max_derivative > 1)
930 out[2 * out_stride] = res2;
935 for (
unsigned int col = 0; col < n_rows; ++col)
938 out[col * in_stride] += shape_values[col] * in[0];
940 out[col * in_stride] = shape_values[col] * in[0];
941 if (max_derivative > 0)
942 out[col * in_stride] +=
943 shape_values[col + n_rows] * in[out_stride];
944 if (max_derivative > 1)
945 out[col * in_stride] +=
946 shape_values[col + 2 * n_rows] * in[2 * out_stride];
953 switch (face_direction)
956 in += contract_onto_face ? n_rows : 1;
957 out += contract_onto_face ? 1 : n_rows;
967 if (contract_onto_face)
981 if (face_direction == 1 && dim == 3)
984 if (contract_onto_face)
986 in += n_rows * (n_rows - 1);
987 out -= n_rows * n_rows - 1;
991 out += n_rows * (n_rows - 1);
992 in -= n_rows * n_rows - 1;
1032 static constexpr
unsigned int n_rows_of_product =
1034 static constexpr
unsigned int n_columns_of_product =
1043 const unsigned int dummy1 = 0,
1044 const unsigned int dummy2 = 0)
1045 : shape_values(shape_values.
begin())
1046 , shape_gradients(shape_gradients.
begin())
1047 , shape_hessians(shape_hessians.
begin())
1050 shape_values.
size() == n_rows * n_columns,
1053 shape_gradients.
size() == n_rows * n_columns,
1056 shape_hessians.
size() == n_rows * n_columns,
1062 template <
int direction,
bool contract_over_rows,
bool add>
1064 values(
const Number in[], Number out[])
const;
1066 template <
int direction,
bool contract_over_rows,
bool add>
1068 gradients(
const Number in[], Number out[])
const;
1070 template <
int direction,
bool contract_over_rows,
bool add>
1072 hessians(
const Number in[], Number out[])
const;
1105 template <
int direction,
bool contract_over_rows,
bool add>
1112 Number2>
::values(
const Number in[], Number out[])
const
1116 constexpr
int mm = contract_over_rows ? n_rows : n_columns,
1117 nn = contract_over_rows ? n_columns : n_rows;
1118 constexpr
int n_cols = nn / 2;
1119 constexpr
int mid = mm / 2;
1122 constexpr
int n_blocks1 = stride;
1123 constexpr
int n_blocks2 =
1124 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1126 for (
int i2 = 0; i2 < n_blocks2; ++i2)
1128 for (
int i1 = 0; i1 < n_blocks1; ++i1)
1130 for (
int col = 0; col < n_cols; ++col)
1133 Number in0, in1, res0, res1;
1134 if (contract_over_rows ==
true)
1136 val0 = shape_values[col];
1137 val1 = shape_values[nn - 1 - col];
1141 val0 = shape_values[col * n_columns];
1142 val1 = shape_values[(col + 1) * n_columns - 1];
1147 in1 = in[stride * (mm - 1)];
1152 for (
int ind = 1; ind < mid; ++ind)
1154 if (contract_over_rows ==
true)
1156 val0 = shape_values[ind * n_columns + col];
1157 val1 = shape_values[ind * n_columns + nn - 1 - col];
1161 val0 = shape_values[col * n_columns + ind];
1163 shape_values[(col + 1) * n_columns - 1 - ind];
1165 in0 = in[stride * ind];
1166 in1 = in[stride * (mm - 1 - ind)];
1174 res0 = res1 = Number();
1175 if (contract_over_rows ==
true)
1179 val0 = shape_values[mid * n_columns + col];
1180 in1 = val0 * in[stride * mid];
1187 if (mm % 2 == 1 && nn % 2 == 0)
1189 val0 = shape_values[col * n_columns + mid];
1190 in1 = val0 * in[stride * mid];
1197 out[stride * col] += res0;
1198 out[stride * (nn - 1 - col)] += res1;
1202 out[stride * col] = res0;
1203 out[stride * (nn - 1 - col)] = res1;
1206 if (contract_over_rows ==
true && nn % 2 == 1 && mm % 2 == 1)
1209 out[stride * n_cols] += in[stride * mid];
1211 out[stride * n_cols] = in[stride * mid];
1213 else if (contract_over_rows ==
true && nn % 2 == 1)
1216 Number2 val0 = shape_values[n_cols];
1219 res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1220 for (
int ind = 1; ind < mid; ++ind)
1222 val0 = shape_values[ind * n_columns + n_cols];
1223 res0 += val0 * (in[stride * ind] +
1224 in[stride * (mm - 1 - ind)]);
1230 out[stride * n_cols] += res0;
1232 out[stride * n_cols] = res0;
1234 else if (contract_over_rows ==
false && nn % 2 == 1)
1239 Number2 val0 = shape_values[n_cols * n_columns];
1240 res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1241 for (
int ind = 1; ind < mid; ++ind)
1243 val0 = shape_values[n_cols * n_columns + ind];
1244 Number in1 = val0 * (in[stride * ind] +
1245 in[stride * (mm - 1 - ind)]);
1249 res0 += in[stride * mid];
1254 out[stride * n_cols] += res0;
1256 out[stride * n_cols] = res0;
1262 in += stride * (mm - 1);
1263 out += stride * (nn - 1);
1293 template <
int direction,
bool contract_over_rows,
bool add>
1305 constexpr
int mm = contract_over_rows ? n_rows : n_columns,
1306 nn = contract_over_rows ? n_columns : n_rows;
1307 constexpr
int n_cols = nn / 2;
1308 constexpr
int mid = mm / 2;
1311 constexpr
int n_blocks1 = stride;
1312 constexpr
int n_blocks2 =
1313 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1315 for (
int i2 = 0; i2 < n_blocks2; ++i2)
1317 for (
int i1 = 0; i1 < n_blocks1; ++i1)
1319 for (
int col = 0; col < n_cols; ++col)
1322 Number in0, in1, res0, res1;
1323 if (contract_over_rows ==
true)
1325 val0 = shape_gradients[col];
1326 val1 = shape_gradients[nn - 1 - col];
1330 val0 = shape_gradients[col * n_columns];
1331 val1 = shape_gradients[(nn - col - 1) * n_columns];
1336 in1 = in[stride * (mm - 1)];
1341 for (
int ind = 1; ind < mid; ++ind)
1343 if (contract_over_rows ==
true)
1345 val0 = shape_gradients[ind * n_columns + col];
1347 shape_gradients[ind * n_columns + nn - 1 - col];
1351 val0 = shape_gradients[col * n_columns + ind];
1353 shape_gradients[(nn - col - 1) * n_columns + ind];
1355 in0 = in[stride * ind];
1356 in1 = in[stride * (mm - 1 - ind)];
1364 res0 = res1 = Number();
1367 if (contract_over_rows ==
true)
1368 val0 = shape_gradients[mid * n_columns + col];
1370 val0 = shape_gradients[col * n_columns + mid];
1371 in1 = val0 * in[stride * mid];
1377 out[stride * col] += res0;
1378 out[stride * (nn - 1 - col)] += res1;
1382 out[stride * col] = res0;
1383 out[stride * (nn - 1 - col)] = res1;
1390 if (contract_over_rows ==
true)
1391 val0 = shape_gradients[n_cols];
1393 val0 = shape_gradients[n_cols * n_columns];
1394 res0 = val0 * (in[0] - in[stride * (mm - 1)]);
1395 for (
int ind = 1; ind < mid; ++ind)
1397 if (contract_over_rows ==
true)
1398 val0 = shape_gradients[ind * n_columns + n_cols];
1400 val0 = shape_gradients[n_cols * n_columns + ind];
1402 val0 * (in[stride * ind] - in[stride * (mm - 1 - ind)]);
1406 out[stride * n_cols] += res0;
1408 out[stride * n_cols] = res0;
1414 in += stride * (mm - 1);
1415 out += stride * (nn - 1);
1429 template <
int direction,
bool contract_over_rows,
bool add>
1441 constexpr
int mm = contract_over_rows ? n_rows : n_columns;
1442 constexpr
int nn = contract_over_rows ? n_columns : n_rows;
1443 constexpr
int n_cols = nn / 2;
1444 constexpr
int mid = mm / 2;
1447 constexpr
int n_blocks1 = stride;
1448 constexpr
int n_blocks2 =
1449 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1451 for (
int i2 = 0; i2 < n_blocks2; ++i2)
1453 for (
int i1 = 0; i1 < n_blocks1; ++i1)
1455 for (
int col = 0; col < n_cols; ++col)
1458 Number in0, in1, res0, res1;
1459 if (contract_over_rows ==
true)
1461 val0 = shape_hessians[col];
1462 val1 = shape_hessians[nn - 1 - col];
1466 val0 = shape_hessians[col * n_columns];
1467 val1 = shape_hessians[(col + 1) * n_columns - 1];
1472 in1 = in[stride * (mm - 1)];
1477 for (
int ind = 1; ind < mid; ++ind)
1479 if (contract_over_rows ==
true)
1481 val0 = shape_hessians[ind * n_columns + col];
1483 shape_hessians[ind * n_columns + nn - 1 - col];
1487 val0 = shape_hessians[col * n_columns + ind];
1489 shape_hessians[(col + 1) * n_columns - 1 - ind];
1491 in0 = in[stride * ind];
1492 in1 = in[stride * (mm - 1 - ind)];
1500 res0 = res1 = Number();
1503 if (contract_over_rows ==
true)
1504 val0 = shape_hessians[mid * n_columns + col];
1506 val0 = shape_hessians[col * n_columns + mid];
1507 in1 = val0 * in[stride * mid];
1513 out[stride * col] += res0;
1514 out[stride * (nn - 1 - col)] += res1;
1518 out[stride * col] = res0;
1519 out[stride * (nn - 1 - col)] = res1;
1526 if (contract_over_rows ==
true)
1527 val0 = shape_hessians[n_cols];
1529 val0 = shape_hessians[n_cols * n_columns];
1532 res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1533 for (
int ind = 1; ind < mid; ++ind)
1535 if (contract_over_rows ==
true)
1536 val0 = shape_hessians[ind * n_columns + n_cols];
1538 val0 = shape_hessians[n_cols * n_columns + ind];
1539 Number in1 = val0 * (in[stride * ind] +
1540 in[stride * (mm - 1 - ind)]);
1548 if (contract_over_rows ==
true)
1549 val0 = shape_hessians[mid * n_columns + n_cols];
1551 val0 = shape_hessians[n_cols * n_columns + mid];
1552 res0 += val0 * in[stride * mid];
1555 out[stride * n_cols] += res0;
1557 out[stride * n_cols] = res0;
1563 in += stride * (mm - 1);
1564 out += stride * (nn - 1);
1572 int n_columns_static,
1576 bool contract_over_rows,
1582 const int n_columns_in,
1587 static_assert(type < 3,
"Only three variants type=0,1,2 implemented");
1588 static_assert(one_line ==
false || direction == dim - 1,
1589 "Single-line evaluation only works for direction=dim-1.");
1591 const int n_rows = n_rows_static == -1 ? n_rows_in : n_rows_static;
1592 const int n_columns =
1593 n_columns_static == -1 ? n_columns_in : n_columns_static;
1595 Assert(dim == direction + 1 || one_line ==
true || n_rows == n_columns ||
1597 ExcMessage(
"In-place operation only supported for "
1598 "n_rows==n_columns or single-line interpolation"));
1604 const int nn = contract_over_rows ? n_columns : n_rows;
1605 const int mm = contract_over_rows ? n_rows : n_columns;
1606 constexpr
int mm_static =
1607 contract_over_rows ? n_rows_static : n_columns_static;
1608 const int n_cols = nn / 2;
1609 const int mid = mm / 2;
1610 constexpr
int mid_static = mm_static / 2;
1611 constexpr
int max_mid = 15;
1613 Assert((n_rows_static != -1 && n_columns_static != -1) || mid <= max_mid,
1617 const int n_blocks1 = one_line ? 1 : stride;
1618 const int n_blocks2 =
1619 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1621 const int offset = (n_columns + 1) / 2;
1627 for (
int i2 = 0; i2 < n_blocks2; ++i2)
1629 for (
int i1 = 0; i1 < n_blocks1; ++i1)
1631 constexpr
unsigned int mid_size =
1632 (n_rows_static == -1 || n_columns_static == -1) ?
1634 (mid_static > 0 ? mid_static : 1);
1635 Number xp[mid_size], xm[mid_size];
1636 for (
int i = 0; i < mid; ++i)
1638 if (contract_over_rows ==
true && type == 1)
1640 xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
1641 xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
1645 xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
1646 xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
1649 Number xmid = in[stride * mid];
1650 for (
int col = 0; col < n_cols; ++col)
1655 if (contract_over_rows ==
true)
1657 r0 = shapes[col] * xp[0];
1658 r1 = shapes[(n_rows - 1) * offset + col] * xm[0];
1662 r0 = shapes[col * offset] * xp[0];
1663 r1 = shapes[(n_rows - 1 - col) * offset] * xm[0];
1665 for (
int ind = 1; ind < mid; ++ind)
1667 if (contract_over_rows ==
true)
1669 r0 += shapes[ind * offset + col] * xp[ind];
1670 r1 += shapes[(n_rows - 1 - ind) * offset + col] *
1675 r0 += shapes[col * offset + ind] * xp[ind];
1676 r1 += shapes[(n_rows - 1 - col) * offset + ind] *
1683 if (mm % 2 == 1 && contract_over_rows ==
true)
1686 r1 += shapes[mid * offset + col] * xmid;
1688 r0 += shapes[mid * offset + col] * xmid;
1690 else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0 || mm == 3))
1691 r0 += shapes[col * offset + mid] * xmid;
1695 out[stride * col] += r0 + r1;
1696 if (type == 1 && contract_over_rows ==
false)
1697 out[stride * (nn - 1 - col)] += r1 - r0;
1699 out[stride * (nn - 1 - col)] += r0 - r1;
1703 out[stride * col] = r0 + r1;
1704 if (type == 1 && contract_over_rows ==
false)
1705 out[stride * (nn - 1 - col)] = r1 - r0;
1707 out[stride * (nn - 1 - col)] = r0 - r1;
1710 if (type == 0 && contract_over_rows ==
true && nn % 2 == 1 &&
1711 mm % 2 == 1 && mm > 3)
1714 out[stride * n_cols] += shapes[mid * offset + n_cols] * xmid;
1716 out[stride * n_cols] = shapes[mid * offset + n_cols] * xmid;
1718 else if (contract_over_rows ==
true && nn % 2 == 1)
1723 r0 = shapes[n_cols] * xp[0];
1724 for (
int ind = 1; ind < mid; ++ind)
1725 r0 += shapes[ind * offset + n_cols] * xp[ind];
1729 if (type != 1 && mm % 2 == 1)
1730 r0 += shapes[mid * offset + n_cols] * xmid;
1733 out[stride * n_cols] += r0;
1735 out[stride * n_cols] = r0;
1737 else if (contract_over_rows ==
false && nn % 2 == 1)
1744 r0 = shapes[n_cols * offset] * xm[0];
1745 for (
int ind = 1; ind < mid; ++ind)
1746 r0 += shapes[n_cols * offset + ind] * xm[ind];
1750 r0 = shapes[n_cols * offset] * xp[0];
1751 for (
int ind = 1; ind < mid; ++ind)
1752 r0 += shapes[n_cols * offset + ind] * xp[ind];
1758 if ((type == 0 || type == 2) && mm % 2 == 1)
1759 r0 += shapes[n_cols * offset + mid] * xmid;
1762 out[stride * n_cols] += r0;
1764 out[stride * n_cols] = r0;
1766 if (one_line ==
false)
1772 if (one_line ==
false)
1774 in += stride * (mm - 1);
1775 out += stride * (nn - 1);
1825 static constexpr
unsigned int n_rows_of_product =
1827 static constexpr
unsigned int n_columns_of_product =
1836 : shape_values(nullptr)
1837 , shape_gradients(nullptr)
1838 , shape_hessians(nullptr)
1846 : shape_values(shape_values.
begin())
1847 , shape_gradients(nullptr)
1848 , shape_hessians(nullptr)
1860 const unsigned int dummy1 = 0,
1861 const unsigned int dummy2 = 0)
1862 : shape_values(shape_values.
begin())
1863 , shape_gradients(shape_gradients.
begin())
1864 , shape_hessians(shape_hessians.
begin())
1868 if (!shape_values.
empty())
1870 if (!shape_gradients.
empty())
1872 if (!shape_hessians.
empty())
1878 template <
int direction,
bool contract_over_rows,
bool add>
1880 values(
const Number in[], Number out[])
const
1883 apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
1886 template <
int direction,
bool contract_over_rows,
bool add>
1891 apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
1894 template <
int direction,
bool contract_over_rows,
bool add>
1899 apply<direction, contract_over_rows, add, 2>(shape_hessians, in, out);
1902 template <
int direction,
bool contract_over_rows,
bool add>
1907 apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
1910 template <
int direction,
bool contract_over_rows,
bool add>
1915 apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
1920 template <
int direction,
bool contract_over_rows,
bool add>
1925 apply<direction, contract_over_rows, add, 2, true>(shape_hessians,
1958 template <
int direction,
1959 bool contract_over_rows,
1962 bool one_line =
false>
1977 one_line>(n_rows, n_columns, shape_data, in, out);
1993 template <
int dim,
typename Number,
typename Number2>
1997 : shape_values(nullptr)
1998 , shape_gradients(nullptr)
1999 , shape_hessians(nullptr)
2005 const unsigned int n_rows = 0,
2006 const unsigned int n_columns = 0)
2007 : shape_values(shape_values.
begin())
2008 , shape_gradients(nullptr)
2009 , shape_hessians(nullptr)
2011 , n_columns(n_columns)
2019 const unsigned int n_rows = 0,
2020 const unsigned int n_columns = 0)
2021 : shape_values(shape_values.
begin())
2022 , shape_gradients(shape_gradients.
begin())
2023 , shape_hessians(shape_hessians.
begin())
2025 , n_columns(n_columns)
2027 if (!shape_values.
empty())
2029 if (!shape_gradients.
empty())
2031 if (!shape_hessians.
empty())
2035 template <
int direction,
bool contract_over_rows,
bool add>
2037 values(
const Number in[], Number out[])
const
2040 apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
2043 template <
int direction,
bool contract_over_rows,
bool add>
2048 apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
2051 template <
int direction,
bool contract_over_rows,
bool add>
2056 apply<direction, contract_over_rows, add, 2>(shape_hessians, in, out);
2059 template <
int direction,
bool contract_over_rows,
bool add>
2064 apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
2067 template <
int direction,
bool contract_over_rows,
bool add>
2072 apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
2077 template <
int direction,
bool contract_over_rows,
bool add>
2082 apply<direction, contract_over_rows, add, 2, true>(shape_hessians,
2087 template <
int direction,
2088 bool contract_over_rows,
2091 bool one_line =
false>
2106 one_line>(n_rows, n_columns, shape_data, in, out);
2159 static constexpr
unsigned int n_rows_of_product =
2161 static constexpr
unsigned int n_columns_of_product =
2170 : shape_values(nullptr)
2171 , shape_gradients(nullptr)
2172 , shape_hessians(nullptr)
2180 : shape_values(shape_values.
begin())
2181 , shape_gradients(nullptr)
2182 , shape_hessians(nullptr)
2192 const unsigned int dummy1 = 0,
2193 const unsigned int dummy2 = 0)
2194 : shape_values(shape_values.
begin())
2195 , shape_gradients(shape_gradients.
begin())
2196 , shape_hessians(shape_hessians.
begin())
2202 template <
int direction,
bool contract_over_rows,
bool add>
2204 values(
const Number in[], Number out[])
const
2207 apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
2210 template <
int direction,
bool contract_over_rows,
bool add>
2215 apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
2218 template <
int direction,
bool contract_over_rows,
bool add>
2223 apply<direction, contract_over_rows, add, 0>(shape_hessians, in, out);
2226 template <
int direction,
bool contract_over_rows,
bool add>
2231 apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
2234 template <
int direction,
bool contract_over_rows,
bool add>
2239 apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
2244 template <
int direction,
bool contract_over_rows,
bool add>
2249 apply<direction, contract_over_rows, add, 0, true>(shape_hessians,
2281 template <
int direction,
2282 bool contract_over_rows,
2285 bool one_line =
false>
2304 template <
int direction,
2305 bool contract_over_rows,
2319 static_assert(one_line ==
false || direction == dim - 1,
2320 "Single-line evaluation only works for direction=dim-1.");
2322 type == 0 || type == 1,
2323 "Only types 0 and 1 implemented for evaluate_symmetric_hierarchical.");
2324 Assert(dim == direction + 1 || one_line ==
true || n_rows == n_columns ||
2326 ExcMessage(
"In-place operation only supported for "
2327 "n_rows==n_columns or single-line interpolation"));
2333 constexpr
int nn = contract_over_rows ? n_columns : n_rows;
2334 constexpr
int mm = contract_over_rows ? n_rows : n_columns;
2335 constexpr
int n_cols = nn / 2;
2336 constexpr
int mid = mm / 2;
2339 constexpr
int n_blocks1 = one_line ? 1 : stride;
2340 constexpr
int n_blocks2 =
2341 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
2347 for (
int i2 = 0; i2 < n_blocks2; ++i2)
2349 for (
int i1 = 0; i1 < n_blocks1; ++i1)
2351 if (contract_over_rows)
2354 for (
unsigned int i = 0; i < mm; ++i)
2355 x[i] = in[stride * i];
2356 for (
unsigned int col = 0; col < n_cols; ++col)
2361 r0 = shapes[col] * x[0];
2362 r1 = shapes[col + n_columns] * x[1];
2363 for (
unsigned int ind = 1; ind < mid; ++ind)
2366 shapes[col + 2 * ind * n_columns] * x[2 * ind];
2367 r1 += shapes[col + (2 * ind + 1) * n_columns] *
2374 r0 += shapes[col + (mm - 1) * n_columns] * x[mm - 1];
2377 out[stride * col] += r0 + r1;
2379 out[stride * (nn - 1 - col)] += r1 - r0;
2381 out[stride * (nn - 1 - col)] += r0 - r1;
2385 out[stride * col] = r0 + r1;
2387 out[stride * (nn - 1 - col)] = r1 - r0;
2389 out[stride * (nn - 1 - col)] = r0 - r1;
2395 const unsigned int shift = type == 1 ? 1 : 0;
2398 r0 = shapes[n_cols +
shift * n_columns] * x[
shift];
2399 for (
unsigned int ind = 1; ind < mid; ++ind)
2400 r0 += shapes[n_cols + (2 * ind +
shift) * n_columns] *
2405 if (type != 1 && mm % 2 == 1)
2406 r0 += shapes[n_cols + (mm - 1) * n_columns] * x[mm - 1];
2408 out[stride * n_cols] += r0;
2410 out[stride * n_cols] = r0;
2415 Number xp[mid + 1], xm[mid > 0 ? mid : 1];
2416 for (
int i = 0; i < mid; ++i)
2419 xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
2420 xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
2424 xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
2425 xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
2428 xp[mid] = in[stride * mid];
2429 for (
unsigned int col = 0; col < n_cols; ++col)
2434 r0 = shapes[2 * col * n_columns] * xp[0];
2435 r1 = shapes[(2 * col + 1) * n_columns] * xm[0];
2436 for (
unsigned int ind = 1; ind < mid; ++ind)
2438 r0 += shapes[2 * col * n_columns + ind] * xp[ind];
2440 shapes[(2 * col + 1) * n_columns + ind] * xm[ind];
2449 shapes[(2 * col + 1) * n_columns + mid] * xp[mid];
2451 r0 += shapes[2 * col * n_columns + mid] * xp[mid];
2455 out[stride * (2 * col)] += r0;
2456 out[stride * (2 * col + 1)] += r1;
2460 out[stride * (2 * col)] = r0;
2461 out[stride * (2 * col + 1)] = r1;
2469 r0 = shapes[(nn - 1) * n_columns] * xp[0];
2470 for (
unsigned int ind = 1; ind < mid; ++ind)
2471 r0 += shapes[(nn - 1) * n_columns + ind] * xp[ind];
2475 if (mm % 2 == 1 && type == 0)
2476 r0 += shapes[(nn - 1) * n_columns + mid] * xp[mid];
2478 out[stride * (nn - 1)] += r0;
2480 out[stride * (nn - 1)] = r0;
2483 if (one_line ==
false)
2489 if (one_line ==
false)
2491 in += stride * (mm - 1);
2492 out += stride * (nn - 1);
2532 static constexpr
unsigned int n_rows_of_product =
2534 static constexpr
unsigned int n_columns_of_product =
2542 : shape_values(nullptr)
2543 , shape_gradients(nullptr)
2544 , shape_hessians(nullptr)
2554 const unsigned int dummy1 = 0,
2555 const unsigned int dummy2 = 0)
2556 : shape_values(shape_values.
begin())
2557 , shape_gradients(shape_gradients.
begin())
2558 , shape_hessians(shape_hessians.
begin())
2565 shape_values.
size() == n_rows * n_columns ||
2566 shape_values.
size() == 3 * n_rows,
2569 shape_gradients.
size() == n_rows * n_columns,
2572 shape_hessians.
size() == n_rows * n_columns,
2578 template <
int direction,
bool contract_over_rows,
bool add>
2580 values(
const Number in[], Number out[])
const
2582 apply<direction, contract_over_rows, add>(shape_values, in, out);
2585 template <
int direction,
bool contract_over_rows,
bool add>
2589 apply<direction, contract_over_rows, add>(shape_gradients, in, out);
2592 template <
int direction,
bool contract_over_rows,
bool add>
2596 apply<direction, contract_over_rows, add>(shape_hessians, in, out);
2627 template <
int direction,
2628 bool contract_over_rows,
2630 bool one_line =
false>
2636 template <
int face_direction,
2637 bool contract_onto_face,
2656 template <
int direction,
bool contract_over_rows,
bool add,
bool one_line>
2669 static_assert(one_line ==
false || direction == dim - 1,
2670 "Single-line evaluation only works for direction=dim-1.");
2671 Assert(shape_data !=
nullptr,
2673 "The given array shape_data must not be the null pointer!"));
2674 Assert(dim == direction + 1 || one_line ==
true || n_rows == n_columns ||
2676 ExcMessage(
"In-place operation only supported for "
2677 "n_rows==n_columns or single-line interpolation"));
2679 constexpr
int mm = contract_over_rows ? n_rows : n_columns,
2680 nn = contract_over_rows ? n_columns : n_rows;
2683 constexpr
int n_blocks1 = one_line ? 1 : stride;
2686 constexpr
int n_blocks2 =
2687 (dim - direction - 1 == 0) ?
2689 ((direction == normal_dir) ?
2691 (direction >= dim) ? 0 : dim - direction - 1) :
2692 (((direction < normal_dir) ? (n_rows + 1) : n_rows) *
2693 ((dim - direction == 3) ? n_rows : 1)));
2695 for (
int i2 = 0; i2 < n_blocks2; ++i2)
2697 for (
int i1 = 0; i1 < n_blocks1; ++i1)
2700 for (
int i = 0; i < mm; ++i)
2701 x[i] = in[stride * i];
2703 for (
int col = 0; col < nn; ++col)
2707 if (contract_over_rows)
2708 val0 = shape_data[col];
2710 val0 = shape_data[col * n_columns];
2712 Number res0 = val0 * x[0];
2713 for (
int i = 1; i < mm; ++i)
2715 if (contract_over_rows)
2716 val0 = shape_data[i * n_columns + col];
2718 val0 = shape_data[col * n_columns + i];
2720 res0 += val0 * x[i];
2723 out[stride * col] += res0;
2726 out[stride * col] = res0;
2729 if (one_line ==
false)
2735 if (one_line ==
false)
2737 in += stride * (mm - 1);
2738 out += stride * (nn - 1);
2749 template <
int face_direction,
2750 bool contract_onto_face,
2765 static_assert(max_derivative >= 0 && max_derivative < 3,
2766 "Only derivative orders 0-2 implemented");
2767 Assert(shape_values !=
nullptr,
2769 "The given array shape_values must not be the null pointer."));
2773 constexpr
int n_blocks1 = (face_direction == normal_dir) ? (n_rows - 1) :
2774 ((face_direction == 0 && normal_dir == 2) ||
2775 (face_direction == 1 && normal_dir == 2) ||
2776 (face_direction == 2 && normal_dir == 1)) ?
2779 constexpr
int n_blocks2 = (dim == 2) ?
2781 ((face_direction == normal_dir) ?
2783 (((face_direction == 0 && normal_dir == 1) ||
2784 (face_direction == 1 && normal_dir == 0) ||
2785 (face_direction == 2 && normal_dir == 0)) ?
2791 constexpr
int in_stride =
2792 (face_direction == normal_dir) ?
2794 ((face_direction == 0) ?
2796 ((face_direction == 2) ?
2797 n_rows * (n_rows + 1) :
2798 ((face_direction == 1 && normal_dir == 0) ? (n_rows + 1) :
2800 constexpr
int out_stride = n_blocks1 * n_blocks2;
2804 for (
int i2 = 0; i2 < n_blocks2; ++i2)
2806 for (
int i1 = 0; i1 < n_blocks1; ++i1)
2808 if (contract_onto_face ==
true)
2810 Number res0 = shape_values[0] * in[0];
2813 if (max_derivative > 0)
2814 res1 = shape_values[n_rows] * in[0];
2816 if (max_derivative > 1)
2817 res2 = shape_values[2 * n_rows] * in[0];
2819 for (
int ind = 1; ind < n_rows; ++ind)
2821 res0 += shape_values[ind] * in[in_stride * ind];
2822 if (max_derivative > 0)
2823 res1 += shape_values[ind + n_rows] * in[in_stride * ind];
2825 if (max_derivative > 1)
2827 shape_values[ind + 2 * n_rows] * in[in_stride * ind];
2833 if (max_derivative > 0)
2834 out[out_stride] += res1;
2836 if (max_derivative > 1)
2837 out[2 * out_stride] += res2;
2843 if (max_derivative > 0)
2844 out[out_stride] = res1;
2846 if (max_derivative > 1)
2847 out[2 * out_stride] = res2;
2852 for (
int col = 0; col < n_rows; ++col)
2855 out[col * in_stride] += shape_values[col] * in[0];
2857 out[col * in_stride] = shape_values[col] * in[0];
2859 if (max_derivative > 0)
2860 out[col * in_stride] +=
2861 shape_values[col + n_rows] * in[out_stride];
2863 if (max_derivative > 1)
2864 out[col * in_stride] +=
2865 shape_values[col + 2 * n_rows] * in[2 * out_stride];
2872 switch (face_direction)
2875 in += contract_onto_face ? n_rows : 1;
2876 out += contract_onto_face ? 1 : n_rows;
2887 if (normal_dir == 0)
2889 if (contract_onto_face)
2894 if (normal_dir == 1)
2896 if (contract_onto_face)
2901 if (normal_dir == 2)
2903 if (contract_onto_face)
2920 if (face_direction == 1 && dim == 3)
2923 if (contract_onto_face)
2925 if (normal_dir == 0)
2927 in += (n_rows + 1) * (n_rows - 1);
2928 out -= n_rows * (n_rows + 1) - 1;
2930 if (normal_dir == 1)
2932 in += (n_rows - 1) * (n_rows - 1);
2933 out -= (n_rows - 1) * (n_rows - 1) - 1;
2935 if (normal_dir == 2)
2937 in += (n_rows - 1) * (n_rows);
2938 out -= (n_rows) * (n_rows + 1) - 1;
2943 if (normal_dir == 0)
2945 out += (n_rows + 1) * (n_rows - 1);
2946 in -= n_rows * (n_rows + 1) - 1;
2948 if (normal_dir == 1)
2950 out += (n_rows - 1) * (n_rows - 1);
2951 in -= (n_rows - 1) * (n_rows - 1) - 1;
2953 if (normal_dir == 2)
2955 out += (n_rows - 1) * (n_rows);
2956 in -= (n_rows) * (n_rows + 1) - 1;
2971 template <
typename Number,
typename Number2>
2977 template <
int dim,
typename Number,
typename Number2>
3019 template <
int dim,
typename Number,
typename Number2>
3025 const std::vector<Number> &
values,
3027 const bool d_linear =
false,
3028 const std::vector<unsigned int> & renumber = {})
3030 static_assert(dim >= 1 && dim <= 3,
"Only dim=1,2,3 implemented");
3038 const int n_shapes = poly.size();
3040 Assert(renumber.empty() || renumber.size() ==
values.size(),
3047 for (
unsigned int i = 0; i < renumber.size(); ++i)
3054 return std::make_pair((1. - p[0]) *
values[0] + p[0] *
values[1],
3059 const Number2 x0 = 1. - p[0], x1 = p[0];
3062 const Number3 mapped = (1. - p[1]) * tmp0 + p[1] * tmp1;
3064 derivative[0] = (1. - p[1]) * (
values[1] -
values[0]) +
3066 derivative[1] = tmp1 - tmp0;
3067 return std::make_pair(mapped, derivative);
3071 const Number2 x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
3072 z0 = 1. - p[2], z1 = p[2];
3075 const Number3 tmpy0 = y0 * tmp0 + y1 * tmp1;
3078 const Number3 tmpy1 = y0 * tmp2 + y1 * tmp3;
3079 const Number3 mapped = z0 * tmpy0 + z1 * tmpy1;
3081 derivative[2] = tmpy1 - tmpy0;
3082 derivative[1] = z0 * (tmp1 - tmp0) + z1 * (tmp3 - tmp2);
3088 return std::make_pair(mapped, derivative);
3096 std::array<Number2, dim>
point;
3097 for (
unsigned int d = 0;
d < dim; ++
d)
3099 for (
int i = 0; i < n_shapes; ++i)
3100 poly[i].values_of_array(
point, 1, &shapes[i][0]);
3104 std::pair<Number3, Tensor<1, dim, Number3>> result = {};
3105 for (
int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
3107 Number3 value_y = {}, deriv_x = {}, deriv_y = {};
3108 for (
int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
3111 Number3
value = {}, deriv = {};
3115 if (renumber.empty())
3116 for (
int i0 = 0; i0 < n_shapes; ++i0, ++i)
3119 deriv += shapes[i0][1][0] *
values[i];
3122 for (
int i0 = 0; i0 < n_shapes; ++i0, ++i)
3125 deriv += shapes[i0][1][0] *
values[renumber[i]];
3131 value_y +=
value * shapes[i1][0][1];
3132 deriv_x += deriv * shapes[i1][0][1];
3133 deriv_y +=
value * shapes[i1][1][1];
3137 result.first =
value;
3138 result.second[0] = deriv;
3144 result.first += value_y * shapes[i2][0][2];
3145 result.second[0] += deriv_x * shapes[i2][0][2];
3146 result.second[1] += deriv_y * shapes[i2][0][2];
3147 result.second[2] += value_y * shapes[i2][1][2];
3151 result.first = value_y;
3152 result.second[0] = deriv_x;
3153 result.second[1] = deriv_y;
3162 template <
int dim,
typename Number,
typename Number2>
3166 const std::vector<Number> &
values,
3168 const std::vector<unsigned int> & renumber = {})
3170 static_assert(dim >= 1 && dim <= 3,
"Only dim=1,2,3 implemented");
3178 const int n_shapes = poly.size();
3180 Assert(renumber.empty() || renumber.size() ==
values.size(),
3187 std::array<Number2, dim>
point;
3188 for (
unsigned int d = 0;
d < dim; ++
d)
3190 for (
int i = 0; i < n_shapes; ++i)
3191 poly[i].values_of_array(
point, 2, &shapes[i][0]);
3196 for (
int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
3198 Number3 value_y = {}, deriv_x = {}, deriv_y = {}, deriv_xx = {},
3199 deriv_xy = {}, deriv_yy = {};
3200 for (
int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
3203 Number3
value = {}, deriv_1 = {}, deriv_2 = {};
3207 if (renumber.empty())
3208 for (
int i0 = 0; i0 < n_shapes; ++i0, ++i)
3211 deriv_1 += shapes[i0][1][0] *
values[i];
3212 deriv_2 += shapes[i0][2][0] *
values[i];
3215 for (
int i0 = 0; i0 < n_shapes; ++i0, ++i)
3218 deriv_1 += shapes[i0][1][0] *
values[renumber[i]];
3219 deriv_2 += shapes[i0][2][0] *
values[renumber[i]];
3227 value_y +=
value * shapes[i1][0][1];
3228 deriv_x += deriv_1 * shapes[i1][0][1];
3229 deriv_y +=
value * shapes[i1][1][1];
3231 deriv_xx += deriv_2 * shapes[i1][0][1];
3232 deriv_xy += deriv_1 * shapes[i1][1][1];
3233 deriv_yy +=
value * shapes[i1][2][1];
3237 result[0][0] = deriv_2;
3243 result[0][0] += deriv_xx * shapes[i2][0][2];
3244 result[0][1] += deriv_xy * shapes[i2][0][2];
3245 result[0][2] += deriv_x * shapes[i2][1][2];
3246 result[1][1] += deriv_yy * shapes[i2][0][2];
3247 result[1][2] += deriv_y * shapes[i2][1][2];
3248 result[2][2] += value_y * shapes[i2][2][2];
3252 result[0][0] = deriv_xx;
3253 result[1][0] = deriv_xy;
3254 result[1][1] = deriv_yy;
3266 template <
int dim,
typename Number,
typename Number2>
3270 const Number2 & value,
3274 const std::vector<unsigned int> & renumber = {})
3276 static_assert(dim >= 1 && dim <= 3,
"Only dim=1,2,3 implemented");
3279 const int n_shapes = poly.size();
3281 Assert(renumber.empty() || renumber.size() ==
values.size(),
3288 std::array<Number, dim>
point;
3289 for (
unsigned int d = 0;
d < dim; ++
d)
3291 for (
int i = 0; i < n_shapes; ++i)
3292 poly[i].values_of_array(
point, 1, &shapes[i][0]);
3295 for (
int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
3297 const Number2 test_value_z =
3299 (value * shapes[i2][0][2] + gradient[2] * shapes[i2][1][2]) :
3301 const Number2 test_grad_x =
3302 dim > 2 ? gradient[0] * shapes[i2][0][2] : gradient[0];
3303 const Number2 test_grad_y = dim > 2 ?
3304 gradient[1] * shapes[i2][0][2] :
3305 (dim > 1 ? gradient[1] : Number2());
3306 for (
int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
3308 const Number2 test_value_y = dim > 1 ?
3309 (test_value_z * shapes[i1][0][1] +
3310 test_grad_y * shapes[i1][1][1]) :
3312 const Number2 test_grad_xy =
3313 dim > 1 ? test_grad_x * shapes[i1][0][1] : test_grad_x;
3314 if (renumber.empty())
3315 for (
int i0 = 0; i0 < n_shapes; ++i0, ++i)
3316 values[i] += shapes[i0][0][0] * test_value_y +
3317 shapes[i0][1][0] * test_grad_xy;
3319 for (
int i0 = 0; i0 < n_shapes; ++i0, ++i)
3320 values[renumber[i]] += shapes[i0][0][0] * test_value_y +
3321 shapes[i0][1][0] * test_grad_xy;
3327 template <
int dim,
int n_po
ints_1d_
template,
typename Number>
3330 const unsigned int n_components,
3331 const int n_points_1d_non_template,
3334 const int n_points_1d = n_points_1d_template != -1 ?
3335 n_points_1d_template :
3336 n_points_1d_non_template;
3341 unsigned int compressed_index[100];
3342 compressed_index[0] = 0;
3343 for (
int i = 1; i < n_points_1d - 1; ++i)
3344 compressed_index[i] = 1;
3345 compressed_index[n_points_1d - 1] = 2;
3347 for (
unsigned int c = 0; c < n_components; ++c)
3348 for (
int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3349 for (
int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
3351 const unsigned int shift =
3352 9 * compressed_index[k] + 3 * compressed_index[j];
3353 data[0] *= weights[
shift];
3356 for (
int i = 1; i < n_points_1d - 1; ++i)
3357 data[i] *= weights[
shift + 1];
3358 data[n_points_1d - 1] *= weights[
shift + 2];
3359 data += n_points_1d;
3364 template <
int dim,
int n_po
ints_1d_
template,
typename Number>
3367 const unsigned int n_components,
3368 const int n_points_1d_non_template,
3371 const int n_points_1d = n_points_1d_template != -1 ?
3372 n_points_1d_template :
3373 n_points_1d_non_template;
3375 Assert((n_points_1d % 2) == 1,
3376 ExcMessage(
"The function can only with add number of points"));
3380 const unsigned int n_inside_1d = n_points_1d / 2;
3382 unsigned int compressed_index[100];
3385 for (
int i = 0; i < n_inside_1d; ++i)
3386 compressed_index[c++] = 0;
3387 compressed_index[c++] = 1;
3388 for (
int i = 0; i < n_inside_1d; ++i)
3389 compressed_index[c++] = 2;
3391 for (
unsigned int c = 0; c < n_components; ++c)
3392 for (
int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3393 for (
int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
3395 const unsigned int shift =
3396 9 * compressed_index[k] + 3 * compressed_index[j];
3399 for (
int i = 0; i < n_inside_1d; ++i)
3400 data[c++] *= weights[
shift];
3401 data[c++] *= weights[
shift + 1];
3402 for (
int i = 0; i < n_inside_1d; ++i)
3403 data[c++] *= weights[
shift + 2];
3404 data += n_points_1d;
3409 template <
int dim,
int n_po
ints_1d_
template,
typename Number>
3412 const unsigned int n_components,
3413 const int n_points_1d_non_template,
3416 const int n_points_1d = n_points_1d_template != -1 ?
3417 n_points_1d_template :
3418 n_points_1d_non_template;
3423 unsigned int compressed_index[100];
3424 compressed_index[0] = 0;
3425 for (
int i = 1; i < n_points_1d - 1; ++i)
3426 compressed_index[i] = 1;
3427 compressed_index[n_points_1d - 1] = 2;
3434 const auto check_and_set = [](Number &weight,
const Number &data) {
3435 if (weight == Number(-1.0) || weight == data)
3444 for (
unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
3445 weights[c] = Number(-1.0);
3447 for (
unsigned int c = 0; c < n_components; ++c)
3448 for (
int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3449 for (
int j = 0; j < (dim > 1 ? n_points_1d : 1);
3450 ++j, data += n_points_1d)
3452 const unsigned int shift =
3453 9 * compressed_index[k] + 3 * compressed_index[j];
3455 if (!check_and_set(weights[
shift], data[0]))
3458 for (
int i = 1; i < n_points_1d - 1; ++i)
3459 if (!check_and_set(weights[
shift + 1], data[i]))
3462 if (!check_and_set(weights[
shift + 2], data[n_points_1d - 1]))
3470 template <
int dim,
int n_po
ints_1d_
template,
typename Number>
3473 const Number * data,
3474 const unsigned int n_components,
3475 const int n_points_1d_non_template,
3478 const int n_points_1d = n_points_1d_template != -1 ?
3479 n_points_1d_template :
3480 n_points_1d_non_template;
3482 Assert((n_points_1d % 2) == 1,
3483 ExcMessage(
"The function can only with add number of points"));
3487 const unsigned int n_inside_1d = n_points_1d / 2;
3489 unsigned int compressed_index[100];
3492 for (
int i = 0; i < n_inside_1d; ++i)
3493 compressed_index[c++] = 0;
3494 compressed_index[c++] = 1;
3495 for (
int i = 0; i < n_inside_1d; ++i)
3496 compressed_index[c++] = 2;
3503 const auto check_and_set = [](Number &weight,
const Number &data) {
3504 if (weight == Number(-1.0) || weight == data)
3513 for (
unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
3514 weights[c] = Number(-1.0);
3516 for (
unsigned int comp = 0; comp < n_components; ++comp)
3517 for (
int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3518 for (
int j = 0; j < (dim > 1 ? n_points_1d : 1);
3519 ++j, data += n_points_1d)
3521 const unsigned int shift =
3522 9 * compressed_index[k] + 3 * compressed_index[j];
3526 for (
int i = 0; i < n_inside_1d; ++i)
3527 if (!check_and_set(weights[
shift], data[c++]))
3530 if (!check_and_set(weights[
shift + 1], data[c++]))
3533 for (
int i = 0; i < n_inside_1d; ++i)
3534 if (!check_and_set(weights[
shift + 2], data[c++]))
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
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)
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)
VectorType::value_type * begin(VectorType &V)
constexpr T pow(const T base, const int iexp)
void integrate_add_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double >> &poly, const Number2 &value, const Tensor< 1, dim, Number2 > &gradient, const Point< dim, Number > &p, AlignedVector< Number2 > &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)
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={})
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 even_odd_apply(const int n_rows_in, const int n_columns_in, const Number2 *DEAL_II_RESTRICT shapes, const Number *in, Number *out)
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)
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)
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)
@ evaluate_raviart_thomas
@ evaluate_symmetric_hierarchical
static const unsigned int invalid_unsigned_int
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
void gradients(const Number in[], Number out[]) const
const Number2 * shape_hessians
void hessians(const Number in[], Number out[]) const
EvaluatorTensorProductAnisotropic()
void values(const Number in[], Number out[]) const
const Number2 * shape_gradients
EvaluatorTensorProductAnisotropic(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
const Number2 * shape_values
void hessians(const Number in[], Number out[]) const
void gradients_one_line(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 unsigned int n_columns
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const unsigned int n_rows=0, const unsigned int n_columns=0)
void gradients(const Number in[], Number out[]) const
void hessians_one_line(const Number in[], Number out[]) const
const Number2 * shape_values
void values(const Number in[], Number out[]) const
const Number2 * shape_gradients
const unsigned int n_rows
void values_one_line(const Number in[], Number out[]) const
void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out) const
const Number2 * shape_hessians
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
void hessians_one_line(const Number in[], Number out[]) const
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values)
void gradients_one_line(const Number in[], Number out[]) const
const Number2 * shape_values
const Number2 * shape_hessians
const Number2 * shape_gradients
void values_one_line(const Number in[], Number out[]) const
static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out)
void values(const Number in[], Number out[]) const
void gradients(const Number in[], Number out[]) const
void hessians(const Number in[], Number out[]) const
EvaluatorTensorProduct(const Number2 *shape_values, const Number2 *shape_gradients, const Number2 *shape_hessians, const unsigned int n_rows, const unsigned int n_columns)
const unsigned int n_columns
const Number2 * shape_gradients
void hessians_one_line(const Number in[], Number out[]) const
void gradients(const Number *in, Number *out) const
void values_one_line(const Number in[], Number out[]) const
void gradients_one_line(const Number in[], Number out[]) const
const Number2 * shape_values
void values(const Number *in, Number *out) const
const Number2 * shape_hessians
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int n_rows, const unsigned int n_columns)
const unsigned int n_rows
void hessians(const Number *in, Number *out) const
void gradients(const Number in[], Number out[]) const
const Number2 * shape_values
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
const Number2 * shape_gradients
void values(const Number in[], Number out[]) const
void hessians(const Number in[], Number out[]) const
void hessians_one_line(const Number in[], Number out[]) const
void values_one_line(const Number in[], Number out[]) const
const Number2 * shape_hessians
void gradients_one_line(const Number in[], Number out[]) const
const Number2 * shape_gradients
const Number2 * shape_hessians
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
const Number2 * shape_values
const Number2 * shape_hessians
void gradients(const Number in[], Number out[]) const
void values(const Number in[], Number out[]) const
void hessians_one_line(const Number in[], Number out[]) const
EvaluatorTensorProduct(const AlignedVector< Number > &shape_values)
void values_one_line(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 dummy1=0, const unsigned int dummy2=0)
void gradients_one_line(const Number in[], Number out[]) const
void hessians(const Number in[], Number out[]) const
const Number2 * shape_gradients
const Number2 * shape_values
typename ProductType< Tensor< 1, dim, Number >, Number2 >::type type
typename ProductType< Number, Number2 >::type type