17 #ifndef dealii_matrix_free_evaluation_kernels_h
18 #define dealii_matrix_free_evaluation_kernels_h
40 template <MatrixFreeFunctions::ElementType element,
bool is_
long>
44 template <
bool is_
long>
62 template <
bool is_
long>
81 template <
bool is_
long>
88 template <
bool is_
long>
129 evaluate(
const unsigned int n_components,
131 const Number * values_dofs_actual,
135 integrate(
const unsigned int n_components,
137 Number * values_dofs_actual,
139 const bool add_into_values_array);
144 *univariate_shape_data)
167 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
175 evaluate(
const unsigned int n_components,
177 const Number * values_dofs_actual,
181 integrate(
const unsigned int n_components,
183 Number * values_dofs_actual,
185 const bool add_into_values_array);
192 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
199 template <
bool integrate>
201 evaluate_or_integrate(
203 Number * values_dofs_actual,
205 const bool add_into_values_array =
false);
208 template <
typename EvalType>
218 template <
int normal_dir>
220 evaluate_tensor_product_per_component(
222 Number * values_dofs_actual,
224 const bool add_into_values_array,
225 std::integral_constant<bool, false>);
227 template <
int normal_dir>
229 evaluate_tensor_product_per_component(
231 Number * values_dofs_actual,
233 const bool add_into_values_array,
234 std::integral_constant<bool, true>);
246 const unsigned int n_components,
248 const Number * values_dofs_actual,
254 std::array<const MatrixFreeFunctions::UnivariateShapeData<Number> *, 3>
255 univariate_shape_data;
259 univariate_shape_data.fill(&shape_data.front());
261 if (shape_data.size() == dim)
262 for (
int i = 1; i < dim; ++i)
263 univariate_shape_data[i] = &shape_data[i];
265 Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
266 Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
267 Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);
269 const unsigned int temp_size =
272 (Eval::n_rows_of_product > Eval::n_columns_of_product ?
273 Eval::n_rows_of_product :
274 Eval::n_columns_of_product);
279 temp2 = temp1 +
std::max(Utilities::fixed_power<dim>(
280 shape_data.front().fe_degree + 1),
281 Utilities::fixed_power<dim>(
282 shape_data.front().n_q_points_1d));
286 temp2 = temp1 + temp_size;
289 const std::size_t n_q_points = temp_size == 0 ?
291 Eval::n_columns_of_product;
292 const std::size_t dofs_per_comp =
296 const Number *values_dofs = values_dofs_actual;
299 const std::size_t n_dofs_per_comp =
301 Number *values_dofs_tmp =
302 temp1 + 2 * (
std::max(n_dofs_per_comp, n_q_points));
304 fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
305 for (
unsigned int c = 0; c < n_components; ++c)
306 for (
int i = 0, count_p = 0, count_q = 0;
307 i < (dim > 2 ? degree + 1 : 1);
310 for (
int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
312 for (
int k = 0; k < degree + 1 - j - i;
313 ++k, ++count_p, ++count_q)
314 values_dofs_tmp[c * dofs_per_comp + count_q] =
315 values_dofs_actual[c * n_dofs_per_comp + count_p];
316 for (
int k = degree + 1 - j - i; k < degree + 1;
318 values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
320 for (
int j = degree + 1 - i; j < degree + 1; ++j)
321 for (
int k = 0; k < degree + 1; ++k, ++count_q)
322 values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
324 values_dofs = values_dofs_tmp;
334 for (
unsigned int c = 0; c < n_components; ++c)
337 eval0.template values<0, true, false>(values_dofs, values_quad);
339 eval0.template gradients<0, true, false>(values_dofs,
342 eval0.template hessians<0, true, false>(values_dofs,
346 values_dofs += dofs_per_comp;
347 values_quad += n_q_points;
348 gradients_quad += n_q_points;
349 hessians_quad += n_q_points;
354 for (
unsigned int c = 0; c < n_components; ++c)
359 eval0.template gradients<0, true, false>(values_dofs, temp1);
360 eval1.template values<1, true, false>(temp1, gradients_quad);
366 eval0.template gradients<0, true, false>(values_dofs,
368 eval1.template gradients<1, true, false>(temp1,
373 eval0.template hessians<0, true, false>(values_dofs, temp1);
374 eval1.template values<1, true, false>(temp1, hessians_quad);
378 eval0.template values<0, true, false>(values_dofs, temp1);
380 eval1.template gradients<1, true, false>(temp1,
386 eval1.template hessians<1, true, false>(temp1,
392 eval1.template values<1, true, false>(temp1, values_quad);
395 values_dofs += dofs_per_comp;
396 values_quad += n_q_points;
397 gradients_quad += 2 * n_q_points;
398 hessians_quad += 3 * n_q_points;
403 for (
unsigned int c = 0; c < n_components; ++c)
408 eval0.template gradients<0, true, false>(values_dofs, temp1);
409 eval1.template values<1, true, false>(temp1, temp2);
410 eval2.template values<2, true, false>(temp2, gradients_quad);
418 eval0.template gradients<0, true, false>(values_dofs,
420 eval1.template values<1, true, false>(temp1, temp2);
422 eval2.template gradients<2, true, false>(temp2,
427 eval1.template gradients<1, true, false>(temp1, temp2);
428 eval2.template values<2, true, false>(temp2,
433 eval0.template hessians<0, true, false>(values_dofs, temp1);
434 eval1.template values<1, true, false>(temp1, temp2);
435 eval2.template values<2, true, false>(temp2, hessians_quad);
439 eval0.template values<0, true, false>(values_dofs, temp1);
442 eval1.template gradients<1, true, false>(temp1, temp2);
443 eval2.template values<2, true, false>(temp2,
452 eval1.template gradients<1, true, false>(temp1, temp2);
453 eval2.template gradients<2, true, false>(temp2,
458 eval1.template hessians<1, true, false>(temp1, temp2);
459 eval2.template values<2, true, false>(temp2,
466 eval1.template values<1, true, false>(temp1, temp2);
468 eval2.template gradients<2, true, false>(temp2,
475 eval2.template hessians<2, true, false>(temp2,
482 eval2.template values<2, true, false>(temp2, values_quad);
485 values_dofs += dofs_per_comp;
486 values_quad += n_q_points;
487 gradients_quad += 3 * n_q_points;
488 hessians_quad += 6 * n_q_points;
501 values_quad -= n_components * n_q_points;
502 values_dofs -= n_components * dofs_per_comp;
503 for (std::size_t c = 0; c < n_components; ++c)
504 for (std::size_t q = 0; q < n_q_points; ++q)
505 values_quad[c * n_q_points + q] +=
506 values_dofs[(c + 1) * dofs_per_comp - 1];
519 const unsigned int n_components,
521 Number * values_dofs_actual,
523 const bool add_into_values_array)
525 std::array<const MatrixFreeFunctions::UnivariateShapeData<Number> *, 3>
526 univariate_shape_data;
529 univariate_shape_data.fill(&shape_data.front());
531 if (shape_data.size() == dim)
532 for (
int i = 1; i < dim; ++i)
533 univariate_shape_data[i] = &shape_data[i];
535 Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
536 Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
537 Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);
539 const unsigned int temp_size =
542 (Eval::n_rows_of_product > Eval::n_columns_of_product ?
543 Eval::n_rows_of_product :
544 Eval::n_columns_of_product);
549 temp2 = temp1 +
std::max(Utilities::fixed_power<dim>(
550 shape_data.front().fe_degree + 1),
551 Utilities::fixed_power<dim>(
552 shape_data.front().n_q_points_1d));
556 temp2 = temp1 + temp_size;
559 const std::size_t n_q_points = temp_size == 0 ?
561 Eval::n_columns_of_product;
562 const unsigned int dofs_per_comp =
564 Utilities::fixed_power<dim>(shape_data.front().fe_degree + 1) :
567 Number *values_dofs =
569 temp1 + 2 * (std::max<std::size_t>(
581 for (
unsigned int c = 0; c < n_components; ++c)
585 if (add_into_values_array ==
false)
586 eval0.template values<0, false, false>(values_quad,
589 eval0.template values<0, false, true>(values_quad,
595 add_into_values_array ==
true)
596 eval0.template gradients<0, false, true>(gradients_quad,
599 eval0.template gradients<0, false, false>(gradients_quad,
606 add_into_values_array ==
true)
607 eval0.template hessians<0, false, true>(hessians_quad,
610 eval0.template hessians<0, false, false>(hessians_quad,
615 values_dofs += dofs_per_comp;
616 values_quad += n_q_points;
617 gradients_quad += n_q_points;
618 hessians_quad += n_q_points;
623 for (
unsigned int c = 0; c < n_components; ++c)
628 eval1.template values<1, false, false>(values_quad, temp1);
629 if (add_into_values_array ==
false)
630 eval0.template values<0, false, false>(temp1, values_dofs);
632 eval0.template values<0, false, true>(temp1, values_dofs);
636 eval1.template gradients<1, false, false>(gradients_quad +
640 eval1.template values<1, false, true>(values_quad, temp1);
641 if (add_into_values_array ==
false)
642 eval0.template values<0, false, false>(temp1, values_dofs);
644 eval0.template values<0, false, true>(temp1, values_dofs);
645 eval1.template values<1, false, false>(gradients_quad, temp1);
646 eval0.template gradients<0, false, true>(temp1, values_dofs);
651 eval1.template values<1, false, false>(hessians_quad, temp1);
655 add_into_values_array ==
true)
656 eval0.template hessians<0, false, true>(temp1, values_dofs);
658 eval0.template hessians<0, false, false>(temp1,
662 eval1.template hessians<1, false, false>(hessians_quad +
665 eval0.template values<0, false, true>(temp1, values_dofs);
668 eval1.template gradients<1, false, false>(hessians_quad +
671 eval0.template gradients<0, false, true>(temp1, values_dofs);
675 values_dofs += dofs_per_comp;
676 values_quad += n_q_points;
677 gradients_quad += 2 * n_q_points;
678 hessians_quad += 3 * n_q_points;
683 for (
unsigned int c = 0; c < n_components; ++c)
688 eval2.template values<2, false, false>(values_quad, temp1);
689 eval1.template values<1, false, false>(temp1, temp2);
690 if (add_into_values_array ==
false)
691 eval0.template values<0, false, false>(temp2, values_dofs);
693 eval0.template values<0, false, true>(temp2, values_dofs);
697 eval2.template gradients<2, false, false>(gradients_quad +
701 eval2.template values<2, false, true>(values_quad, temp1);
702 eval1.template values<1, false, false>(temp1, temp2);
703 eval2.template values<2, false, false>(gradients_quad +
706 eval1.template gradients<1, false, true>(temp1, temp2);
707 if (add_into_values_array ==
false)
708 eval0.template values<0, false, false>(temp2, values_dofs);
710 eval0.template values<0, false, true>(temp2, values_dofs);
711 eval2.template values<2, false, false>(gradients_quad, temp1);
712 eval1.template values<1, false, false>(temp1, temp2);
713 eval0.template gradients<0, false, true>(temp2, values_dofs);
718 eval2.template values<2, false, false>(hessians_quad, temp1);
719 eval1.template values<1, false, false>(temp1, temp2);
723 add_into_values_array ==
true)
724 eval0.template hessians<0, false, true>(temp2, values_dofs);
726 eval0.template hessians<0, false, false>(temp2,
730 eval2.template values<2, false, false>(hessians_quad +
733 eval1.template hessians<1, false, false>(temp1, temp2);
734 eval0.template values<0, false, true>(temp2, values_dofs);
737 eval2.template hessians<2, false, false>(hessians_quad +
740 eval1.template values<1, false, false>(temp1, temp2);
741 eval0.template values<0, false, true>(temp2, values_dofs);
744 eval2.template values<2, false, false>(hessians_quad +
747 eval1.template gradients<1, false, false>(temp1, temp2);
748 eval0.template gradients<0, false, true>(temp2, values_dofs);
751 eval2.template gradients<2, false, false>(hessians_quad +
754 eval1.template values<1, false, false>(temp1, temp2);
755 eval0.template gradients<0, false, true>(temp2, values_dofs);
758 eval2.template gradients<2, false, false>(hessians_quad +
761 eval1.template gradients<1, false, false>(temp1, temp2);
762 eval0.template values<0, false, true>(temp2, values_dofs);
766 values_dofs += dofs_per_comp;
767 values_quad += n_q_points;
768 gradients_quad += 3 * n_q_points;
769 hessians_quad += 6 * n_q_points;
780 values_dofs -= n_components * dofs_per_comp - dofs_per_comp + 1;
781 values_quad -= n_components * n_q_points;
783 for (
unsigned int c = 0; c < n_components; ++c)
785 values_dofs[0] = values_quad[0];
786 for (
unsigned int q = 1; q < n_q_points; ++q)
787 values_dofs[0] += values_quad[q];
788 values_dofs += dofs_per_comp;
789 values_quad += n_q_points;
793 for (
unsigned int c = 0; c < n_components; ++c)
794 values_dofs[c * dofs_per_comp] = Number();
795 values_dofs += n_components * dofs_per_comp;
801 const std::size_t n_dofs_per_comp =
803 values_dofs -= dofs_per_comp * n_components;
805 fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
806 for (
unsigned int c = 0; c < n_components; ++c)
807 for (
int i = 0, count_p = 0, count_q = 0;
808 i < (dim > 2 ? degree + 1 : 1);
811 for (
int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
813 for (
int k = 0; k < degree + 1 - j - i;
814 ++k, ++count_p, ++count_q)
815 values_dofs_actual[c * n_dofs_per_comp + count_p] =
816 values_dofs[c * dofs_per_comp + count_q];
819 count_q += i * (degree + 1);
826 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
833 Number>::evaluate(
const unsigned int n_components,
835 const Number * values_dofs_actual,
840 const std::size_t n_dofs =
851 const auto shape_values = shape_data.front().shape_values.data();
853 auto values_dofs_actual_ptr = values_dofs_actual;
855 Eval eval(shape_values,
nullptr,
nullptr, n_dofs, n_q_points);
856 for (
unsigned int c = 0; c < n_components; ++c)
858 eval.template values<0, true, false>(values_dofs_actual_ptr,
861 values_quad_ptr += n_q_points;
862 values_dofs_actual_ptr += n_dofs;
868 const auto shape_gradients = shape_data.front().shape_gradients.data();
870 auto values_dofs_actual_ptr = values_dofs_actual;
872 for (
unsigned int c = 0; c < n_components; ++c)
874 for (
unsigned int d = 0;
d < dim; ++
d)
877 shape_gradients + n_q_points * n_dofs *
d,
882 eval.template gradients<0, true, false>(values_dofs_actual_ptr,
885 gradients_quad_ptr += n_q_points;
887 values_dofs_actual_ptr += n_dofs;
894 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
901 Number>::integrate(
const unsigned int n_components,
903 Number * values_dofs_actual,
905 const bool add_into_values_array)
910 const std::size_t n_dofs =
921 const auto shape_values = shape_data.front().shape_values.data();
923 auto values_dofs_actual_ptr = values_dofs_actual;
925 Eval eval(shape_values,
nullptr,
nullptr, n_dofs, n_q_points);
926 for (
unsigned int c = 0; c < n_components; ++c)
928 if (add_into_values_array ==
false)
929 eval.template values<0, false, false>(values_quad_ptr,
930 values_dofs_actual_ptr);
932 eval.template values<0, false, true>(values_quad_ptr,
933 values_dofs_actual_ptr);
935 values_quad_ptr += n_q_points;
936 values_dofs_actual_ptr += n_dofs;
942 const auto shape_gradients = shape_data.front().shape_gradients.data();
944 auto values_dofs_actual_ptr = values_dofs_actual;
946 for (
unsigned int c = 0; c < n_components; ++c)
948 for (
unsigned int d = 0;
d < dim; ++
d)
951 shape_gradients + n_q_points * n_dofs *
d,
956 if ((add_into_values_array ==
false &&
959 eval.template gradients<0, false, false>(
960 gradients_quad_ptr, values_dofs_actual_ptr);
962 eval.template gradients<0, false, true>(
963 gradients_quad_ptr, values_dofs_actual_ptr);
965 gradients_quad_ptr += n_q_points;
967 values_dofs_actual_ptr += n_dofs;
973 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
974 template <
bool integrate>
981 evaluate_or_integrate(
983 Number * values_dofs_actual,
985 const bool add_into_values_array)
992 evaluate_tensor_product_per_component<0>(
996 add_into_values_array,
997 std::integral_constant<bool, integrate>());
999 evaluate_tensor_product_per_component<1>(
1003 add_into_values_array,
1004 std::integral_constant<bool, integrate>());
1008 evaluate_tensor_product_per_component<2>(
1012 add_into_values_array,
1013 std::integral_constant<bool, integrate>());
1020 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
1021 template <
int normal_dir>
1028 evaluate_tensor_product_per_component(
1030 Number * values_dofs_actual,
1032 const bool add_into_values_array,
1033 std::integral_constant<bool, false>)
1035 (void)add_into_values_array;
1040 (fe_degree == -1) ? 1 : fe_degree + 1,
1048 (fe_degree == -1) ? 1 : fe_degree,
1053 typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
1055 typename std::conditional<normal_dir == 1, EvalNormal, EvalTangent>::type;
1057 typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
1061 Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
1062 ((normal_dir == 0) ? shape_info.
data[0] : shape_info.
data[1]));
1063 Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
1064 ((normal_dir == 1) ? shape_info.
data[0] : shape_info.
data[1]));
1065 Eval2 eval2 = create_evaluator_tensor_product<Eval2>(
1066 ((normal_dir == 2) ? shape_info.
data[0] : shape_info.
data[1]));
1073 std::max(Utilities::fixed_power<dim>(shape_info.
data[0].fe_degree + 1),
1074 Utilities::fixed_power<dim>(shape_info.
data[0].n_q_points_1d));
1076 const std::size_t n_q_points = shape_info.
n_q_points;
1080 Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
1081 Number *values_quad = fe_eval.
begin_values() + n_q_points * normal_dir;
1082 Number *gradients_quad =
1084 Number *hessians_quad =
1085 (dim == 2) ? fe_eval.
begin_hessians() + 3 * n_q_points * normal_dir :
1093 eval0.template gradients<0, true, false>(values_dofs, temp1);
1094 eval1.template values<1, true, false>(temp1, gradients_quad);
1103 eval0.template gradients<0, true, false>(values_dofs, temp1);
1104 eval1.template gradients<1, true, false>(temp1,
1109 eval0.template hessians<0, true, false>(values_dofs, temp1);
1110 eval1.template values<1, true, false>(temp1, hessians_quad);
1114 eval0.template values<0, true, false>(values_dofs, temp1);
1116 eval1.template gradients<1, true, false>(temp1,
1122 eval1.template hessians<1, true, false>(temp1,
1123 hessians_quad + n_q_points);
1127 eval1.template values<1, true, false>(temp1, values_quad);
1133 eval0.template gradients<0, true, false>(values_dofs, temp1);
1134 eval1.template values<1, true, false>(temp1, temp2);
1135 eval2.template values<2, true, false>(temp2, gradients_quad);
1146 eval0.template gradients<0, true, false>(values_dofs, temp1);
1147 eval1.template values<1, true, false>(temp1, temp2);
1149 eval2.template gradients<2, true, false>(temp2,
1154 eval1.template gradients<1, true, false>(temp1, temp2);
1155 eval2.template values<2, true, false>(temp2,
1160 eval0.template hessians<0, true, false>(values_dofs, temp1);
1161 eval1.template values<1, true, false>(temp1, temp2);
1162 eval2.template values<2, true, false>(temp2, hessians_quad);
1166 eval0.template values<0, true, false>(values_dofs, temp1);
1169 eval1.template gradients<1, true, false>(temp1, temp2);
1170 eval2.template values<2, true, false>(temp2,
1179 eval1.template gradients<1, true, false>(temp1, temp2);
1180 eval2.template gradients<2, true, false>(temp2,
1185 eval1.template hessians<1, true, false>(temp1, temp2);
1186 eval2.template values<2, true, false>(temp2,
1187 hessians_quad + n_q_points);
1192 eval1.template values<1, true, false>(temp1, temp2);
1194 eval2.template gradients<2, true, false>(temp2,
1201 eval2.template hessians<2, true, false>(temp2,
1208 eval2.template values<2, true, false>(temp2, values_quad);
1215 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
1216 template <
int normal_dir>
1223 evaluate_tensor_product_per_component(
1225 Number * values_dofs_actual,
1227 const bool add_into_values_array,
1228 std::integral_constant<bool, true>)
1233 (fe_degree == -1) ? 1 : fe_degree + 1,
1241 (fe_degree == -1) ? 1 : fe_degree,
1246 typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
1248 typename std::conditional<normal_dir == 1, EvalNormal, EvalTangent>::type;
1250 typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
1254 Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
1255 ((normal_dir == 0) ? shape_info.
data[0] : shape_info.
data[1]));
1256 Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
1257 ((normal_dir == 1) ? shape_info.
data[0] : shape_info.
data[1]));
1258 Eval2 eval2 = create_evaluator_tensor_product<Eval2>(
1259 ((normal_dir == 2) ? shape_info.
data[0] : shape_info.
data[1]));
1266 std::max(Utilities::fixed_power<dim>(shape_info.
data[0].fe_degree + 1),
1267 Utilities::fixed_power<dim>(shape_info.
data[0].n_q_points_1d));
1269 const std::size_t n_q_points = shape_info.
n_q_points;
1273 Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
1274 Number *values_quad = fe_eval.
begin_values() + n_q_points * normal_dir;
1275 Number *gradients_quad =
1277 Number *hessians_quad =
1278 (dim == 2) ? fe_eval.
begin_hessians() + 3 * n_q_points * normal_dir :
1288 eval1.template values<1, false, false>(values_quad, temp1);
1289 if (add_into_values_array ==
false)
1290 eval0.template values<0, false, false>(temp1, values_dofs);
1292 eval0.template values<0, false, true>(temp1, values_dofs);
1296 eval1.template gradients<1, false, false>(gradients_quad +
1300 eval1.template values<1, false, true>(values_quad, temp1);
1301 if (add_into_values_array ==
false)
1302 eval0.template values<0, false, false>(temp1, values_dofs);
1304 eval0.template values<0, false, true>(temp1, values_dofs);
1305 eval1.template values<1, false, false>(gradients_quad, temp1);
1306 eval0.template gradients<0, false, true>(temp1, values_dofs);
1311 eval1.template values<1, false, false>(hessians_quad, temp1);
1315 add_into_values_array ==
true)
1316 eval0.template hessians<0, false, true>(temp1, values_dofs);
1318 eval0.template hessians<0, false, false>(temp1, values_dofs);
1321 eval1.template hessians<1, false, false>(hessians_quad +
1324 eval0.template values<0, false, true>(temp1, values_dofs);
1327 eval1.template gradients<1, false, false>(hessians_quad +
1330 eval0.template gradients<0, false, true>(temp1, values_dofs);
1338 eval2.template values<2, false, false>(values_quad, temp1);
1339 eval1.template values<1, false, false>(temp1, temp2);
1340 if (add_into_values_array ==
false)
1341 eval0.template values<0, false, false>(temp2, values_dofs);
1343 eval0.template values<0, false, true>(temp2, values_dofs);
1347 eval2.template gradients<2, false, false>(gradients_quad +
1351 eval2.template values<2, false, true>(values_quad, temp1);
1352 eval1.template values<1, false, false>(temp1, temp2);
1353 eval2.template values<2, false, false>(gradients_quad +
1356 eval1.template gradients<1, false, true>(temp1, temp2);
1357 if (add_into_values_array ==
false)
1358 eval0.template values<0, false, false>(temp2, values_dofs);
1360 eval0.template values<0, false, true>(temp2, values_dofs);
1361 eval2.template values<2, false, false>(gradients_quad, temp1);
1362 eval1.template values<1, false, false>(temp1, temp2);
1363 eval0.template gradients<0, false, true>(temp2, values_dofs);
1368 eval2.template values<2, false, false>(hessians_quad, temp1);
1369 eval1.template values<1, false, false>(temp1, temp2);
1373 add_into_values_array ==
true)
1374 eval0.template hessians<0, false, true>(temp2, values_dofs);
1376 eval0.template hessians<0, false, false>(temp2, values_dofs);
1379 eval2.template values<2, false, false>(hessians_quad + n_q_points,
1381 eval1.template hessians<1, false, false>(temp1, temp2);
1382 eval0.template values<0, false, true>(temp2, values_dofs);
1385 eval2.template hessians<2, false, false>(hessians_quad +
1388 eval1.template values<1, false, false>(temp1, temp2);
1389 eval0.template values<0, false, true>(temp2, values_dofs);
1392 eval2.template values<2, false, false>(hessians_quad +
1395 eval1.template gradients<1, false, false>(temp1, temp2);
1396 eval0.template gradients<0, false, true>(temp2, values_dofs);
1399 eval2.template gradients<2, false, false>(hessians_quad +
1402 eval1.template values<1, false, false>(temp1, temp2);
1403 eval0.template gradients<0, false, true>(temp2, values_dofs);
1406 eval2.template gradients<2, false, false>(hessians_quad +
1409 eval1.template gradients<1, false, false>(temp1, temp2);
1410 eval0.template values<0, false, true>(temp2, values_dofs);
1437 static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
1438 "The second dimension must not be smaller than the first");
1467 const unsigned int n_components,
1469 const Number * values_in,
1470 Number * values_out,
1475 basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
1476 ExcMessage(
"The second dimension must not be smaller than the first"));
1484 constexpr
int next_dim =
1486 ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
1493 (basis_size_1 == 0 ? 0 : basis_size_2),
1496 eval_val(transformation_matrix,
1499 basis_size_1_variable,
1500 basis_size_2_variable);
1501 const unsigned int np_1 =
1502 basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
1503 const unsigned int np_2 =
1504 basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
1506 ExcMessage(
"Cannot transform with 0-point basis"));
1508 ExcMessage(
"Cannot transform with 0-point basis"));
1512 values_in = values_in + n_components * Utilities::fixed_power<dim>(np_1);
1514 values_out + n_components * Utilities::fixed_power<dim>(np_2);
1515 for (
unsigned int c = n_components; c != 0; --c)
1517 values_in -= Utilities::fixed_power<dim>(np_1);
1518 values_out -= Utilities::fixed_power<dim>(np_2);
1520 for (
unsigned int q = np_1; q != 0; --q)
1529 transformation_matrix,
1532 Utilities::fixed_power<next_dim>(np_1),
1535 Utilities::fixed_power<next_dim>(np_2),
1536 basis_size_1_variable,
1537 basis_size_2_variable);
1542 if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
1544 eval_val.template values<0, true, false>(values_in, values_out);
1545 eval_val.template values<1, true, false>(values_out, values_out);
1548 eval_val.template
values<dim - 1,
true,
false>(values_in,
1551 eval_val.template
values<dim - 1,
true,
false>(values_out,
1591 const unsigned int n_components,
1593 const bool add_into_result,
1595 Number * values_out,
1600 basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
1601 ExcMessage(
"The second dimension must not be smaller than the first"));
1602 Assert(add_into_result ==
false || values_in != values_out,
1604 "Input and output cannot alias with each other when "
1605 "adding the result of the basis change to existing data"));
1611 constexpr
int next_dim =
1613 ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
1619 (basis_size_1 == 0 ? 0 : basis_size_2),
1622 eval_val(transformation_matrix,
1623 transformation_matrix,
1624 transformation_matrix,
1625 basis_size_1_variable,
1626 basis_size_2_variable);
1627 const unsigned int np_1 =
1628 basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
1629 const unsigned int np_2 =
1630 basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
1632 ExcMessage(
"Cannot transform with 0-point basis"));
1634 ExcMessage(
"Cannot transform with 0-point basis"));
1636 for (
unsigned int c = 0; c < n_components; ++c)
1638 if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
1641 eval_val.template values<1, false, false>(values_in, values_in);
1643 eval_val.template hessians<1, false, false>(values_in,
1646 if (add_into_result)
1649 eval_val.template values<0, false, true>(values_in,
1652 eval_val.template hessians<0, false, true>(values_in,
1658 eval_val.template values<0, false, false>(values_in,
1661 eval_val.template hessians<0, false, false>(values_in,
1667 if (dim == 1 && add_into_result)
1670 eval_val.template values<0, false, true>(values_in,
1673 eval_val.template hessians<0, false, true>(values_in,
1679 eval_val.template values<0, false, false>(values_in,
1682 eval_val.template hessians<0, false, false>(values_in,
1688 eval_val.template
values<dim - 1,
false,
false>(values_in,
1691 eval_val.template
hessians<dim - 1,
false,
false>(
1692 values_in, values_in);
1696 for (
unsigned int q = 0; q < np_1; ++q)
1705 transformation_matrix,
1708 q * Utilities::fixed_power<next_dim>(np_2),
1710 q * Utilities::fixed_power<next_dim>(np_1),
1711 basis_size_1_variable,
1712 basis_size_2_variable);
1714 values_in += Utilities::fixed_power<dim>(np_2);
1715 values_out += Utilities::fixed_power<dim>(np_1);
1743 const Number * values_in,
1744 Number * scratch_data,
1745 Number * values_out)
1747 constexpr
int next_dim = dim > 1 ? dim - 1 : dim;
1748 Number * my_scratch =
1749 basis_size_1 != basis_size_2 ? scratch_data : values_out;
1751 const unsigned int size_per_component =
Utilities::pow(basis_size_2, dim);
1752 Assert(coefficients.
size() == size_per_component ||
1753 coefficients.
size() == n_components * size_per_component,
1755 const unsigned int stride =
1756 coefficients.
size() == size_per_component ? 0 : 1;
1758 for (
unsigned int q = basis_size_1; q != 0; --q)
1767 transformation_matrix,
1780 eval_val(transformation_matrix);
1781 const unsigned int n_inner_blocks =
1782 (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
1784 for (
unsigned int ii = 0; ii <
n_blocks; ii += n_inner_blocks)
1785 for (
unsigned int c = 0; c < n_components; ++c)
1787 for (
unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1788 eval_val.template values_one_line<dim - 1, true, false>(
1789 my_scratch + i, my_scratch + i);
1790 for (
unsigned int q = 0; q < basis_size_2; ++q)
1791 for (
unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1792 my_scratch[i + q *
n_blocks + c * size_per_component] *=
1794 c * stride * size_per_component];
1795 for (
unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1796 eval_val.template values_one_line<dim - 1, false, false>(
1797 my_scratch + i, my_scratch + i);
1799 for (
unsigned int q = 0; q < basis_size_1; ++q)
1808 transformation_matrix,
1831 template <
int dim,
int fe_degree,
typename Number>
1835 evaluate(
const unsigned int n_components,
1837 const Number * values_dofs,
1843 const Number * values_dofs,
1844 Number * gradients_quad,
1845 Number * hessians_quad);
1848 integrate(
const unsigned int n_components,
1850 Number * values_dofs,
1852 const bool add_into_values_array);
1857 Number * values_dofs,
1858 Number * gradients_quad,
1859 const Number * hessians_quad,
1860 const bool add_into_values_array);
1865 template <
int dim,
int fe_degree,
typename Number>
1868 const unsigned int n_components,
1870 const Number * values_dofs,
1873 constexpr std::size_t n_points =
Utilities::pow(fe_degree + 1, dim);
1875 for (
unsigned int c = 0; c < n_components; ++c)
1878 for (
unsigned int i = 0; i < n_points; ++i)
1880 values_dofs[n_points * c + i];
1884 values_dofs + c * n_points,
1887 c * dim * (dim + 1) / 2 * n_points);
1893 template <
int dim,
int fe_degree,
typename Number>
1898 const Number * values_dofs,
1899 Number * gradients_quad,
1900 Number * hessians_quad)
1903 (fe_degree + 2) / 2 * (fe_degree + 1));
1904 constexpr std::size_t n_points =
Utilities::pow(fe_degree + 1, dim);
1914 if ((evaluation_flag &
1917 eval.template gradients<0, true, false>(values_dofs, gradients_quad);
1919 eval.template gradients<1, true, false>(values_dofs,
1920 gradients_quad + n_points);
1922 eval.template gradients<2, true, false>(values_dofs,
1928 eval.template hessians<0, true, false>(values_dofs, hessians_quad);
1931 eval.template gradients<1, true, false>(gradients_quad,
1934 eval.template hessians<1, true, false>(values_dofs,
1935 hessians_quad + n_points);
1939 eval.template gradients<2, true, false>(gradients_quad,
1942 eval.template gradients<2, true, false>(gradients_quad + n_points,
1945 eval.template hessians<2, true, false>(values_dofs,
1954 template <
int dim,
int fe_degree,
typename Number>
1957 const unsigned int n_components,
1959 Number * values_dofs,
1961 const bool add_into_values_array)
1963 constexpr std::size_t n_points =
Utilities::pow(fe_degree + 1, dim);
1965 for (
unsigned int c = 0; c < n_components; ++c)
1969 if (add_into_values_array)
1970 for (
unsigned int i = 0; i < n_points; ++i)
1971 values_dofs[n_points * c + i] +=
1974 for (
unsigned int i = 0; i < n_points; ++i)
1975 values_dofs[n_points * c + i] =
1981 values_dofs + c * n_points,
1984 c * dim * (dim + 1) / 2 * n_points,
1985 add_into_values_array ||
1992 template <
int dim,
int fe_degree,
typename Number>
1997 Number * values_dofs,
1998 Number * gradients_quad,
1999 const Number * hessians_quad,
2000 const bool add_into_values_array)
2003 (fe_degree + 2) / 2 * (fe_degree + 1));
2013 constexpr std::size_t n_points =
Utilities::pow(fe_degree + 1, dim);
2019 if (add_into_values_array ==
true)
2020 eval.template hessians<0, false, true>(hessians_quad, values_dofs);
2022 eval.template hessians<0, false, false>(hessians_quad, values_dofs);
2025 eval.template hessians<1, false, true>(hessians_quad + n_points,
2029 eval.template hessians<2, false, true>(hessians_quad + 2 * n_points,
2036 eval.template gradients<1, false, true>(hessians_quad +
2040 eval.template gradients<1, false, false>(hessians_quad +
2048 eval.template gradients<1, false, true>(hessians_quad +
2052 eval.template gradients<1, false, false>(hessians_quad +
2057 eval.template gradients<2, false, true>(hessians_quad +
2063 eval.template gradients<2, false, true>(
2064 hessians_quad + 5 * n_points, gradients_quad + n_points);
2066 eval.template gradients<2, false, false>(
2067 hessians_quad + 5 * n_points, gradients_quad + n_points);
2074 for (
unsigned int q = 0; q < n_points; ++q)
2075 gradients_quad[(dim - 1) * n_points + q] = Number();
2078 if ((integration_flag &
2081 if (add_into_values_array ||
2083 eval.template gradients<0, false, true>(gradients_quad, values_dofs);
2085 eval.template gradients<0, false, false>(gradients_quad, values_dofs);
2087 eval.template gradients<1, false, true>(gradients_quad + n_points,
2090 eval.template gradients<2, false, true>(gradients_quad + 2 * n_points,
2107 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
2111 evaluate(
const unsigned int n_components,
2113 const Number * values_dofs,
2117 integrate(
const unsigned int n_components,
2119 Number * values_dofs,
2121 const bool add_into_values_array);
2126 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
2132 Number>::evaluate(
const unsigned int n_components,
2134 const Number * values_dofs,
2139 Assert(n_q_points_1d > fe_degree,
2140 ExcMessage(
"You lose information when going to a collocation space "
2141 "of lower degree, so the evaluation results would be "
2142 "wrong. Thus, this class does not permit the desired "
2144 constexpr std::size_t n_dofs =
Utilities::pow(fe_degree + 1, dim);
2145 constexpr std::size_t n_q_points =
Utilities::pow(n_q_points_1d, dim);
2147 for (
unsigned int c = 0; c < n_components; ++c)
2153 (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
2156 Number>::do_forward(1,
2157 shape_data.shape_values_eo,
2158 values_dofs + c * n_dofs,
2162 if (evaluation_flag &
2165 do_evaluate(shape_data,
2171 c * dim * (dim + 1) / 2 * n_q_points);
2177 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
2183 Number>::integrate(
const unsigned int n_components,
2185 Number * values_dofs,
2187 const bool add_into_values_array)
2191 Assert(n_q_points_1d > fe_degree,
2192 ExcMessage(
"You lose information when going to a collocation space "
2193 "of lower degree, so the evaluation results would be "
2194 "wrong. Thus, this class does not permit the desired "
2196 constexpr std::size_t n_q_points =
Utilities::pow(n_q_points_1d, dim);
2198 for (
unsigned int c = 0; c < n_components; ++c)
2201 if (integration_flag &
2210 c * dim * (dim + 1) / 2 * n_q_points,
2219 (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
2222 Number>::do_backward(1,
2223 shape_data.shape_values_eo,
2224 add_into_values_array,
2242 const unsigned int n_q_points_1d)
2244 return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
2245 (n_q_points_1d <= 3 * fe_degree / 2 + 1);
2264 template <
int dim,
typename Number>
2267 template <
int fe_degree,
int n_q_po
ints_1d>
2269 run(
const unsigned int n_components,
2271 const Number * values_dofs,
2283 if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
2287 n_components, evaluation_flag, values_dofs, fe_eval);
2291 else if (fe_degree >= 0 &&
2299 Number>::evaluate(n_components,
2304 else if (fe_degree >= 0 &&
2311 Number>::evaluate(n_components,
2322 Number>::evaluate(n_components,
2333 Number>::evaluate(n_components,
2344 Number>::evaluate(n_components,
2354 (fe_degree == -1) ? 1 : fe_degree,
2355 (n_q_points_1d < 1) ? 1 : n_q_points_1d,
2356 Number>::template evaluate_or_integrate<false>(evaluation_flag,
2357 const_cast<Number *
>(
2367 Number>::evaluate(n_components,
2394 template <
int dim,
typename Number>
2397 template <
int fe_degree,
int n_q_po
ints_1d>
2399 run(
const unsigned int n_components,
2401 Number * values_dofs,
2403 const bool sum_into_values_array)
2414 if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
2422 sum_into_values_array);
2426 else if (fe_degree >= 0 &&
2434 Number>::integrate(n_components,
2438 sum_into_values_array);
2440 else if (fe_degree >= 0 &&
2447 Number>::integrate(n_components,
2451 sum_into_values_array);
2459 Number>::integrate(n_components,
2463 sum_into_values_array);
2471 Number>::integrate(n_components,
2475 sum_into_values_array);
2483 Number>::integrate(n_components,
2487 sum_into_values_array);
2493 (fe_degree == -1) ? 1 : fe_degree,
2494 (n_q_points_1d < 1) ? 1 : n_q_points_1d,
2496 template evaluate_or_integrate<true>(integration_flag,
2497 const_cast<Number *
>(
2500 sum_into_values_array);
2508 Number>::integrate(n_components,
2512 sum_into_values_array);
2521 template <
bool symmetric_evaluate,
2543 const unsigned int subface_index,
2544 const unsigned int direction)
2546 if (symmetric_evaluate)
2560 const unsigned int index =
2561 direction == 0 ? subface_index % 2 : subface_index / 2;
2572 const unsigned int n_components,
2575 Number * values_dofs,
2576 Number * values_quad,
2577 Number * gradients_quad,
2578 Number * hessians_quad,
2579 Number * scratch_data,
2580 const unsigned int subface_index)
2585 const std::size_t n_dofs = fe_degree > -1 ?
2588 const std::size_t n_q_points =
2593 Number *values_dofs_ptr = values_dofs;
2597 for (
unsigned int c = 0; c < n_components; ++c)
2602 eval0.template values<0, true, false>(values_dofs,
2604 eval1.template values<1, true, false>(values_quad,
2608 eval0.template values<0, true, false>(values_dofs,
2612 values_quad[0] = values_dofs[0];
2619 values_dofs += 3 * n_dofs;
2620 values_quad += n_q_points;
2623 for (
unsigned int c = 0; c < n_components; ++c)
2628 if (symmetric_evaluate &&
2631 eval0.template values<0, true, false>(values_dofs,
2633 eval0.template values<1, true, false>(values_quad,
2643 eval_grad.template gradients<0, true, false>(
2644 values_quad, gradients_quad);
2645 eval_grad.template gradients<1, true, false>(
2646 values_quad, gradients_quad + n_q_points);
2651 eval0.template gradients<0, true, false>(values_dofs,
2653 eval1.template values<1, true, false>(scratch_data,
2657 eval0.template values<0, true, false>(values_dofs,
2659 eval1.template gradients<1, true, false>(scratch_data,
2664 eval1.template values<1, true, false>(scratch_data,
2668 eval0.template values<0, true, false>(values_dofs + n_dofs,
2670 eval1.template values<1, true, false>(
2671 scratch_data, gradients_quad + (dim - 1) * n_q_points);
2675 eval0.template values<0, true, false>(values_dofs + n_dofs,
2678 eval0.template gradients<0, true, false>(values_dofs,
2681 eval0.template values<0, true, false>(values_dofs,
2685 values_quad[0] = values_dofs[0];
2686 gradients_quad[0] = values_dofs[1];
2691 values_dofs += 3 * n_dofs;
2692 values_quad += n_q_points;
2693 gradients_quad += dim * n_q_points;
2698 values_dofs = values_dofs_ptr;
2699 for (
unsigned int c = 0; c < n_components; ++c)
2705 eval0.template hessians<0, true, false>(values_dofs,
2707 eval1.template values<1, true, false>(scratch_data,
2711 eval0.template values<0, true, false>(values_dofs,
2713 eval1.template hessians<1, true, false>(scratch_data,
2718 eval0.template values<0, true, false>(values_dofs +
2721 eval1.template values<1, true, false>(scratch_data,
2726 eval0.template gradients<0, true, false>(values_dofs,
2728 eval1.template gradients<1, true, false>(scratch_data,
2733 eval0.template gradients<0, true, false>(values_dofs +
2736 eval1.template values<1, true, false>(scratch_data,
2741 eval0.template values<0, true, false>(values_dofs + n_dofs,
2743 eval1.template gradients<1, true, false>(scratch_data,
2750 eval0.template hessians<0, true, false>(values_dofs,
2753 eval0.template values<0, true, false>(
2754 values_dofs + 2 * n_dofs, hessians_quad + n_q_points);
2756 eval0.template gradients<0, true, false>(
2757 values_dofs + n_dofs, hessians_quad + 2 * n_q_points);
2760 hessians_quad[0] = values_dofs[2];
2765 values_dofs += 3 * n_dofs;
2766 hessians_quad += dim * (dim + 1) / 2 * n_q_points;
2773 const unsigned int n_components,
2776 Number * values_dofs,
2777 Number * values_quad,
2778 Number * gradients_quad,
2779 Number * hessians_quad,
2780 Number * scratch_data,
2781 const unsigned int subface_index)
2786 const std::size_t n_dofs =
2790 const std::size_t n_q_points =
2795 Number *values_dofs_ptr = values_dofs;
2799 for (
unsigned int c = 0; c < n_components; ++c)
2804 eval1.template values<1, false, false>(values_quad,
2806 eval0.template values<0, false, false>(values_quad,
2810 eval0.template values<0, false, false>(values_quad,
2814 values_dofs[0] = values_quad[0];
2819 values_dofs += 3 * n_dofs;
2820 values_quad += n_q_points;
2823 for (
unsigned int c = 0; c < n_components; ++c)
2829 eval1.template values<1, false, false>(gradients_quad +
2833 eval0.template values<0, false, false>(gradients_quad +
2835 values_dofs + n_dofs);
2836 if (symmetric_evaluate &&
2848 eval_grad.template gradients<1, false, true>(
2849 gradients_quad + n_q_points, values_quad);
2851 eval_grad.template gradients<1, false, false>(
2852 gradients_quad + n_q_points, values_quad);
2853 eval_grad.template gradients<0, false, true>(
2854 gradients_quad, values_quad);
2855 eval0.template values<1, false, false>(values_quad,
2857 eval0.template values<0, false, false>(values_quad,
2864 eval1.template values<1, false, false>(values_quad,
2866 eval1.template gradients<1, false, true>(
2867 gradients_quad + n_q_points, scratch_data);
2870 eval1.template gradients<1, false, false>(
2871 gradients_quad + n_q_points, scratch_data);
2874 eval0.template values<0, false, false>(scratch_data,
2878 eval1.template values<1, false, false>(gradients_quad,
2880 eval0.template gradients<0, false, true>(scratch_data,
2885 eval0.template values<0, false, false>(gradients_quad +
2887 values_dofs + n_dofs);
2888 eval0.template gradients<0, false, false>(gradients_quad,
2891 eval0.template values<0, false, true>(values_quad,
2895 values_dofs[0] = values_quad[0];
2896 values_dofs[1] = gradients_quad[0];
2901 values_dofs += 3 * n_dofs;
2902 values_quad += n_q_points;
2903 gradients_quad += dim * n_q_points;
2908 values_dofs = values_dofs_ptr;
2909 for (
unsigned int c = 0; c < n_components; ++c)
2915 eval1.template values<1, false, false>(hessians_quad,
2919 eval0.template hessians<0, false, true>(scratch_data,
2922 eval0.template hessians<0, false, false>(scratch_data,
2926 eval1.template hessians<1, false, false>(hessians_quad +
2929 eval0.template values<0, false, true>(scratch_data,
2933 eval1.template values<1, false, false>(hessians_quad +
2936 eval0.template values<0, false, false>(scratch_data,
2941 eval1.template gradients<1, false, false>(hessians_quad +
2944 eval0.template gradients<0, false, true>(scratch_data,
2948 eval1.template values<1, false, false>(hessians_quad +
2952 eval0.template gradients<0, false, true>(scratch_data,
2956 eval0.template gradients<0, false, false>(scratch_data,
2961 eval1.template gradients<1, false, false>(hessians_quad +
2964 eval0.template values<0, false, true>(scratch_data,
2965 values_dofs + n_dofs);
2972 eval0.template hessians<0, false, true>(hessians_quad,
2975 eval0.template hessians<0, false, false>(hessians_quad,
2979 eval0.template values<0, false, false>(
2980 hessians_quad + n_q_points, values_dofs + 2 * n_dofs);
2983 eval0.template gradients<0, false, true>(
2984 hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
2986 eval0.template gradients<0, false, false>(
2987 hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
2990 values_dofs[2] = hessians_quad[0];
2999 values_dofs += 3 * n_dofs;
3000 hessians_quad += dim * (dim + 1) / 2 * n_q_points;
3006 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
typename Number>
3014 template <
typename EvalType>
3018 const unsigned int subface_index,
3019 const unsigned int direction)
3027 const unsigned int index =
3028 direction == 0 ? subface_index % 2 : subface_index / 2;
3035 template <
bool integrate>
3039 Number * values_dofs,
3041 Number * scratch_data,
3042 const unsigned int subface_index,
3043 const unsigned int face_no)
3045 const unsigned int face_direction = face_no / 2;
3065 evaluate_in_face_apply<0>(values_dofs,
3071 std::integral_constant<bool, integrate>());
3074 evaluate_in_face_apply<1>(values_dofs,
3080 std::integral_constant<bool, integrate>());
3082 evaluate_in_face_apply<2>(values_dofs,
3088 std::integral_constant<bool, integrate>());
3099 template <
int normal_dir>
3102 Number * values_dofs,
3104 Number * scratch_data,
3106 const unsigned int face_direction,
3107 const unsigned int subface_index,
3108 std::integral_constant<bool, false>)
3113 (fe_degree == -1) ? 1 : fe_degree + 1,
3120 (fe_degree == -1) ? 1 : fe_degree,
3125 using TempEval0 =
typename std::
3126 conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
3127 using TempEval1 =
typename std::
3128 conditional<normal_dir == 0, EvalTangent, EvalNormal>::type;
3129 using Eval0 =
typename std::
3130 conditional<normal_dir == 2, EvalGeneral, TempEval0>::type;
3131 using Eval1 =
typename std::
3132 conditional<normal_dir == 2, EvalGeneral, TempEval1>::type;
3136 Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
3137 ((normal_dir == 0) ? shape_info.
data[0] : shape_info.
data[1]),
3140 Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
3141 ((normal_dir == 1) ? shape_info.
data[0] : shape_info.
data[1]),
3145 constexpr std::size_t n_q_points =
Utilities::pow(n_q_points_1d, dim - 1);
3147 const std::size_t n_dofs_normal =
3149 const std::size_t dofs_stride =
3150 (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
3153 static constexpr ::ndarray<unsigned int, 3, 3> component_table = {
3154 {{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
3155 const unsigned int component =
3156 (dim == 2 && normal_dir == 0 && face_direction == 1) ?
3158 component_table[face_direction][normal_dir];
3162 3 * ((component == 0) ?
3165 ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) :
3166 ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent :
3167 n_dofs_normal + n_dofs_tangent)));
3168 const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir;
3170 Number *gradients_quad =
3172 Number *hessians_quad =
3182 eval0.template values<0, true, false>(values_dofs, values_quad);
3183 eval1.template values<1, true, false>(values_quad, values_quad);
3186 eval0.template values<0, true, false>(values_dofs, values_quad);
3198 eval0.template gradients<0, true, false>(values_dofs,
3200 eval1.template values<1, true, false>(scratch_data,
3204 eval0.template values<0, true, false>(values_dofs,
3206 eval1.template gradients<1, true, false>(scratch_data,
3211 eval1.template values<1, true, false>(scratch_data,
3215 eval0.template values<0, true, false>(values_dofs + dofs_stride,
3217 eval1.template values<1, true, false>(scratch_data,
3223 eval0.template values<0, true, false>(values_dofs + dofs_stride,
3226 eval0.template gradients<0, true, false>(values_dofs,
3229 eval0.template values<0, true, false>(values_dofs,
3243 eval0.template hessians<0, true, false>(values_dofs,
3245 eval1.template values<1, true, false>(scratch_data,
3249 eval0.template values<0, true, false>(values_dofs,
3251 eval1.template hessians<1, true, false>(scratch_data,
3256 eval0.template values<0, true, false>(values_dofs +
3259 eval1.template values<1, true, false>(scratch_data,
3264 eval0.template gradients<0, true, false>(values_dofs,
3266 eval1.template gradients<1, true, false>(scratch_data,
3271 eval0.template gradients<0, true, false>(values_dofs +
3274 eval1.template values<1, true, false>(scratch_data,
3279 eval0.template values<0, true, false>(values_dofs + dofs_stride,
3281 eval1.template gradients<1, true, false>(scratch_data,
3288 eval0.template hessians<0, true, false>(values_dofs,
3291 eval0.template values<0, true, false>(
3292 values_dofs + 2 * dofs_stride, hessians_quad + n_q_points);
3294 eval0.template gradients<0, true, false>(
3295 values_dofs + dofs_stride, hessians_quad + 2 * n_q_points);
3303 template <
int normal_dir>
3306 Number * values_dofs,
3308 Number * scratch_data,
3310 const unsigned int face_direction,
3311 const unsigned int subface_index,
3312 std::integral_constant<bool, true>)
3317 (fe_degree == -1) ? 1 : fe_degree + 1,
3324 (fe_degree == -1) ? 1 : fe_degree,
3329 using TempEval0 =
typename std::
3330 conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
3331 using TempEval1 =
typename std::
3332 conditional<normal_dir == 0, EvalTangent, EvalNormal>::type;
3333 using Eval0 =
typename std::
3334 conditional<normal_dir == 2, EvalGeneral, TempEval0>::type;
3335 using Eval1 =
typename std::
3336 conditional<normal_dir == 2, EvalGeneral, TempEval1>::type;
3340 Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
3341 ((normal_dir == 0) ? shape_info.
data[0] : shape_info.
data[1]),
3344 Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
3345 ((normal_dir == 1) ? shape_info.
data[0] : shape_info.
data[1]),
3349 constexpr std::size_t n_q_points =
Utilities::pow(n_q_points_1d, dim - 1);
3351 const std::size_t n_dofs_normal =
3353 const std::size_t dofs_stride =
3354 (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
3357 static constexpr ::ndarray<unsigned int, 3, 3> component_table = {
3358 {{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
3359 const unsigned int component =
3360 (dim == 2 && normal_dir == 0 && face_direction == 1) ?
3362 component_table[face_direction][normal_dir];
3366 3 * ((component == 0) ?
3369 ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) :
3370 ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent :
3371 n_dofs_normal + n_dofs_tangent)));
3372 const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir;
3374 Number *gradients_quad =
3376 Number *hessians_quad =
3386 eval1.template values<1, false, false>(values_quad,
3388 eval0.template values<0, false, false>(values_quad,
3392 eval0.template values<0, false, false>(values_quad,
3405 eval1.template values<1, false, false>(gradients_quad +
3409 eval0.template values<0, false, false>(
3410 gradients_quad + 2 * n_q_points, values_dofs + dofs_stride);
3414 eval1.template values<1, false, false>(values_quad,
3416 eval1.template gradients<1, false, true>(gradients_quad +
3421 eval1.template gradients<1, false, false>(gradients_quad +
3426 eval0.template values<0, false, false>(scratch_data,
3430 eval1.template values<1, false, false>(gradients_quad,
3432 eval0.template gradients<0, false, true>(scratch_data,
3437 eval0.template values<0, false, false>(
3438 gradients_quad + n_q_points, values_dofs + dofs_stride);
3439 eval0.template gradients<0, false, false>(gradients_quad,
3442 eval0.template values<0, false, true>(values_quad,
3456 eval1.template values<1, false, false>(hessians_quad,
3458 if ((evaluation_flag &
3460 eval0.template hessians<0, false, true>(scratch_data,
3463 eval0.template hessians<0, false, false>(scratch_data,
3467 eval1.template hessians<1, false, false>(hessians_quad +
3470 eval0.template values<0, false, true>(scratch_data,
3474 eval1.template values<1, false, false>(hessians_quad +
3477 eval0.template values<0, false, false>(scratch_data,
3482 eval1.template gradients<1, false, false>(hessians_quad +
3485 eval0.template gradients<0, false, true>(scratch_data,
3489 eval1.template values<1, false, false>(hessians_quad +
3493 eval0.template gradients<0, false, true>(scratch_data,
3497 eval0.template gradients<0, false, false>(scratch_data,
3502 eval1.template gradients<1, false, false>(hessians_quad +
3505 eval0.template values<0, false, true>(scratch_data,
3512 if (evaluation_flag &
3514 eval0.template hessians<0, false, true>(hessians_quad,
3517 eval0.template hessians<0, false, false>(hessians_quad,
3521 eval0.template values<0, false, false>(
3522 hessians_quad + n_q_points, values_dofs + 2 * dofs_stride);
3525 eval0.template gradients<0, false, true>(
3526 hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
3528 eval0.template gradients<0, false, false>(
3529 hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
3539 template <
int dim,
int fe_degree,
typename Number>
3542 template <
bool do_evaluate,
bool add_
into_output>
3547 const Number * input,
3549 const unsigned int face_no)
3551 Assert(
static_cast<unsigned int>(fe_degree) ==
3552 shape_info.
data.front().fe_degree ||
3556 interpolate_generic_raviart_thomas<do_evaluate, add_into_output>(
3557 n_components, input, output, flags, face_no, shape_info);
3559 interpolate_generic<do_evaluate, add_into_output>(
3565 shape_info.
data.front().fe_degree + 1,
3566 shape_info.
data.front().shape_data_on_face,
3574 template <
bool do_evaluate,
bool add_
into_output>
3577 const unsigned int n_components,
3580 const Number * input,
3582 const unsigned int face_no)
3584 Assert(
static_cast<unsigned int>(fe_degree + 1) ==
3585 shape_info.
data.front().n_q_points_1d ||
3589 interpolate_generic<do_evaluate, add_into_output>(
3595 shape_info.
data.front().quadrature.size(),
3596 shape_info.
data.front().quadrature_data_on_face,
3602 template <
bool do_evaluate,
bool add_
into_output,
int face_direction = 0>
3605 const Number * input,
3608 const unsigned int face_no,
3609 const unsigned int n_points_1d,
3611 const unsigned int dofs_per_component_on_cell,
3612 const unsigned int dofs_per_component_on_face)
3614 if (face_direction == face_no / 2)
3621 evalf(shape_data[face_no % 2],
3627 const unsigned int in_stride = do_evaluate ?
3628 dofs_per_component_on_cell :
3629 dofs_per_component_on_face;
3630 const unsigned int out_stride = do_evaluate ?
3631 dofs_per_component_on_face :
3632 dofs_per_component_on_cell;
3634 for (
unsigned int c = 0; c < n_components; ++c)
3637 evalf.template apply_face<face_direction,
3642 evalf.template apply_face<face_direction,
3647 evalf.template apply_face<face_direction,
3652 output += out_stride;
3655 else if (face_direction < dim)
3659 std::min(face_direction + 1, dim - 1)>(
3667 dofs_per_component_on_cell,
3668 dofs_per_component_on_face);
3672 template <
typename EvalType>
3676 const unsigned int face_no)
3683 template <
bool do_evaluate,
3684 bool add_into_output,
3685 int face_direction = 0,
3686 int max_derivative = 0>
3689 const unsigned int n_components,
3690 const Number * input,
3693 const unsigned int face_no,
3705 bool increase_max_der =
false;
3708 increase_max_der =
true;
3710 if (face_direction == face_no / 2 && !increase_max_der)
3716 shape_info, face_no, input, output);
3718 else if (face_direction == face_no / 2)
3725 n_components, input, output, flag, face_no, shape_info);
3727 else if (face_direction < dim)
3729 if (increase_max_der)
3734 std::min(face_direction + 1, dim - 1),
3736 n_components, input, output, flag, face_no, shape_info);
3745 n_components, input, output, flag, face_no, shape_info);
3751 template <
bool do_evaluate,
3752 bool add_into_output,
3758 const unsigned int face_no,
3759 const Number * input,
3765 using Evalf0 =
typename std::conditional<
3766 face_direction == 0,
3769 (fe_degree == -1) ? 1 : fe_degree + 1,
3775 (fe_degree == -1) ? 1 : fe_degree,
3779 using Evalf1 =
typename std::conditional<
3780 face_direction == 1,
3783 (fe_degree == -1) ? 1 : fe_degree + 1,
3789 (fe_degree == -1) ? 1 : fe_degree,
3793 using Evalf2 =
typename std::conditional<
3794 face_direction == 2,
3797 (fe_degree == -1) ? 1 : fe_degree + 1,
3803 (fe_degree == -1) ? 1 : fe_degree,
3809 create_evaluator_tensor_product<Evalf0>((face_direction == 0) ?
3810 shape_info.
data[0] :
3814 create_evaluator_tensor_product<Evalf1>((face_direction == 1) ?
3815 shape_info.
data[0] :
3819 create_evaluator_tensor_product<Evalf2>((face_direction == 2) ?
3820 shape_info.
data[0] :
3824 const unsigned int dofs_per_component_on_cell =
3826 const unsigned int dofs_per_component_on_face =
3832 const unsigned int in_stride =
3833 do_evaluate ? dofs_per_component_on_cell : dofs_per_component_on_face;
3834 const unsigned int out_stride =
3835 do_evaluate ? dofs_per_component_on_face : dofs_per_component_on_cell;
3837 const unsigned int in_stride_after_normal =
3839 dofs_per_component_on_cell :
3840 dofs_per_component_on_face - 3 *
Utilities::pow(fe_degree, dim - 2);
3841 const unsigned int out_stride_after_normal =
3843 dofs_per_component_on_face - 3 *
Utilities::pow(fe_degree, dim - 2) :
3844 dofs_per_component_on_cell;
3846 evalf0.template apply_face<face_direction,
3849 max_derivative>(input, output);
3851 input += (face_direction == 0) ? in_stride_after_normal : in_stride;
3852 output += (face_direction == 0) ? out_stride_after_normal : out_stride;
3854 evalf1.template apply_face<face_direction,
3857 max_derivative>(input, output);
3862 input += (face_direction == 1) ? in_stride_after_normal : in_stride;
3864 (face_direction == 1) ? out_stride_after_normal : out_stride;
3866 evalf2.template apply_face<face_direction,
3869 max_derivative>(input, output);
3877 template <
typename VectorizedArrayType,
typename Number2>
3881 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3882 dst[v] = src_ptr[v];
3889 template <
typename Number, std::
size_t w
idth>
3899 template <
typename VectorizedArrayType,
typename Number2>
3902 const unsigned int * indices,
3903 VectorizedArrayType &dst)
3905 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3906 dst[v] = src_ptr[indices[v]];
3913 template <
typename Number, std::
size_t w
idth>
3916 const unsigned int * indices,
3919 dst.
gather(src_ptr, indices);
3925 template <
typename VectorizedArrayType,
typename Number2>
3929 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3930 dst_ptr[v] += src[v];
3937 template <
typename Number, std::
size_t w
idth>
3943 (tmp + src).store(dst_ptr);
3949 template <
typename VectorizedArrayType,
typename Number2>
3952 const unsigned int * indices,
3955 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3956 dst_ptr[indices[v]] += src[v];
3963 template <
typename Number, std::
size_t w
idth>
3966 const unsigned int * indices,
3969 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
3970 for (
unsigned int v = 0; v < width; ++v)
3971 dst_ptr[indices[v]] += src[v];
3974 tmp.
gather(dst_ptr, indices);
3975 (tmp + src).
scatter(indices, dst_ptr);
3981 template <
typename Number>
3984 const unsigned int n_components,
3986 const unsigned int *orientation,
3987 const bool integrate,
3988 const std::size_t n_q_points,
3989 Number * tmp_values,
3990 Number * values_quad,
3991 Number * gradients_quad,
3992 Number * hessians_quad)
3994 for (
unsigned int c = 0; c < n_components; ++c)
3999 for (
unsigned int q = 0; q < n_q_points; ++q)
4000 tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
4002 for (
unsigned int q = 0; q < n_q_points; ++q)
4003 tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
4004 for (
unsigned int q = 0; q < n_q_points; ++q)
4005 values_quad[c * n_q_points + q] = tmp_values[q];
4008 for (
unsigned int d = 0;
d < dim; ++
d)
4011 for (
unsigned int q = 0; q < n_q_points; ++q)
4013 gradients_quad[(c * dim +
d) * n_q_points + orientation[q]];
4015 for (
unsigned int q = 0; q < n_q_points; ++q)
4016 tmp_values[orientation[q]] =
4017 gradients_quad[(c * dim +
d) * n_q_points + q];
4018 for (
unsigned int q = 0; q < n_q_points; ++q)
4019 gradients_quad[(c * dim +
d) * n_q_points + q] = tmp_values[q];
4023 const unsigned int hdim = (dim * (dim + 1)) / 2;
4024 for (
unsigned int d = 0;
d < hdim; ++
d)
4027 for (
unsigned int q = 0; q < n_q_points; ++q)
4028 tmp_values[q] = hessians_quad[(c * hdim +
d) * n_q_points +
4031 for (
unsigned int q = 0; q < n_q_points; ++q)
4032 tmp_values[orientation[q]] =
4033 hessians_quad[(c * hdim +
d) * n_q_points + q];
4034 for (
unsigned int q = 0; q < n_q_points; ++q)
4035 hessians_quad[(c * hdim +
d) * n_q_points + q] =
4044 template <
typename Number,
typename VectorizedArrayType>
4047 const unsigned int dim,
4048 const unsigned int n_components,
4049 const unsigned int v,
4051 const unsigned int * orientation,
4052 const bool integrate,
4053 const std::size_t n_q_points,
4054 Number * tmp_values,
4055 VectorizedArrayType * values_quad,
4056 VectorizedArrayType * gradients_quad =
nullptr,
4057 VectorizedArrayType * hessians_quad =
nullptr)
4059 for (
unsigned int c = 0; c < n_components; ++c)
4064 for (
unsigned int q = 0; q < n_q_points; ++q)
4065 tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v];
4067 for (
unsigned int q = 0; q < n_q_points; ++q)
4068 tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v];
4069 for (
unsigned int q = 0; q < n_q_points; ++q)
4070 values_quad[c * n_q_points + q][v] = tmp_values[q];
4073 for (
unsigned int d = 0;
d < dim; ++
d)
4077 for (
unsigned int q = 0; q < n_q_points; ++q)
4078 tmp_values[q] = gradients_quad[(c * dim +
d) * n_q_points +
4081 for (
unsigned int q = 0; q < n_q_points; ++q)
4082 tmp_values[orientation[q]] =
4083 gradients_quad[(c * dim +
d) * n_q_points + q][v];
4084 for (
unsigned int q = 0; q < n_q_points; ++q)
4085 gradients_quad[(c * dim +
d) * n_q_points + q][v] =
4091 const unsigned int hdim = (dim * (dim + 1)) / 2;
4092 for (
unsigned int d = 0;
d < hdim; ++
d)
4095 for (
unsigned int q = 0; q < n_q_points; ++q)
4096 tmp_values[q] = hessians_quad[(c * hdim +
d) * n_q_points +
4099 for (
unsigned int q = 0; q < n_q_points; ++q)
4100 tmp_values[orientation[q]] =
4101 hessians_quad[(c * hdim +
d) * n_q_points + q][v];
4102 for (
unsigned int q = 0; q < n_q_points; ++q)
4103 hessians_quad[(c * hdim +
d) * n_q_points + q][v] =
4112 template <
int dim,
typename Number>
4115 template <
int fe_degree,
int n_q_po
ints_1d>
4117 run(
const unsigned int n_components,
4119 const Number * values_dofs,
4123 const auto &shape_data = shape_info.
data.front();
4132 const unsigned int face_no = fe_eval.
get_face_no();
4134 const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
4135 const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
4142 const auto shape_values =
4143 &shape_data.shape_values_face(face_no, face_orientation, 0);
4146 auto values_dofs_actual_ptr = values_dofs;
4148 Eval eval(shape_values,
nullptr,
nullptr, n_dofs, n_q_points);
4149 for (
unsigned int c = 0; c < n_components; ++c)
4151 eval.template values<0, true, false>(values_dofs_actual_ptr,
4154 values_quad_ptr += n_q_points;
4155 values_dofs_actual_ptr += n_dofs;
4162 auto values_dofs_actual_ptr = values_dofs;
4164 std::array<const Number *, dim> shape_gradients;
4165 for (
unsigned int d = 0;
d < dim; ++
d)
4166 shape_gradients[
d] = &shape_data.shape_gradients_face(
4167 face_no, face_orientation,
d, 0);
4169 for (
unsigned int c = 0; c < n_components; ++c)
4171 for (
unsigned int d = 0;
d < dim; ++
d)
4179 eval.template gradients<0, true, false>(
4180 values_dofs_actual_ptr, gradients_quad_ptr);
4182 gradients_quad_ptr += n_q_points;
4184 values_dofs_actual_ptr += n_dofs;
4194 const unsigned int dofs_per_face =
4201 Number *scratch_data = temp + 3 * n_components * dofs_per_face;
4203 bool use_vectorization =
true;
4208 for (
unsigned int v = 0; v < Number::size(); ++v)
4211 use_vectorization =
false;
4213 if (use_vectorization ==
false)
4215 for (
unsigned int v = 0; v < Number::size(); ++v)
4222 for (
unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4229 template interpolate<true, false>(n_components,
4236 for (
unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4238 temp[i][v] = scratch_data[i][v];
4243 template interpolate<true, false>(n_components,
4251 constexpr
unsigned int n_q_points_1d_actual =
4252 fe_degree > -1 ? n_q_points_1d : 0;
4254 if (fe_degree >= 1 &&
4258 (fe_degree == -1) ? 1 : fe_degree,
4259 (n_q_points_1d < 1) ? 1 :
4262 template evaluate_or_integrate_in_face<false>(
4270 else if (fe_degree > -1 &&
4276 n_q_points_1d_actual,
4277 Number>::evaluate_in_face(n_components,
4291 n_q_points_1d_actual,
4292 Number>::evaluate_in_face(n_components,
4303 if (use_vectorization ==
false)
4305 for (
unsigned int v = 0; v < Number::size(); ++v)
4322 shape_info.n_q_points_face,
4337 shape_info.n_q_points_face,
4349 template <
int dim,
typename Number>
4352 template <
int fe_degree,
int n_q_po
ints_1d>
4354 run(
const unsigned int n_components,
4356 Number * values_dofs,
4360 const auto &shape_data = shape_info.
data.front();
4369 const unsigned int face_no = fe_eval.
get_face_no();
4371 const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
4372 const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
4379 const auto shape_values =
4380 &shape_data.shape_values_face(face_no, face_orientation, 0);
4383 auto values_dofs_actual_ptr = values_dofs;
4385 Eval eval(shape_values,
nullptr,
nullptr, n_dofs, n_q_points);
4386 for (
unsigned int c = 0; c < n_components; ++c)
4388 eval.template values<0, false, false>(values_quad_ptr,
4389 values_dofs_actual_ptr);
4391 values_quad_ptr += n_q_points;
4392 values_dofs_actual_ptr += n_dofs;
4399 auto values_dofs_actual_ptr = values_dofs;
4401 std::array<const Number *, dim> shape_gradients;
4402 for (
unsigned int d = 0;
d < dim; ++
d)
4403 shape_gradients[
d] = &shape_data.shape_gradients_face(
4404 face_no, face_orientation,
d, 0);
4406 for (
unsigned int c = 0; c < n_components; ++c)
4408 for (
unsigned int d = 0;
d < dim; ++
d)
4418 eval.template gradients<0, false, false>(
4419 gradients_quad_ptr, values_dofs_actual_ptr);
4421 eval.template gradients<0, false, true>(
4422 gradients_quad_ptr, values_dofs_actual_ptr);
4424 gradients_quad_ptr += n_q_points;
4426 values_dofs_actual_ptr += n_dofs;
4436 const unsigned int dofs_per_face =
4441 Number *scratch_data = temp + 3 * n_components * dofs_per_face;
4443 bool use_vectorization =
true;
4452 [&](
const auto &v) {
4453 return v == fe_eval.get_cell_ids()[0] ||
4454 v == numbers::invalid_unsigned_int;
4457 if (use_vectorization ==
false)
4459 for (
unsigned int v = 0; v < Number::size(); ++v)
4476 shape_info.n_q_points_face,
4491 shape_info.n_q_points_face,
4497 const unsigned int n_q_points_1d_actual =
4498 fe_degree > -1 ? n_q_points_1d : 0;
4501 if (fe_degree >= 1 &&
4505 (fe_degree == -1) ? 1 : fe_degree,
4506 (n_q_points_1d < 1) ? 1 :
4509 template evaluate_or_integrate_in_face<true>(integration_flag,
4516 else if (fe_degree > -1 &&
4524 n_q_points_1d_actual,
4525 Number>::integrate_in_face(n_components,
4539 n_q_points_1d_actual,
4540 Number>::integrate_in_face(n_components,
4550 if (use_vectorization ==
false)
4552 for (
unsigned int v = 0; v < Number::size(); ++v)
4561 template interpolate<false, false>(n_components,
4568 for (
unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4570 temp[i][v] = scratch_data[i][v];
4575 template interpolate<false, false>(n_components,
4587 template <
int n_face_orientations,
4589 typename EvaluationData,
4590 const bool check_face_orientations =
false>
4594 const unsigned int n_components,
4596 typename Processor::Number2_ * global_vector_ptr,
4598 const EvaluationData & fe_eval,
4599 typename Processor::VectorizedArrayType_ * temp1)
4601 constexpr
int dim = Processor::dim_;
4602 constexpr
int fe_degree = Processor::fe_degree_;
4603 using VectorizedArrayType =
typename Processor::VectorizedArrayType_;
4604 constexpr
int n_lanes = VectorizedArrayType::size();
4606 using Number =
typename Processor::Number_;
4607 using Number2_ =
typename Processor::Number2_;
4609 const auto & shape_data = fe_eval.get_shape_info().data.front();
4610 constexpr
bool integrate = Processor::do_integrate;
4611 const unsigned int face_no = fe_eval.get_face_no();
4612 const auto & dof_info = fe_eval.get_dof_info();
4613 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
4615 fe_eval.get_dof_access_index();
4617 dof_info.index_storage_variants[dof_access_index].size());
4618 constexpr
unsigned int dofs_per_face =
4620 const unsigned int subface_index = fe_eval.get_subface_index();
4622 const unsigned int n_filled_lanes =
4623 dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
4625 bool all_faces_are_same = n_filled_lanes == n_lanes;
4626 if (n_face_orientations == n_lanes)
4627 for (
unsigned int v = 1; v < n_lanes; ++v)
4628 if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) ||
4629 fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0))
4631 all_faces_are_same =
false;
4636 std::array<const unsigned int *, n_face_orientations> orientation = {};
4638 if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same &&
4639 fe_eval.is_interior_face() == 0)
4640 for (
unsigned int v = 0; v < n_lanes; ++v)
4648 if (shape_data.nodal_at_cell_boundaries &&
4649 fe_eval.get_face_orientation(v) != 0)
4654 check_face_orientations ==
false)
4668 orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
4669 fe_eval.get_face_orientation(v), 0);
4672 else if (dim == 3 && fe_eval.get_face_orientation() != 0)
4676 check_face_orientations ==
false)
4690 for (
unsigned int v = 0; v < n_face_orientations; ++v)
4691 orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
4692 fe_eval.get_face_orientation(), 0);
4698 VectorizedArrayType grad_weight =
4700 .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) :
4701 (1 + face_no % 2))];
4704 std::array<const unsigned int *, n_face_orientations> index_array_hermite =
4708 if (n_face_orientations == 1)
4709 index_array_hermite[0] =
4710 &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0);
4713 for (
unsigned int v = 0; v < n_lanes; ++v)
4718 const auto face_no = fe_eval.get_face_no(v);
4721 shape_data.shape_data_on_face[0][fe_degree +
4723 (2 - (face_no % 2)) :
4724 (1 + (face_no % 2)))][0];
4726 index_array_hermite[v] =
4727 &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no,
4734 std::array<const unsigned int *, n_face_orientations> index_array_nodal =
4736 if (shape_data.nodal_at_cell_boundaries ==
true)
4738 if (n_face_orientations == 1)
4739 index_array_nodal[0] =
4740 &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0);
4743 for (
unsigned int v = 0; v < n_lanes; ++v)
4748 const auto face_no = fe_eval.get_face_no(v);
4750 index_array_nodal[v] =
4751 &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no,
4758 const auto reorientate = [&](
const unsigned int v,
const unsigned int i) {
4759 return (!check_face_orientations || orientation[v] ==
nullptr) ?
4766 fe_eval.get_cell_ids()[0] :
4768 const unsigned int *dof_indices =
4769 &dof_info.dof_indices_contiguous[dof_access_index][
cell_index];
4771 for (
unsigned int comp = 0; comp < n_components; ++comp)
4773 const std::size_t index_offset =
4774 dof_info.component_dof_indices_offset
4775 [fe_eval.get_active_fe_index()]
4776 [fe_eval.get_first_selected_component()] +
4780 if (n_face_orientations == 1 &&
4781 dof_info.index_storage_variants[dof_access_index][cell] ==
4783 interleaved_contiguous)
4786 dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
4788 Number2_ *vector_ptr =
4789 global_vector_ptr + dof_indices[0] + index_offset * n_lanes;
4793 for (
unsigned int i = 0; i < dofs_per_face; ++i)
4797 const unsigned int ind1 = index_array_hermite[0][2 * i];
4798 const unsigned int ind2 = index_array_hermite[0][2 * i + 1];
4799 const unsigned int i_ = reorientate(0, i);
4800 proc.hermite_grad_vectorized(temp1[i_],
4801 temp1[i_ + dofs_per_face],
4802 vector_ptr + ind1 * n_lanes,
4803 vector_ptr + ind2 * n_lanes,
4809 for (
unsigned int i = 0; i < dofs_per_face; ++i)
4813 const unsigned int i_ = reorientate(0, i);
4814 const unsigned int ind = index_array_nodal[0][i];
4815 proc.value_vectorized(temp1[i_],
4816 vector_ptr + ind * n_lanes);
4822 else if (n_face_orientations == 1 &&
4823 dof_info.index_storage_variants[dof_access_index][cell] ==
4825 interleaved_contiguous_strided)
4828 dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
4830 Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes;
4833 for (
unsigned int i = 0; i < dofs_per_face; ++i)
4837 const unsigned int i_ = reorientate(0, i);
4838 const unsigned int ind1 =
4839 index_array_hermite[0][2 * i] * n_lanes;
4840 const unsigned int ind2 =
4841 index_array_hermite[0][2 * i + 1] * n_lanes;
4842 proc.hermite_grad_vectorized_indexed(
4844 temp1[i_ + dofs_per_face],
4854 for (
unsigned int i = 0; i < dofs_per_face; ++i)
4858 const unsigned int i_ = reorientate(0, i);
4859 const unsigned int ind = index_array_nodal[0][i] * n_lanes;
4860 proc.value_vectorized_indexed(temp1[i_],
4868 else if (n_face_orientations == 1 &&
4869 dof_info.index_storage_variants[dof_access_index][cell] ==
4871 interleaved_contiguous_mixed_strides)
4873 const unsigned int *strides =
4874 &dof_info.dof_indices_interleave_strides[dof_access_index]
4876 unsigned int indices[n_lanes];
4877 for (
unsigned int v = 0; v < n_lanes; ++v)
4878 indices[v] = dof_indices[v] + index_offset * strides[v];
4879 const unsigned int n_filled_lanes =
4880 dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
4884 if (n_filled_lanes == n_lanes)
4885 for (
unsigned int i = 0; i < dofs_per_face; ++i)
4889 const unsigned int i_ = reorientate(0, i);
4890 unsigned int ind1[n_lanes];
4892 for (
unsigned int v = 0; v < n_lanes; ++v)
4893 ind1[v] = indices[v] +
4894 index_array_hermite[0][2 * i] * strides[v];
4895 unsigned int ind2[n_lanes];
4897 for (
unsigned int v = 0; v < n_lanes; ++v)
4901 index_array_hermite[0][2 * i + 1] * strides[v];
4902 proc.hermite_grad_vectorized_indexed(
4904 temp1[i_ + dofs_per_face],
4913 if (integrate ==
false)
4914 for (
unsigned int i = 0; i < 2 * dofs_per_face; ++i)
4915 temp1[i] = VectorizedArrayType();
4917 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4918 for (
unsigned int i = 0; i < dofs_per_face; ++i)
4920 const unsigned int i_ =
4921 reorientate(n_face_orientations == 1 ? 0 : v, i);
4924 temp1[i_ + dofs_per_face][v],
4928 [n_face_orientations == 1 ? 0 : v][2 * i] *
4932 index_array_hermite[n_face_orientations == 1 ?
4936 grad_weight[n_face_orientations == 1 ? 0 : v]);
4942 if (n_filled_lanes == n_lanes)
4943 for (
unsigned int i = 0; i < dofs_per_face; ++i)
4946 unsigned int ind[n_lanes];
4948 for (
unsigned int v = 0; v < n_lanes; ++v)
4950 indices[v] + index_array_nodal[0][i] * strides[v];
4951 const unsigned int i_ = reorientate(0, i);
4952 proc.value_vectorized_indexed(temp1[i_],
4958 if (integrate ==
false)
4959 for (
unsigned int i = 0; i < dofs_per_face; ++i)
4960 temp1[i] = VectorizedArrayType();
4962 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4963 for (
unsigned int i = 0; i < dofs_per_face; ++i)
4965 temp1[reorientate(n_face_orientations == 1 ? 0 : v,
4969 index_array_nodal[n_face_orientations == 1 ? 0 : v]
4977 else if (n_face_orientations > 1 ||
4978 dof_info.index_storage_variants[dof_access_index][cell] ==
4982 Number2_ *vector_ptr = global_vector_ptr + index_offset;
4984 const bool vectorization_possible =
4985 all_faces_are_same && (sm_ptr ==
nullptr);
4987 std::array<Number2_ *, n_lanes> vector_ptrs;
4988 std::array<unsigned int, n_lanes> reordered_indices;
4990 if (vectorization_possible ==
false)
4993 if (n_face_orientations == 1)
4995 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4996 if (sm_ptr ==
nullptr)
4998 vector_ptrs[v] = vector_ptr + dof_indices[v];
5004 .dof_indices_contiguous_sm[dof_access_index]
5005 [cell * n_lanes + v];
5006 vector_ptrs[v] =
const_cast<Number2_ *
>(
5007 sm_ptr->operator[](temp.first).data() +
5008 temp.second + index_offset);
5011 else if (n_face_orientations == n_lanes)
5013 const auto &cells = fe_eval.get_cell_ids();
5014 for (
unsigned int v = 0; v < n_lanes; ++v)
5017 if (sm_ptr ==
nullptr)
5022 .dof_indices_contiguous[dof_access_index]
5029 .dof_indices_contiguous_sm[dof_access_index]
5031 vector_ptrs[v] =
const_cast<Number2_ *
>(
5032 sm_ptr->operator[](temp.first).data() +
5033 temp.second + index_offset);
5042 else if (n_face_orientations == n_lanes)
5044 for (
unsigned int v = 0; v < n_lanes; ++v)
5045 reordered_indices[v] =
5046 dof_info.dof_indices_contiguous[dof_access_index]
5047 [fe_eval.get_cell_ids()[v]];
5048 dof_indices = reordered_indices.data();
5053 if (vectorization_possible)
5054 for (
unsigned int i = 0; i < dofs_per_face; ++i)
5056 const unsigned int ind1 = index_array_hermite[0][2 * i];
5057 const unsigned int ind2 =
5058 index_array_hermite[0][2 * i + 1];
5059 const unsigned int i_ = reorientate(0, i);
5061 proc.hermite_grad_vectorized_indexed(
5063 temp1[i_ + dofs_per_face],
5070 else if (n_face_orientations == 1)
5071 for (
unsigned int i = 0; i < dofs_per_face; ++i)
5073 const unsigned int ind1 = index_array_hermite[0][2 * i];
5074 const unsigned int ind2 =
5075 index_array_hermite[0][2 * i + 1];
5076 const unsigned int i_ = reorientate(0, i);
5078 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
5079 proc.hermite_grad(temp1[i_][v],
5080 temp1[i_ + dofs_per_face][v],
5081 vector_ptrs[v][ind1],
5082 vector_ptrs[v][ind2],
5085 if (integrate ==
false)
5086 for (
unsigned int v = n_filled_lanes; v < n_lanes; ++v)
5089 temp1[i + dofs_per_face][v] = 0.0;
5094 if (integrate ==
false && n_filled_lanes < n_lanes)
5095 for (
unsigned int i = 0; i < dofs_per_face; ++i)
5096 temp1[i] = temp1[i + dofs_per_face] = Number();
5098 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
5099 for (
unsigned int i = 0; i < dofs_per_face; ++i)
5101 temp1[reorientate(v, i)][v],
5102 temp1[reorientate(v, i) + dofs_per_face][v],
5103 vector_ptrs[v][index_array_hermite[v][2 * i]],
5104 vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
5110 if (vectorization_possible)
5111 for (
unsigned int i = 0; i < dofs_per_face; ++i)
5113 const unsigned int ind = index_array_nodal[0][i];
5114 const unsigned int i_ = reorientate(0, i);
5116 proc.value_vectorized_indexed(temp1[i_],
5120 else if (n_face_orientations == 1)
5121 for (
unsigned int i = 0; i < dofs_per_face; ++i)
5123 const unsigned int ind = index_array_nodal[0][i];
5124 const unsigned int i_ = reorientate(0, i);
5126 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
5127 proc.value(temp1[i_][v], vector_ptrs[v][ind]);
5129 if (integrate ==
false)
5130 for (
unsigned int v = n_filled_lanes; v < n_lanes; ++v)
5135 if (integrate ==
false && n_filled_lanes < n_lanes)
5136 for (
unsigned int i = 0; i < dofs_per_face; ++i)
5137 temp1[i] = Number();
5139 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
5140 for (
unsigned int i = 0; i < dofs_per_face; ++i)
5141 proc.value(temp1[reorientate(v, i)][v],
5142 vector_ptrs[v][index_array_nodal[v][i]]);
5152 temp1 += 3 * dofs_per_face;
5158 template <
int dim,
typename Number2,
typename VectorizedArrayType>
5161 using Number =
typename VectorizedArrayType::value_type;
5163 template <
int fe_degree,
int n_q_po
ints_1d>
5165 run(
const unsigned int n_components,
5167 const Number2 * src_ptr,
5176 const unsigned int dofs_per_face =
Utilities::pow(fe_degree + 1, dim - 1);
5179 VectorizedArrayType *scratch_data =
5180 temp + 3 * n_components * dofs_per_face;
5187 fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
5188 p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
5190 fe_face_evaluation_process_and_io<1>(
5191 p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
5200 VectorizedArrayType>::
5201 evaluate_in_face(n_components,
5215 VectorizedArrayType>::
5216 evaluate_in_face(n_components,
5233 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5280 const Number2 * vector_ptr,
5283 const unsigned int fe_degree = shape_info.
data.front().fe_degree;
5284 if (fe_degree < 1 || !shape_info.
data.front().nodal_at_cell_boundaries ||
5287 shape_info.
data.front().element_type !=
5290 vector_ptr ==
nullptr ||
5291 shape_info.
data.front().element_type >
5301 template <
int fe_degree>
5311 template <
typename T0,
typename T1,
typename T2>
5317 const T2 &grad_weight)
5321 temp_2 = grad_weight * (temp_1 - temp_2);
5324 template <
typename T1,
typename T2>
5331 template <
typename T0,
typename T1,
typename T2,
typename T3>
5337 const T2 &grad_weight,
5338 const T3 &indices_1,
5339 const T3 &indices_2)
5343 temp_2 = grad_weight * (temp_1 - temp_2);
5346 template <
typename T0,
typename T1,
typename T2>
5353 template <
typename T0,
typename T1,
typename T2>
5357 const T1 &src_ptr_1,
5358 const T1 &src_ptr_2,
5359 const T2 &grad_weight)
5363 temp_2 = grad_weight * (temp_1 - src_ptr_2);
5366 template <
typename T1,
typename T2>
5378 template <
int dim,
typename Number2,
typename VectorizedArrayType>
5381 using Number =
typename VectorizedArrayType::value_type;
5383 template <
int fe_degree,
int n_q_po
ints_1d>
5385 run(
const unsigned int n_components,
5396 const unsigned int dofs_per_face =
Utilities::pow(fe_degree + 1, dim - 1);
5399 VectorizedArrayType *scratch_data =
5400 temp + 3 * n_components * dofs_per_face;
5410 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5454 VectorizedArrayType>::
5455 integrate_in_face(n_components,
5469 VectorizedArrayType>::
5470 integrate_in_face(n_components,
5485 fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
5486 p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
5488 fe_face_evaluation_process_and_io<1>(
5489 p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
5495 template <
int fe_degree>
5505 template <
typename T0,
typename T1,
typename T2,
typename T3,
typename T4>
5511 const T4 &grad_weight)
5514 const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
5515 const VectorizedArrayType grad = grad_weight * temp_2;
5520 template <
typename T0,
typename T1>
5528 template <
typename T0,
typename T1,
typename T2,
typename T3>
5534 const T2 &grad_weight,
5535 const T3 &indices_1,
5536 const T3 &indices_2)
5539 const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
5540 const VectorizedArrayType grad = grad_weight * temp_2;
5545 template <
typename T0,
typename T1,
typename T2>
5553 template <
typename T0,
typename T1,
typename T2>
5559 const T2 &grad_weight)
5562 const Number val = temp_1 - grad_weight * temp_2;
5563 const Number grad = grad_weight * temp_2;
5568 template <
typename T0,
typename T1>
5584 template <
int dim,
typename Number>
5587 template <
int fe_degree,
int = 0>
5589 run(
const unsigned int n_components,
5591 const Number * in_array,
5593 std::enable_if_t<fe_degree != -1> * =
nullptr)
5595 constexpr
unsigned int dofs_per_component =
5613 for (
unsigned int d = 0;
d < n_components; ++
d)
5615 const Number *in = in_array +
d * dofs_per_component;
5616 Number * out = out_array +
d * dofs_per_component;
5619 evaluator.template hessians<0, true, false>(in, out);
5621 evaluator.template hessians<1, true, false>(out, out);
5623 evaluator.template hessians<2, true, false>(out, out);
5625 for (
unsigned int q = 0; q < dofs_per_component; ++q)
5627 const Number inverse_JxW_q = Number(1.) / fe_eval.
JxW(q);
5628 for (
unsigned int d = 0;
d < n_components; ++
d)
5629 out_array[q +
d * dofs_per_component] *= inverse_JxW_q;
5631 for (
unsigned int d = 0;
d < n_components; ++
d)
5633 Number *out = out_array +
d * dofs_per_component;
5635 evaluator.template hessians<2, false, false>(out, out);
5637 evaluator.template hessians<1, false, false>(out, out);
5638 evaluator.template hessians<0, false, false>(out, out);
5643 template <
int fe_degree,
int = 0>
5645 run(
const unsigned int n_components,
5647 const Number * in_array,
5649 std::enable_if_t<fe_degree == -1> * =
nullptr)
5651 static_assert(fe_degree == -1,
"Only usable for degree -1");
5652 const unsigned int dofs_per_component =
5664 for (
unsigned int d = 0;
d < n_components; ++
d)
5666 const Number *in = in_array +
d * dofs_per_component;
5667 Number * out = out_array +
d * dofs_per_component;
5670 evaluator.template values<0, true, false>(in, out);
5672 evaluator.template values<1, true, false>(out, out);
5674 evaluator.template values<2, true, false>(out, out);
5676 for (
unsigned int q = 0; q < dofs_per_component; ++q)
5678 const Number inverse_JxW_q = Number(1.) / fe_eval.
JxW(q);
5679 for (
unsigned int d = 0;
d < n_components; ++
d)
5680 out_array[q +
d * dofs_per_component] *= inverse_JxW_q;
5682 for (
unsigned int d = 0;
d < n_components; ++
d)
5684 Number *out = out_array +
d * dofs_per_component;
5686 evaluator.template values<2, false, false>(out, out);
5688 evaluator.template values<1, false, false>(out, out);
5689 evaluator.template values<0, false, false>(out, out);
5703 template <
int dim,
typename Number>
5706 template <
int fe_degree,
int = 0>
5708 run(
const unsigned int n_desired_components,
5711 const bool dyadic_coefficients,
5712 const Number * in_array,
5715 const unsigned int given_degree =
5716 (fe_degree > -1) ? fe_degree :
5719 const unsigned int dofs_per_component =
5723 inverse_coefficients.
size() % dofs_per_component == 0,
5725 "Expected diagonal to be a multiple of scalar dof per cells"));
5727 if (!dyadic_coefficients)
5729 if (inverse_coefficients.
size() != dofs_per_component)
5731 inverse_coefficients.
size());
5737 inverse_coefficients.
size());
5756 const Number *in = in_array;
5757 Number * out = out_array;
5759 const Number *inv_coefficient = inverse_coefficients.
data();
5761 const unsigned int shift_coefficient =
5762 inverse_coefficients.
size() > dofs_per_component ? dofs_per_component :
5765 const auto n_comp_outer = dyadic_coefficients ? 1 : n_desired_components;
5766 const auto n_comp_inner = dyadic_coefficients ? n_desired_components : 1;
5768 for (
unsigned int d = 0;
d < n_comp_outer; ++
d)
5770 for (
unsigned int di = 0; di < n_comp_inner; ++di)
5772 const Number *in_ = in + di * dofs_per_component;
5773 Number * out_ = out + di * dofs_per_component;
5774 evaluator.template hessians<0, true, false>(in_, out_);
5776 evaluator.template hessians<1, true, false>(out_, out_);
5778 evaluator.template hessians<2, true, false>(out_, out_);
5780 if (dyadic_coefficients)
5782 const auto n_coeff_components =
5783 n_desired_components * n_desired_components;
5784 if (n_desired_components == dim)
5786 for (
unsigned int q = 0; q < dofs_per_component; ++q)
5787 vmult<dim>(&inv_coefficient[q * n_coeff_components],
5790 dofs_per_component);
5794 for (
unsigned int q = 0; q < dofs_per_component; ++q)
5795 vmult<-1>(&inv_coefficient[q * n_coeff_components],
5799 n_desired_components);
5803 for (
unsigned int q = 0; q < dofs_per_component; ++q)
5804 out[q] *= inv_coefficient[q];
5806 for (
unsigned int di = 0; di < n_comp_inner; ++di)
5808 Number *out_ = out + di * dofs_per_component;
5810 evaluator.template hessians<2, false, false>(out_, out_);
5812 evaluator.template hessians<1, false, false>(out_, out_);
5813 evaluator.template hessians<0, false, false>(out_, out_);
5816 in += dofs_per_component;
5817 out += dofs_per_component;
5818 inv_coefficient += shift_coefficient;
5825 template <
int n_components>
5827 vmult(
const Number * inverse_coefficients,
5830 const unsigned int dofs_per_component,
5831 const unsigned int n_given_components = 0)
5833 const unsigned int n_desired_components =
5834 (n_components > -1) ? n_components : n_given_components;
5836 std::array<Number, dim + 2> tmp = {};
5837 Assert(n_desired_components <= dim + 2,
5839 "Number of components larger than dim+2 not supported."));
5841 for (
unsigned int d = 0;
d < n_desired_components; ++
d)
5842 tmp[
d] = src[
d * dofs_per_component];
5844 for (
unsigned int d1 = 0; d1 < n_desired_components; ++d1)
5846 const Number *inv_coeff_row =
5847 &inverse_coefficients[d1 * n_desired_components];
5848 Number
sum = inv_coeff_row[0] * tmp[0];
5849 for (
unsigned int d2 = 1; d2 < n_desired_components; ++d2)
5850 sum += inv_coeff_row[d2] * tmp[d2];
5851 dst[d1 * dofs_per_component] =
sum;
5864 template <
int dim,
typename Number>
5867 template <
int fe_degree,
int n_q_po
ints_1d>
5869 run(
const unsigned int n_desired_components,
5871 const Number * in_array,
5874 static const bool do_inplace =
5875 fe_degree > -1 && (fe_degree + 1 == n_q_points_1d);
5881 const auto &inverse_shape =
5886 const std::size_t dofs_per_component =
5889 const std::size_t n_q_points = do_inplace ?
5904 for (
unsigned int d = 0;
d < n_desired_components; ++
d)
5906 const Number *in = in_array +
d * n_q_points;
5907 Number * out = out_array +
d * dofs_per_component;
5910 auto temp_2 = do_inplace ?
5912 (temp_1 +
std::max(n_q_points, dofs_per_component));
5916 evaluator.template hessians<2, false, false>(in, temp_1);
5917 evaluator.template hessians<1, false, false>(temp_1, temp_2);
5918 evaluator.template hessians<0, false, false>(temp_2, out);
5922 evaluator.template hessians<1, false, false>(in, temp_1);
5923 evaluator.template hessians<0, false, false>(temp_1, out);
5926 evaluator.template hessians<0, false, false>(in, out);
value_type * data() const noexcept
std::uint8_t get_face_no(const unsigned int v=0) const
internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex get_dof_access_index() const
Number JxW(const unsigned int q_point) const
const Number * begin_values() const
const Number * begin_hessians() const
const std::array< unsigned int, n_lanes > & get_cell_ids() const
unsigned int get_subface_index() const
const Number * begin_gradients() const
bool is_interior_face() const
ArrayView< Number > get_scratch_data() const
const ShapeInfoType & get_shape_info() const
std::uint8_t get_face_orientation(const unsigned int v=0) const
void gather(const Number *base_ptr, const unsigned int *offsets)
void load(const OtherNumber *ptr)
#define DEAL_II_ALWAYS_INLINE
#define DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcInternalError()
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)
#define AssertThrow(cond, exc)
@ tensor_symmetric_no_collocation
@ tensor_symmetric_collocation
@ tensor_symmetric_plus_dg0
@ tensor_symmetric_hermite
EvaluationFlags
The EvaluationFlags enum.
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm &mpi_communicator)
T scatter(const MPI_Comm &comm, const std::vector< T > &objects_to_send, const unsigned int root_process=0)
constexpr T pow(const T base, const int iexp)
void do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
void adjust_for_face_orientation_per_lane(const unsigned int dim, const unsigned int n_components, const unsigned int v, const EvaluationFlags::EvaluationFlags flag, const unsigned int *orientation, const bool integrate, const std::size_t n_q_points, Number *tmp_values, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad=nullptr, VectorizedArrayType *hessians_quad=nullptr)
void do_vectorized_scatter_add(const VectorizedArrayType src, const unsigned int *indices, Number2 *dst_ptr)
void do_vectorized_gather(const Number2 *src_ptr, const unsigned int *indices, VectorizedArrayType &dst)
void do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
void fe_face_evaluation_process_and_io(Processor &proc, const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, typename Processor::Number2_ *global_vector_ptr, const std::vector< ArrayView< const typename Processor::Number2_ >> *sm_ptr, const EvaluationData &fe_eval, typename Processor::VectorizedArrayType_ *temp1)
void adjust_for_face_orientation(const unsigned int dim, const unsigned int n_components, const EvaluationFlags::EvaluationFlags flag, const unsigned int *orientation, const bool integrate, const std::size_t n_q_points, Number *tmp_values, Number *values_quad, Number *gradients_quad, Number *hessians_quad)
@ evaluate_raviart_thomas
static const unsigned int invalid_unsigned_int
static bool run(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array, std::enable_if_t< fe_degree==-1 > *=nullptr)
static bool run(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array, std::enable_if_t< fe_degree !=-1 > *=nullptr)
static bool run(const unsigned int n_desired_components, const FEEvaluationData< dim, Number, false > &fe_eval, const ArrayView< const Number > &inverse_coefficients, const bool dyadic_coefficients, const Number *in_array, Number *out_array)
static void vmult(const Number *inverse_coefficients, const Number *src, Number *dst, const unsigned int dofs_per_component, const unsigned int n_given_components=0)
static void do_mass(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const AlignedVector< Number > &coefficients, const Number *values_in, Number *scratch_data, Number *values_out)
static void do_backward(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const bool add_into_result, Number *values_in, Number *values_out, const unsigned int basis_size_1_variable=numbers::invalid_unsigned_int, const unsigned int basis_size_2_variable=numbers::invalid_unsigned_int)
static void do_forward(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const Number *values_in, Number *values_out, const unsigned int basis_size_1_variable=numbers::invalid_unsigned_int, const unsigned int basis_size_2_variable=numbers::invalid_unsigned_int)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
static void do_integrate(const MatrixFreeFunctions::UnivariateShapeData< Number > &shape, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, Number *gradients_quad, const Number *hessians_quad, const bool add_into_values_array)
static void do_evaluate(const MatrixFreeFunctions::UnivariateShapeData< Number > &shape, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, Number *gradients_quad, Number *hessians_quad)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array)
static EvalType create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number > &shape_data)