Reference documentation for deal.II version GIT relicensing-224-gc660c0d696 2024-03-28 18:40:02+00:00
\(\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
23
25
26
28
29
30
31namespace internal
32{
68
69
70
75 {
79 value,
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];
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];
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 typename Number,
197 typename Number2>
198 std::enable_if_t<(variant == evaluate_general), void>
199 apply_matrix_vector_product(const Number2 *matrix,
200 const Number *in,
201 Number *out,
202 const int n_rows,
203 const int n_columns,
204 const int stride_in,
205 const int stride_out)
206 {
207 const int mm = transpose_matrix ? n_rows : n_columns,
208 nn = transpose_matrix ? n_columns : n_rows;
209 Assert(n_rows > 0 && n_columns > 0,
210 ExcInternalError("Empty evaluation task!"));
211 Assert(n_rows > 0 && n_columns > 0,
212 ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
213 std::to_string(n_rows) + ", " +
214 std::to_string(n_columns) + " was passed!"));
215
216 static_assert(quantity == EvaluatorQuantity::value,
217 "This function should only use EvaluatorQuantity::value");
218
219 // specialization for n_rows = 2 that manually unrolls the innermost loop
220 // to make the operation perform better (not completely as good as the
221 // templated one, but much better than the generic version down below,
222 // because the loop over col can be more effectively unrolled by the
223 // compiler)
224 if (transpose_matrix && n_rows == 2)
225 {
226 const Number2 *matrix_1 = matrix + n_columns;
227 const Number x0 = in[0], x1 = in[stride_in];
228 for (int col = 0; col < nn; ++col)
229 {
230 const Number result = matrix[col] * x0 + matrix_1[col] * x1;
231 if (add)
232 out[stride_out * col] += result;
233 else
234 out[stride_out * col] = result;
235 }
236 }
237 else if (transpose_matrix && n_rows == 3)
238 {
239 const Number2 *matrix_1 = matrix + n_columns;
240 const Number2 *matrix_2 = matrix_1 + n_columns;
241 const Number x0 = in[0], x1 = in[stride_in], x2 = in[2 * stride_in];
242 for (int col = 0; col < nn; ++col)
243 {
244 const Number result =
245 matrix[col] * x0 + matrix_1[col] * x1 + matrix_2[col] * x2;
246 if (add)
247 out[stride_out * col] += result;
248 else
249 out[stride_out * col] = result;
250 }
251 }
252 else if (mm <= 128)
253 {
254 std::array<Number, 129> x;
255 for (int i = 0; i < mm; ++i)
256 x[i] = in[stride_in * i];
257
258 for (int col = 0; col < nn; ++col)
259 {
260 Number res0;
261 if (transpose_matrix == true)
262 {
263 res0 = matrix[col] * x[0];
264 for (int i = 1; i < mm; ++i)
265 res0 += matrix[i * n_columns + col] * x[i];
266 }
267 else
268 {
269 res0 = matrix[col * n_columns] * x[0];
270 for (int i = 1; i < mm; ++i)
271 res0 += matrix[col * n_columns + i] * x[i];
272 }
273 if (add)
274 out[stride_out * col] += res0;
275 else
276 out[stride_out * col] = res0;
277 }
278 }
279 else
280 {
281 Assert(in != out,
282 ExcNotImplemented("For large sizes, arrays may not overlap"));
283 for (int col = 0; col < nn; ++col)
284 {
285 Number res0;
286 if (transpose_matrix == true)
287 {
288 res0 = matrix[col] * in[0];
289 for (int i = 1; i < mm; ++i)
290 res0 += matrix[i * n_columns + col] * in[stride_in * i];
291 }
292 else
293 {
294 res0 = matrix[col * n_columns] * in[0];
295 for (int i = 1; i < mm; ++i)
296 res0 += matrix[col * n_columns + i] * in[stride_in * i];
297 }
298 if (add)
299 out[stride_out * col] += res0;
300 else
301 out[stride_out * col] = res0;
302 }
303 }
304 }
305
306
307
314 template <EvaluatorVariant variant,
315 EvaluatorQuantity quantity,
316 int n_rows,
317 int n_columns,
318 int stride_in,
319 int stride_out,
320 bool transpose_matrix,
321 bool add,
322 typename Number,
323 typename Number2>
324 std::enable_if_t<(variant == evaluate_symmetric), void>
325 apply_matrix_vector_product(const Number2 *matrix,
326 const Number *in,
327 Number *out)
328 {
329 // We can only statically assert that one argument is non-zero because
330 // face evaluation might instantiate some functions, so we need to use the
331 // run-time assert to verify that we do not end up involuntarily.
332 static_assert(n_rows > 0 || n_columns > 0,
333 "Specialization only for n_rows, n_columns > 0");
334 Assert(n_rows > 0 && n_columns > 0,
335 ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
336 std::to_string(n_rows) + ", " +
337 std::to_string(n_columns) + " was passed!"));
338
339 constexpr int mm = transpose_matrix ? n_rows : n_columns,
340 nn = transpose_matrix ? n_columns : n_rows;
341 constexpr int n_cols = nn / 2;
342 constexpr int mid = mm / 2;
343
344 std::array<Number, mm> x;
345 for (int i = 0; i < mm; ++i)
346 x[i] = in[stride_in * i];
347
348 if (quantity == EvaluatorQuantity::value)
349 {
350 // In this case, the 1d shape values read (sorted lexicographically,
351 // rows run over 1d dofs, columns over quadrature points):
352 // Q2 --> [ 0.687 0 -0.087 ]
353 // [ 0.4 1 0.4 ]
354 // [-0.087 0 0.687 ]
355 // Q3 --> [ 0.66 0.003 0.002 0.049 ]
356 // [ 0.521 1.005 -0.01 -0.230 ]
357 // [-0.230 -0.01 1.005 0.521 ]
358 // [ 0.049 0.002 0.003 0.66 ]
359 // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ]
360 // [ 0.608 1.059 0 0.039 0.176 ]
361 // [-0.409 -0.113 1 -0.113 -0.409 ]
362 // [ 0.176 0.039 0 1.059 0.608 ]
363 // [-0.032 -0.007 0 0.022 0.658 ]
364 //
365 // In these matrices, we want to use avoid computations involving
366 // zeros and ones and use the symmetry in entries starting from (1,1)
367 // forward and (N,N) backward, respectively to reduce the number of
368 // read operations.
369 for (int col = 0; col < n_cols; ++col)
370 {
371 Number2 val0, val1;
372 Number res0, res1;
373 if (transpose_matrix == true)
374 {
375 val0 = matrix[col];
376 val1 = matrix[nn - 1 - col];
377 }
378 else
379 {
380 val0 = matrix[col * n_columns];
381 val1 = matrix[(col + 1) * n_columns - 1];
382 }
383 if (mid > 0)
384 {
385 res0 = val0 * x[0];
386 res1 = val1 * x[0];
387 res0 += val1 * x[mm - 1];
388 res1 += val0 * x[mm - 1];
389 for (int ind = 1; ind < mid; ++ind)
390 {
391 if (transpose_matrix == true)
392 {
393 val0 = matrix[ind * n_columns + col];
394 val1 = matrix[ind * n_columns + nn - 1 - col];
395 }
396 else
397 {
398 val0 = matrix[col * n_columns + ind];
399 val1 = matrix[(col + 1) * n_columns - 1 - ind];
400 }
401 res0 += val0 * x[ind];
402 res1 += val1 * x[ind];
403 res0 += val1 * x[mm - 1 - ind];
404 res1 += val0 * x[mm - 1 - ind];
405 }
406 }
407 else
408 res0 = res1 = Number();
409 if (transpose_matrix == true)
410 {
411 if (mm % 2 == 1)
412 {
413 const Number tmp = matrix[mid * n_columns + col] * x[mid];
414 res0 += tmp;
415 res1 += tmp;
416 }
417 }
418 else
419 {
420 if (mm % 2 == 1 && nn % 2 == 0)
421 {
422 const Number tmp = matrix[col * n_columns + mid] * x[mid];
423 res0 += tmp;
424 res1 += tmp;
425 }
426 }
427 if (add)
428 {
429 out[stride_out * col] += res0;
430 out[stride_out * (nn - 1 - col)] += res1;
431 }
432 else
433 {
434 out[stride_out * col] = res0;
435 out[stride_out * (nn - 1 - col)] = res1;
436 }
437 }
438 if (transpose_matrix == true && nn % 2 == 1 && mm % 2 == 1)
439 {
440 if (add)
441 out[stride_out * n_cols] += x[mid];
442 else
443 out[stride_out * n_cols] = x[mid];
444 }
445 else if (transpose_matrix == true && nn % 2 == 1)
446 {
447 Number res0;
448 if (mid > 0)
449 {
450 res0 = matrix[n_cols] * (x[0] + x[mm - 1]);
451 for (int ind = 1; ind < mid; ++ind)
452 {
453 const Number2 val0 = matrix[ind * n_columns + n_cols];
454 res0 += val0 * (x[ind] + in[mm - 1 - ind]);
455 }
456 }
457 else
458 res0 = Number();
459 if (add)
460 out[stride_out * n_cols] += res0;
461 else
462 out[stride_out * n_cols] = res0;
463 }
464 else if (transpose_matrix == false && nn % 2 == 1)
465 {
466 Number res0;
467 if (mid > 0)
468 {
469 res0 = matrix[n_cols * n_columns] * (x[0] + x[mm - 1]);
470 for (int ind = 1; ind < mid; ++ind)
471 {
472 const Number2 val0 = matrix[n_cols * n_columns + ind];
473 res0 += val0 * (x[ind] + x[mm - 1 - ind]);
474 ;
475 }
476 if (mm % 2)
477 res0 += x[mid];
478 }
479 else
480 res0 = in[0];
481 if (add)
482 out[stride_out * n_cols] += res0;
483 else
484 out[stride_out * n_cols] = res0;
485 }
486 }
487 else if (quantity == EvaluatorQuantity::gradient)
488 {
489 // For the specialized loop used for gradient computations we again
490 // exploit symmetries according to the following entries (sorted
491 // lexicographically, rows run over 1d dofs, columns over quadrature
492 // points):
493 // Q2 --> [-2.549 -1 0.549 ]
494 // [ 3.098 0 -3.098 ]
495 // [-0.549 1 2.549 ]
496 // Q3 --> [-4.315 -1.03 0.5 -0.44 ]
497 // [ 6.07 -1.44 -2.97 2.196 ]
498 // [-2.196 2.97 1.44 -6.07 ]
499 // [ 0.44 -0.5 1.03 4.315 ]
500 // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ]
501 // [10.111 -2.76 -2.667 2.066 -2.306 ]
502 // [-5.688 5.773 0 -5.773 5.688 ]
503 // [ 2.306 -2.066 2.667 2.76 -10.111 ]
504 // [-0.413 0.353 -0.333 -0.353 0.413 ]
505 for (int col = 0; col < n_cols; ++col)
506 {
507 Number2 val0, val1;
508 Number res0, res1;
509 if (transpose_matrix == true)
510 {
511 val0 = matrix[col];
512 val1 = matrix[nn - 1 - col];
513 }
514 else
515 {
516 val0 = matrix[col * n_columns];
517 val1 = matrix[(nn - col - 1) * n_columns];
518 }
519 if (mid > 0)
520 {
521 res0 = val0 * x[0];
522 res1 = val1 * x[0];
523 res0 -= val1 * x[mm - 1];
524 res1 -= val0 * x[mm - 1];
525 for (int ind = 1; ind < mid; ++ind)
526 {
527 if (transpose_matrix == true)
528 {
529 val0 = matrix[ind * n_columns + col];
530 val1 = matrix[ind * n_columns + nn - 1 - col];
531 }
532 else
533 {
534 val0 = matrix[col * n_columns + ind];
535 val1 = matrix[(nn - col - 1) * n_columns + ind];
536 }
537 res0 += val0 * x[ind];
538 res1 += val1 * x[ind];
539 res0 -= val1 * x[mm - 1 - ind];
540 res1 -= val0 * x[mm - 1 - ind];
541 }
542 }
543 else
544 res0 = res1 = Number();
545 if (mm % 2 == 1)
546 {
547 if (transpose_matrix == true)
548 val0 = matrix[mid * n_columns + col];
549 else
550 val0 = matrix[col * n_columns + mid];
551 const Number tmp = val0 * x[mid];
552 res0 += tmp;
553 res1 -= tmp;
554 }
555 if (add)
556 {
557 out[stride_out * col] += res0;
558 out[stride_out * (nn - 1 - col)] += res1;
559 }
560 else
561 {
562 out[stride_out * col] = res0;
563 out[stride_out * (nn - 1 - col)] = res1;
564 }
565 }
566 if (nn % 2 == 1)
567 {
568 Number2 val0;
569 Number res0;
570 if (transpose_matrix == true)
571 val0 = matrix[n_cols];
572 else
573 val0 = matrix[n_cols * n_columns];
574 res0 = val0 * (x[0] - x[mm - 1]);
575 for (int ind = 1; ind < mid; ++ind)
576 {
577 if (transpose_matrix == true)
578 val0 = matrix[ind * n_columns + n_cols];
579 else
580 val0 = matrix[n_cols * n_columns + ind];
581 Number in1 = val0 * (x[ind] - x[mm - 1 - ind]);
582 res0 += in1;
583 }
584 if (add)
585 out[stride_out * n_cols] += res0;
586 else
587 out[stride_out * n_cols] = res0;
588 }
589 }
590 else
591 {
592 // Hessians are almost the same as values, apart from some missing '1'
593 // entries
594 for (int col = 0; col < n_cols; ++col)
595 {
596 Number2 val0, val1;
597 Number res0, res1;
598 if (transpose_matrix == true)
599 {
600 val0 = matrix[col];
601 val1 = matrix[nn - 1 - col];
602 }
603 else
604 {
605 val0 = matrix[col * n_columns];
606 val1 = matrix[(col + 1) * n_columns - 1];
607 }
608 if (mid > 0)
609 {
610 res0 = val0 * x[0];
611 res1 = val1 * x[0];
612 res0 += val1 * x[mm - 1];
613 res1 += val0 * x[mm - 1];
614 for (int ind = 1; ind < mid; ++ind)
615 {
616 if (transpose_matrix == true)
617 {
618 val0 = matrix[ind * n_columns + col];
619 val1 = matrix[ind * n_columns + nn - 1 - col];
620 }
621 else
622 {
623 val0 = matrix[col * n_columns + ind];
624 val1 = matrix[(col + 1) * n_columns - 1 - ind];
625 }
626 res0 += val0 * x[ind];
627 res1 += val1 * x[ind];
628 res0 += val1 * x[mm - 1 - ind];
629 res1 += val0 * x[mm - 1 - ind];
630 }
631 }
632 else
633 res0 = res1 = Number();
634 if (mm % 2 == 1)
635 {
636 if (transpose_matrix == true)
637 val0 = matrix[mid * n_columns + col];
638 else
639 val0 = matrix[col * n_columns + mid];
640 const Number tmp = val0 * x[mid];
641 res0 += tmp;
642 res1 += tmp;
643 }
644 if (add)
645 {
646 out[stride_out * col] += res0;
647 out[stride_out * (nn - 1 - col)] += res1;
648 }
649 else
650 {
651 out[stride_out * col] = res0;
652 out[stride_out * (nn - 1 - col)] = res1;
653 }
654 }
655 if (nn % 2 == 1)
656 {
657 Number2 val0;
658 Number res0;
659 if (transpose_matrix == true)
660 val0 = matrix[n_cols];
661 else
662 val0 = matrix[n_cols * n_columns];
663 if (mid > 0)
664 {
665 res0 = val0 * (x[0] + x[mm - 1]);
666 for (int ind = 1; ind < mid; ++ind)
667 {
668 if (transpose_matrix == true)
669 val0 = matrix[ind * n_columns + n_cols];
670 else
671 val0 = matrix[n_cols * n_columns + ind];
672 Number in1 = val0 * (x[ind] + x[mm - 1 - ind]);
673 res0 += in1;
674 }
675 }
676 else
677 res0 = Number();
678 if (mm % 2 == 1)
679 {
680 if (transpose_matrix == true)
681 val0 = matrix[mid * n_columns + n_cols];
682 else
683 val0 = matrix[n_cols * n_columns + mid];
684 res0 += val0 * x[mid];
685 }
686 if (add)
687 out[stride_out * n_cols] += res0;
688 else
689 out[stride_out * n_cols] = res0;
690 }
691 }
692 }
693
694
695
714 template <EvaluatorVariant variant,
715 EvaluatorQuantity quantity,
716 int n_rows_static,
717 int n_columns_static,
718 int stride_in_static,
719 int stride_out_static,
720 bool transpose_matrix,
721 bool add,
722 typename Number,
723 typename Number2>
724#ifndef DEBUG
726#endif
727 std::enable_if_t<(variant == evaluate_evenodd), void>
729 const Number *in,
730 Number *out,
731 int n_rows_runtime = 0,
732 int n_columns_runtime = 0,
733 int stride_in_runtime = 0,
734 int stride_out_runtime = 0)
735 {
736 static_assert(n_rows_static >= 0 && n_columns_static >= 0,
737 "Negative loop ranges are not allowed!");
738
739 const int n_rows = n_rows_static == 0 ? n_rows_runtime : n_rows_static;
740 const int n_columns =
741 n_rows_static == 0 ? n_columns_runtime : n_columns_static;
742 const int stride_in =
743 n_rows_static == 0 ? stride_in_runtime : stride_in_static;
744 const int stride_out =
745 n_rows_static == 0 ? stride_out_runtime : stride_out_static;
746
747 Assert(n_rows > 0 && n_columns > 0,
748 ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
749 std::to_string(n_rows) + ", " +
750 std::to_string(n_columns) + " was passed!"));
751
752 const int mm = transpose_matrix ? n_rows : n_columns,
753 nn = transpose_matrix ? n_columns : n_rows;
754 const int n_half = nn / 2;
755 const int m_half = mm / 2;
756
757 constexpr int array_length =
758 (n_rows_static == 0) ?
759 16 // for non-templated execution
760 :
761 (1 + (transpose_matrix ? n_rows_static : n_columns_static) / 2);
762 const int offset = (n_columns + 1) / 2;
763
764 Assert(m_half <= array_length, ExcNotImplemented());
765
766 std::array<Number, array_length> xp, xm;
767 for (int i = 0; i < m_half; ++i)
768 {
769 if (transpose_matrix == true && quantity == EvaluatorQuantity::gradient)
770 {
771 xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
772 xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
773 }
774 else
775 {
776 xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
777 xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
778 }
779 }
780 Number xmid = in[stride_in * m_half];
781 for (int col = 0; col < n_half; ++col)
782 {
783 Number r0, r1;
784 if (m_half > 0)
785 {
786 if (transpose_matrix == true)
787 {
788 r0 = matrix[col] * xp[0];
789 r1 = matrix[(n_rows - 1) * offset + col] * xm[0];
790 }
791 else
792 {
793 r0 = matrix[col * offset] * xp[0];
794 r1 = matrix[(n_rows - 1 - col) * offset] * xm[0];
795 }
796 for (int ind = 1; ind < m_half; ++ind)
797 {
798 if (transpose_matrix == true)
799 {
800 r0 += matrix[ind * offset + col] * xp[ind];
801 r1 += matrix[(n_rows - 1 - ind) * offset + col] * xm[ind];
802 }
803 else
804 {
805 r0 += matrix[col * offset + ind] * xp[ind];
806 r1 += matrix[(n_rows - 1 - col) * offset + ind] * xm[ind];
807 }
808 }
809 }
810 else
811 r0 = r1 = Number();
812 if (mm % 2 == 1 && transpose_matrix == true)
813 {
814 if (quantity == EvaluatorQuantity::gradient)
815 r1 += matrix[m_half * offset + col] * xmid;
816 else
817 r0 += matrix[m_half * offset + col] * xmid;
818 }
819 else if (mm % 2 == 1 &&
820 (nn % 2 == 0 || quantity != EvaluatorQuantity::value ||
821 mm == 3))
822 r0 += matrix[col * offset + m_half] * xmid;
823
824 if (add)
825 {
826 out[stride_out * col] += r0 + r1;
827 if (quantity == EvaluatorQuantity::gradient &&
828 transpose_matrix == false)
829 out[stride_out * (nn - 1 - col)] += r1 - r0;
830 else
831 out[stride_out * (nn - 1 - col)] += r0 - r1;
832 }
833 else
834 {
835 out[stride_out * col] = r0 + r1;
836 if (quantity == EvaluatorQuantity::gradient &&
837 transpose_matrix == false)
838 out[stride_out * (nn - 1 - col)] = r1 - r0;
839 else
840 out[stride_out * (nn - 1 - col)] = r0 - r1;
841 }
842 }
843 if (quantity == EvaluatorQuantity::value && transpose_matrix == true &&
844 nn % 2 == 1 && mm % 2 == 1 && mm > 3)
845 {
846 if (add)
847 out[stride_out * n_half] += matrix[m_half * offset + n_half] * xmid;
848 else
849 out[stride_out * n_half] = matrix[m_half * offset + n_half] * xmid;
850 }
851 else if (transpose_matrix == true && nn % 2 == 1)
852 {
853 Number r0;
854 if (m_half > 0)
855 {
856 r0 = matrix[n_half] * xp[0];
857 for (int ind = 1; ind < m_half; ++ind)
858 r0 += matrix[ind * offset + n_half] * xp[ind];
859 }
860 else
861 r0 = Number();
862 if (quantity != EvaluatorQuantity::gradient && mm % 2 == 1)
863 r0 += matrix[m_half * offset + n_half] * xmid;
864
865 if (add)
866 out[stride_out * n_half] += r0;
867 else
868 out[stride_out * n_half] = r0;
869 }
870 else if (transpose_matrix == false && nn % 2 == 1)
871 {
872 Number r0;
873 if (m_half > 0)
874 {
875 if (quantity == EvaluatorQuantity::gradient)
876 {
877 r0 = matrix[n_half * offset] * xm[0];
878 for (int ind = 1; ind < m_half; ++ind)
879 r0 += matrix[n_half * offset + ind] * xm[ind];
880 }
881 else
882 {
883 r0 = matrix[n_half * offset] * xp[0];
884 for (int ind = 1; ind < m_half; ++ind)
885 r0 += matrix[n_half * offset + ind] * xp[ind];
886 }
887 }
888 else
889 r0 = Number();
890
891 if (quantity != EvaluatorQuantity::gradient && mm % 2 == 1)
892 r0 += matrix[n_half * offset + m_half] * xmid;
893
894 if (add)
895 out[stride_out * n_half] += r0;
896 else
897 out[stride_out * n_half] = r0;
898 }
899 }
900
901
902
907 template <EvaluatorVariant variant,
908 EvaluatorQuantity quantity,
909 bool transpose_matrix,
910 bool add,
911 typename Number,
912 typename Number2>
913 std::enable_if_t<(variant == evaluate_evenodd), void>
914 apply_matrix_vector_product(const Number2 *matrix,
915 const Number *in,
916 Number *out,
917 int n_rows,
918 int n_columns,
919 int stride_in,
920 int stride_out)
921 {
923 quantity,
924 0,
925 0,
926 0,
927 0,
928 transpose_matrix,
929 add>(
930 matrix, in, out, n_rows, n_columns, stride_in, stride_out);
931 }
932
933
934
950 template <EvaluatorVariant variant,
951 EvaluatorQuantity quantity,
952 int n_rows,
953 int n_columns,
954 int stride_in,
955 int stride_out,
956 bool transpose_matrix,
957 bool add,
958 typename Number,
959 typename Number2>
960 std::enable_if_t<(variant == evaluate_symmetric_hierarchical), void>
961 apply_matrix_vector_product(const Number2 *matrix,
962 const Number *in,
963 Number *out)
964 {
965 static_assert(n_rows > 0 && n_columns > 0,
966 "Specialization requires n_rows, n_columns > 0");
967
968 constexpr bool evaluate_antisymmetric =
969 (quantity == EvaluatorQuantity::gradient);
970
971 constexpr int mm = transpose_matrix ? n_rows : n_columns,
972 nn = transpose_matrix ? n_columns : n_rows;
973 constexpr int n_half = nn / 2;
974 constexpr int m_half = mm / 2;
975
976 if (transpose_matrix)
977 {
978 std::array<Number, mm> x;
979 for (unsigned int i = 0; i < mm; ++i)
980 x[i] = in[stride_in * i];
981 for (unsigned int col = 0; col < n_half; ++col)
982 {
983 Number r0, r1;
984 if (m_half > 0)
985 {
986 r0 = matrix[col] * x[0];
987 r1 = matrix[col + n_columns] * x[1];
988 for (unsigned int ind = 1; ind < m_half; ++ind)
989 {
990 r0 += matrix[col + 2 * ind * n_columns] * x[2 * ind];
991 r1 +=
992 matrix[col + (2 * ind + 1) * n_columns] * x[2 * ind + 1];
993 }
994 }
995 else
996 r0 = r1 = Number();
997 if (mm % 2 == 1)
998 r0 += matrix[col + (mm - 1) * n_columns] * x[mm - 1];
999 if (add)
1000 {
1001 out[stride_out * col] += r0 + r1;
1002 if (evaluate_antisymmetric)
1003 out[stride_out * (nn - 1 - col)] += r1 - r0;
1004 else
1005 out[stride_out * (nn - 1 - col)] += r0 - r1;
1006 }
1007 else
1008 {
1009 out[stride_out * col] = r0 + r1;
1010 if (evaluate_antisymmetric)
1011 out[stride_out * (nn - 1 - col)] = r1 - r0;
1012 else
1013 out[stride_out * (nn - 1 - col)] = r0 - r1;
1014 }
1015 }
1016 if (nn % 2 == 1)
1017 {
1018 Number r0;
1019 const unsigned int shift = evaluate_antisymmetric ? 1 : 0;
1020 if (m_half > 0)
1021 {
1022 r0 = matrix[n_half + shift * n_columns] * x[shift];
1023 for (unsigned int ind = 1; ind < m_half; ++ind)
1024 r0 += matrix[n_half + (2 * ind + shift) * n_columns] *
1025 x[2 * ind + shift];
1026 }
1027 else
1028 r0 = 0;
1029 if (!evaluate_antisymmetric && mm % 2 == 1)
1030 r0 += matrix[n_half + (mm - 1) * n_columns] * x[mm - 1];
1031 if (add)
1032 out[stride_out * n_half] += r0;
1033 else
1034 out[stride_out * n_half] = r0;
1035 }
1036 }
1037 else
1038 {
1039 std::array<Number, m_half + 1> xp, xm;
1040 for (int i = 0; i < m_half; ++i)
1041 if (!evaluate_antisymmetric)
1042 {
1043 xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1044 xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1045 }
1046 else
1047 {
1048 xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1049 xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1050 }
1051 if (mm % 2 == 1)
1052 xp[m_half] = in[stride_in * m_half];
1053 for (unsigned int col = 0; col < n_half; ++col)
1054 {
1055 Number r0, r1;
1056 if (m_half > 0)
1057 {
1058 r0 = matrix[2 * col * n_columns] * xp[0];
1059 r1 = matrix[(2 * col + 1) * n_columns] * xm[0];
1060 for (unsigned int ind = 1; ind < m_half; ++ind)
1061 {
1062 r0 += matrix[2 * col * n_columns + ind] * xp[ind];
1063 r1 += matrix[(2 * col + 1) * n_columns + ind] * xm[ind];
1064 }
1065 }
1066 else
1067 r0 = r1 = Number();
1068 if (mm % 2 == 1)
1069 {
1070 if (evaluate_antisymmetric)
1071 r1 += matrix[(2 * col + 1) * n_columns + m_half] * xp[m_half];
1072 else
1073 r0 += matrix[2 * col * n_columns + m_half] * xp[m_half];
1074 }
1075 if (add)
1076 {
1077 out[stride_out * (2 * col)] += r0;
1078 out[stride_out * (2 * col + 1)] += r1;
1079 }
1080 else
1081 {
1082 out[stride_out * (2 * col)] = r0;
1083 out[stride_out * (2 * col + 1)] = r1;
1084 }
1085 }
1086 if (nn % 2 == 1)
1087 {
1088 Number r0;
1089 if (m_half > 0)
1090 {
1091 r0 = matrix[(nn - 1) * n_columns] * xp[0];
1092 for (unsigned int ind = 1; ind < m_half; ++ind)
1093 r0 += matrix[(nn - 1) * n_columns + ind] * xp[ind];
1094 }
1095 else
1096 r0 = Number();
1097 if (mm % 2 == 1 && !evaluate_antisymmetric)
1098 r0 += matrix[(nn - 1) * n_columns + m_half] * xp[m_half];
1099 if (add)
1100 out[stride_out * (nn - 1)] += r0;
1101 else
1102 out[stride_out * (nn - 1)] = r0;
1103 }
1104 }
1105 }
1106
1107
1108
1131 template <EvaluatorVariant variant,
1132 int dim,
1133 int n_rows,
1134 int n_columns,
1135 typename Number,
1136 typename Number2 = Number>
1138 {
1139 static constexpr unsigned int n_rows_of_product =
1140 Utilities::pow(n_rows, dim);
1141 static constexpr unsigned int n_columns_of_product =
1142 Utilities::pow(n_columns, dim);
1143
1149 : shape_values(nullptr)
1150 , shape_gradients(nullptr)
1151 , shape_hessians(nullptr)
1152 {}
1153
1160 const unsigned int = 0,
1161 const unsigned int = 0)
1162 : shape_values(shape_values.begin())
1165 {
1166 if (variant == evaluate_evenodd)
1167 {
1168 if (!shape_values.empty())
1170 n_rows * ((n_columns + 1) / 2));
1171 if (!shape_gradients.empty())
1173 n_rows * ((n_columns + 1) / 2));
1174 if (!shape_hessians.empty())
1176 n_rows * ((n_columns + 1) / 2));
1177 }
1178 else
1179 {
1180 Assert(shape_values.empty() ||
1181 shape_values.size() == n_rows * n_columns,
1182 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1183 Assert(shape_gradients.empty() ||
1184 shape_gradients.size() == n_rows * n_columns,
1186 n_rows * n_columns));
1187 Assert(shape_hessians.empty() ||
1188 shape_hessians.size() == n_rows * n_columns,
1190 n_rows * n_columns));
1191 }
1192 }
1193
1198 const Number2 *shape_gradients,
1199 const Number2 *shape_hessians,
1200 const unsigned int dummy1 = 0,
1201 const unsigned int dummy2 = 0)
1205 {
1206 (void)dummy1;
1207 (void)dummy2;
1208 }
1209
1235 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1236 void
1237 values(const Number in[], Number out[]) const
1238 {
1239 constexpr EvaluatorQuantity value_type = EvaluatorQuantity::value;
1240 apply<direction, contract_over_rows, add, false, value_type, stride>(
1241 shape_values, in, out);
1242 }
1243
1249 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1250 void
1251 gradients(const Number in[], Number out[]) const
1252 {
1253 constexpr EvaluatorQuantity gradient_type =
1256 apply<direction, contract_over_rows, add, false, gradient_type, stride>(
1257 shape_gradients, in, out);
1258 }
1259
1265 template <int direction, bool contract_over_rows, bool add>
1266 void
1267 hessians(const Number in[], Number out[]) const
1268 {
1269 constexpr EvaluatorQuantity hessian_type =
1270 (((variant == evaluate_general) |
1271 (variant == evaluate_symmetric_hierarchical)) ?
1274 apply<direction, contract_over_rows, add, false, hessian_type>(
1275 shape_hessians, in, out);
1276 }
1277
1285 template <int direction, bool contract_over_rows, bool add>
1286 void
1287 values_one_line(const Number in[], Number out[]) const
1288 {
1289 Assert(shape_values != nullptr, ExcNotInitialized());
1290 apply<direction, contract_over_rows, add, true, EvaluatorQuantity::value>(
1291 shape_values, in, out);
1292 }
1293
1301 template <int direction, bool contract_over_rows, bool add>
1302 void
1303 gradients_one_line(const Number in[], Number out[]) const
1304 {
1306 constexpr EvaluatorQuantity gradient_type =
1309 apply<direction, contract_over_rows, add, true, gradient_type>(
1310 shape_gradients, in, out);
1311 }
1312
1320 template <int direction, bool contract_over_rows, bool add>
1321 void
1322 hessians_one_line(const Number in[], Number out[]) const
1323 {
1325 constexpr EvaluatorQuantity hessian_type =
1326 (((variant == evaluate_general) |
1327 (variant == evaluate_symmetric_hierarchical)) ?
1330 apply<direction, contract_over_rows, add, true, hessian_type>(
1331 shape_hessians, in, out);
1332 }
1333
1370 template <int direction,
1371 bool contract_over_rows,
1372 bool add,
1373 bool one_line = false,
1375 int stride = 1>
1376 static void
1377 apply(const Number2 *DEAL_II_RESTRICT shape_data,
1378 const Number *in,
1379 Number *out);
1380
1381 private:
1382 const Number2 *shape_values;
1383 const Number2 *shape_gradients;
1384 const Number2 *shape_hessians;
1385 };
1386
1387
1388
1389 template <EvaluatorVariant variant,
1390 int dim,
1391 int n_rows,
1392 int n_columns,
1393 typename Number,
1394 typename Number2>
1395 template <int direction,
1396 bool contract_over_rows,
1397 bool add,
1398 bool one_line,
1399 EvaluatorQuantity quantity,
1400 int stride>
1401 inline void
1403 apply(const Number2 *DEAL_II_RESTRICT shape_data,
1404 const Number *in,
1405 Number *out)
1406 {
1407 static_assert(one_line == false || direction == dim - 1,
1408 "Single-line evaluation only works for direction=dim-1.");
1409 Assert(shape_data != nullptr,
1410 ExcMessage(
1411 "The given array shape_data must not be the null pointer!"));
1412 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
1413 in != out,
1414 ExcMessage("In-place operation only supported for "
1415 "n_rows==n_columns or single-line interpolation"));
1416 AssertIndexRange(direction, dim);
1417 constexpr int mm = contract_over_rows ? n_rows : n_columns,
1418 nn = contract_over_rows ? n_columns : n_rows;
1419
1420 constexpr int stride_operation = Utilities::pow(n_columns, direction);
1421 constexpr int n_blocks1 = one_line ? 1 : stride_operation;
1422 constexpr int n_blocks2 =
1423 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1424
1425 constexpr int stride_in = !contract_over_rows ? stride : 1;
1426 constexpr int stride_out = contract_over_rows ? stride : 1;
1427 for (int i2 = 0; i2 < n_blocks2; ++i2)
1428 {
1429 for (int i1 = 0; i1 < n_blocks1; ++i1)
1430 {
1432 quantity,
1433 n_rows,
1434 n_columns,
1435 stride_operation * stride_in,
1436 stride_operation * stride_out,
1437 contract_over_rows,
1438 add>(shape_data, in, out);
1439
1440 if (one_line == false)
1441 {
1442 in += stride_in;
1443 out += stride_out;
1444 }
1445 }
1446 if (one_line == false)
1447 {
1448 in += stride_operation * (mm - 1) * stride_in;
1449 out += stride_operation * (nn - 1) * stride_out;
1450 }
1451 }
1452 }
1453
1454
1455
1469 template <EvaluatorVariant variant,
1470 int dim,
1471 typename Number,
1472 typename Number2>
1473 struct EvaluatorTensorProduct<variant, dim, 0, 0, Number, Number2>
1474 {
1475 static constexpr unsigned int n_rows_of_product =
1477 static constexpr unsigned int n_columns_of_product =
1479
1485 : shape_values(nullptr)
1486 , shape_gradients(nullptr)
1487 , shape_hessians(nullptr)
1488 , n_rows(numbers::invalid_unsigned_int)
1489 , n_columns(numbers::invalid_unsigned_int)
1490 {}
1491
1498 const unsigned int n_rows = 0,
1499 const unsigned int n_columns = 0)
1500 : shape_values(shape_values.begin())
1503 , n_rows(n_rows)
1504 , n_columns(n_columns)
1505 {
1506 if (variant == evaluate_evenodd)
1507 {
1508 if (!shape_values.empty())
1510 n_rows * ((n_columns + 1) / 2));
1511 if (!shape_gradients.empty())
1513 n_rows * ((n_columns + 1) / 2));
1514 if (!shape_hessians.empty())
1516 n_rows * ((n_columns + 1) / 2));
1517 }
1518 else
1519 {
1520 Assert(shape_values.empty() ||
1521 shape_values.size() == n_rows * n_columns,
1522 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1523 Assert(shape_gradients.empty() ||
1524 shape_gradients.size() == n_rows * n_columns,
1526 n_rows * n_columns));
1527 Assert(shape_hessians.empty() ||
1528 shape_hessians.size() == n_rows * n_columns,
1530 n_rows * n_columns));
1531 }
1532 }
1533
1538 const Number2 *shape_gradients,
1539 const Number2 *shape_hessians,
1540 const unsigned int n_rows = 0,
1541 const unsigned int n_columns = 0)
1545 , n_rows(n_rows)
1546 , n_columns(n_columns)
1547 {}
1548
1549 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1550 void
1551 values(const Number *in, Number *out) const
1552 {
1553 constexpr EvaluatorQuantity value_type = EvaluatorQuantity::value;
1554 apply<direction, contract_over_rows, add, false, value_type, stride>(
1555 shape_values, in, out);
1556 }
1557
1558 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1559 void
1560 gradients(const Number *in, Number *out) const
1561 {
1562 constexpr EvaluatorQuantity gradient_type =
1565 apply<direction, contract_over_rows, add, false, gradient_type, stride>(
1566 shape_gradients, in, out);
1567 }
1568
1569 template <int direction, bool contract_over_rows, bool add>
1570 void
1571 hessians(const Number *in, Number *out) const
1572 {
1573 constexpr EvaluatorQuantity hessian_type =
1576 apply<direction, contract_over_rows, add, false, hessian_type>(
1577 shape_hessians, in, out);
1578 }
1579
1580 template <int direction, bool contract_over_rows, bool add>
1581 void
1582 values_one_line(const Number in[], Number out[]) const
1583 {
1584 Assert(shape_values != nullptr, ExcNotInitialized());
1585 apply<direction, contract_over_rows, add, true, EvaluatorQuantity::value>(
1586 shape_values, in, out);
1587 }
1588
1589 template <int direction, bool contract_over_rows, bool add>
1590 void
1591 gradients_one_line(const Number in[], Number out[]) const
1592 {
1594 constexpr EvaluatorQuantity gradient_type =
1597 apply<direction, contract_over_rows, add, true, gradient_type>(
1598 shape_gradients, in, out);
1599 }
1600
1601 template <int direction, bool contract_over_rows, bool add>
1602 void
1603 hessians_one_line(const Number in[], Number out[]) const
1604 {
1606 constexpr EvaluatorQuantity hessian_type =
1609 apply<direction, contract_over_rows, add, true, hessian_type>(
1610 shape_hessians, in, out);
1611 }
1612
1613 template <int direction,
1614 bool contract_over_rows,
1615 bool add,
1616 bool one_line = false,
1618 int stride = 1>
1619 void
1620 apply(const Number2 *DEAL_II_RESTRICT shape_data,
1621 const Number *in,
1622 Number *out) const;
1623
1624 const Number2 *shape_values;
1625 const Number2 *shape_gradients;
1626 const Number2 *shape_hessians;
1627 const unsigned int n_rows;
1628 const unsigned int n_columns;
1629 };
1630
1631
1632
1633 template <EvaluatorVariant variant,
1634 int dim,
1635 typename Number,
1636 typename Number2>
1637 template <int direction,
1638 bool contract_over_rows,
1639 bool add,
1640 bool one_line,
1641 EvaluatorQuantity quantity,
1642 int stride>
1643 inline void
1645 const Number2 *DEAL_II_RESTRICT shape_data,
1646 const Number *in,
1647 Number *out) const
1648 {
1649 static_assert(one_line == false || direction == dim - 1,
1650 "Single-line evaluation only works for direction=dim-1.");
1651 Assert(shape_data != nullptr,
1652 ExcMessage(
1653 "The given array shape_data must not be the null pointer!"));
1654 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
1655 in != out,
1656 ExcMessage("In-place operation only supported for "
1657 "n_rows==n_columns or single-line interpolation"));
1658 AssertIndexRange(direction, dim);
1659 const int mm = contract_over_rows ? n_rows : n_columns,
1660 nn = contract_over_rows ? n_columns : n_rows;
1661
1662 const int stride_operation =
1663 direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
1664 const int n_blocks1 = one_line ? 1 : stride_operation;
1665 const int n_blocks2 = direction >= dim - 1 ?
1666 1 :
1667 Utilities::fixed_power<dim - direction - 1>(n_rows);
1668 Assert(n_rows <= 128, ExcNotImplemented());
1669
1670 constexpr int stride_in = !contract_over_rows ? stride : 1;
1671 constexpr int stride_out = contract_over_rows ? stride : 1;
1672 for (int i2 = 0; i2 < n_blocks2; ++i2)
1673 {
1674 for (int i1 = 0; i1 < n_blocks1; ++i1)
1675 {
1676 // the empty template case can only run the general evaluator or
1677 // evenodd
1678 constexpr EvaluatorVariant restricted_variant =
1680 apply_matrix_vector_product<restricted_variant,
1681 quantity,
1682 contract_over_rows,
1683 add>(shape_data,
1684 in,
1685 out,
1686 n_rows,
1687 n_columns,
1688 stride_operation * stride_in,
1689 stride_operation * stride_out);
1690
1691 if (one_line == false)
1692 {
1693 in += stride_in;
1694 out += stride_out;
1695 }
1696 }
1697 if (one_line == false)
1698 {
1699 in += stride_operation * (mm - 1) * stride_in;
1700 out += stride_operation * (nn - 1) * stride_out;
1701 }
1702 }
1703 }
1704
1705
1706
1707 template <int dim,
1708 int fe_degree,
1709 int n_q_points_1d,
1710 bool contract_over_rows,
1711 bool symmetric_evaluate = true>
1713 {
1714 template <int direction,
1715 int stride = 1,
1716 typename Number = double,
1717 typename Number2 = double>
1718 static void
1720 const Number *in,
1721 Number *out,
1722 const bool add_into_result = false,
1723 const int subface_index_1d = 0)
1724 {
1725 AssertIndexRange(direction, dim);
1726 AssertDimension(fe_degree, data.fe_degree);
1727 AssertDimension(n_q_points_1d, data.n_q_points_1d);
1728 constexpr int n_rows = fe_degree + 1;
1729 constexpr int n_columns = n_q_points_1d;
1730 constexpr int mm = contract_over_rows ? n_rows : n_columns;
1731 constexpr int nn = contract_over_rows ? n_columns : n_rows;
1732 const Number2 *shape_data =
1733 symmetric_evaluate ?
1734 data.shape_values_eo.data() :
1735 data.values_within_subface[subface_index_1d].data();
1736 Assert(shape_data != nullptr, ExcNotInitialized());
1737 Assert(contract_over_rows == false || !add_into_result,
1738 ExcMessage("Cannot add into result if contract_over_rows = true"));
1739
1740 constexpr int n_blocks1 = Utilities::pow(fe_degree, direction);
1741 constexpr int n_blocks2 = Utilities::pow(fe_degree, dim - direction - 1);
1742 constexpr int stride_in = contract_over_rows ? 1 : stride;
1743 constexpr int stride_out = contract_over_rows ? stride : 1;
1744 constexpr EvaluatorVariant variant =
1745 symmetric_evaluate ? evaluate_evenodd : evaluate_general;
1746
1747 for (int i2 = 0; i2 < n_blocks2; ++i2)
1748 {
1749 for (int i1 = 0; i1 < n_blocks1; ++i1)
1750 {
1751 if (contract_over_rows == false && add_into_result)
1754 n_rows,
1755 n_columns,
1756 n_blocks1 * stride_in,
1757 n_blocks1 * stride_out,
1758 contract_over_rows,
1759 true>(shape_data, in, out);
1760 else
1763 n_rows,
1764 n_columns,
1765 n_blocks1 * stride_in,
1766 n_blocks1 * stride_out,
1767 contract_over_rows,
1768 false>(shape_data, in, out);
1769
1770 in += stride_in;
1771 out += stride_out;
1772 }
1773 in += n_blocks1 * (mm - 1) * stride_in;
1774 out += n_blocks1 * (nn - 1) * stride_out;
1775 }
1776 }
1777
1778 template <int direction,
1779 int normal_direction,
1780 int stride = 1,
1781 typename Number = double,
1782 typename Number2 = double>
1783 static void
1785 const Number *in,
1786 Number *out,
1787 const int subface_index_1d = 0)
1788 {
1789 AssertIndexRange(direction, dim);
1790 AssertDimension(fe_degree - 1, data.fe_degree);
1791 AssertDimension(n_q_points_1d, data.n_q_points_1d);
1792 static_assert(direction != normal_direction,
1793 "Cannot interpolate tangentially in normal direction");
1794
1795 constexpr int n_rows = std::max(fe_degree, 0);
1796 constexpr int n_columns = n_q_points_1d;
1797 const Number2 *shape_data =
1798 symmetric_evaluate ?
1799 data.shape_values_eo.data() :
1800 data.values_within_subface[subface_index_1d].data();
1801 Assert(shape_data != nullptr, ExcNotInitialized());
1802
1803 constexpr int n_blocks1 =
1804 (direction > normal_direction) ?
1805 Utilities::pow(n_q_points_1d, direction) :
1806 (direction > 0 ?
1807 (Utilities::pow(fe_degree, direction - 1) * n_q_points_1d) :
1808 1);
1809 constexpr int n_blocks2 =
1810 (direction > normal_direction) ?
1811 Utilities::pow(fe_degree, dim - 1 - direction) :
1812 ((direction + 1 < dim) ?
1813 (Utilities::pow(fe_degree, dim - 2 - direction) * n_q_points_1d) :
1814 1);
1815
1816 constexpr EvaluatorVariant variant =
1817 symmetric_evaluate ? evaluate_evenodd : evaluate_general;
1818
1819 // Since we may perform an in-place interpolation, we must run the step
1820 // expanding the size of the basis backward ('contract_over_rows' aka
1821 // 'evaluate' case), so shift the pointers and decrement during the loop
1822 if (contract_over_rows)
1823 {
1824 in += (n_blocks2 - 1) * n_blocks1 * n_rows + n_blocks1 - 1;
1825 out +=
1826 stride * ((n_blocks2 - 1) * n_blocks1 * n_columns + n_blocks1 - 1);
1827 for (int i2 = 0; i2 < n_blocks2; ++i2)
1828 {
1829 for (int i1 = 0; i1 < n_blocks1; ++i1)
1830 {
1833 n_rows,
1834 n_columns,
1835 n_blocks1,
1836 n_blocks1 * stride,
1837 true,
1838 false>(shape_data, in, out);
1839
1840 --in;
1841 out -= stride;
1842 }
1843 in -= n_blocks1 * (n_rows - 1);
1844 out -= n_blocks1 * (n_columns - 1) * stride;
1845 }
1846 }
1847 else
1848 {
1849 for (int i2 = 0; i2 < n_blocks2; ++i2)
1850 {
1851 for (int i1 = 0; i1 < n_blocks1; ++i1)
1852 {
1855 n_rows,
1856 n_columns,
1857 n_blocks1 * stride,
1858 n_blocks1,
1859 false,
1860 false>(shape_data, in, out);
1861
1862 in += stride;
1863 ++out;
1864 }
1865 in += n_blocks1 * (n_columns - 1) * stride;
1866 out += n_blocks1 * (n_rows - 1);
1867 }
1868 }
1869 }
1870 };
1871
1872
1873
1919 template <int n_rows_template,
1920 int stride_template,
1921 bool contract_onto_face,
1922 bool add,
1923 int max_derivative,
1924 typename Number,
1925 typename Number2>
1926 inline std::enable_if_t<contract_onto_face, void>
1927 interpolate_to_face(const Number2 *shape_values,
1928 const std::array<int, 2> &n_blocks,
1929 const std::array<int, 2> &steps,
1930 const Number *input,
1931 Number *DEAL_II_RESTRICT output,
1932 const int n_rows_runtime = 0,
1933 const int stride_runtime = 1)
1934 {
1935 const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
1936 const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
1937
1938 Number *output1 = output + n_blocks[0] * n_blocks[1];
1939 Number *output2 = output1 + n_blocks[0] * n_blocks[1];
1940 for (int i2 = 0; i2 < n_blocks[1]; ++i2)
1941 {
1942 for (int i1 = 0; i1 < n_blocks[0]; ++i1)
1943 {
1944 Number res0 = shape_values[0] * input[0];
1945 Number res1, res2;
1946 if (max_derivative > 0)
1947 res1 = shape_values[n_rows] * input[0];
1948 if (max_derivative > 1)
1949 res2 = shape_values[2 * n_rows] * input[0];
1950 for (int ind = 1; ind < n_rows; ++ind)
1951 {
1952 res0 += shape_values[ind] * input[stride * ind];
1953 if (max_derivative > 0)
1954 res1 += shape_values[ind + n_rows] * input[stride * ind];
1955 if (max_derivative > 1)
1956 res2 += shape_values[ind + 2 * n_rows] * input[stride * ind];
1957 }
1958 if (add)
1959 {
1960 output[i1] += res0;
1961 if (max_derivative > 0)
1962 output1[i1] += res1;
1963 if (max_derivative > 1)
1964 output2[i2] += res2;
1965 }
1966 else
1967 {
1968 output[i1] = res0;
1969 if (max_derivative > 0)
1970 output1[i1] = res1;
1971 if (max_derivative > 1)
1972 output2[i1] = res2;
1973 }
1974 input += steps[0];
1975 }
1976 output += n_blocks[0];
1977 if (max_derivative > 0)
1978 output1 += n_blocks[0];
1979 if (max_derivative > 1)
1980 output2 += n_blocks[0];
1981 input += steps[1];
1982 }
1983 }
1984
1985
1986
1994 constexpr bool
1995 use_collocation_evaluation(const unsigned int fe_degree,
1996 const unsigned int n_q_points_1d)
1997 {
1998 return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
1999 (n_q_points_1d <= 3 * fe_degree / 2 + 1);
2000 }
2001
2002
2003
2009 template <int n_rows_template,
2010 int stride_template,
2011 bool contract_onto_face,
2012 bool add,
2013 int max_derivative,
2014 typename Number,
2015 typename Number2>
2016 inline std::enable_if_t<!contract_onto_face, void>
2017 interpolate_to_face(const Number2 *shape_values,
2018 const std::array<int, 2> &n_blocks,
2019 const std::array<int, 2> &steps,
2020 const Number *input,
2021 Number *DEAL_II_RESTRICT output,
2022 const int n_rows_runtime = 0,
2023 const int stride_runtime = 1)
2024 {
2025 const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
2026 const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
2027
2028 const Number *input1 = input + n_blocks[0] * n_blocks[1];
2029 const Number *input2 = input1 + n_blocks[0] * n_blocks[1];
2030 for (int i2 = 0; i2 < n_blocks[1]; ++i2)
2031 {
2032 for (int i1 = 0; i1 < n_blocks[0]; ++i1)
2033 {
2034 const Number in = input[i1];
2035 Number in1, in2;
2036 if (max_derivative > 0)
2037 in1 = input1[i1];
2038 if (max_derivative > 1)
2039 in2 = input2[i1];
2040 for (int col = 0; col < n_rows; ++col)
2041 {
2042 Number result =
2043 add ? (output[col * stride] + shape_values[col] * in) :
2044 (shape_values[col] * in);
2045 if (max_derivative > 0)
2046 result += shape_values[col + n_rows] * in1;
2047 if (max_derivative > 1)
2048 result += shape_values[col + 2 * n_rows] * in2;
2049
2050 output[col * stride] = result;
2051 }
2052 output += steps[0];
2053 }
2054 input += n_blocks[0];
2055 if (max_derivative > 0)
2056 input1 += n_blocks[0];
2057 if (max_derivative > 1)
2058 input2 += n_blocks[0];
2059 output += steps[1];
2060 }
2061 }
2062
2063 template <int dim, int n_points_1d_template, typename Number>
2064 inline void
2065 weight_fe_q_dofs_by_entity(const Number *weights,
2066 const unsigned int n_components,
2067 const int n_points_1d_non_template,
2068 Number *data)
2069 {
2070 const int n_points_1d = n_points_1d_template != -1 ?
2071 n_points_1d_template :
2072 n_points_1d_non_template;
2073
2074 Assert(n_points_1d > 0, ExcNotImplemented());
2075 Assert(n_points_1d < 100, ExcNotImplemented());
2076
2077 unsigned int compressed_index[100];
2078 compressed_index[0] = 0;
2079 for (int i = 1; i < n_points_1d - 1; ++i)
2080 compressed_index[i] = 1;
2081 compressed_index[n_points_1d - 1] = 2;
2082
2083 for (unsigned int c = 0; c < n_components; ++c)
2084 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2085 for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
2086 {
2087 const unsigned int shift =
2088 9 * compressed_index[k] + 3 * compressed_index[j];
2089 data[0] *= weights[shift];
2090 // loop bound as int avoids compiler warnings in case n_points_1d
2091 // == 1 (polynomial degree 0)
2092 const Number weight = weights[shift + 1];
2093 for (int i = 1; i < n_points_1d - 1; ++i)
2094 data[i] *= weight;
2095 data[n_points_1d - 1] *= weights[shift + 2];
2096 data += n_points_1d;
2097 }
2098 }
2099
2100
2101 template <int dim, int n_points_1d_template, typename Number>
2102 inline void
2104 const unsigned int n_components,
2105 const int n_points_1d_non_template,
2106 Number *data)
2107 {
2108 const int n_points_1d = n_points_1d_template != -1 ?
2109 n_points_1d_template :
2110 n_points_1d_non_template;
2111
2112 Assert((n_points_1d % 2) == 1,
2113 ExcMessage("The function can only with add number of points"));
2114 Assert(n_points_1d > 0, ExcNotImplemented());
2115 Assert(n_points_1d < 100, ExcNotImplemented());
2116
2117 const unsigned int n_inside_1d = n_points_1d / 2;
2118
2119 unsigned int compressed_index[100];
2120
2121 unsigned int c = 0;
2122 for (int i = 0; i < n_inside_1d; ++i)
2123 compressed_index[c++] = 0;
2124 compressed_index[c++] = 1;
2125 for (int i = 0; i < n_inside_1d; ++i)
2126 compressed_index[c++] = 2;
2127
2128 for (unsigned int c = 0; c < n_components; ++c)
2129 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2130 for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
2131 {
2132 const unsigned int shift =
2133 9 * compressed_index[k] + 3 * compressed_index[j];
2134
2135 unsigned int c = 0;
2136 const Number weight1 = weights[shift];
2137 for (int i = 0; i < n_inside_1d; ++i)
2138 data[c++] *= weight1;
2139 data[c++] *= weights[shift + 1];
2140 const Number weight2 = weights[shift + 2];
2141 for (int i = 0; i < n_inside_1d; ++i)
2142 data[c++] *= weight2;
2143 data += n_points_1d;
2144 }
2145 }
2146
2147
2148 template <int dim, int n_points_1d_template, typename Number>
2149 inline bool
2151 const unsigned int n_components,
2152 const int n_points_1d_non_template,
2153 Number *weights)
2154 {
2155 const int n_points_1d = n_points_1d_template != -1 ?
2156 n_points_1d_template :
2157 n_points_1d_non_template;
2158
2159 Assert(n_points_1d > 0, ExcNotImplemented());
2160 Assert(n_points_1d < 100, ExcNotImplemented());
2161
2162 unsigned int compressed_index[100];
2163 compressed_index[0] = 0;
2164 for (int i = 1; i < n_points_1d - 1; ++i)
2165 compressed_index[i] = 1;
2166 compressed_index[n_points_1d - 1] = 2;
2167
2168 // Insert the number data into a storage position for weight,
2169 // ensuring that the array has either not been touched before
2170 // or the previous content is the same. In case the previous
2171 // content has a different value, we exit this function and
2172 // signal to outer functions that the compression was not possible.
2173 const auto check_and_set = [](Number &weight, const Number &data) {
2174 if (weight == Number(-1.0) || weight == data)
2175 {
2176 weight = data;
2177 return true; // success for the entry
2178 }
2179
2180 return false; // failure for the entry
2181 };
2182
2183 for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
2184 weights[c] = Number(-1.0);
2185
2186 for (unsigned int c = 0; c < n_components; ++c)
2187 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2188 for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
2189 ++j, data += n_points_1d)
2190 {
2191 const unsigned int shift =
2192 9 * compressed_index[k] + 3 * compressed_index[j];
2193
2194 if (!check_and_set(weights[shift], data[0]))
2195 return false; // failure
2196
2197 for (int i = 1; i < n_points_1d - 1; ++i)
2198 if (!check_and_set(weights[shift + 1], data[i]))
2199 return false; // failure
2200
2201 if (!check_and_set(weights[shift + 2], data[n_points_1d - 1]))
2202 return false; // failure
2203 }
2204
2205 return true; // success
2206 }
2207
2208
2209 template <int dim, int n_points_1d_template, typename Number>
2210 inline bool
2212 const Number *data,
2213 const unsigned int n_components,
2214 const int n_points_1d_non_template,
2215 Number *weights)
2216 {
2217 const int n_points_1d = n_points_1d_template != -1 ?
2218 n_points_1d_template :
2219 n_points_1d_non_template;
2220
2221 Assert((n_points_1d % 2) == 1,
2222 ExcMessage("The function can only with add number of points"));
2223 Assert(n_points_1d > 0, ExcNotImplemented());
2224 Assert(n_points_1d < 100, ExcNotImplemented());
2225
2226 const unsigned int n_inside_1d = n_points_1d / 2;
2227
2228 unsigned int compressed_index[100];
2229
2230 unsigned int c = 0;
2231 for (int i = 0; i < n_inside_1d; ++i)
2232 compressed_index[c++] = 0;
2233 compressed_index[c++] = 1;
2234 for (int i = 0; i < n_inside_1d; ++i)
2235 compressed_index[c++] = 2;
2236
2237 // Insert the number data into a storage position for weight,
2238 // ensuring that the array has either not been touched before
2239 // or the previous content is the same. In case the previous
2240 // content has a different value, we exit this function and
2241 // signal to outer functions that the compression was not possible.
2242 const auto check_and_set = [](Number &weight, const Number &data) {
2243 if (weight == Number(-1.0) || weight == data)
2244 {
2245 weight = data;
2246 return true; // success for the entry
2247 }
2248
2249 return false; // failure for the entry
2250 };
2251
2252 for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
2253 weights[c] = Number(-1.0);
2254
2255 for (unsigned int comp = 0; comp < n_components; ++comp)
2256 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2257 for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
2258 ++j, data += n_points_1d)
2259 {
2260 const unsigned int shift =
2261 9 * compressed_index[k] + 3 * compressed_index[j];
2262
2263 unsigned int c = 0;
2264
2265 for (int i = 0; i < n_inside_1d; ++i)
2266 if (!check_and_set(weights[shift], data[c++]))
2267 return false; // failure
2268
2269 if (!check_and_set(weights[shift + 1], data[c++]))
2270 return false; // failure
2271
2272 for (int i = 0; i < n_inside_1d; ++i)
2273 if (!check_and_set(weights[shift + 2], data[c++]))
2274 return false; // failure
2275 }
2276
2277 return true; // success
2278 }
2279
2280
2281} // end of namespace internal
2282
2283
2285
2286#endif
pointer data()
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_RESTRICT
Definition config.h:110
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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)
constexpr T fixed_power(const T t)
Definition utilities.h:939
constexpr T pow(const T base, const int iexp)
Definition utilities.h:963
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)
static const unsigned int invalid_unsigned_int
Definition types.h:220
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, 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)
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)
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)
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)
static constexpr unsigned int n_rows_of_product
void hessians(const Number in[], Number out[]) const
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
std::array< AlignedVector< Number >, 2 > values_within_subface
Definition shape_info.h:314