Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
tensor_product_kernels.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_matrix_free_tensor_product_kernels_h
17#define dealii_matrix_free_tensor_product_kernels_h
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/base/aligned_vector.h>
22#include <deal.II/base/utilities.h>
23
24#include <deal.II/matrix_free/shape_info.h>
25
26
27DEAL_II_NAMESPACE_OPEN
28
29
30
31namespace internal
32{
68
69
70
74 enum class EvaluatorQuantity
75 {
79 value,
83 gradient,
87 hessian
88 };
89
90
91
107 template <EvaluatorVariant variant,
108 EvaluatorQuantity quantity,
109 int n_rows,
110 int n_columns,
111 int stride_in,
112 int stride_out,
113 bool transpose_matrix,
114 bool add,
115 typename Number,
116 typename Number2>
117 std::enable_if_t<(variant == evaluate_general), void>
118 apply_matrix_vector_product(const Number2 *matrix,
119 const Number *in,
120 Number *out)
121 {
122 // We can only statically assert that one argument is non-zero because
123 // face evaluation might instantiate some functions, so we need to use the
124 // run-time assert to verify that we do not end up involuntarily.
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,
128 ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
129 std::to_string(n_rows) + ", " +
130 std::to_string(n_columns) + " was passed!"));
131 static_assert(quantity == EvaluatorQuantity::value,
132 "This function should only use EvaluatorQuantity::value");
133
134 constexpr int mm = transpose_matrix ? n_rows : n_columns,
135 nn = transpose_matrix ? n_columns : n_rows;
136
137 std::array<Number, mm> x;
138 for (int i = 0; i < mm; ++i)
139 x[i] = in[stride_in * i];
140 for (int col = 0; col < nn; ++col)
141 {
142 Number res0;
143 if (transpose_matrix == true)
144 {
145 res0 = matrix[col] * x[0];
146 for (int i = 1; i < mm; ++i)
147 {
148 const Number2 mji = matrix[i * n_columns + col];
149 if constexpr (numbers::NumberTraits<Number>::is_complex &&
150 numbers::NumberTraits<Number2>::is_complex)
151 {
152 res0.real(res0.real() + mji.real() * x[i].real() -
153 mji.imag() * x[i].imag());
154 res0.imag(res0.imag() + mji.imag() * x[i].real() +
155 mji.real() * x[i].imag());
156 }
157 else
158 res0 += mji * x[i];
159 }
160 }
161 else
162 {
163 res0 = matrix[col * n_columns] * x[0];
164 for (int i = 1; i < mm; ++i)
165 {
166 const Number2 mij = matrix[col * n_columns + i];
167 if constexpr (numbers::NumberTraits<Number>::is_complex &&
168 numbers::NumberTraits<Number2>::is_complex)
169 {
170 res0.real(res0.real() + mij.real() * x[i].real() -
171 mij.imag() * x[i].imag());
172 res0.imag(res0.imag() + mij.imag() * x[i].real() +
173 mij.real() * x[i].imag());
174 }
175 else
176 res0 += mij * x[i];
177 }
178 }
179 if (add)
180 out[stride_out * col] += res0;
181 else
182 out[stride_out * col] = res0;
183 }
184 }
185
186
187
192 template <EvaluatorVariant variant,
193 EvaluatorQuantity quantity,
194 bool transpose_matrix,
195 bool add,
196 bool consider_strides,
197 typename Number,
198 typename Number2>
199 std::enable_if_t<(variant == evaluate_general), void>
200 apply_matrix_vector_product(const Number2 *matrix,
201 const Number *in,
202 Number *out,
203 const int n_rows,
204 const int n_columns,
205 const int stride_in_given,
206 const int stride_out_given)
207 {
208 const int mm = transpose_matrix ? n_rows : n_columns,
209 nn = transpose_matrix ? n_columns : n_rows;
210 Assert(n_rows > 0 && n_columns > 0,
211 ExcInternalError("Empty evaluation task!"));
212 Assert(n_rows > 0 && n_columns > 0,
213 ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
214 std::to_string(n_rows) + ", " +
215 std::to_string(n_columns) + " was passed!"));
216
217 static_assert(quantity == EvaluatorQuantity::value,
218 "This function should only use EvaluatorQuantity::value");
219
220 Assert(consider_strides || (stride_in_given == 1 && stride_out_given == 1),
221 ExcInternalError());
222 const int stride_in = consider_strides ? stride_in_given : 1;
223 const int stride_out = consider_strides ? stride_out_given : 1;
224
225 // specialization for n_rows = 2 that manually unrolls the innermost loop
226 // to make the operation perform better (not completely as good as the
227 // templated one, but much better than the generic version down below,
228 // because the loop over col can be more effectively unrolled by the
229 // compiler)
230 if (transpose_matrix && n_rows == 2)
231 {
232 const Number2 *matrix_1 = matrix + n_columns;
233 const Number x0 = in[0], x1 = in[stride_in];
234 for (int col = 0; col < nn; ++col)
235 {
236 const Number result = matrix[col] * x0 + matrix_1[col] * x1;
237 if (add)
238 out[stride_out * col] += result;
239 else
240 out[stride_out * col] = result;
241 }
242 }
243 else if (transpose_matrix && n_rows == 3)
244 {
245 const Number2 *matrix_1 = matrix + n_columns;
246 const Number2 *matrix_2 = matrix_1 + n_columns;
247 const Number x0 = in[0], x1 = in[stride_in], x2 = in[2 * stride_in];
248 for (int col = 0; col < nn; ++col)
249 {
250 const Number result =
251 matrix[col] * x0 + matrix_1[col] * x1 + matrix_2[col] * x2;
252 if (add)
253 out[stride_out * col] += result;
254 else
255 out[stride_out * col] = result;
256 }
257 }
258 else if (std::abs(in - out) < std::min(stride_out * nn, stride_in * mm))
259 {
260 Assert(mm <= 128,
261 ExcNotImplemented("For large sizes, arrays may not overlap"));
262 std::array<Number, 129> x;
263 for (int i = 0; i < mm; ++i)
264 x[i] = in[stride_in * i];
265
266 for (int col = 0; col < nn; ++col)
267 {
268 Number res0;
269 if (transpose_matrix == true)
270 {
271 res0 = matrix[col] * x[0];
272 for (int i = 1; i < mm; ++i)
273 res0 += matrix[i * n_columns + col] * x[i];
274 }
275 else
276 {
277 res0 = matrix[col * n_columns] * x[0];
278 for (int i = 1; i < mm; ++i)
279 res0 += matrix[col * n_columns + i] * x[i];
280 }
281 if (add)
282 out[stride_out * col] += res0;
283 else
284 out[stride_out * col] = res0;
285 }
286 }
287 else
288 {
289 int nn_regular = (nn / 4) * 4;
290 for (int col = 0; col < nn_regular; col += 4)
291 {
292 Number res0, res1, res2, res3;
293 if (transpose_matrix == true)
294 {
295 const Number2 *matrix_ptr = matrix + col;
296 res0 = matrix_ptr[0] * in[0];
297 res1 = matrix_ptr[1] * in[0];
298 res2 = matrix_ptr[2] * in[0];
299 res3 = matrix_ptr[3] * in[0];
300 matrix_ptr += n_columns;
301 for (int i = 1; i < mm; ++i, matrix_ptr += n_columns)
302 {
303 res0 += matrix_ptr[0] * in[stride_in * i];
304 res1 += matrix_ptr[1] * in[stride_in * i];
305 res2 += matrix_ptr[2] * in[stride_in * i];
306 res3 += matrix_ptr[3] * in[stride_in * i];
307 }
308 }
309 else
310 {
311 const Number2 *matrix_0 = matrix + col * n_columns;
312 const Number2 *matrix_1 = matrix + (col + 1) * n_columns;
313 const Number2 *matrix_2 = matrix + (col + 2) * n_columns;
314 const Number2 *matrix_3 = matrix + (col + 3) * n_columns;
315
316 res0 = matrix_0[0] * in[0];
317 res1 = matrix_1[0] * in[0];
318 res2 = matrix_2[0] * in[0];
319 res3 = matrix_3[0] * in[0];
320 for (int i = 1; i < mm; ++i)
321 {
322 res0 += matrix_0[i] * in[stride_in * i];
323 res1 += matrix_1[i] * in[stride_in * i];
324 res2 += matrix_2[i] * in[stride_in * i];
325 res3 += matrix_3[i] * in[stride_in * i];
326 }
327 }
328 if (add)
329 {
330 out[0] += res0;
331 out[stride_out] += res1;
332 out[2 * stride_out] += res2;
333 out[3 * stride_out] += res3;
334 }
335 else
336 {
337 out[0] = res0;
338 out[stride_out] = res1;
339 out[2 * stride_out] = res2;
340 out[3 * stride_out] = res3;
341 }
342 out += 4 * stride_out;
343 }
344 if (nn - nn_regular == 3)
345 {
346 Number res0, res1, res2;
347 if (transpose_matrix == true)
348 {
349 const Number2 *matrix_ptr = matrix + nn_regular;
350 res0 = matrix_ptr[0] * in[0];
351 res1 = matrix_ptr[1] * in[0];
352 res2 = matrix_ptr[2] * in[0];
353 matrix_ptr += n_columns;
354 for (int i = 1; i < mm; ++i, matrix_ptr += n_columns)
355 {
356 res0 += matrix_ptr[0] * in[stride_in * i];
357 res1 += matrix_ptr[1] * in[stride_in * i];
358 res2 += matrix_ptr[2] * in[stride_in * i];
359 }
360 }
361 else
362 {
363 const Number2 *matrix_0 = matrix + nn_regular * n_columns;
364 const Number2 *matrix_1 = matrix + (nn_regular + 1) * n_columns;
365 const Number2 *matrix_2 = matrix + (nn_regular + 2) * n_columns;
366
367 res0 = matrix_0[0] * in[0];
368 res1 = matrix_1[0] * in[0];
369 res2 = matrix_2[0] * in[0];
370 for (int i = 1; i < mm; ++i)
371 {
372 res0 += matrix_0[i] * in[stride_in * i];
373 res1 += matrix_1[i] * in[stride_in * i];
374 res2 += matrix_2[i] * in[stride_in * i];
375 }
376 }
377 if (add)
378 {
379 out[0] += res0;
380 out[stride_out] += res1;
381 out[2 * stride_out] += res2;
382 }
383 else
384 {
385 out[0] = res0;
386 out[stride_out] = res1;
387 out[2 * stride_out] = res2;
388 }
389 }
390 else if (nn - nn_regular == 2)
391 {
392 Number res0, res1;
393 if (transpose_matrix == true)
394 {
395 const Number2 *matrix_ptr = matrix + nn_regular;
396 res0 = matrix_ptr[0] * in[0];
397 res1 = matrix_ptr[1] * in[0];
398 matrix_ptr += n_columns;
399 for (int i = 1; i < mm; ++i, matrix_ptr += n_columns)
400 {
401 res0 += matrix_ptr[0] * in[stride_in * i];
402 res1 += matrix_ptr[1] * in[stride_in * i];
403 }
404 }
405 else
406 {
407 const Number2 *matrix_0 = matrix + nn_regular * n_columns;
408 const Number2 *matrix_1 = matrix + (nn_regular + 1) * n_columns;
409
410 res0 = matrix_0[0] * in[0];
411 res1 = matrix_1[0] * in[0];
412 for (int i = 1; i < mm; ++i)
413 {
414 res0 += matrix_0[i] * in[stride_in * i];
415 res1 += matrix_1[i] * in[stride_in * i];
416 }
417 }
418 if (add)
419 {
420 out[0] += res0;
421 out[stride_out] += res1;
422 }
423 else
424 {
425 out[0] = res0;
426 out[stride_out] = res1;
427 }
428 }
429 else if (nn - nn_regular == 1)
430 {
431 Number res0;
432 if (transpose_matrix == true)
433 {
434 const Number2 *matrix_ptr = matrix + nn_regular;
435 res0 = matrix_ptr[0] * in[0];
436 matrix_ptr += n_columns;
437 for (int i = 1; i < mm; ++i, matrix_ptr += n_columns)
438 res0 += matrix_ptr[0] * in[stride_in * i];
439 }
440 else
441 {
442 const Number2 *matrix_ptr = matrix + nn_regular * n_columns;
443 res0 = matrix_ptr[0] * in[0];
444 for (int i = 1; i < mm; ++i)
445 res0 += matrix_ptr[i] * in[stride_in * i];
446 }
447 if (add)
448 out[0] += res0;
449 else
450 out[0] = res0;
451 }
452 }
453 }
454
455
456
463 template <EvaluatorVariant variant,
464 EvaluatorQuantity quantity,
465 int n_rows,
466 int n_columns,
467 int stride_in,
468 int stride_out,
469 bool transpose_matrix,
470 bool add,
471 typename Number,
472 typename Number2>
473 std::enable_if_t<(variant == evaluate_symmetric), void>
474 apply_matrix_vector_product(const Number2 *matrix,
475 const Number *in,
476 Number *out)
477 {
478 // We can only statically assert that one argument is non-zero because
479 // face evaluation might instantiate some functions, so we need to use the
480 // run-time assert to verify that we do not end up involuntarily.
481 static_assert(n_rows > 0 || n_columns > 0,
482 "Specialization only for n_rows, n_columns > 0");
483 Assert(n_rows > 0 && n_columns > 0,
484 ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
485 std::to_string(n_rows) + ", " +
486 std::to_string(n_columns) + " was passed!"));
487
488 constexpr int mm = transpose_matrix ? n_rows : n_columns,
489 nn = transpose_matrix ? n_columns : n_rows;
490 constexpr int n_cols = nn / 2;
491 constexpr int mid = mm / 2;
492
493 std::array<Number, mm> x;
494 for (int i = 0; i < mm; ++i)
495 x[i] = in[stride_in * i];
496
497 if (quantity == EvaluatorQuantity::value)
498 {
499 // In this case, the 1d shape values read (sorted lexicographically,
500 // rows run over 1d dofs, columns over quadrature points):
501 // Q2 --> [ 0.687 0 -0.087 ]
502 // [ 0.4 1 0.4 ]
503 // [-0.087 0 0.687 ]
504 // Q3 --> [ 0.66 0.003 0.002 0.049 ]
505 // [ 0.521 1.005 -0.01 -0.230 ]
506 // [-0.230 -0.01 1.005 0.521 ]
507 // [ 0.049 0.002 0.003 0.66 ]
508 // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ]
509 // [ 0.608 1.059 0 0.039 0.176 ]
510 // [-0.409 -0.113 1 -0.113 -0.409 ]
511 // [ 0.176 0.039 0 1.059 0.608 ]
512 // [-0.032 -0.007 0 0.022 0.658 ]
513 //
514 // In these matrices, we want to use avoid computations involving
515 // zeros and ones and use the symmetry in entries starting from (1,1)
516 // forward and (N,N) backward, respectively to reduce the number of
517 // read operations.
518 for (int col = 0; col < n_cols; ++col)
519 {
520 Number2 val0, val1;
521 Number res0, res1;
522 if (transpose_matrix == true)
523 {
524 val0 = matrix[col];
525 val1 = matrix[nn - 1 - col];
526 }
527 else
528 {
529 val0 = matrix[col * n_columns];
530 val1 = matrix[(col + 1) * n_columns - 1];
531 }
532 if (mid > 0)
533 {
534 res0 = val0 * x[0];
535 res1 = val1 * x[0];
536 res0 += val1 * x[mm - 1];
537 res1 += val0 * x[mm - 1];
538 for (int ind = 1; ind < mid; ++ind)
539 {
540 if (transpose_matrix == true)
541 {
542 val0 = matrix[ind * n_columns + col];
543 val1 = matrix[ind * n_columns + nn - 1 - col];
544 }
545 else
546 {
547 val0 = matrix[col * n_columns + ind];
548 val1 = matrix[(col + 1) * n_columns - 1 - ind];
549 }
550 res0 += val0 * x[ind];
551 res1 += val1 * x[ind];
552 res0 += val1 * x[mm - 1 - ind];
553 res1 += val0 * x[mm - 1 - ind];
554 }
555 }
556 else
557 res0 = res1 = Number();
558 if (transpose_matrix == true)
559 {
560 if (mm % 2 == 1)
561 {
562 const Number tmp = matrix[mid * n_columns + col] * x[mid];
563 res0 += tmp;
564 res1 += tmp;
565 }
566 }
567 else
568 {
569 if (mm % 2 == 1 && nn % 2 == 0)
570 {
571 const Number tmp = matrix[col * n_columns + mid] * x[mid];
572 res0 += tmp;
573 res1 += tmp;
574 }
575 }
576 if (add)
577 {
578 out[stride_out * col] += res0;
579 out[stride_out * (nn - 1 - col)] += res1;
580 }
581 else
582 {
583 out[stride_out * col] = res0;
584 out[stride_out * (nn - 1 - col)] = res1;
585 }
586 }
587 if (transpose_matrix == true && nn % 2 == 1 && mm % 2 == 1)
588 {
589 if (add)
590 out[stride_out * n_cols] += x[mid];
591 else
592 out[stride_out * n_cols] = x[mid];
593 }
594 else if (transpose_matrix == true && nn % 2 == 1)
595 {
596 Number res0;
597 if (mid > 0)
598 {
599 res0 = matrix[n_cols] * (x[0] + x[mm - 1]);
600 for (int ind = 1; ind < mid; ++ind)
601 {
602 const Number2 val0 = matrix[ind * n_columns + n_cols];
603 res0 += val0 * (x[ind] + in[mm - 1 - ind]);
604 }
605 }
606 else
607 res0 = Number();
608 if (add)
609 out[stride_out * n_cols] += res0;
610 else
611 out[stride_out * n_cols] = res0;
612 }
613 else if (transpose_matrix == false && nn % 2 == 1)
614 {
615 Number res0;
616 if (mid > 0)
617 {
618 res0 = matrix[n_cols * n_columns] * (x[0] + x[mm - 1]);
619 for (int ind = 1; ind < mid; ++ind)
620 {
621 const Number2 val0 = matrix[n_cols * n_columns + ind];
622 res0 += val0 * (x[ind] + x[mm - 1 - ind]);
623 ;
624 }
625 if (mm % 2)
626 res0 += x[mid];
627 }
628 else
629 res0 = in[0];
630 if (add)
631 out[stride_out * n_cols] += res0;
632 else
633 out[stride_out * n_cols] = res0;
634 }
635 }
636 else if (quantity == EvaluatorQuantity::gradient)
637 {
638 // For the specialized loop used for gradient computations we again
639 // exploit symmetries according to the following entries (sorted
640 // lexicographically, rows run over 1d dofs, columns over quadrature
641 // points):
642 // Q2 --> [-2.549 -1 0.549 ]
643 // [ 3.098 0 -3.098 ]
644 // [-0.549 1 2.549 ]
645 // Q3 --> [-4.315 -1.03 0.5 -0.44 ]
646 // [ 6.07 -1.44 -2.97 2.196 ]
647 // [-2.196 2.97 1.44 -6.07 ]
648 // [ 0.44 -0.5 1.03 4.315 ]
649 // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ]
650 // [10.111 -2.76 -2.667 2.066 -2.306 ]
651 // [-5.688 5.773 0 -5.773 5.688 ]
652 // [ 2.306 -2.066 2.667 2.76 -10.111 ]
653 // [-0.413 0.353 -0.333 -0.353 0.413 ]
654 for (int col = 0; col < n_cols; ++col)
655 {
656 Number2 val0, val1;
657 Number res0, res1;
658 if (transpose_matrix == true)
659 {
660 val0 = matrix[col];
661 val1 = matrix[nn - 1 - col];
662 }
663 else
664 {
665 val0 = matrix[col * n_columns];
666 val1 = matrix[(nn - col - 1) * n_columns];
667 }
668 if (mid > 0)
669 {
670 res0 = val0 * x[0];
671 res1 = val1 * x[0];
672 res0 -= val1 * x[mm - 1];
673 res1 -= val0 * x[mm - 1];
674 for (int ind = 1; ind < mid; ++ind)
675 {
676 if (transpose_matrix == true)
677 {
678 val0 = matrix[ind * n_columns + col];
679 val1 = matrix[ind * n_columns + nn - 1 - col];
680 }
681 else
682 {
683 val0 = matrix[col * n_columns + ind];
684 val1 = matrix[(nn - col - 1) * n_columns + ind];
685 }
686 res0 += val0 * x[ind];
687 res1 += val1 * x[ind];
688 res0 -= val1 * x[mm - 1 - ind];
689 res1 -= val0 * x[mm - 1 - ind];
690 }
691 }
692 else
693 res0 = res1 = Number();
694 if (mm % 2 == 1)
695 {
696 if (transpose_matrix == true)
697 val0 = matrix[mid * n_columns + col];
698 else
699 val0 = matrix[col * n_columns + mid];
700 const Number tmp = val0 * x[mid];
701 res0 += tmp;
702 res1 -= tmp;
703 }
704 if (add)
705 {
706 out[stride_out * col] += res0;
707 out[stride_out * (nn - 1 - col)] += res1;
708 }
709 else
710 {
711 out[stride_out * col] = res0;
712 out[stride_out * (nn - 1 - col)] = res1;
713 }
714 }
715 if (nn % 2 == 1)
716 {
717 Number2 val0;
718 Number res0;
719 if (transpose_matrix == true)
720 val0 = matrix[n_cols];
721 else
722 val0 = matrix[n_cols * n_columns];
723 res0 = val0 * (x[0] - x[mm - 1]);
724 for (int ind = 1; ind < mid; ++ind)
725 {
726 if (transpose_matrix == true)
727 val0 = matrix[ind * n_columns + n_cols];
728 else
729 val0 = matrix[n_cols * n_columns + ind];
730 Number in1 = val0 * (x[ind] - x[mm - 1 - ind]);
731 res0 += in1;
732 }
733 if (add)
734 out[stride_out * n_cols] += res0;
735 else
736 out[stride_out * n_cols] = res0;
737 }
738 }
739 else
740 {
741 // Hessians are almost the same as values, apart from some missing '1'
742 // entries
743 for (int col = 0; col < n_cols; ++col)
744 {
745 Number2 val0, val1;
746 Number res0, res1;
747 if (transpose_matrix == true)
748 {
749 val0 = matrix[col];
750 val1 = matrix[nn - 1 - col];
751 }
752 else
753 {
754 val0 = matrix[col * n_columns];
755 val1 = matrix[(col + 1) * n_columns - 1];
756 }
757 if (mid > 0)
758 {
759 res0 = val0 * x[0];
760 res1 = val1 * x[0];
761 res0 += val1 * x[mm - 1];
762 res1 += val0 * x[mm - 1];
763 for (int ind = 1; ind < mid; ++ind)
764 {
765 if (transpose_matrix == true)
766 {
767 val0 = matrix[ind * n_columns + col];
768 val1 = matrix[ind * n_columns + nn - 1 - col];
769 }
770 else
771 {
772 val0 = matrix[col * n_columns + ind];
773 val1 = matrix[(col + 1) * n_columns - 1 - ind];
774 }
775 res0 += val0 * x[ind];
776 res1 += val1 * x[ind];
777 res0 += val1 * x[mm - 1 - ind];
778 res1 += val0 * x[mm - 1 - ind];
779 }
780 }
781 else
782 res0 = res1 = Number();
783 if (mm % 2 == 1)
784 {
785 if (transpose_matrix == true)
786 val0 = matrix[mid * n_columns + col];
787 else
788 val0 = matrix[col * n_columns + mid];
789 const Number tmp = val0 * x[mid];
790 res0 += tmp;
791 res1 += tmp;
792 }
793 if (add)
794 {
795 out[stride_out * col] += res0;
796 out[stride_out * (nn - 1 - col)] += res1;
797 }
798 else
799 {
800 out[stride_out * col] = res0;
801 out[stride_out * (nn - 1 - col)] = res1;
802 }
803 }
804 if (nn % 2 == 1)
805 {
806 Number2 val0;
807 Number res0;
808 if (transpose_matrix == true)
809 val0 = matrix[n_cols];
810 else
811 val0 = matrix[n_cols * n_columns];
812 if (mid > 0)
813 {
814 res0 = val0 * (x[0] + x[mm - 1]);
815 for (int ind = 1; ind < mid; ++ind)
816 {
817 if (transpose_matrix == true)
818 val0 = matrix[ind * n_columns + n_cols];
819 else
820 val0 = matrix[n_cols * n_columns + ind];
821 Number in1 = val0 * (x[ind] + x[mm - 1 - ind]);
822 res0 += in1;
823 }
824 }
825 else
826 res0 = Number();
827 if (mm % 2 == 1)
828 {
829 if (transpose_matrix == true)
830 val0 = matrix[mid * n_columns + n_cols];
831 else
832 val0 = matrix[n_cols * n_columns + mid];
833 res0 += val0 * x[mid];
834 }
835 if (add)
836 out[stride_out * n_cols] += res0;
837 else
838 out[stride_out * n_cols] = res0;
839 }
840 }
841 }
842
843
844
863 template <EvaluatorVariant variant,
864 EvaluatorQuantity quantity,
865 int n_rows_static,
866 int n_columns_static,
867 int stride_in_static,
868 int stride_out_static,
869 bool transpose_matrix,
870 bool add,
871 typename Number,
872 typename Number2>
873#ifndef DEBUG
874 inline DEAL_II_ALWAYS_INLINE
875#endif
876 std::enable_if_t<(variant == evaluate_evenodd), void>
877 apply_matrix_vector_product(const Number2 *DEAL_II_RESTRICT matrix,
878 const Number *in,
879 Number *out,
880 int n_rows_runtime = 0,
881 int n_columns_runtime = 0,
882 int stride_in_runtime = 0,
883 int stride_out_runtime = 0)
884 {
885 static_assert(n_rows_static >= 0 && n_columns_static >= 0,
886 "Negative loop ranges are not allowed!");
887
888 const int n_rows = n_rows_static == 0 ? n_rows_runtime : n_rows_static;
889 const int n_columns =
890 n_rows_static == 0 ? n_columns_runtime : n_columns_static;
891 const int stride_in =
892 stride_in_static == 0 ? stride_in_runtime : stride_in_static;
893 const int stride_out =
894 stride_out_static == 0 ? stride_out_runtime : stride_out_static;
895
896 Assert(n_rows > 0 && n_columns > 0,
897 ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
898 std::to_string(n_rows) + ", " +
899 std::to_string(n_columns) + " was passed!"));
900
901 const int mm = transpose_matrix ? n_rows : n_columns,
902 nn = transpose_matrix ? n_columns : n_rows;
903 const int n_half = nn / 2;
904 const int m_half = mm / 2;
905
906 constexpr int array_length =
907 (n_rows_static == 0) ?
908 16 // for non-templated execution
909 :
910 (1 + (transpose_matrix ? n_rows_static : n_columns_static) / 2);
911 const int offset = (n_columns + 1) / 2;
912
913 Assert(m_half <= array_length, ExcNotImplemented());
914
915 std::array<Number, array_length> xp, xm;
916 for (int i = 0; i < m_half; ++i)
917 {
918 if (transpose_matrix == true && quantity == EvaluatorQuantity::gradient)
919 {
920 xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
921 xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
922 }
923 else
924 {
925 xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
926 xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
927 }
928 }
929 Number xmid = in[stride_in * m_half];
930 for (int col = 0; col < n_half; ++col)
931 {
932 Number r0, r1;
933 if (m_half > 0)
934 {
935 if (transpose_matrix == true)
936 {
937 r0 = matrix[col] * xp[0];
938 r1 = matrix[(n_rows - 1) * offset + col] * xm[0];
939 }
940 else
941 {
942 r0 = matrix[col * offset] * xp[0];
943 r1 = matrix[(n_rows - 1 - col) * offset] * xm[0];
944 }
945 for (int ind = 1; ind < m_half; ++ind)
946 {
947 if (transpose_matrix == true)
948 {
949 r0 += matrix[ind * offset + col] * xp[ind];
950 r1 += matrix[(n_rows - 1 - ind) * offset + col] * xm[ind];
951 }
952 else
953 {
954 r0 += matrix[col * offset + ind] * xp[ind];
955 r1 += matrix[(n_rows - 1 - col) * offset + ind] * xm[ind];
956 }
957 }
958 }
959 else
960 r0 = r1 = Number();
961 if (mm % 2 == 1 && transpose_matrix == true)
962 {
963 if (quantity == EvaluatorQuantity::gradient)
964 r1 += matrix[m_half * offset + col] * xmid;
965 else
966 r0 += matrix[m_half * offset + col] * xmid;
967 }
968 else if (mm % 2 == 1 &&
969 (nn % 2 == 0 || quantity != EvaluatorQuantity::value ||
970 mm == 3))
971 r0 += matrix[col * offset + m_half] * xmid;
972
973 if (add)
974 {
975 out[stride_out * col] += r0 + r1;
976 if (quantity == EvaluatorQuantity::gradient &&
977 transpose_matrix == false)
978 out[stride_out * (nn - 1 - col)] += r1 - r0;
979 else
980 out[stride_out * (nn - 1 - col)] += r0 - r1;
981 }
982 else
983 {
984 out[stride_out * col] = r0 + r1;
985 if (quantity == EvaluatorQuantity::gradient &&
986 transpose_matrix == false)
987 out[stride_out * (nn - 1 - col)] = r1 - r0;
988 else
989 out[stride_out * (nn - 1 - col)] = r0 - r1;
990 }
991 }
992 if (quantity == EvaluatorQuantity::value && transpose_matrix == true &&
993 nn % 2 == 1 && mm % 2 == 1 && mm > 3)
994 {
995 if (add)
996 out[stride_out * n_half] += matrix[m_half * offset + n_half] * xmid;
997 else
998 out[stride_out * n_half] = matrix[m_half * offset + n_half] * xmid;
999 }
1000 else if (transpose_matrix == true && nn % 2 == 1)
1001 {
1002 Number r0;
1003 if (m_half > 0)
1004 {
1005 r0 = matrix[n_half] * xp[0];
1006 for (int ind = 1; ind < m_half; ++ind)
1007 r0 += matrix[ind * offset + n_half] * xp[ind];
1008 }
1009 else
1010 r0 = Number();
1011 if (quantity != EvaluatorQuantity::gradient && mm % 2 == 1)
1012 r0 += matrix[m_half * offset + n_half] * xmid;
1013
1014 if (add)
1015 out[stride_out * n_half] += r0;
1016 else
1017 out[stride_out * n_half] = r0;
1018 }
1019 else if (transpose_matrix == false && nn % 2 == 1)
1020 {
1021 Number r0;
1022 if (m_half > 0)
1023 {
1024 if (quantity == EvaluatorQuantity::gradient)
1025 {
1026 r0 = matrix[n_half * offset] * xm[0];
1027 for (int ind = 1; ind < m_half; ++ind)
1028 r0 += matrix[n_half * offset + ind] * xm[ind];
1029 }
1030 else
1031 {
1032 r0 = matrix[n_half * offset] * xp[0];
1033 for (int ind = 1; ind < m_half; ++ind)
1034 r0 += matrix[n_half * offset + ind] * xp[ind];
1035 }
1036 }
1037 else
1038 r0 = Number();
1039
1040 if (quantity != EvaluatorQuantity::gradient && mm % 2 == 1)
1041 r0 += matrix[n_half * offset + m_half] * xmid;
1042
1043 if (add)
1044 out[stride_out * n_half] += r0;
1045 else
1046 out[stride_out * n_half] = r0;
1047 }
1048 }
1049
1050
1051
1056 template <EvaluatorVariant variant,
1057 EvaluatorQuantity quantity,
1058 bool transpose_matrix,
1059 bool add,
1060 bool consider_strides,
1061 typename Number,
1062 typename Number2>
1063 std::enable_if_t<(variant == evaluate_evenodd), void>
1064 apply_matrix_vector_product(const Number2 *matrix,
1065 const Number *in,
1066 Number *out,
1067 int n_rows,
1068 int n_columns,
1069 int stride_in,
1070 int stride_out)
1071 {
1072 apply_matrix_vector_product<evaluate_evenodd,
1073 quantity,
1074 0,
1075 0,
1076 consider_strides ? 0 : 1,
1077 consider_strides ? 0 : 1,
1078 transpose_matrix,
1079 add>(
1080 matrix, in, out, n_rows, n_columns, stride_in, stride_out);
1081 }
1082
1083
1084
1100 template <EvaluatorVariant variant,
1101 EvaluatorQuantity quantity,
1102 int n_rows,
1103 int n_columns,
1104 int stride_in,
1105 int stride_out,
1106 bool transpose_matrix,
1107 bool add,
1108 typename Number,
1109 typename Number2>
1110 std::enable_if_t<(variant == evaluate_symmetric_hierarchical), void>
1111 apply_matrix_vector_product(const Number2 *matrix,
1112 const Number *in,
1113 Number *out)
1114 {
1115 static_assert(n_rows > 0 && n_columns > 0,
1116 "Specialization requires n_rows, n_columns > 0");
1117
1118 constexpr bool evaluate_antisymmetric =
1119 (quantity == EvaluatorQuantity::gradient);
1120
1121 constexpr int mm = transpose_matrix ? n_rows : n_columns,
1122 nn = transpose_matrix ? n_columns : n_rows;
1123 constexpr int n_half = nn / 2;
1124 constexpr int m_half = mm / 2;
1125
1126 if (transpose_matrix)
1127 {
1128 std::array<Number, mm> x;
1129 for (unsigned int i = 0; i < mm; ++i)
1130 x[i] = in[stride_in * i];
1131 for (unsigned int col = 0; col < n_half; ++col)
1132 {
1133 Number r0, r1;
1134 if (m_half > 0)
1135 {
1136 r0 = matrix[col] * x[0];
1137 r1 = matrix[col + n_columns] * x[1];
1138 for (unsigned int ind = 1; ind < m_half; ++ind)
1139 {
1140 r0 += matrix[col + 2 * ind * n_columns] * x[2 * ind];
1141 r1 +=
1142 matrix[col + (2 * ind + 1) * n_columns] * x[2 * ind + 1];
1143 }
1144 }
1145 else
1146 r0 = r1 = Number();
1147 if (mm % 2 == 1)
1148 r0 += matrix[col + (mm - 1) * n_columns] * x[mm - 1];
1149 if (add)
1150 {
1151 out[stride_out * col] += r0 + r1;
1152 if (evaluate_antisymmetric)
1153 out[stride_out * (nn - 1 - col)] += r1 - r0;
1154 else
1155 out[stride_out * (nn - 1 - col)] += r0 - r1;
1156 }
1157 else
1158 {
1159 out[stride_out * col] = r0 + r1;
1160 if (evaluate_antisymmetric)
1161 out[stride_out * (nn - 1 - col)] = r1 - r0;
1162 else
1163 out[stride_out * (nn - 1 - col)] = r0 - r1;
1164 }
1165 }
1166 if (nn % 2 == 1)
1167 {
1168 Number r0;
1169 const unsigned int shift = evaluate_antisymmetric ? 1 : 0;
1170 if (m_half > 0)
1171 {
1172 r0 = matrix[n_half + shift * n_columns] * x[shift];
1173 for (unsigned int ind = 1; ind < m_half; ++ind)
1174 r0 += matrix[n_half + (2 * ind + shift) * n_columns] *
1175 x[2 * ind + shift];
1176 }
1177 else
1178 r0 = 0;
1179 if (!evaluate_antisymmetric && mm % 2 == 1)
1180 r0 += matrix[n_half + (mm - 1) * n_columns] * x[mm - 1];
1181 if (add)
1182 out[stride_out * n_half] += r0;
1183 else
1184 out[stride_out * n_half] = r0;
1185 }
1186 }
1187 else
1188 {
1189 std::array<Number, m_half + 1> xp, xm;
1190 for (int i = 0; i < m_half; ++i)
1191 if (!evaluate_antisymmetric)
1192 {
1193 xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1194 xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1195 }
1196 else
1197 {
1198 xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1199 xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1200 }
1201 if (mm % 2 == 1)
1202 xp[m_half] = in[stride_in * m_half];
1203 for (unsigned int col = 0; col < n_half; ++col)
1204 {
1205 Number r0, r1;
1206 if (m_half > 0)
1207 {
1208 r0 = matrix[2 * col * n_columns] * xp[0];
1209 r1 = matrix[(2 * col + 1) * n_columns] * xm[0];
1210 for (unsigned int ind = 1; ind < m_half; ++ind)
1211 {
1212 r0 += matrix[2 * col * n_columns + ind] * xp[ind];
1213 r1 += matrix[(2 * col + 1) * n_columns + ind] * xm[ind];
1214 }
1215 }
1216 else
1217 r0 = r1 = Number();
1218 if (mm % 2 == 1)
1219 {
1220 if (evaluate_antisymmetric)
1221 r1 += matrix[(2 * col + 1) * n_columns + m_half] * xp[m_half];
1222 else
1223 r0 += matrix[2 * col * n_columns + m_half] * xp[m_half];
1224 }
1225 if (add)
1226 {
1227 out[stride_out * (2 * col)] += r0;
1228 out[stride_out * (2 * col + 1)] += r1;
1229 }
1230 else
1231 {
1232 out[stride_out * (2 * col)] = r0;
1233 out[stride_out * (2 * col + 1)] = r1;
1234 }
1235 }
1236 if (nn % 2 == 1)
1237 {
1238 Number r0;
1239 if (m_half > 0)
1240 {
1241 r0 = matrix[(nn - 1) * n_columns] * xp[0];
1242 for (unsigned int ind = 1; ind < m_half; ++ind)
1243 r0 += matrix[(nn - 1) * n_columns + ind] * xp[ind];
1244 }
1245 else
1246 r0 = Number();
1247 if (mm % 2 == 1 && !evaluate_antisymmetric)
1248 r0 += matrix[(nn - 1) * n_columns + m_half] * xp[m_half];
1249 if (add)
1250 out[stride_out * (nn - 1)] += r0;
1251 else
1252 out[stride_out * (nn - 1)] = r0;
1253 }
1254 }
1255 }
1256
1257
1258
1281 template <EvaluatorVariant variant,
1282 int dim,
1283 int n_rows,
1284 int n_columns,
1285 typename Number,
1286 typename Number2 = Number>
1287 struct EvaluatorTensorProduct
1288 {
1289 static constexpr unsigned int n_rows_of_product =
1290 Utilities::pow(n_rows, dim);
1291 static constexpr unsigned int n_columns_of_product =
1292 Utilities::pow(n_columns, dim);
1293
1298 EvaluatorTensorProduct()
1299 : shape_values(nullptr)
1300 , shape_gradients(nullptr)
1301 , shape_hessians(nullptr)
1302 {}
1303
1307 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
1308 const AlignedVector<Number2> &shape_gradients,
1309 const AlignedVector<Number2> &shape_hessians,
1310 const unsigned int = 0,
1311 const unsigned int = 0)
1312 : shape_values(shape_values.begin())
1313 , shape_gradients(shape_gradients.begin())
1314 , shape_hessians(shape_hessians.begin())
1315 {
1316 if (variant == evaluate_evenodd)
1317 {
1318 if (!shape_values.empty())
1319 AssertDimension(shape_values.size(),
1320 n_rows * ((n_columns + 1) / 2));
1321 if (!shape_gradients.empty())
1322 AssertDimension(shape_gradients.size(),
1323 n_rows * ((n_columns + 1) / 2));
1324 if (!shape_hessians.empty())
1325 AssertDimension(shape_hessians.size(),
1326 n_rows * ((n_columns + 1) / 2));
1327 }
1328 else
1329 {
1330 Assert(shape_values.empty() ||
1331 shape_values.size() == n_rows * n_columns,
1332 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1333 Assert(shape_gradients.empty() ||
1334 shape_gradients.size() == n_rows * n_columns,
1335 ExcDimensionMismatch(shape_gradients.size(),
1336 n_rows * n_columns));
1337 Assert(shape_hessians.empty() ||
1338 shape_hessians.size() == n_rows * n_columns,
1339 ExcDimensionMismatch(shape_hessians.size(),
1340 n_rows * n_columns));
1341 }
1342 }
1343
1347 EvaluatorTensorProduct(const Number2 *shape_values,
1348 const Number2 *shape_gradients,
1349 const Number2 *shape_hessians,
1350 const unsigned int dummy1 = 0,
1351 const unsigned int dummy2 = 0)
1352 : shape_values(shape_values)
1353 , shape_gradients(shape_gradients)
1354 , shape_hessians(shape_hessians)
1355 {
1356 (void)dummy1;
1357 (void)dummy2;
1358 }
1359
1385 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1386 void
1387 values(const Number in[], Number out[]) const
1388 {
1389 constexpr EvaluatorQuantity value_type = EvaluatorQuantity::value;
1390 apply<direction, contract_over_rows, add, false, value_type, stride>(
1391 shape_values, in, out);
1392 }
1393
1399 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1400 void
1401 gradients(const Number in[], Number out[]) const
1402 {
1403 constexpr EvaluatorQuantity gradient_type =
1404 (variant == evaluate_general ? EvaluatorQuantity::value :
1405 EvaluatorQuantity::gradient);
1406 apply<direction, contract_over_rows, add, false, gradient_type, stride>(
1407 shape_gradients, in, out);
1408 }
1409
1415 template <int direction, bool contract_over_rows, bool add>
1416 void
1417 hessians(const Number in[], Number out[]) const
1418 {
1419 constexpr EvaluatorQuantity hessian_type =
1420 (((variant == evaluate_general) |
1421 (variant == evaluate_symmetric_hierarchical)) ?
1422 EvaluatorQuantity::value :
1423 EvaluatorQuantity::hessian);
1424 apply<direction, contract_over_rows, add, false, hessian_type>(
1425 shape_hessians, in, out);
1426 }
1427
1435 template <int direction, bool contract_over_rows, bool add>
1436 void
1437 values_one_line(const Number in[], Number out[]) const
1438 {
1439 Assert(shape_values != nullptr, ExcNotInitialized());
1440 apply<direction, contract_over_rows, add, true, EvaluatorQuantity::value>(
1441 shape_values, in, out);
1442 }
1443
1451 template <int direction, bool contract_over_rows, bool add>
1452 void
1453 gradients_one_line(const Number in[], Number out[]) const
1454 {
1455 Assert(shape_gradients != nullptr, ExcNotInitialized());
1456 constexpr EvaluatorQuantity gradient_type =
1457 (variant == evaluate_general ? EvaluatorQuantity::value :
1458 EvaluatorQuantity::gradient);
1459 apply<direction, contract_over_rows, add, true, gradient_type>(
1460 shape_gradients, in, out);
1461 }
1462
1470 template <int direction, bool contract_over_rows, bool add>
1471 void
1472 hessians_one_line(const Number in[], Number out[]) const
1473 {
1474 Assert(shape_hessians != nullptr, ExcNotInitialized());
1475 constexpr EvaluatorQuantity hessian_type =
1476 (((variant == evaluate_general) |
1477 (variant == evaluate_symmetric_hierarchical)) ?
1478 EvaluatorQuantity::value :
1479 EvaluatorQuantity::hessian);
1480 apply<direction, contract_over_rows, add, true, hessian_type>(
1481 shape_hessians, in, out);
1482 }
1483
1520 template <int direction,
1521 bool contract_over_rows,
1522 bool add,
1523 bool one_line = false,
1524 EvaluatorQuantity quantity = EvaluatorQuantity::value,
1525 int stride = 1>
1526 static void
1527 apply(const Number2 *DEAL_II_RESTRICT shape_data,
1528 const Number *in,
1529 Number *out);
1530
1531 private:
1532 const Number2 *shape_values;
1533 const Number2 *shape_gradients;
1534 const Number2 *shape_hessians;
1535 };
1536
1537
1538
1539 template <EvaluatorVariant variant,
1540 int dim,
1541 int n_rows,
1542 int n_columns,
1543 typename Number,
1544 typename Number2>
1545 template <int direction,
1546 bool contract_over_rows,
1547 bool add,
1548 bool one_line,
1549 EvaluatorQuantity quantity,
1550 int stride>
1551 inline void
1552 EvaluatorTensorProduct<variant, dim, n_rows, n_columns, Number, Number2>::
1553 apply(const Number2 *DEAL_II_RESTRICT shape_data,
1554 const Number *in,
1555 Number *out)
1556 {
1557 static_assert(one_line == false || direction == dim - 1,
1558 "Single-line evaluation only works for direction=dim-1.");
1559 Assert(shape_data != nullptr,
1560 ExcMessage(
1561 "The given array shape_data must not be the null pointer!"));
1562 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
1563 in != out,
1564 ExcMessage("In-place operation only supported for "
1565 "n_rows==n_columns or single-line interpolation"));
1566 AssertIndexRange(direction, dim);
1567 constexpr int mm = contract_over_rows ? n_rows : n_columns,
1568 nn = contract_over_rows ? n_columns : n_rows;
1569
1570 constexpr int stride_operation = Utilities::pow(n_columns, direction);
1571 constexpr int n_blocks1 = one_line ? 1 : stride_operation;
1572 constexpr int n_blocks2 =
1573 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1574
1575 constexpr int stride_in = !contract_over_rows ? stride : 1;
1576 constexpr int stride_out = contract_over_rows ? stride : 1;
1577 for (int i2 = 0; i2 < n_blocks2; ++i2)
1578 {
1579 for (int i1 = 0; i1 < n_blocks1; ++i1)
1580 {
1581 apply_matrix_vector_product<variant,
1582 quantity,
1583 n_rows,
1584 n_columns,
1585 stride_operation * stride_in,
1586 stride_operation * stride_out,
1587 contract_over_rows,
1588 add>(shape_data, in, out);
1589
1590 if (one_line == false)
1591 {
1592 in += stride_in;
1593 out += stride_out;
1594 }
1595 }
1596 if (one_line == false)
1597 {
1598 in += stride_operation * (mm - 1) * stride_in;
1599 out += stride_operation * (nn - 1) * stride_out;
1600 }
1601 }
1602 }
1603
1604
1605
1619 template <EvaluatorVariant variant,
1620 int dim,
1621 typename Number,
1622 typename Number2>
1623 struct EvaluatorTensorProduct<variant, dim, 0, 0, Number, Number2>
1624 {
1625 static constexpr unsigned int n_rows_of_product =
1626 numbers::invalid_unsigned_int;
1627 static constexpr unsigned int n_columns_of_product =
1628 numbers::invalid_unsigned_int;
1629
1634 EvaluatorTensorProduct()
1635 : shape_values(nullptr)
1636 , shape_gradients(nullptr)
1637 , shape_hessians(nullptr)
1638 , n_rows(numbers::invalid_unsigned_int)
1639 , n_columns(numbers::invalid_unsigned_int)
1640 {}
1641
1645 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
1646 const AlignedVector<Number2> &shape_gradients,
1647 const AlignedVector<Number2> &shape_hessians,
1648 const unsigned int n_rows = 0,
1649 const unsigned int n_columns = 0)
1650 : shape_values(shape_values.begin())
1651 , shape_gradients(shape_gradients.begin())
1652 , shape_hessians(shape_hessians.begin())
1653 , n_rows(n_rows)
1654 , n_columns(n_columns)
1655 {
1656 if (variant == evaluate_evenodd)
1657 {
1658 if (!shape_values.empty())
1659 AssertDimension(shape_values.size(),
1660 n_rows * ((n_columns + 1) / 2));
1661 if (!shape_gradients.empty())
1662 AssertDimension(shape_gradients.size(),
1663 n_rows * ((n_columns + 1) / 2));
1664 if (!shape_hessians.empty())
1665 AssertDimension(shape_hessians.size(),
1666 n_rows * ((n_columns + 1) / 2));
1667 }
1668 else
1669 {
1670 Assert(shape_values.empty() ||
1671 shape_values.size() == n_rows * n_columns,
1672 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1673 Assert(shape_gradients.empty() ||
1674 shape_gradients.size() == n_rows * n_columns,
1675 ExcDimensionMismatch(shape_gradients.size(),
1676 n_rows * n_columns));
1677 Assert(shape_hessians.empty() ||
1678 shape_hessians.size() == n_rows * n_columns,
1679 ExcDimensionMismatch(shape_hessians.size(),
1680 n_rows * n_columns));
1681 }
1682 }
1683
1687 EvaluatorTensorProduct(const Number2 *shape_values,
1688 const Number2 *shape_gradients,
1689 const Number2 *shape_hessians,
1690 const unsigned int n_rows = 0,
1691 const unsigned int n_columns = 0)
1692 : shape_values(shape_values)
1693 , shape_gradients(shape_gradients)
1694 , shape_hessians(shape_hessians)
1695 , n_rows(n_rows)
1696 , n_columns(n_columns)
1697 {}
1698
1699 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1700 void
1701 values(const Number *in, Number *out) const
1702 {
1703 constexpr EvaluatorQuantity value_type = EvaluatorQuantity::value;
1704 apply<direction, contract_over_rows, add, false, value_type, stride>(
1705 shape_values, in, out);
1706 }
1707
1708 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1709 void
1710 gradients(const Number *in, Number *out) const
1711 {
1712 constexpr EvaluatorQuantity gradient_type =
1713 (variant != evaluate_evenodd ? EvaluatorQuantity::value :
1714 EvaluatorQuantity::gradient);
1715 apply<direction, contract_over_rows, add, false, gradient_type, stride>(
1716 shape_gradients, in, out);
1717 }
1718
1719 template <int direction, bool contract_over_rows, bool add>
1720 void
1721 hessians(const Number *in, Number *out) const
1722 {
1723 constexpr EvaluatorQuantity hessian_type =
1724 (variant != evaluate_evenodd ? EvaluatorQuantity::value :
1725 EvaluatorQuantity::hessian);
1726 apply<direction, contract_over_rows, add, false, hessian_type>(
1727 shape_hessians, in, out);
1728 }
1729
1730 template <int direction, bool contract_over_rows, bool add>
1731 void
1732 values_one_line(const Number in[], Number out[]) const
1733 {
1734 Assert(shape_values != nullptr, ExcNotInitialized());
1735 apply<direction, contract_over_rows, add, true, EvaluatorQuantity::value>(
1736 shape_values, in, out);
1737 }
1738
1739 template <int direction, bool contract_over_rows, bool add>
1740 void
1741 gradients_one_line(const Number in[], Number out[]) const
1742 {
1743 Assert(shape_gradients != nullptr, ExcNotInitialized());
1744 constexpr EvaluatorQuantity gradient_type =
1745 (variant != evaluate_evenodd ? EvaluatorQuantity::value :
1746 EvaluatorQuantity::gradient);
1747 apply<direction, contract_over_rows, add, true, gradient_type>(
1748 shape_gradients, in, out);
1749 }
1750
1751 template <int direction, bool contract_over_rows, bool add>
1752 void
1753 hessians_one_line(const Number in[], Number out[]) const
1754 {
1755 Assert(shape_hessians != nullptr, ExcNotInitialized());
1756 constexpr EvaluatorQuantity hessian_type =
1757 (variant != evaluate_evenodd ? EvaluatorQuantity::value :
1758 EvaluatorQuantity::hessian);
1759 apply<direction, contract_over_rows, add, true, hessian_type>(
1760 shape_hessians, in, out);
1761 }
1762
1763 template <int direction,
1764 bool contract_over_rows,
1765 bool add,
1766 bool one_line = false,
1767 EvaluatorQuantity quantity = EvaluatorQuantity::value,
1768 int stride = 1>
1769 void
1770 apply(const Number2 *DEAL_II_RESTRICT shape_data,
1771 const Number *in,
1772 Number *out) const;
1773
1774 const Number2 *shape_values;
1775 const Number2 *shape_gradients;
1776 const Number2 *shape_hessians;
1777 const unsigned int n_rows;
1778 const unsigned int n_columns;
1779 };
1780
1781
1782
1783 template <EvaluatorVariant variant,
1784 int dim,
1785 typename Number,
1786 typename Number2>
1787 template <int direction,
1788 bool contract_over_rows,
1789 bool add,
1790 bool one_line,
1791 EvaluatorQuantity quantity,
1792 int stride>
1793 inline void
1794 EvaluatorTensorProduct<variant, dim, 0, 0, Number, Number2>::apply(
1795 const Number2 *DEAL_II_RESTRICT shape_data,
1796 const Number *in,
1797 Number *out) const
1798 {
1799 static_assert(one_line == false || direction == dim - 1,
1800 "Single-line evaluation only works for direction=dim-1.");
1801 Assert(shape_data != nullptr,
1802 ExcMessage(
1803 "The given array shape_data must not be the null pointer!"));
1804 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
1805 in != out,
1806 ExcMessage("In-place operation only supported for "
1807 "n_rows==n_columns or single-line interpolation"));
1808 AssertIndexRange(direction, dim);
1809 const int mm = contract_over_rows ? n_rows : n_columns,
1810 nn = contract_over_rows ? n_columns : n_rows;
1811
1812 const int stride_operation =
1813 direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
1814 const int n_blocks1 = one_line ? 1 : stride_operation;
1815 const int n_blocks2 = direction >= dim - 1 ?
1816 1 :
1817 Utilities::fixed_power<dim - direction - 1>(n_rows);
1818 Assert(n_rows <= 128, ExcNotImplemented());
1819
1820 constexpr int stride_in = !contract_over_rows ? stride : 1;
1821 constexpr int stride_out = contract_over_rows ? stride : 1;
1822 for (int i2 = 0; i2 < n_blocks2; ++i2)
1823 {
1824 for (int i1 = 0; i1 < n_blocks1; ++i1)
1825 {
1826 // the empty template case can only run the general evaluator or
1827 // evenodd
1828 constexpr EvaluatorVariant restricted_variant =
1829 variant == evaluate_evenodd ? evaluate_evenodd : evaluate_general;
1830 apply_matrix_vector_product<restricted_variant,
1831 quantity,
1832 contract_over_rows,
1833 add,
1834 (direction != 0 || stride != 1)>(
1835 shape_data,
1836 in,
1837 out,
1838 n_rows,
1839 n_columns,
1840 stride_operation * stride_in,
1841 stride_operation * stride_out);
1842
1843 if (one_line == false)
1844 {
1845 in += stride_in;
1846 out += stride_out;
1847 }
1848 }
1849 if (one_line == false)
1850 {
1851 in += stride_operation * (mm - 1) * stride_in;
1852 out += stride_operation * (nn - 1) * stride_out;
1853 }
1854 }
1855 }
1856
1857
1858
1859 template <int dim,
1860 int fe_degree,
1861 int n_q_points_1d,
1862 bool contract_over_rows,
1863 bool symmetric_evaluate = true>
1864 struct EvaluatorTensorProductAnisotropic
1865 {
1866 template <int direction,
1867 int stride = 1,
1868 typename Number = double,
1869 typename Number2 = double>
1870 static void
1871 normal(const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
1872 const Number *in,
1873 Number *out,
1874 const bool add_into_result = false,
1875 const int subface_index_1d = 0)
1876 {
1877 AssertIndexRange(direction, dim);
1878 AssertDimension(fe_degree, data.fe_degree);
1879 AssertDimension(n_q_points_1d, data.n_q_points_1d);
1880 constexpr int n_rows = fe_degree + 1;
1881 constexpr int n_columns = n_q_points_1d;
1882 constexpr int mm = contract_over_rows ? n_rows : n_columns;
1883 constexpr int nn = contract_over_rows ? n_columns : n_rows;
1884 const Number2 *shape_data =
1885 symmetric_evaluate ?
1886 data.shape_values_eo.data() :
1887 data.values_within_subface[subface_index_1d].data();
1888 Assert(shape_data != nullptr, ExcNotInitialized());
1889 Assert(contract_over_rows == false || !add_into_result,
1890 ExcMessage("Cannot add into result if contract_over_rows = true"));
1891
1892 constexpr int n_blocks1 = Utilities::pow(fe_degree, direction);
1893 constexpr int n_blocks2 = Utilities::pow(fe_degree, dim - direction - 1);
1894 constexpr int stride_in = contract_over_rows ? 1 : stride;
1895 constexpr int stride_out = contract_over_rows ? stride : 1;
1896 constexpr EvaluatorVariant variant =
1897 symmetric_evaluate ? evaluate_evenodd : evaluate_general;
1898
1899 for (int i2 = 0; i2 < n_blocks2; ++i2)
1900 {
1901 for (int i1 = 0; i1 < n_blocks1; ++i1)
1902 {
1903 if (contract_over_rows == false && add_into_result)
1904 apply_matrix_vector_product<variant,
1905 EvaluatorQuantity::value,
1906 n_rows,
1907 n_columns,
1908 n_blocks1 * stride_in,
1909 n_blocks1 * stride_out,
1910 contract_over_rows,
1911 true>(shape_data, in, out);
1912 else
1913 apply_matrix_vector_product<variant,
1914 EvaluatorQuantity::value,
1915 n_rows,
1916 n_columns,
1917 n_blocks1 * stride_in,
1918 n_blocks1 * stride_out,
1919 contract_over_rows,
1920 false>(shape_data, in, out);
1921
1922 in += stride_in;
1923 out += stride_out;
1924 }
1925 in += n_blocks1 * (mm - 1) * stride_in;
1926 out += n_blocks1 * (nn - 1) * stride_out;
1927 }
1928 }
1929
1930 template <int direction,
1931 int normal_direction,
1932 int stride = 1,
1933 typename Number = double,
1934 typename Number2 = double>
1935 static void
1936 tangential(const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
1937 const Number *in,
1938 Number *out,
1939 const int subface_index_1d = 0)
1940 {
1941 AssertIndexRange(direction, dim);
1942 AssertDimension(fe_degree - 1, data.fe_degree);
1943 AssertDimension(n_q_points_1d, data.n_q_points_1d);
1944 static_assert(direction != normal_direction,
1945 "Cannot interpolate tangentially in normal direction");
1946
1947 constexpr int n_rows = std::max(fe_degree, 0);
1948 constexpr int n_columns = n_q_points_1d;
1949 const Number2 *shape_data =
1950 symmetric_evaluate ?
1951 data.shape_values_eo.data() :
1952 data.values_within_subface[subface_index_1d].data();
1953 Assert(shape_data != nullptr, ExcNotInitialized());
1954
1955 constexpr int n_blocks1 =
1956 (direction > normal_direction) ?
1957 Utilities::pow(n_q_points_1d, direction) :
1958 (direction > 0 ?
1959 (Utilities::pow(fe_degree, direction - 1) * n_q_points_1d) :
1960 1);
1961 constexpr int n_blocks2 =
1962 (direction > normal_direction) ?
1963 Utilities::pow(fe_degree, dim - 1 - direction) :
1964 ((direction + 1 < dim) ?
1965 (Utilities::pow(fe_degree, dim - 2 - direction) * n_q_points_1d) :
1966 1);
1967
1968 constexpr EvaluatorVariant variant =
1969 symmetric_evaluate ? evaluate_evenodd : evaluate_general;
1970
1971 // Since we may perform an in-place interpolation, we must run the step
1972 // expanding the size of the basis backward ('contract_over_rows' aka
1973 // 'evaluate' case), so shift the pointers and decrement during the loop
1974 if (contract_over_rows)
1975 {
1976 in += (n_blocks2 - 1) * n_blocks1 * n_rows + n_blocks1 - 1;
1977 out +=
1978 stride * ((n_blocks2 - 1) * n_blocks1 * n_columns + n_blocks1 - 1);
1979 for (int i2 = 0; i2 < n_blocks2; ++i2)
1980 {
1981 for (int i1 = 0; i1 < n_blocks1; ++i1)
1982 {
1983 apply_matrix_vector_product<variant,
1984 EvaluatorQuantity::value,
1985 n_rows,
1986 n_columns,
1987 n_blocks1,
1988 n_blocks1 * stride,
1989 true,
1990 false>(shape_data, in, out);
1991
1992 --in;
1993 out -= stride;
1994 }
1995 in -= n_blocks1 * (n_rows - 1);
1996 out -= n_blocks1 * (n_columns - 1) * stride;
1997 }
1998 }
1999 else
2000 {
2001 for (int i2 = 0; i2 < n_blocks2; ++i2)
2002 {
2003 for (int i1 = 0; i1 < n_blocks1; ++i1)
2004 {
2005 apply_matrix_vector_product<variant,
2006 EvaluatorQuantity::value,
2007 n_rows,
2008 n_columns,
2009 n_blocks1 * stride,
2010 n_blocks1,
2011 false,
2012 false>(shape_data, in, out);
2013
2014 in += stride;
2015 ++out;
2016 }
2017 in += n_blocks1 * (n_columns - 1) * stride;
2018 out += n_blocks1 * (n_rows - 1);
2019 }
2020 }
2021 }
2022 };
2023
2024
2025
2071 template <int n_rows_template,
2072 int stride_template,
2073 bool contract_onto_face,
2074 bool add,
2075 int max_derivative,
2076 typename Number,
2077 typename Number2>
2078 inline std::enable_if_t<contract_onto_face, void>
2079 interpolate_to_face(const Number2 *shape_values,
2080 const std::array<int, 2> &n_blocks,
2081 const std::array<int, 2> &steps,
2082 const Number *input,
2083 Number *DEAL_II_RESTRICT output,
2084 const int n_rows_runtime = 0,
2085 const int stride_runtime = 1)
2086 {
2087 const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
2088 const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
2089
2090 Number *output1 = output + n_blocks[0] * n_blocks[1];
2091 Number *output2 = output1 + n_blocks[0] * n_blocks[1];
2092 for (int i2 = 0; i2 < n_blocks[1]; ++i2)
2093 {
2094 for (int i1 = 0; i1 < n_blocks[0]; ++i1)
2095 {
2096 Number res0 = shape_values[0] * input[0];
2097 Number res1, res2;
2098 if (max_derivative > 0)
2099 res1 = shape_values[n_rows] * input[0];
2100 if (max_derivative > 1)
2101 res2 = shape_values[2 * n_rows] * input[0];
2102 for (int ind = 1; ind < n_rows; ++ind)
2103 {
2104 res0 += shape_values[ind] * input[stride * ind];
2105 if (max_derivative > 0)
2106 res1 += shape_values[ind + n_rows] * input[stride * ind];
2107 if (max_derivative > 1)
2108 res2 += shape_values[ind + 2 * n_rows] * input[stride * ind];
2109 }
2110 if (add)
2111 {
2112 output[i1] += res0;
2113 if (max_derivative > 0)
2114 output1[i1] += res1;
2115 if (max_derivative > 1)
2116 output2[i2] += res2;
2117 }
2118 else
2119 {
2120 output[i1] = res0;
2121 if (max_derivative > 0)
2122 output1[i1] = res1;
2123 if (max_derivative > 1)
2124 output2[i1] = res2;
2125 }
2126 input += steps[0];
2127 }
2128 output += n_blocks[0];
2129 if (max_derivative > 0)
2130 output1 += n_blocks[0];
2131 if (max_derivative > 1)
2132 output2 += n_blocks[0];
2133 input += steps[1];
2134 }
2135 }
2136
2137
2138
2146 constexpr bool
2147 use_collocation_evaluation(const unsigned int fe_degree,
2148 const unsigned int n_q_points_1d)
2149 {
2150 return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
2151 (n_q_points_1d <= 3 * fe_degree / 2 + 1);
2152 }
2153
2154
2155
2161 template <int n_rows_template,
2162 int stride_template,
2163 bool contract_onto_face,
2164 bool add,
2165 int max_derivative,
2166 typename Number,
2167 typename Number2>
2168 inline std::enable_if_t<!contract_onto_face, void>
2169 interpolate_to_face(const Number2 *shape_values,
2170 const std::array<int, 2> &n_blocks,
2171 const std::array<int, 2> &steps,
2172 const Number *input,
2173 Number *DEAL_II_RESTRICT output,
2174 const int n_rows_runtime = 0,
2175 const int stride_runtime = 1)
2176 {
2177 const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
2178 const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
2179
2180 const Number *input1 = input + n_blocks[0] * n_blocks[1];
2181 const Number *input2 = input1 + n_blocks[0] * n_blocks[1];
2182 for (int i2 = 0; i2 < n_blocks[1]; ++i2)
2183 {
2184 for (int i1 = 0; i1 < n_blocks[0]; ++i1)
2185 {
2186 const Number in = input[i1];
2187 Number in1, in2;
2188 if (max_derivative > 0)
2189 in1 = input1[i1];
2190 if (max_derivative > 1)
2191 in2 = input2[i1];
2192 for (int col = 0; col < n_rows; ++col)
2193 {
2194 Number result =
2195 add ? (output[col * stride] + shape_values[col] * in) :
2196 (shape_values[col] * in);
2197 if (max_derivative > 0)
2198 result += shape_values[col + n_rows] * in1;
2199 if (max_derivative > 1)
2200 result += shape_values[col + 2 * n_rows] * in2;
2201
2202 output[col * stride] = result;
2203 }
2204 output += steps[0];
2205 }
2206 input += n_blocks[0];
2207 if (max_derivative > 0)
2208 input1 += n_blocks[0];
2209 if (max_derivative > 1)
2210 input2 += n_blocks[0];
2211 output += steps[1];
2212 }
2213 }
2214
2215 template <int dim, int n_points_1d_template, typename Number>
2216 inline void
2217 weight_fe_q_dofs_by_entity(const Number *weights,
2218 const unsigned int n_components,
2219 const int n_points_1d_non_template,
2220 Number *data)
2221 {
2222 const int n_points_1d = n_points_1d_template != -1 ?
2223 n_points_1d_template :
2224 n_points_1d_non_template;
2225
2226 Assert(n_points_1d > 0, ExcNotImplemented());
2227 Assert(n_points_1d < 100, ExcNotImplemented());
2228
2229 unsigned int compressed_index[100];
2230 compressed_index[0] = 0;
2231 for (int i = 1; i < n_points_1d - 1; ++i)
2232 compressed_index[i] = 1;
2233 compressed_index[n_points_1d - 1] = 2;
2234
2235 for (unsigned int c = 0; c < n_components; ++c)
2236 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2237 for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
2238 {
2239 const unsigned int shift =
2240 9 * compressed_index[k] + 3 * compressed_index[j];
2241 data[0] *= weights[shift];
2242 // loop bound as int avoids compiler warnings in case n_points_1d
2243 // == 1 (polynomial degree 0)
2244 const Number weight = weights[shift + 1];
2245 for (int i = 1; i < n_points_1d - 1; ++i)
2246 data[i] *= weight;
2247 data[n_points_1d - 1] *= weights[shift + 2];
2248 data += n_points_1d;
2249 }
2250 }
2251
2252
2253 template <int dim, int n_points_1d_template, typename Number>
2254 inline void
2255 weight_fe_q_dofs_by_entity_shifted(const Number *weights,
2256 const unsigned int n_components,
2257 const int n_points_1d_non_template,
2258 Number *data)
2259 {
2260 const int n_points_1d = n_points_1d_template != -1 ?
2261 n_points_1d_template :
2262 n_points_1d_non_template;
2263
2264 Assert((n_points_1d % 2) == 1,
2265 ExcMessage("The function can only with add number of points"));
2266 Assert(n_points_1d > 0, ExcNotImplemented());
2267 Assert(n_points_1d < 100, ExcNotImplemented());
2268
2269 const unsigned int n_inside_1d = n_points_1d / 2;
2270
2271 unsigned int compressed_index[100];
2272
2273 unsigned int c = 0;
2274 for (int i = 0; i < n_inside_1d; ++i)
2275 compressed_index[c++] = 0;
2276 compressed_index[c++] = 1;
2277 for (int i = 0; i < n_inside_1d; ++i)
2278 compressed_index[c++] = 2;
2279
2280 for (unsigned int c = 0; c < n_components; ++c)
2281 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2282 for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
2283 {
2284 const unsigned int shift =
2285 9 * compressed_index[k] + 3 * compressed_index[j];
2286
2287 unsigned int c = 0;
2288 const Number weight1 = weights[shift];
2289 for (int i = 0; i < n_inside_1d; ++i)
2290 data[c++] *= weight1;
2291 data[c++] *= weights[shift + 1];
2292 const Number weight2 = weights[shift + 2];
2293 for (int i = 0; i < n_inside_1d; ++i)
2294 data[c++] *= weight2;
2295 data += n_points_1d;
2296 }
2297 }
2298
2299
2300 template <int dim, int n_points_1d_template, typename Number>
2301 inline bool
2302 compute_weights_fe_q_dofs_by_entity(const Number *data,
2303 const unsigned int n_components,
2304 const int n_points_1d_non_template,
2305 Number *weights)
2306 {
2307 const int n_points_1d = n_points_1d_template != -1 ?
2308 n_points_1d_template :
2309 n_points_1d_non_template;
2310
2311 Assert(n_points_1d > 0, ExcNotImplemented());
2312 Assert(n_points_1d < 100, ExcNotImplemented());
2313
2314 unsigned int compressed_index[100];
2315 compressed_index[0] = 0;
2316 for (int i = 1; i < n_points_1d - 1; ++i)
2317 compressed_index[i] = 1;
2318 compressed_index[n_points_1d - 1] = 2;
2319
2320 // Insert the number data into a storage position for weight,
2321 // ensuring that the array has either not been touched before
2322 // or the previous content is the same. In case the previous
2323 // content has a different value, we exit this function and
2324 // signal to outer functions that the compression was not possible.
2325 const auto check_and_set = [](Number &weight, const Number &data) {
2326 if (weight == Number(-1.0) || weight == data)
2327 {
2328 weight = data;
2329 return true; // success for the entry
2330 }
2331
2332 return false; // failure for the entry
2333 };
2334
2335 for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
2336 weights[c] = Number(-1.0);
2337
2338 for (unsigned int c = 0; c < n_components; ++c)
2339 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2340 for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
2341 ++j, data += n_points_1d)
2342 {
2343 const unsigned int shift =
2344 9 * compressed_index[k] + 3 * compressed_index[j];
2345
2346 if (!check_and_set(weights[shift], data[0]))
2347 return false; // failure
2348
2349 for (int i = 1; i < n_points_1d - 1; ++i)
2350 if (!check_and_set(weights[shift + 1], data[i]))
2351 return false; // failure
2352
2353 if (!check_and_set(weights[shift + 2], data[n_points_1d - 1]))
2354 return false; // failure
2355 }
2356
2357 return true; // success
2358 }
2359
2360
2361 template <int dim, int n_points_1d_template, typename Number>
2362 inline bool
2363 compute_weights_fe_q_dofs_by_entity_shifted(
2364 const Number *data,
2365 const unsigned int n_components,
2366 const int n_points_1d_non_template,
2367 Number *weights)
2368 {
2369 const int n_points_1d = n_points_1d_template != -1 ?
2370 n_points_1d_template :
2371 n_points_1d_non_template;
2372
2373 Assert((n_points_1d % 2) == 1,
2374 ExcMessage("The function can only with add number of points"));
2375 Assert(n_points_1d > 0, ExcNotImplemented());
2376 Assert(n_points_1d < 100, ExcNotImplemented());
2377
2378 const unsigned int n_inside_1d = n_points_1d / 2;
2379
2380 unsigned int compressed_index[100];
2381
2382 unsigned int c = 0;
2383 for (int i = 0; i < n_inside_1d; ++i)
2384 compressed_index[c++] = 0;
2385 compressed_index[c++] = 1;
2386 for (int i = 0; i < n_inside_1d; ++i)
2387 compressed_index[c++] = 2;
2388
2389 // Insert the number data into a storage position for weight,
2390 // ensuring that the array has either not been touched before
2391 // or the previous content is the same. In case the previous
2392 // content has a different value, we exit this function and
2393 // signal to outer functions that the compression was not possible.
2394 const auto check_and_set = [](Number &weight, const Number &data) {
2395 if (weight == Number(-1.0) || weight == data)
2396 {
2397 weight = data;
2398 return true; // success for the entry
2399 }
2400
2401 return false; // failure for the entry
2402 };
2403
2404 for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
2405 weights[c] = Number(-1.0);
2406
2407 for (unsigned int comp = 0; comp < n_components; ++comp)
2408 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2409 for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
2410 ++j, data += n_points_1d)
2411 {
2412 const unsigned int shift =
2413 9 * compressed_index[k] + 3 * compressed_index[j];
2414
2415 unsigned int c = 0;
2416
2417 for (int i = 0; i < n_inside_1d; ++i)
2418 if (!check_and_set(weights[shift], data[c++]))
2419 return false; // failure
2420
2421 if (!check_and_set(weights[shift + 1], data[c++]))
2422 return false; // failure
2423
2424 for (int i = 0; i < n_inside_1d; ++i)
2425 if (!check_and_set(weights[shift + 2], data[c++]))
2426 return false; // failure
2427 }
2428
2429 return true; // success
2430 }
2431
2432
2433} // end of namespace internal
2434
2435
2436DEAL_II_NAMESPACE_CLOSE
2437
2438#endif
aligned_vector.h
AlignedVector::data
pointer data()
DEAL_II_ALWAYS_INLINE
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
DEAL_II_RESTRICT
#define DEAL_II_RESTRICT
Definition config.h:110
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1638
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1985
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:2053
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
Utilities::fixed_power
constexpr T fixed_power(const T t)
Definition utilities.h:942
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
internal::use_collocation_evaluation
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
Definition tensor_product_kernels.h:2147
internal::weight_fe_q_dofs_by_entity
void weight_fe_q_dofs_by_entity(const Number *weights, const unsigned int n_components, const int n_points_1d_non_template, Number *data)
Definition tensor_product_kernels.h:2217
internal::apply_matrix_vector_product
std::enable_if_t<(variant==evaluate_general), void > apply_matrix_vector_product(const Number2 *matrix, const Number *in, Number *out)
Definition tensor_product_kernels.h:118
internal::weight_fe_q_dofs_by_entity_shifted
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)
Definition tensor_product_kernels.h:2255
internal::interpolate_to_face
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)
Definition tensor_product_kernels.h:2079
internal::compute_weights_fe_q_dofs_by_entity
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)
Definition tensor_product_kernels.h:2302
internal::compute_weights_fe_q_dofs_by_entity_shifted
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)
Definition tensor_product_kernels.h:2363
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:220
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6960
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6943
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6927
shape_info.h
internal::EvaluatorTensorProductAnisotropic::normal
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)
Definition tensor_product_kernels.h:1871
internal::EvaluatorTensorProductAnisotropic::tangential
static void tangential(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const Number *in, Number *out, const int subface_index_1d=0)
Definition tensor_product_kernels.h:1936
internal::EvaluatorTensorProduct< variant, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
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)
Definition tensor_product_kernels.h:1645
internal::EvaluatorTensorProduct< variant, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
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)
Definition tensor_product_kernels.h:1687
internal::EvaluatorTensorProduct::values
void values(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1387
internal::EvaluatorTensorProduct::EvaluatorTensorProduct
EvaluatorTensorProduct(const Number2 *shape_values, const Number2 *shape_gradients, const Number2 *shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition tensor_product_kernels.h:1347
internal::EvaluatorTensorProduct::values_one_line
void values_one_line(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1437
internal::EvaluatorTensorProduct::gradients
void gradients(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1401
internal::EvaluatorTensorProduct::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int=0, const unsigned int=0)
Definition tensor_product_kernels.h:1307
internal::EvaluatorTensorProduct::n_rows_of_product
static constexpr unsigned int n_rows_of_product
Definition tensor_product_kernels.h:1289
internal::EvaluatorTensorProduct::hessians
void hessians(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1417
internal::EvaluatorTensorProduct::apply
static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out)
Definition tensor_product_kernels.h:1553
internal::EvaluatorTensorProduct::n_columns_of_product
static constexpr unsigned int n_columns_of_product
Definition tensor_product_kernels.h:1291
internal::EvaluatorTensorProduct::gradients_one_line
void gradients_one_line(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1453
internal::EvaluatorTensorProduct::hessians_one_line
void hessians_one_line(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1472
internal::MatrixFreeFunctions::UnivariateShapeData::values_within_subface
std::array< AlignedVector< Number >, 2 > values_within_subface
Definition shape_info.h:314
Generated by doxygen 1.11.0