16#ifndef dealii_matrix_free_tensor_product_kernels_h
17#define dealii_matrix_free_tensor_product_kernels_h
125 static_assert(n_rows > 0 || n_columns > 0,
126 "Specialization only for n_rows, n_columns > 0");
127 Assert(n_rows > 0 && n_columns > 0,
129 std::to_string(n_rows) +
", " +
130 std::to_string(n_columns) +
" was passed!"));
132 "This function should only use EvaluatorQuantity::value");
137 std::array<Number, mm>
x;
138 for (
int i = 0; i <
mm; ++i)
140 for (
int col = 0; col <
nn; ++col)
145 res0 = matrix[col] *
x[0];
146 for (
int i = 1; i <
mm; ++i)
148 const Number2
mji = matrix[i * n_columns + col];
153 mji.imag() *
x[i].imag());
155 mji.real() *
x[i].imag());
163 res0 = matrix[col * n_columns] *
x[0];
164 for (
int i = 1; i <
mm; ++i)
166 const Number2
mij = matrix[col * n_columns + i];
171 mij.imag() *
x[i].imag());
173 mij.real() *
x[i].imag());
199 int n_components = 1>
211 Assert(n_rows > 0 && n_columns > 0,
213 Assert(n_rows > 0 && n_columns > 0,
215 std::to_string(n_rows) +
", " +
216 std::to_string(n_columns) +
" was passed!"));
219 "This function should only use EvaluatorQuantity::value");
226 static_assert(n_components > 0 && n_components < 4,
227 "Invalid number of components");
236 const Number2 *
matrix_1 = matrix + n_columns;
238 for (
int col = 0; col <
nn; ++col)
249 const Number2 *
matrix_1 = matrix + n_columns;
252 for (
int col = 0; col <
nn; ++col)
267 std::array<Number, 129>
x;
268 for (
int i = 0; i <
mm; ++i)
271 for (
int col = 0; col <
nn; ++col)
276 res0 = matrix[col] *
x[0];
277 for (
int i = 1; i <
mm; ++i)
278 res0 += matrix[i * n_columns + col] *
x[i];
282 res0 = matrix[col * n_columns] *
x[0];
283 for (
int i = 1; i <
mm; ++i)
284 res0 += matrix[col * n_columns + i] *
x[i];
294 const Number *
in0 = in;
295 const Number *
in1 = n_components > 1 ? in +
mm :
nullptr;
296 const Number *
in2 = n_components > 2 ? in + 2 *
mm :
nullptr;
299 Number *
out1 = n_components > 1 ? out +
nn :
nullptr;
300 Number *
out2 = n_components > 2 ? out + 2 *
nn :
nullptr;
308 const Number2 *matrix_ptr = matrix + col;
309 const Number a =
in0[0];
310 res[0] = matrix_ptr[0] * a;
311 res[1] = matrix_ptr[1] * a;
312 res[2] = matrix_ptr[2] * a;
313 res[3] = matrix_ptr[3] * a;
315 if (n_components > 1)
317 const Number b =
in1[0];
318 res[4] = matrix_ptr[0] * b;
319 res[5] = matrix_ptr[1] * b;
320 res[6] = matrix_ptr[2] * b;
321 res[7] = matrix_ptr[3] * b;
324 if (n_components > 2)
326 const Number c =
in2[0];
327 res[8] = matrix_ptr[0] * c;
328 res[9] = matrix_ptr[1] * c;
329 res[10] = matrix_ptr[2] * c;
330 res[11] = matrix_ptr[3] * c;
333 matrix_ptr += n_columns;
334 for (
int i = 1; i <
mm; ++i, matrix_ptr += n_columns)
337 res[0] += matrix_ptr[0] * a;
338 res[1] += matrix_ptr[1] * a;
339 res[2] += matrix_ptr[2] * a;
340 res[3] += matrix_ptr[3] * a;
342 if (n_components > 1)
345 res[4] += matrix_ptr[0] * b;
346 res[5] += matrix_ptr[1] * b;
347 res[6] += matrix_ptr[2] * b;
348 res[7] += matrix_ptr[3] * b;
350 if (n_components > 2)
353 res[8] += matrix_ptr[0] * c;
354 res[9] += matrix_ptr[1] * c;
355 res[10] += matrix_ptr[2] * c;
356 res[11] += matrix_ptr[3] * c;
362 const Number2 *
matrix_0 = matrix + col * n_columns;
363 const Number2 *
matrix_1 = matrix + (col + 1) * n_columns;
364 const Number2 *
matrix_2 = matrix + (col + 2) * n_columns;
365 const Number2 *
matrix_3 = matrix + (col + 3) * n_columns;
367 const Number a =
in0[0];
373 if (n_components > 1)
375 const Number b =
in1[0];
382 if (n_components > 2)
384 const Number c =
in2[0];
391 for (
int i = 1; i <
mm; ++i)
399 if (n_components > 1)
408 if (n_components > 2)
424 if (n_components > 1)
431 if (n_components > 2)
445 if (n_components > 1)
452 if (n_components > 2)
461 if (n_components > 1)
463 if (n_components > 2)
471 const Number2 *matrix_ptr = matrix +
nn_regular;
474 res2 = matrix_ptr[2] *
in0[0];
475 if (n_components > 1)
481 if (n_components > 2)
487 matrix_ptr += n_columns;
488 for (
int i = 1; i <
mm; ++i, matrix_ptr += n_columns)
493 if (n_components > 1)
499 if (n_components > 2)
516 if (n_components > 1)
522 if (n_components > 2)
528 for (
int i = 1; i <
mm; ++i)
533 if (n_components > 1)
539 if (n_components > 2)
552 if (n_components > 1)
558 if (n_components > 2)
570 if (n_components > 1)
576 if (n_components > 2)
589 const Number2 *matrix_ptr = matrix +
nn_regular;
592 if (n_components > 1)
594 res2 = matrix_ptr[0] *
in1[0];
597 if (n_components > 2)
602 matrix_ptr += n_columns;
603 for (
int i = 1; i <
mm; ++i, matrix_ptr += n_columns)
607 if (n_components > 1)
612 if (n_components > 2)
626 if (n_components > 1)
631 if (n_components > 2)
636 for (
int i = 1; i <
mm; ++i)
640 if (n_components > 1)
645 if (n_components > 2)
656 if (n_components > 1)
661 if (n_components > 2)
671 if (n_components > 1)
676 if (n_components > 2)
688 const Number2 *matrix_ptr = matrix +
nn_regular;
690 if (n_components > 1)
692 if (n_components > 2)
693 res2 = matrix_ptr[0] *
in2[0];
694 matrix_ptr += n_columns;
695 for (
int i = 1; i <
mm; ++i, matrix_ptr += n_columns)
698 if (n_components > 1)
700 if (n_components > 2)
706 const Number2 *matrix_ptr = matrix +
nn_regular * n_columns;
708 if (n_components > 1)
710 if (n_components > 2)
711 res2 = matrix_ptr[0] *
in2[0];
712 for (
int i = 1; i <
mm; ++i)
715 if (n_components > 1)
717 if (n_components > 2)
724 if (n_components > 1)
726 if (n_components > 2)
732 if (n_components > 1)
734 if (n_components > 2)
767 static_assert(n_rows > 0 || n_columns > 0,
768 "Specialization only for n_rows, n_columns > 0");
769 Assert(n_rows > 0 && n_columns > 0,
771 std::to_string(n_rows) +
", " +
772 std::to_string(n_columns) +
" was passed!"));
776 constexpr int n_cols =
nn / 2;
777 constexpr int mid =
mm / 2;
779 std::array<Number, mm>
x;
780 for (
int i = 0; i <
mm; ++i)
804 for (
int col = 0; col < n_cols; ++col)
811 val1 = matrix[
nn - 1 - col];
815 val0 = matrix[col * n_columns];
816 val1 = matrix[(col + 1) * n_columns - 1];
828 val0 = matrix[
ind * n_columns + col];
829 val1 = matrix[
ind * n_columns +
nn - 1 - col];
833 val0 = matrix[col * n_columns +
ind];
834 val1 = matrix[(col + 1) * n_columns - 1 -
ind];
848 const Number tmp = matrix[
mid * n_columns + col] *
x[
mid];
855 if (
mm % 2 == 1 &&
nn % 2 == 0)
857 const Number tmp = matrix[col * n_columns +
mid] *
x[
mid];
885 res0 = matrix[n_cols] * (
x[0] +
x[
mm - 1]);
888 const Number2
val0 = matrix[
ind * n_columns + n_cols];
904 res0 = matrix[n_cols * n_columns] * (
x[0] +
x[
mm - 1]);
907 const Number2
val0 = matrix[n_cols * n_columns +
ind];
940 for (
int col = 0; col < n_cols; ++col)
947 val1 = matrix[
nn - 1 - col];
951 val0 = matrix[col * n_columns];
952 val1 = matrix[(
nn - col - 1) * n_columns];
964 val0 = matrix[
ind * n_columns + col];
965 val1 = matrix[
ind * n_columns +
nn - 1 - col];
969 val0 = matrix[col * n_columns +
ind];
970 val1 = matrix[(
nn - col - 1) * n_columns +
ind];
983 val0 = matrix[
mid * n_columns + col];
985 val0 = matrix[col * n_columns +
mid];
1006 val0 = matrix[n_cols];
1008 val0 = matrix[n_cols * n_columns];
1013 val0 = matrix[
ind * n_columns + n_cols];
1015 val0 = matrix[n_cols * n_columns +
ind];
1029 for (
int col = 0; col < n_cols; ++col)
1036 val1 = matrix[
nn - 1 - col];
1040 val0 = matrix[col * n_columns];
1041 val1 = matrix[(col + 1) * n_columns - 1];
1053 val0 = matrix[
ind * n_columns + col];
1054 val1 = matrix[
ind * n_columns +
nn - 1 - col];
1058 val0 = matrix[col * n_columns +
ind];
1059 val1 = matrix[(col + 1) * n_columns - 1 -
ind];
1072 val0 = matrix[
mid * n_columns + col];
1074 val0 = matrix[col * n_columns +
mid];
1095 val0 = matrix[n_cols];
1097 val0 = matrix[n_cols * n_columns];
1104 val0 = matrix[
ind * n_columns + n_cols];
1106 val0 = matrix[n_cols * n_columns +
ind];
1116 val0 = matrix[
mid * n_columns + n_cols];
1118 val0 = matrix[n_cols * n_columns +
mid];
1172 "Negative loop ranges are not allowed!");
1175 const int n_columns =
1182 Assert(n_rows > 0 && n_columns > 0,
1184 std::to_string(n_rows) +
", " +
1185 std::to_string(n_columns) +
" was passed!"));
1197 const int offset = (n_columns + 1) / 2;
1201 std::array<Number, array_length>
xp,
xm;
1202 for (
int i = 0; i <
m_half; ++i)
1216 for (
int col = 0; col <
n_half; ++col)
1223 r0 = matrix[col] *
xp[0];
1224 r1 = matrix[(n_rows - 1) * offset + col] *
xm[0];
1228 r0 = matrix[col * offset] *
xp[0];
1229 r1 = matrix[(n_rows - 1 - col) * offset] *
xm[0];
1236 r1 += matrix[(n_rows - 1 -
ind) * offset + col] *
xm[
ind];
1241 r1 += matrix[(n_rows - 1 - col) * offset +
ind] *
xm[
ind];
1254 else if (
mm % 2 == 1 &&
1279 nn % 2 == 1 &&
mm % 2 == 1 &&
mm > 3)
1401 static_assert(n_rows > 0 && n_columns > 0,
1402 "Specialization requires n_rows, n_columns > 0");
1414 std::array<Number, mm>
x;
1415 for (
unsigned int i = 0; i <
mm; ++i)
1417 for (
unsigned int col = 0; col <
n_half; ++col)
1422 r0 = matrix[col] *
x[0];
1423 r1 = matrix[col + n_columns] *
x[1];
1426 r0 += matrix[col + 2 *
ind * n_columns] *
x[2 *
ind];
1428 matrix[col + (2 *
ind + 1) * n_columns] *
x[2 *
ind + 1];
1434 r0 += matrix[col + (
mm - 1) * n_columns] *
x[
mm - 1];
1458 r0 = matrix[
n_half + shift * n_columns] *
x[shift];
1460 r0 += matrix[
n_half + (2 *
ind + shift) * n_columns] *
1475 std::array<Number, m_half + 1>
xp,
xm;
1476 for (
int i = 0; i <
m_half; ++i)
1489 for (
unsigned int col = 0; col <
n_half; ++col)
1494 r0 = matrix[2 * col * n_columns] *
xp[0];
1495 r1 = matrix[(2 * col + 1) * n_columns] *
xm[0];
1498 r0 += matrix[2 * col * n_columns +
ind] *
xp[
ind];
1499 r1 += matrix[(2 * col + 1) * n_columns +
ind] *
xm[
ind];
1527 r0 = matrix[(
nn - 1) * n_columns] *
xp[0];
1572 typename Number2 = Number>
1596 const unsigned int = 0,
1597 const unsigned int = 0)
1606 n_rows * ((n_columns + 1) / 2));
1609 n_rows * ((n_columns + 1) / 2));
1612 n_rows * ((n_columns + 1) / 2));
1622 n_rows * n_columns));
1626 n_rows * n_columns));
1636 const unsigned int dummy1 = 0,
1637 const unsigned int dummy2 = 0)
1671 template <
int direction,
bool contract_over_rows,
bool add,
int str
ide = 1>
1673 values(
const Number in[], Number out[])
const
1685 template <
int direction,
bool contract_over_rows,
bool add,
int str
ide = 1>
1701 template <
int direction,
bool contract_over_rows,
bool add>
1721 template <
int direction,
bool contract_over_rows,
bool add>
1737 template <
int direction,
bool contract_over_rows,
bool add>
1756 template <
int direction,
bool contract_over_rows,
bool add>
1806 template <
int direction,
1831 template <
int direction,
1843 static_assert(
one_line ==
false || direction == dim - 1,
1844 "Single-line evaluation only works for direction=dim-1.");
1847 "The given array shape_data must not be the null pointer!"));
1848 Assert(dim == direction + 1 ||
one_line ==
true || n_rows == n_columns ||
1850 ExcMessage(
"In-place operation only supported for "
1851 "n_rows==n_columns or single-line interpolation"));
1859 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1924 , n_rows(
numbers::invalid_unsigned_int)
1925 , n_columns(
numbers::invalid_unsigned_int)
1934 const unsigned int n_rows = 0,
1935 const unsigned int n_columns = 0)
1940 , n_columns(n_columns)
1946 n_rows * ((n_columns + 1) / 2));
1949 n_rows * ((n_columns + 1) / 2));
1952 n_rows * ((n_columns + 1) / 2));
1962 n_rows * n_columns));
1966 n_rows * n_columns));
1976 const unsigned int n_rows = 0,
1977 const unsigned int n_columns = 0)
1982 , n_columns(n_columns)
1985 template <
int direction,
bool contract_over_rows,
bool add,
int str
ide = 1>
1994 template <
int direction,
bool contract_over_rows,
bool add,
int str
ide = 1>
2005 template <
int direction,
bool contract_over_rows,
bool add>
2016 template <
int direction,
bool contract_over_rows,
bool add>
2025 template <
int direction,
bool contract_over_rows,
bool add>
2037 template <
int direction,
bool contract_over_rows,
bool add>
2049 template <
int direction,
2073 template <
int direction,
2085 static_assert(
one_line ==
false || direction == dim - 1,
2086 "Single-line evaluation only works for direction=dim-1.");
2089 "The given array shape_data must not be the null pointer!"));
2090 Assert(dim == direction + 1 ||
one_line ==
true || n_rows == n_columns ||
2092 ExcMessage(
"In-place operation only supported for "
2093 "n_rows==n_columns or single-line interpolation"));
2099 direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
2101 const int n_blocks2 = direction >= dim - 1 ?
2120 (direction != 0 || stride != 1)>(
2152 template <
int direction,
2154 typename Number = double,
2155 typename Number2 =
double>
2166 constexpr int n_rows = fe_degree + 1;
2167 constexpr int n_columns = n_q_points_1d;
2172 data.shape_values_eo.data() :
2176 ExcMessage(
"Cannot add into result if contract_over_rows = true"));
2216 template <
int direction,
2217 int normal_direction,
2219 typename Number = double,
2220 typename Number2 =
double>
2230 static_assert(direction != normal_direction,
2231 "Cannot interpolate tangentially in normal direction");
2233 constexpr int n_rows =
std::max(fe_degree, 0);
2234 constexpr int n_columns = n_q_points_1d;
2237 data.shape_values_eo.data() :
2242 (direction > normal_direction) ?
2248 (direction > normal_direction) ?
2250 ((direction + 1 < dim) ?
2251 (
Utilities::pow(fe_degree, dim - 2 - direction) * n_q_points_1d) :
2282 out -=
n_blocks1 * (n_columns - 1) * stride;
2303 in +=
n_blocks1 * (n_columns - 1) * stride;
2364 inline std::enable_if_t<contract_onto_face, void>
2366 const std::array<int, 2> &n_blocks,
2367 const std::array<int, 2> &steps,
2368 const Number *input,
2376 Number *
output1 = output + n_blocks[0] * n_blocks[1];
2378 for (
int i2 = 0;
i2 < n_blocks[1]; ++
i2)
2380 for (
int i1 = 0;
i1 < n_blocks[0]; ++
i1)
2382 Number
res0 = shape_values[0] * input[0];
2385 res1 = shape_values[n_rows] * input[0];
2387 res2 = shape_values[2 * n_rows] * input[0];
2390 res0 += shape_values[
ind] * input[stride *
ind];
2392 res1 += shape_values[
ind + n_rows] * input[stride *
ind];
2394 res2 += shape_values[
ind + 2 * n_rows] * input[stride *
ind];
2414 output += n_blocks[0];
2434 const unsigned int n_q_points_1d)
2436 return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
2437 (n_q_points_1d <= 3 * fe_degree / 2 + 1);
2454 inline std::enable_if_t<!contract_onto_face, void>
2456 const std::array<int, 2> &n_blocks,
2457 const std::array<int, 2> &steps,
2458 const Number *input,
2466 const Number *
input1 = input + n_blocks[0] * n_blocks[1];
2467 const Number *
input2 =
input1 + n_blocks[0] * n_blocks[1];
2468 for (
int i2 = 0;
i2 < n_blocks[1]; ++
i2)
2470 for (
int i1 = 0;
i1 < n_blocks[0]; ++
i1)
2472 const Number in = input[
i1];
2478 for (
int col = 0; col < n_rows; ++col)
2481 add ? (output[col * stride] + shape_values[col] * in) :
2482 (shape_values[col] * in);
2484 result += shape_values[col + n_rows] *
in1;
2486 result += shape_values[col + 2 * n_rows] *
in2;
2488 output[col * stride] =
result;
2492 input += n_blocks[0];
2501 template <
int dim,
int n_po
ints_1d_
template,
typename Number>
2504 const unsigned int n_components,
2521 for (
unsigned int c = 0; c < n_components; ++c)
2525 const unsigned int shift =
2527 data[0] *= weights[shift];
2530 const Number weight = weights[shift + 1];
2539 template <
int dim,
int n_po
ints_1d_
template,
typename Number>
2542 const unsigned int n_components,
2551 ExcMessage(
"The function can only with add number of points"));
2566 for (
unsigned int c = 0; c < n_components; ++c)
2570 const unsigned int shift =
2574 const Number
weight1 = weights[shift];
2577 data[c++] *= weights[shift + 1];
2578 const Number
weight2 = weights[shift + 2];
2586 template <
int dim,
int n_po
ints_1d_
template,
typename Number>
2589 const unsigned int n_components,
2612 if (weight == Number(-1.0) || weight ==
data)
2622 weights[c] = Number(-1.0);
2624 for (
unsigned int c = 0; c < n_components; ++c)
2629 const unsigned int shift =
2647 template <
int dim,
int n_po
ints_1d_
template,
typename Number>
2651 const unsigned int n_components,
2660 ExcMessage(
"The function can only with add number of points"));
2681 if (weight == Number(-1.0) || weight ==
data)
2691 weights[c] = Number(-1.0);
2693 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
2698 const unsigned int shift =
#define DEAL_II_ALWAYS_INLINE
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
std::vector< index_type > data
constexpr T fixed_power(const T t)
constexpr T pow(const T base, const int iexp)
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
void weight_fe_q_dofs_by_entity(const Number *weights, const unsigned int n_components, const int n_points_1d_non_template, Number *data)
std::enable_if_t<(variant==evaluate_general), void > apply_matrix_vector_product(const Number2 *matrix, const Number *in, Number *out)
void weight_fe_q_dofs_by_entity_shifted(const Number *weights, const unsigned int n_components, const int n_points_1d_non_template, Number *data)
std::enable_if_t< contract_onto_face, void > interpolate_to_face(const Number2 *shape_values, const std::array< int, 2 > &n_blocks, const std::array< int, 2 > &steps, const Number *input, Number *DEAL_II_RESTRICT output, const int n_rows_runtime=0, const int stride_runtime=1)
bool compute_weights_fe_q_dofs_by_entity(const Number *data, const unsigned int n_components, const int n_points_1d_non_template, Number *weights)
bool compute_weights_fe_q_dofs_by_entity_shifted(const Number *data, const unsigned int n_components, const int n_points_1d_non_template, Number *weights)
@ evaluate_symmetric_hierarchical
constexpr unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static void normal(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const Number *in, Number *out, const bool add_into_result=false, const int subface_index_1d=0)
static void tangential(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const Number *in, Number *out, const int subface_index_1d=0)
const Number2 * shape_hessians
const Number2 * shape_gradients
const unsigned int n_rows
void hessians_one_line(const Number in[], Number out[]) const
void gradients(const Number *in, Number *out) const
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int n_rows=0, const unsigned int n_columns=0)
const Number2 * shape_values
void values(const Number *in, Number *out) const
void hessians(const Number *in, Number *out) const
void gradients_one_line(const Number in[], Number out[]) const
void values_one_line(const Number in[], Number out[]) const
EvaluatorTensorProduct(const Number2 *shape_values, const Number2 *shape_gradients, const Number2 *shape_hessians, const unsigned int n_rows=0, const unsigned int n_columns=0)
const unsigned int n_columns
void values(const Number in[], Number out[]) const
EvaluatorTensorProduct(const Number2 *shape_values, const Number2 *shape_gradients, const Number2 *shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
void values_one_line(const Number in[], Number out[]) const
void gradients(const Number in[], Number out[]) const
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int=0, const unsigned int=0)
const Number2 * shape_gradients
static constexpr unsigned int n_rows_of_product
void hessians(const Number in[], Number out[]) const
const Number2 * shape_values
const Number2 * shape_hessians
static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out)
static constexpr unsigned int n_columns_of_product
void gradients_one_line(const Number in[], Number out[]) const
void hessians_one_line(const Number in[], Number out[]) const