deal.II version GIT relicensing-6834-g5b78e6bcdf 2026-10-01 11:20:01+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: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2017 - 2025 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13
14#ifndef dealii_matrix_free_tensor_product_kernels_h
15#define dealii_matrix_free_tensor_product_kernels_h
16
17#include <deal.II/base/config.h>
18
19#include <deal.II/base/aligned_vector.h>
20#include <deal.II/base/utilities.h>
21
22#include <deal.II/matrix_free/shape_info.h>
23
24
25DEAL_II_NAMESPACE_OPEN
26
27
28
29namespace internal
30{
66
67
68
72 enum class EvaluatorQuantity
73 {
77 value,
81 gradient,
85 hessian
86 };
87
88
89
105 template <EvaluatorVariant variant,
106 EvaluatorQuantity quantity,
107 int n_rows,
108 int n_columns,
109 int stride_in,
110 int stride_out,
111 bool transpose_matrix,
112 bool add,
113 typename Number,
114 typename Number2>
115 std::enable_if_t<(variant == evaluate_general), void>
116 apply_matrix_vector_product(const Number2 *matrix,
117 const Number *in,
118 Number *out)
119 {
120 // We can only statically assert that one argument is non-zero because
121 // face evaluation might instantiate some functions, so we need to use the
122 // run-time assert to verify that we do not end up involuntarily.
123 static_assert(n_rows > 0 || n_columns > 0,
124 "Specialization only for n_rows, n_columns > 0");
125 Assert(n_rows > 0 && n_columns > 0,
126 ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
127 std::to_string(n_rows) + ", " +
128 std::to_string(n_columns) + " was passed!"));
129 static_assert(quantity == EvaluatorQuantity::value,
130 "This function should only use EvaluatorQuantity::value");
131
132 constexpr int mm = transpose_matrix ? n_rows : n_columns,
133 nn = transpose_matrix ? n_columns : n_rows;
134
135 std::array<Number, mm> x;
136 for (int i = 0; i < mm; ++i)
137 x[i] = in[stride_in * i];
138 for (int col = 0; col < nn; ++col)
139 {
140 Number res0;
141 if (transpose_matrix == true)
142 {
143 res0 = matrix[col] * x[0];
144 for (int i = 1; i < mm; ++i)
145 {
146 const Number2 mji = matrix[i * n_columns + col];
147 if constexpr (numbers::NumberTraits<Number>::is_complex &&
148 numbers::NumberTraits<Number2>::is_complex)
149 {
150 res0.real(res0.real() + mji.real() * x[i].real() -
151 mji.imag() * x[i].imag());
152 res0.imag(res0.imag() + mji.imag() * x[i].real() +
153 mji.real() * x[i].imag());
154 }
155 else
156 res0 += mji * x[i];
157 }
158 }
159 else
160 {
161 res0 = matrix[col * n_columns] * x[0];
162 for (int i = 1; i < mm; ++i)
163 {
164 const Number2 mij = matrix[col * n_columns + i];
165 if constexpr (numbers::NumberTraits<Number>::is_complex &&
166 numbers::NumberTraits<Number2>::is_complex)
167 {
168 res0.real(res0.real() + mij.real() * x[i].real() -
169 mij.imag() * x[i].imag());
170 res0.imag(res0.imag() + mij.imag() * x[i].real() +
171 mij.real() * x[i].imag());
172 }
173 else
174 res0 += mij * x[i];
175 }
176 }
177 if (add)
178 out[stride_out * col] += res0;
179 else
180 out[stride_out * col] = res0;
181 }
182 }
183
184
185
190 template <EvaluatorVariant variant,
191 EvaluatorQuantity quantity,
192 bool transpose_matrix,
193 bool add,
194 bool consider_strides,
195 typename Number,
196 typename Number2,
197 int n_components = 1>
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_given,
205 const int stride_out_given)
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 Assert(consider_strides || (stride_in_given == 1 && stride_out_given == 1),
220 ExcInternalError());
221 const int stride_in = consider_strides ? stride_in_given : 1;
222 const int stride_out = consider_strides ? stride_out_given : 1;
223
224 static_assert(n_components > 0 && n_components < 4,
225 "Invalid number of components");
226
227 // specialization for n_rows = 2 that manually unrolls the innermost loop
228 // to make the operation perform better (not completely as good as the
229 // templated one, but much better than the generic version down below,
230 // because the loop over col can be more effectively unrolled by the
231 // compiler)
232 if (transpose_matrix && n_rows == 2 && n_components == 1)
233 {
234 const Number2 *matrix_1 = matrix + n_columns;
235 const Number x0 = in[0], x1 = in[stride_in];
236 for (int col = 0; col < nn; ++col)
237 {
238 const Number result = matrix[col] * x0 + matrix_1[col] * x1;
239 if (add)
240 out[stride_out * col] += result;
241 else
242 out[stride_out * col] = result;
243 }
244 }
245 else if (transpose_matrix && n_rows == 3 && n_components == 1)
246 {
247 const Number2 *matrix_1 = matrix + n_columns;
248 const Number2 *matrix_2 = matrix_1 + n_columns;
249 const Number x0 = in[0], x1 = in[stride_in], x2 = in[2 * stride_in];
250 for (int col = 0; col < nn; ++col)
251 {
252 const Number result =
253 matrix[col] * x0 + matrix_1[col] * x1 + matrix_2[col] * x2;
254 if (add)
255 out[stride_out * col] += result;
256 else
257 out[stride_out * col] = result;
258 }
259 }
260 else if (std::abs(in - out) < std::min(stride_out * nn, stride_in * mm) &&
261 n_components == 1)
262 {
263 Assert(mm <= 128,
264 ExcNotImplemented("For large sizes, arrays may not overlap"));
265 std::array<Number, 129> x;
266 for (int i = 0; i < mm; ++i)
267 x[i] = in[stride_in * i];
268
269 for (int col = 0; col < nn; ++col)
270 {
271 Number res0;
272 if (transpose_matrix == true)
273 {
274 res0 = matrix[col] * x[0];
275 for (int i = 1; i < mm; ++i)
276 res0 += matrix[i * n_columns + col] * x[i];
277 }
278 else
279 {
280 res0 = matrix[col * n_columns] * x[0];
281 for (int i = 1; i < mm; ++i)
282 res0 += matrix[col * n_columns + i] * x[i];
283 }
284 if (add)
285 out[stride_out * col] += res0;
286 else
287 out[stride_out * col] = res0;
288 }
289 }
290 else
291 {
292 const Number *in0 = in;
293 const Number *in1 = n_components > 1 ? in + mm : nullptr;
294 const Number *in2 = n_components > 2 ? in + 2 * mm : nullptr;
295
296 Number *out0 = out;
297 Number *out1 = n_components > 1 ? out + nn : nullptr;
298 Number *out2 = n_components > 2 ? out + 2 * nn : nullptr;
299
300 int nn_regular = (nn / 4) * 4;
301 for (int col = 0; col < nn_regular; col += 4)
302 {
303 Number res[12];
304 if (transpose_matrix == true)
305 {
306 const Number2 *matrix_ptr = matrix + col;
307 const Number a = in0[0];
308 res[0] = matrix_ptr[0] * a;
309 res[1] = matrix_ptr[1] * a;
310 res[2] = matrix_ptr[2] * a;
311 res[3] = matrix_ptr[3] * a;
312
313 if (n_components > 1)
314 {
315 const Number b = in1[0];
316 res[4] = matrix_ptr[0] * b;
317 res[5] = matrix_ptr[1] * b;
318 res[6] = matrix_ptr[2] * b;
319 res[7] = matrix_ptr[3] * b;
320 }
321
322 if (n_components > 2)
323 {
324 const Number c = in2[0];
325 res[8] = matrix_ptr[0] * c;
326 res[9] = matrix_ptr[1] * c;
327 res[10] = matrix_ptr[2] * c;
328 res[11] = matrix_ptr[3] * c;
329 }
330
331 matrix_ptr += n_columns;
332 for (int i = 1; i < mm; ++i, matrix_ptr += n_columns)
333 {
334 const Number a = in0[stride_in * i];
335 res[0] += matrix_ptr[0] * a;
336 res[1] += matrix_ptr[1] * a;
337 res[2] += matrix_ptr[2] * a;
338 res[3] += matrix_ptr[3] * a;
339
340 if (n_components > 1)
341 {
342 const Number b = in1[stride_in * i];
343 res[4] += matrix_ptr[0] * b;
344 res[5] += matrix_ptr[1] * b;
345 res[6] += matrix_ptr[2] * b;
346 res[7] += matrix_ptr[3] * b;
347 }
348 if (n_components > 2)
349 {
350 const Number c = in2[stride_in * i];
351 res[8] += matrix_ptr[0] * c;
352 res[9] += matrix_ptr[1] * c;
353 res[10] += matrix_ptr[2] * c;
354 res[11] += matrix_ptr[3] * c;
355 }
356 }
357 }
358 else
359 {
360 const Number2 *matrix_0 = matrix + col * n_columns;
361 const Number2 *matrix_1 = matrix + (col + 1) * n_columns;
362 const Number2 *matrix_2 = matrix + (col + 2) * n_columns;
363 const Number2 *matrix_3 = matrix + (col + 3) * n_columns;
364
365 const Number a = in0[0];
366 res[0] = matrix_0[0] * a;
367 res[1] = matrix_1[0] * a;
368 res[2] = matrix_2[0] * a;
369 res[3] = matrix_3[0] * a;
370
371 if (n_components > 1)
372 {
373 const Number b = in1[0];
374 res[4] = matrix_0[0] * b;
375 res[5] = matrix_1[0] * b;
376 res[6] = matrix_2[0] * b;
377 res[7] = matrix_3[0] * b;
378 }
379
380 if (n_components > 2)
381 {
382 const Number c = in2[0];
383 res[8] = matrix_0[0] * c;
384 res[9] = matrix_1[0] * c;
385 res[10] = matrix_2[0] * c;
386 res[11] = matrix_3[0] * c;
387 }
388
389 for (int i = 1; i < mm; ++i)
390 {
391 const Number a = in0[stride_in * i];
392 res[0] += matrix_0[i] * a;
393 res[1] += matrix_1[i] * a;
394 res[2] += matrix_2[i] * a;
395 res[3] += matrix_3[i] * a;
396
397 if (n_components > 1)
398 {
399 const Number b = in1[stride_in * i];
400 res[4] += matrix_0[i] * b;
401 res[5] += matrix_1[i] * b;
402 res[6] += matrix_2[i] * b;
403 res[7] += matrix_3[i] * b;
404 }
405
406 if (n_components > 2)
407 {
408 const Number c = in2[stride_in * i];
409 res[8] += matrix_0[i] * c;
410 res[9] += matrix_1[i] * c;
411 res[10] += matrix_2[i] * c;
412 res[11] += matrix_3[i] * c;
413 }
414 }
415 }
416 if (add)
417 {
418 out0[0] += res[0];
419 out0[stride_out] += res[1];
420 out0[2 * stride_out] += res[2];
421 out0[3 * stride_out] += res[3];
422 if (n_components > 1)
423 {
424 out1[0] += res[4];
425 out1[stride_out] += res[5];
426 out1[2 * stride_out] += res[6];
427 out1[3 * stride_out] += res[7];
428 }
429 if (n_components > 2)
430 {
431 out2[0] += res[8];
432 out2[stride_out] += res[9];
433 out2[2 * stride_out] += res[10];
434 out2[3 * stride_out] += res[11];
435 }
436 }
437 else
438 {
439 out0[0] = res[0];
440 out0[stride_out] = res[1];
441 out0[2 * stride_out] = res[2];
442 out0[3 * stride_out] = res[3];
443 if (n_components > 1)
444 {
445 out1[0] = res[4];
446 out1[stride_out] = res[5];
447 out1[2 * stride_out] = res[6];
448 out1[3 * stride_out] = res[7];
449 }
450 if (n_components > 2)
451 {
452 out2[0] = res[8];
453 out2[stride_out] = res[9];
454 out2[2 * stride_out] = res[10];
455 out2[3 * stride_out] = res[11];
456 }
457 }
458 out0 += 4 * stride_out;
459 if (n_components > 1)
460 out1 += 4 * stride_out;
461 if (n_components > 2)
462 out2 += 4 * stride_out;
463 }
464 if (nn - nn_regular == 3)
465 {
466 Number res0, res1, res2, res3, res4, res5, res6, res7, res8;
467 if (transpose_matrix == true)
468 {
469 const Number2 *matrix_ptr = matrix + nn_regular;
470 res0 = matrix_ptr[0] * in0[0];
471 res1 = matrix_ptr[1] * in0[0];
472 res2 = matrix_ptr[2] * in0[0];
473 if (n_components > 1)
474 {
475 res3 = matrix_ptr[0] * in1[0];
476 res4 = matrix_ptr[1] * in1[0];
477 res5 = matrix_ptr[2] * in1[0];
478 }
479 if (n_components > 2)
480 {
481 res6 = matrix_ptr[0] * in2[0];
482 res7 = matrix_ptr[1] * in2[0];
483 res8 = matrix_ptr[2] * in2[0];
484 }
485 matrix_ptr += n_columns;
486 for (int i = 1; i < mm; ++i, matrix_ptr += n_columns)
487 {
488 res0 += matrix_ptr[0] * in0[stride_in * i];
489 res1 += matrix_ptr[1] * in0[stride_in * i];
490 res2 += matrix_ptr[2] * in0[stride_in * i];
491 if (n_components > 1)
492 {
493 res3 += matrix_ptr[0] * in1[stride_in * i];
494 res4 += matrix_ptr[1] * in1[stride_in * i];
495 res5 += matrix_ptr[2] * in1[stride_in * i];
496 }
497 if (n_components > 2)
498 {
499 res6 += matrix_ptr[0] * in2[stride_in * i];
500 res7 += matrix_ptr[1] * in2[stride_in * i];
501 res8 += matrix_ptr[2] * in2[stride_in * i];
502 }
503 }
504 }
505 else
506 {
507 const Number2 *matrix_0 = matrix + nn_regular * n_columns;
508 const Number2 *matrix_1 = matrix + (nn_regular + 1) * n_columns;
509 const Number2 *matrix_2 = matrix + (nn_regular + 2) * n_columns;
510
511 res0 = matrix_0[0] * in0[0];
512 res1 = matrix_1[0] * in0[0];
513 res2 = matrix_2[0] * in0[0];
514 if (n_components > 1)
515 {
516 res3 = matrix_0[0] * in1[0];
517 res4 = matrix_1[0] * in1[0];
518 res5 = matrix_2[0] * in1[0];
519 }
520 if (n_components > 2)
521 {
522 res6 = matrix_0[0] * in2[0];
523 res7 = matrix_1[0] * in2[0];
524 res8 = matrix_2[0] * in2[0];
525 }
526 for (int i = 1; i < mm; ++i)
527 {
528 res0 += matrix_0[i] * in0[stride_in * i];
529 res1 += matrix_1[i] * in0[stride_in * i];
530 res2 += matrix_2[i] * in0[stride_in * i];
531 if (n_components > 1)
532 {
533 res3 += matrix_0[i] * in1[stride_in * i];
534 res4 += matrix_1[i] * in1[stride_in * i];
535 res5 += matrix_2[i] * in1[stride_in * i];
536 }
537 if (n_components > 2)
538 {
539 res6 += matrix_0[i] * in2[stride_in * i];
540 res7 += matrix_1[i] * in2[stride_in * i];
541 res8 += matrix_2[i] * in2[stride_in * i];
542 }
543 }
544 }
545 if (add)
546 {
547 out0[0] += res0;
548 out0[stride_out] += res1;
549 out0[2 * stride_out] += res2;
550 if (n_components > 1)
551 {
552 out1[0] += res3;
553 out1[stride_out] += res4;
554 out1[2 * stride_out] += res5;
555 }
556 if (n_components > 2)
557 {
558 out2[0] += res6;
559 out2[stride_out] += res7;
560 out2[2 * stride_out] += res8;
561 }
562 }
563 else
564 {
565 out0[0] = res0;
566 out0[stride_out] = res1;
567 out0[2 * stride_out] = res2;
568 if (n_components > 1)
569 {
570 out1[0] = res3;
571 out1[stride_out] = res4;
572 out1[2 * stride_out] = res5;
573 }
574 if (n_components > 2)
575 {
576 out2[0] = res6;
577 out2[stride_out] = res7;
578 out2[2 * stride_out] = res8;
579 }
580 }
581 }
582 else if (nn - nn_regular == 2)
583 {
584 Number res0, res1, res2, res3, res4, res5;
585 if (transpose_matrix == true)
586 {
587 const Number2 *matrix_ptr = matrix + nn_regular;
588 res0 = matrix_ptr[0] * in0[0];
589 res1 = matrix_ptr[1] * in0[0];
590 if (n_components > 1)
591 {
592 res2 = matrix_ptr[0] * in1[0];
593 res3 = matrix_ptr[1] * in1[0];
594 }
595 if (n_components > 2)
596 {
597 res4 = matrix_ptr[0] * in2[0];
598 res5 = matrix_ptr[1] * in2[0];
599 }
600 matrix_ptr += n_columns;
601 for (int i = 1; i < mm; ++i, matrix_ptr += n_columns)
602 {
603 res0 += matrix_ptr[0] * in0[stride_in * i];
604 res1 += matrix_ptr[1] * in0[stride_in * i];
605 if (n_components > 1)
606 {
607 res2 += matrix_ptr[0] * in1[stride_in * i];
608 res3 += matrix_ptr[1] * in1[stride_in * i];
609 }
610 if (n_components > 2)
611 {
612 res4 += matrix_ptr[0] * in2[stride_in * i];
613 res5 += matrix_ptr[1] * in2[stride_in * i];
614 }
615 }
616 }
617 else
618 {
619 const Number2 *matrix_0 = matrix + nn_regular * n_columns;
620 const Number2 *matrix_1 = matrix + (nn_regular + 1) * n_columns;
621
622 res0 = matrix_0[0] * in0[0];
623 res1 = matrix_1[0] * in0[0];
624 if (n_components > 1)
625 {
626 res2 = matrix_0[0] * in1[0];
627 res3 = matrix_1[0] * in1[0];
628 }
629 if (n_components > 2)
630 {
631 res4 = matrix_0[0] * in2[0];
632 res5 = matrix_1[0] * in2[0];
633 }
634 for (int i = 1; i < mm; ++i)
635 {
636 res0 += matrix_0[i] * in0[stride_in * i];
637 res1 += matrix_1[i] * in0[stride_in * i];
638 if (n_components > 1)
639 {
640 res2 += matrix_0[i] * in1[stride_in * i];
641 res3 += matrix_1[i] * in1[stride_in * i];
642 }
643 if (n_components > 2)
644 {
645 res4 += matrix_0[i] * in2[stride_in * i];
646 res5 += matrix_1[i] * in2[stride_in * i];
647 }
648 }
649 }
650 if (add)
651 {
652 out0[0] += res0;
653 out0[stride_out] += res1;
654 if (n_components > 1)
655 {
656 out1[0] += res2;
657 out1[stride_out] += res3;
658 }
659 if (n_components > 2)
660 {
661 out2[0] += res4;
662 out2[stride_out] += res5;
663 }
664 }
665 else
666 {
667 out0[0] = res0;
668 out0[stride_out] = res1;
669 if (n_components > 1)
670 {
671 out1[0] = res2;
672 out1[stride_out] = res3;
673 }
674 if (n_components > 2)
675 {
676 out2[0] = res4;
677 out2[stride_out] = res5;
678 }
679 }
680 }
681 else if (nn - nn_regular == 1)
682 {
683 Number res0, res1, res2;
684 if (transpose_matrix == true)
685 {
686 const Number2 *matrix_ptr = matrix + nn_regular;
687 res0 = matrix_ptr[0] * in0[0];
688 if (n_components > 1)
689 res1 = matrix_ptr[0] * in1[0];
690 if (n_components > 2)
691 res2 = matrix_ptr[0] * in2[0];
692 matrix_ptr += n_columns;
693 for (int i = 1; i < mm; ++i, matrix_ptr += n_columns)
694 {
695 res0 += matrix_ptr[0] * in0[stride_in * i];
696 if (n_components > 1)
697 res1 += matrix_ptr[0] * in1[stride_in * i];
698 if (n_components > 2)
699 res2 += matrix_ptr[0] * in2[stride_in * i];
700 }
701 }
702 else
703 {
704 const Number2 *matrix_ptr = matrix + nn_regular * n_columns;
705 res0 = matrix_ptr[0] * in0[0];
706 if (n_components > 1)
707 res1 = matrix_ptr[0] * in1[0];
708 if (n_components > 2)
709 res2 = matrix_ptr[0] * in2[0];
710 for (int i = 1; i < mm; ++i)
711 {
712 res0 += matrix_ptr[i] * in0[stride_in * i];
713 if (n_components > 1)
714 res1 += matrix_ptr[i] * in1[stride_in * i];
715 if (n_components > 2)
716 res2 += matrix_ptr[i] * in2[stride_in * i];
717 }
718 }
719 if (add)
720 {
721 out0[0] += res0;
722 if (n_components > 1)
723 out1[0] += res1;
724 if (n_components > 2)
725 out2[0] += res2;
726 }
727 else
728 {
729 out0[0] = res0;
730 if (n_components > 1)
731 out1[0] = res1;
732 if (n_components > 2)
733 out2[0] = res2;
734 }
735 }
736 }
737 }
738
739
740
747 template <EvaluatorVariant variant,
748 EvaluatorQuantity quantity,
749 int n_rows,
750 int n_columns,
751 int stride_in,
752 int stride_out,
753 bool transpose_matrix,
754 bool add,
755 typename Number,
756 typename Number2>
757 std::enable_if_t<(variant == evaluate_symmetric), void>
758 apply_matrix_vector_product(const Number2 *matrix,
759 const Number *in,
760 Number *out)
761 {
762 // We can only statically assert that one argument is non-zero because
763 // face evaluation might instantiate some functions, so we need to use the
764 // run-time assert to verify that we do not end up involuntarily.
765 static_assert(n_rows > 0 || n_columns > 0,
766 "Specialization only for n_rows, n_columns > 0");
767 Assert(n_rows > 0 && n_columns > 0,
768 ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
769 std::to_string(n_rows) + ", " +
770 std::to_string(n_columns) + " was passed!"));
771
772 constexpr int mm = transpose_matrix ? n_rows : n_columns,
773 nn = transpose_matrix ? n_columns : n_rows;
774 constexpr int n_cols = nn / 2;
775 constexpr int mid = mm / 2;
776
777 std::array<Number, mm> x;
778 for (int i = 0; i < mm; ++i)
779 x[i] = in[stride_in * i];
780
781 if (quantity == EvaluatorQuantity::value)
782 {
783 // In this case, the 1d shape values read (sorted lexicographically,
784 // rows run over 1d dofs, columns over quadrature points):
785 // Q2 --> [ 0.687 0 -0.087 ]
786 // [ 0.4 1 0.4 ]
787 // [-0.087 0 0.687 ]
788 // Q3 --> [ 0.66 0.003 0.002 0.049 ]
789 // [ 0.521 1.005 -0.01 -0.230 ]
790 // [-0.230 -0.01 1.005 0.521 ]
791 // [ 0.049 0.002 0.003 0.66 ]
792 // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ]
793 // [ 0.608 1.059 0 0.039 0.176 ]
794 // [-0.409 -0.113 1 -0.113 -0.409 ]
795 // [ 0.176 0.039 0 1.059 0.608 ]
796 // [-0.032 -0.007 0 0.022 0.658 ]
797 //
798 // In these matrices, we want to use avoid computations involving
799 // zeros and ones and use the symmetry in entries starting from (1,1)
800 // forward and (N,N) backward, respectively to reduce the number of
801 // read operations.
802 for (int col = 0; col < n_cols; ++col)
803 {
804 Number2 val0, val1;
805 Number res0, res1;
806 if (transpose_matrix == true)
807 {
808 val0 = matrix[col];
809 val1 = matrix[nn - 1 - col];
810 }
811 else
812 {
813 val0 = matrix[col * n_columns];
814 val1 = matrix[(col + 1) * n_columns - 1];
815 }
816 if (mid > 0)
817 {
818 res0 = val0 * x[0];
819 res1 = val1 * x[0];
820 res0 += val1 * x[mm - 1];
821 res1 += val0 * x[mm - 1];
822 for (int ind = 1; ind < mid; ++ind)
823 {
824 if (transpose_matrix == true)
825 {
826 val0 = matrix[ind * n_columns + col];
827 val1 = matrix[ind * n_columns + nn - 1 - col];
828 }
829 else
830 {
831 val0 = matrix[col * n_columns + ind];
832 val1 = matrix[(col + 1) * n_columns - 1 - ind];
833 }
834 res0 += val0 * x[ind];
835 res1 += val1 * x[ind];
836 res0 += val1 * x[mm - 1 - ind];
837 res1 += val0 * x[mm - 1 - ind];
838 }
839 }
840 else
841 res0 = res1 = Number();
842 if (transpose_matrix == true)
843 {
844 if (mm % 2 == 1)
845 {
846 const Number tmp = matrix[mid * n_columns + col] * x[mid];
847 res0 += tmp;
848 res1 += tmp;
849 }
850 }
851 else
852 {
853 if (mm % 2 == 1 && nn % 2 == 0)
854 {
855 const Number tmp = matrix[col * n_columns + mid] * x[mid];
856 res0 += tmp;
857 res1 += tmp;
858 }
859 }
860 if (add)
861 {
862 out[stride_out * col] += res0;
863 out[stride_out * (nn - 1 - col)] += res1;
864 }
865 else
866 {
867 out[stride_out * col] = res0;
868 out[stride_out * (nn - 1 - col)] = res1;
869 }
870 }
871 if (transpose_matrix == true && nn % 2 == 1 && mm % 2 == 1)
872 {
873 if (add)
874 out[stride_out * n_cols] += x[mid];
875 else
876 out[stride_out * n_cols] = x[mid];
877 }
878 else if (transpose_matrix == true && nn % 2 == 1)
879 {
880 Number res0;
881 if (mid > 0)
882 {
883 res0 = matrix[n_cols] * (x[0] + x[mm - 1]);
884 for (int ind = 1; ind < mid; ++ind)
885 {
886 const Number2 val0 = matrix[ind * n_columns + n_cols];
887 res0 += val0 * (x[ind] + in[mm - 1 - ind]);
888 }
889 }
890 else
891 res0 = Number();
892 if (add)
893 out[stride_out * n_cols] += res0;
894 else
895 out[stride_out * n_cols] = res0;
896 }
897 else if (transpose_matrix == false && nn % 2 == 1)
898 {
899 Number res0;
900 if (mid > 0)
901 {
902 res0 = matrix[n_cols * n_columns] * (x[0] + x[mm - 1]);
903 for (int ind = 1; ind < mid; ++ind)
904 {
905 const Number2 val0 = matrix[n_cols * n_columns + ind];
906 res0 += val0 * (x[ind] + x[mm - 1 - ind]);
907 }
908 if (mm % 2)
909 res0 += x[mid];
910 }
911 else
912 res0 = in[0];
913 if (add)
914 out[stride_out * n_cols] += res0;
915 else
916 out[stride_out * n_cols] = res0;
917 }
918 }
919 else if (quantity == EvaluatorQuantity::gradient)
920 {
921 // For the specialized loop used for gradient computations we again
922 // exploit symmetries according to the following entries (sorted
923 // lexicographically, rows run over 1d dofs, columns over quadrature
924 // points):
925 // Q2 --> [-2.549 -1 0.549 ]
926 // [ 3.098 0 -3.098 ]
927 // [-0.549 1 2.549 ]
928 // Q3 --> [-4.315 -1.03 0.5 -0.44 ]
929 // [ 6.07 -1.44 -2.97 2.196 ]
930 // [-2.196 2.97 1.44 -6.07 ]
931 // [ 0.44 -0.5 1.03 4.315 ]
932 // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ]
933 // [10.111 -2.76 -2.667 2.066 -2.306 ]
934 // [-5.688 5.773 0 -5.773 5.688 ]
935 // [ 2.306 -2.066 2.667 2.76 -10.111 ]
936 // [-0.413 0.353 -0.333 -0.353 0.413 ]
937 for (int col = 0; col < n_cols; ++col)
938 {
939 Number2 val0, val1;
940 Number res0, res1;
941 if (transpose_matrix == true)
942 {
943 val0 = matrix[col];
944 val1 = matrix[nn - 1 - col];
945 }
946 else
947 {
948 val0 = matrix[col * n_columns];
949 val1 = matrix[(nn - col - 1) * n_columns];
950 }
951 if (mid > 0)
952 {
953 res0 = val0 * x[0];
954 res1 = val1 * x[0];
955 res0 -= val1 * x[mm - 1];
956 res1 -= val0 * x[mm - 1];
957 for (int ind = 1; ind < mid; ++ind)
958 {
959 if (transpose_matrix == true)
960 {
961 val0 = matrix[ind * n_columns + col];
962 val1 = matrix[ind * n_columns + nn - 1 - col];
963 }
964 else
965 {
966 val0 = matrix[col * n_columns + ind];
967 val1 = matrix[(nn - col - 1) * n_columns + ind];
968 }
969 res0 += val0 * x[ind];
970 res1 += val1 * x[ind];
971 res0 -= val1 * x[mm - 1 - ind];
972 res1 -= val0 * x[mm - 1 - ind];
973 }
974 }
975 else
976 res0 = res1 = Number();
977 if (mm % 2 == 1)
978 {
979 if (transpose_matrix == true)
980 val0 = matrix[mid * n_columns + col];
981 else
982 val0 = matrix[col * n_columns + mid];
983 const Number tmp = val0 * x[mid];
984 res0 += tmp;
985 res1 -= tmp;
986 }
987 if (add)
988 {
989 out[stride_out * col] += res0;
990 out[stride_out * (nn - 1 - col)] += res1;
991 }
992 else
993 {
994 out[stride_out * col] = res0;
995 out[stride_out * (nn - 1 - col)] = res1;
996 }
997 }
998 if (nn % 2 == 1)
999 {
1000 Number2 val0;
1001 Number res0;
1002 if (transpose_matrix == true)
1003 val0 = matrix[n_cols];
1004 else
1005 val0 = matrix[n_cols * n_columns];
1006 res0 = val0 * (x[0] - x[mm - 1]);
1007 for (int ind = 1; ind < mid; ++ind)
1008 {
1009 if (transpose_matrix == true)
1010 val0 = matrix[ind * n_columns + n_cols];
1011 else
1012 val0 = matrix[n_cols * n_columns + ind];
1013 Number in1 = val0 * (x[ind] - x[mm - 1 - ind]);
1014 res0 += in1;
1015 }
1016 if (add)
1017 out[stride_out * n_cols] += res0;
1018 else
1019 out[stride_out * n_cols] = res0;
1020 }
1021 }
1022 else
1023 {
1024 // Hessians are almost the same as values, apart from some missing '1'
1025 // entries
1026 for (int col = 0; col < n_cols; ++col)
1027 {
1028 Number2 val0, val1;
1029 Number res0, res1;
1030 if (transpose_matrix == true)
1031 {
1032 val0 = matrix[col];
1033 val1 = matrix[nn - 1 - col];
1034 }
1035 else
1036 {
1037 val0 = matrix[col * n_columns];
1038 val1 = matrix[(col + 1) * n_columns - 1];
1039 }
1040 if (mid > 0)
1041 {
1042 res0 = val0 * x[0];
1043 res1 = val1 * x[0];
1044 res0 += val1 * x[mm - 1];
1045 res1 += val0 * x[mm - 1];
1046 for (int ind = 1; ind < mid; ++ind)
1047 {
1048 if (transpose_matrix == true)
1049 {
1050 val0 = matrix[ind * n_columns + col];
1051 val1 = matrix[ind * n_columns + nn - 1 - col];
1052 }
1053 else
1054 {
1055 val0 = matrix[col * n_columns + ind];
1056 val1 = matrix[(col + 1) * n_columns - 1 - ind];
1057 }
1058 res0 += val0 * x[ind];
1059 res1 += val1 * x[ind];
1060 res0 += val1 * x[mm - 1 - ind];
1061 res1 += val0 * x[mm - 1 - ind];
1062 }
1063 }
1064 else
1065 res0 = res1 = Number();
1066 if (mm % 2 == 1)
1067 {
1068 if (transpose_matrix == true)
1069 val0 = matrix[mid * n_columns + col];
1070 else
1071 val0 = matrix[col * n_columns + mid];
1072 const Number tmp = val0 * x[mid];
1073 res0 += tmp;
1074 res1 += tmp;
1075 }
1076 if (add)
1077 {
1078 out[stride_out * col] += res0;
1079 out[stride_out * (nn - 1 - col)] += res1;
1080 }
1081 else
1082 {
1083 out[stride_out * col] = res0;
1084 out[stride_out * (nn - 1 - col)] = res1;
1085 }
1086 }
1087 if (nn % 2 == 1)
1088 {
1089 Number2 val0;
1090 Number res0;
1091 if (transpose_matrix == true)
1092 val0 = matrix[n_cols];
1093 else
1094 val0 = matrix[n_cols * n_columns];
1095 if (mid > 0)
1096 {
1097 res0 = val0 * (x[0] + x[mm - 1]);
1098 for (int ind = 1; ind < mid; ++ind)
1099 {
1100 if (transpose_matrix == true)
1101 val0 = matrix[ind * n_columns + n_cols];
1102 else
1103 val0 = matrix[n_cols * n_columns + ind];
1104 Number in1 = val0 * (x[ind] + x[mm - 1 - ind]);
1105 res0 += in1;
1106 }
1107 }
1108 else
1109 res0 = Number();
1110 if (mm % 2 == 1)
1111 {
1112 if (transpose_matrix == true)
1113 val0 = matrix[mid * n_columns + n_cols];
1114 else
1115 val0 = matrix[n_cols * n_columns + mid];
1116 res0 += val0 * x[mid];
1117 }
1118 if (add)
1119 out[stride_out * n_cols] += res0;
1120 else
1121 out[stride_out * n_cols] = res0;
1122 }
1123 }
1124 }
1125
1126
1127
1146 template <EvaluatorVariant variant,
1147 EvaluatorQuantity quantity,
1148 int n_rows_static,
1149 int n_columns_static,
1150 int stride_in_static,
1151 int stride_out_static,
1152 bool transpose_matrix,
1153 bool add,
1154 typename Number,
1155 typename Number2>
1156#ifndef DEBUG
1157 inline DEAL_II_ALWAYS_INLINE
1158#endif
1159 std::enable_if_t<(variant == evaluate_evenodd), void>
1160 apply_matrix_vector_product(const Number2 *DEAL_II_RESTRICT matrix,
1161 const Number *in,
1162 Number *out,
1163 int n_rows_runtime = 0,
1164 int n_columns_runtime = 0,
1165 int stride_in_runtime = 0,
1166 int stride_out_runtime = 0)
1167 {
1168 static_assert(n_rows_static >= 0 && n_columns_static >= 0,
1169 "Negative loop ranges are not allowed!");
1170
1171 const int n_rows = n_rows_static == 0 ? n_rows_runtime : n_rows_static;
1172 const int n_columns =
1173 n_rows_static == 0 ? n_columns_runtime : n_columns_static;
1174 const int stride_in =
1175 stride_in_static == 0 ? stride_in_runtime : stride_in_static;
1176 const int stride_out =
1177 stride_out_static == 0 ? stride_out_runtime : stride_out_static;
1178
1179 Assert(n_rows > 0 && n_columns > 0,
1180 ExcInternalError("The evaluation needs n_rows, n_columns > 0, but " +
1181 std::to_string(n_rows) + ", " +
1182 std::to_string(n_columns) + " was passed!"));
1183
1184 const int mm = transpose_matrix ? n_rows : n_columns,
1185 nn = transpose_matrix ? n_columns : n_rows;
1186 const int n_half = nn / 2;
1187 const int m_half = mm / 2;
1188
1189 constexpr int array_length =
1190 (n_rows_static == 0) ?
1191 16 // for non-templated execution
1192 :
1193 (1 + (transpose_matrix ? n_rows_static : n_columns_static) / 2);
1194 const int offset = (n_columns + 1) / 2;
1195
1196 Assert(m_half <= array_length, ExcNotImplemented());
1197
1198 std::array<Number, array_length> xp, xm;
1199 for (int i = 0; i < m_half; ++i)
1200 {
1201 if (transpose_matrix == true && quantity == EvaluatorQuantity::gradient)
1202 {
1203 xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1204 xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1205 }
1206 else
1207 {
1208 xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1209 xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1210 }
1211 }
1212 Number xmid = in[stride_in * m_half];
1213 for (int col = 0; col < n_half; ++col)
1214 {
1215 Number r0, r1;
1216 if (m_half > 0)
1217 {
1218 if (transpose_matrix == true)
1219 {
1220 r0 = matrix[col] * xp[0];
1221 r1 = matrix[(n_rows - 1) * offset + col] * xm[0];
1222 }
1223 else
1224 {
1225 r0 = matrix[col * offset] * xp[0];
1226 r1 = matrix[(n_rows - 1 - col) * offset] * xm[0];
1227 }
1228 for (int ind = 1; ind < m_half; ++ind)
1229 {
1230 if (transpose_matrix == true)
1231 {
1232 r0 += matrix[ind * offset + col] * xp[ind];
1233 r1 += matrix[(n_rows - 1 - ind) * offset + col] * xm[ind];
1234 }
1235 else
1236 {
1237 r0 += matrix[col * offset + ind] * xp[ind];
1238 r1 += matrix[(n_rows - 1 - col) * offset + ind] * xm[ind];
1239 }
1240 }
1241 }
1242 else
1243 r0 = r1 = Number();
1244 if (mm % 2 == 1 && transpose_matrix == true)
1245 {
1246 if (quantity == EvaluatorQuantity::gradient)
1247 r1 += matrix[m_half * offset + col] * xmid;
1248 else
1249 r0 += matrix[m_half * offset + col] * xmid;
1250 }
1251 else if (mm % 2 == 1 &&
1252 (nn % 2 == 0 || quantity != EvaluatorQuantity::value ||
1253 mm == 3))
1254 r0 += matrix[col * offset + m_half] * xmid;
1255
1256 if (add)
1257 {
1258 out[stride_out * col] += r0 + r1;
1259 if (quantity == EvaluatorQuantity::gradient &&
1260 transpose_matrix == false)
1261 out[stride_out * (nn - 1 - col)] += r1 - r0;
1262 else
1263 out[stride_out * (nn - 1 - col)] += r0 - r1;
1264 }
1265 else
1266 {
1267 out[stride_out * col] = r0 + r1;
1268 if (quantity == EvaluatorQuantity::gradient &&
1269 transpose_matrix == false)
1270 out[stride_out * (nn - 1 - col)] = r1 - r0;
1271 else
1272 out[stride_out * (nn - 1 - col)] = r0 - r1;
1273 }
1274 }
1275 if (quantity == EvaluatorQuantity::value && transpose_matrix == true &&
1276 nn % 2 == 1 && mm % 2 == 1 && mm > 3)
1277 {
1278 if (add)
1279 out[stride_out * n_half] += matrix[m_half * offset + n_half] * xmid;
1280 else
1281 out[stride_out * n_half] = matrix[m_half * offset + n_half] * xmid;
1282 }
1283 else if (transpose_matrix == true && nn % 2 == 1)
1284 {
1285 Number r0;
1286 if (m_half > 0)
1287 {
1288 r0 = matrix[n_half] * xp[0];
1289 for (int ind = 1; ind < m_half; ++ind)
1290 r0 += matrix[ind * offset + n_half] * xp[ind];
1291 }
1292 else
1293 r0 = Number();
1294 if (quantity != EvaluatorQuantity::gradient && mm % 2 == 1)
1295 r0 += matrix[m_half * offset + n_half] * xmid;
1296
1297 if (add)
1298 out[stride_out * n_half] += r0;
1299 else
1300 out[stride_out * n_half] = r0;
1301 }
1302 else if (transpose_matrix == false && nn % 2 == 1)
1303 {
1304 Number r0;
1305 if (m_half > 0)
1306 {
1307 if (quantity == EvaluatorQuantity::gradient)
1308 {
1309 r0 = matrix[n_half * offset] * xm[0];
1310 for (int ind = 1; ind < m_half; ++ind)
1311 r0 += matrix[n_half * offset + ind] * xm[ind];
1312 }
1313 else
1314 {
1315 r0 = matrix[n_half * offset] * xp[0];
1316 for (int ind = 1; ind < m_half; ++ind)
1317 r0 += matrix[n_half * offset + ind] * xp[ind];
1318 }
1319 }
1320 else
1321 r0 = Number();
1322
1323 if (quantity != EvaluatorQuantity::gradient && mm % 2 == 1)
1324 r0 += matrix[n_half * offset + m_half] * xmid;
1325
1326 if (add)
1327 out[stride_out * n_half] += r0;
1328 else
1329 out[stride_out * n_half] = r0;
1330 }
1331 }
1332
1333
1334
1339 template <EvaluatorVariant variant,
1340 EvaluatorQuantity quantity,
1341 bool transpose_matrix,
1342 bool add,
1343 bool consider_strides,
1344 typename Number,
1345 typename Number2>
1346 std::enable_if_t<(variant == evaluate_evenodd), void>
1347 apply_matrix_vector_product(const Number2 *matrix,
1348 const Number *in,
1349 Number *out,
1350 int n_rows,
1351 int n_columns,
1352 int stride_in,
1353 int stride_out)
1354 {
1355 apply_matrix_vector_product<evaluate_evenodd,
1356 quantity,
1357 0,
1358 0,
1359 consider_strides ? 0 : 1,
1360 consider_strides ? 0 : 1,
1361 transpose_matrix,
1362 add>(
1363 matrix, in, out, n_rows, n_columns, stride_in, stride_out);
1364 }
1365
1366
1367
1383 template <EvaluatorVariant variant,
1384 EvaluatorQuantity quantity,
1385 int n_rows,
1386 int n_columns,
1387 int stride_in,
1388 int stride_out,
1389 bool transpose_matrix,
1390 bool add,
1391 typename Number,
1392 typename Number2>
1393 std::enable_if_t<(variant == evaluate_symmetric_hierarchical), void>
1394 apply_matrix_vector_product(const Number2 *matrix,
1395 const Number *in,
1396 Number *out)
1397 {
1398 static_assert(n_rows > 0 && n_columns > 0,
1399 "Specialization requires n_rows, n_columns > 0");
1400
1401 constexpr bool evaluate_antisymmetric =
1402 (quantity == EvaluatorQuantity::gradient);
1403
1404 constexpr int mm = transpose_matrix ? n_rows : n_columns,
1405 nn = transpose_matrix ? n_columns : n_rows;
1406 constexpr int n_half = nn / 2;
1407 constexpr int m_half = mm / 2;
1408
1409 if (transpose_matrix)
1410 {
1411 std::array<Number, mm> x;
1412 for (unsigned int i = 0; i < mm; ++i)
1413 x[i] = in[stride_in * i];
1414 for (unsigned int col = 0; col < n_half; ++col)
1415 {
1416 Number r0, r1;
1417 if (m_half > 0)
1418 {
1419 r0 = matrix[col] * x[0];
1420 r1 = matrix[col + n_columns] * x[1];
1421 for (unsigned int ind = 1; ind < m_half; ++ind)
1422 {
1423 r0 += matrix[col + 2 * ind * n_columns] * x[2 * ind];
1424 r1 +=
1425 matrix[col + (2 * ind + 1) * n_columns] * x[2 * ind + 1];
1426 }
1427 }
1428 else
1429 r0 = r1 = Number();
1430 if (mm % 2 == 1)
1431 r0 += matrix[col + (mm - 1) * n_columns] * x[mm - 1];
1432 if (add)
1433 {
1434 out[stride_out * col] += r0 + r1;
1435 if (evaluate_antisymmetric)
1436 out[stride_out * (nn - 1 - col)] += r1 - r0;
1437 else
1438 out[stride_out * (nn - 1 - col)] += r0 - r1;
1439 }
1440 else
1441 {
1442 out[stride_out * col] = r0 + r1;
1443 if (evaluate_antisymmetric)
1444 out[stride_out * (nn - 1 - col)] = r1 - r0;
1445 else
1446 out[stride_out * (nn - 1 - col)] = r0 - r1;
1447 }
1448 }
1449 if (nn % 2 == 1)
1450 {
1451 Number r0;
1452 const unsigned int shift = evaluate_antisymmetric ? 1 : 0;
1453 if (m_half > 0)
1454 {
1455 r0 = matrix[n_half + shift * n_columns] * x[shift];
1456 for (unsigned int ind = 1; ind < m_half; ++ind)
1457 r0 += matrix[n_half + (2 * ind + shift) * n_columns] *
1458 x[2 * ind + shift];
1459 }
1460 else
1461 r0 = 0;
1462 if (!evaluate_antisymmetric && mm % 2 == 1)
1463 r0 += matrix[n_half + (mm - 1) * n_columns] * x[mm - 1];
1464 if (add)
1465 out[stride_out * n_half] += r0;
1466 else
1467 out[stride_out * n_half] = r0;
1468 }
1469 }
1470 else
1471 {
1472 std::array<Number, m_half + 1> xp, xm;
1473 for (int i = 0; i < m_half; ++i)
1474 if (!evaluate_antisymmetric)
1475 {
1476 xp[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1477 xm[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1478 }
1479 else
1480 {
1481 xp[i] = in[stride_in * i] - in[stride_in * (mm - 1 - i)];
1482 xm[i] = in[stride_in * i] + in[stride_in * (mm - 1 - i)];
1483 }
1484 if (mm % 2 == 1)
1485 xp[m_half] = in[stride_in * m_half];
1486 for (unsigned int col = 0; col < n_half; ++col)
1487 {
1488 Number r0, r1;
1489 if (m_half > 0)
1490 {
1491 r0 = matrix[2 * col * n_columns] * xp[0];
1492 r1 = matrix[(2 * col + 1) * n_columns] * xm[0];
1493 for (unsigned int ind = 1; ind < m_half; ++ind)
1494 {
1495 r0 += matrix[2 * col * n_columns + ind] * xp[ind];
1496 r1 += matrix[(2 * col + 1) * n_columns + ind] * xm[ind];
1497 }
1498 }
1499 else
1500 r0 = r1 = Number();
1501 if (mm % 2 == 1)
1502 {
1503 if (evaluate_antisymmetric)
1504 r1 += matrix[(2 * col + 1) * n_columns + m_half] * xp[m_half];
1505 else
1506 r0 += matrix[2 * col * n_columns + m_half] * xp[m_half];
1507 }
1508 if (add)
1509 {
1510 out[stride_out * (2 * col)] += r0;
1511 out[stride_out * (2 * col + 1)] += r1;
1512 }
1513 else
1514 {
1515 out[stride_out * (2 * col)] = r0;
1516 out[stride_out * (2 * col + 1)] = r1;
1517 }
1518 }
1519 if (nn % 2 == 1)
1520 {
1521 Number r0;
1522 if (m_half > 0)
1523 {
1524 r0 = matrix[(nn - 1) * n_columns] * xp[0];
1525 for (unsigned int ind = 1; ind < m_half; ++ind)
1526 r0 += matrix[(nn - 1) * n_columns + ind] * xp[ind];
1527 }
1528 else
1529 r0 = Number();
1530 if (mm % 2 == 1 && !evaluate_antisymmetric)
1531 r0 += matrix[(nn - 1) * n_columns + m_half] * xp[m_half];
1532 if (add)
1533 out[stride_out * (nn - 1)] += r0;
1534 else
1535 out[stride_out * (nn - 1)] = r0;
1536 }
1537 }
1538 }
1539
1540
1541
1564 template <EvaluatorVariant variant,
1565 int dim,
1566 int n_rows,
1567 int n_columns,
1568 typename Number,
1569 typename Number2 = Number>
1570 struct EvaluatorTensorProduct
1571 {
1572 static constexpr unsigned int n_rows_of_product =
1573 Utilities::pow(n_rows, dim);
1574 static constexpr unsigned int n_columns_of_product =
1575 Utilities::pow(n_columns, dim);
1576
1581 EvaluatorTensorProduct()
1582 : shape_values(nullptr)
1583 , shape_gradients(nullptr)
1584 , shape_hessians(nullptr)
1585 {}
1586
1590 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
1591 const AlignedVector<Number2> &shape_gradients,
1592 const AlignedVector<Number2> &shape_hessians,
1593 const unsigned int = 0,
1594 const unsigned int = 0)
1595 : shape_values(shape_values.begin())
1596 , shape_gradients(shape_gradients.begin())
1597 , shape_hessians(shape_hessians.begin())
1598 {
1599 if (variant == evaluate_evenodd)
1600 {
1601 if (!shape_values.empty())
1602 AssertDimension(shape_values.size(),
1603 n_rows * ((n_columns + 1) / 2));
1604 if (!shape_gradients.empty())
1605 AssertDimension(shape_gradients.size(),
1606 n_rows * ((n_columns + 1) / 2));
1607 if (!shape_hessians.empty())
1608 AssertDimension(shape_hessians.size(),
1609 n_rows * ((n_columns + 1) / 2));
1610 }
1611 else
1612 {
1613 Assert(shape_values.empty() ||
1614 shape_values.size() == n_rows * n_columns,
1615 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1616 Assert(shape_gradients.empty() ||
1617 shape_gradients.size() == n_rows * n_columns,
1618 ExcDimensionMismatch(shape_gradients.size(),
1619 n_rows * n_columns));
1620 Assert(shape_hessians.empty() ||
1621 shape_hessians.size() == n_rows * n_columns,
1622 ExcDimensionMismatch(shape_hessians.size(),
1623 n_rows * n_columns));
1624 }
1625 }
1626
1630 EvaluatorTensorProduct(const Number2 *shape_values,
1631 const Number2 *shape_gradients,
1632 const Number2 *shape_hessians,
1633 const unsigned int dummy1 = 0,
1634 const unsigned int dummy2 = 0)
1635 : shape_values(shape_values)
1636 , shape_gradients(shape_gradients)
1637 , shape_hessians(shape_hessians)
1638 {
1639 (void)dummy1;
1640 (void)dummy2;
1641 }
1642
1668 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1669 void
1670 values(const Number in[], Number out[]) const
1671 {
1672 constexpr EvaluatorQuantity value_type = EvaluatorQuantity::value;
1673 apply<direction, contract_over_rows, add, false, value_type, stride>(
1674 shape_values, in, out);
1675 }
1676
1682 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1683 void
1684 gradients(const Number in[], Number out[]) const
1685 {
1686 constexpr EvaluatorQuantity gradient_type =
1687 (variant == evaluate_general ? EvaluatorQuantity::value :
1688 EvaluatorQuantity::gradient);
1689 apply<direction, contract_over_rows, add, false, gradient_type, stride>(
1690 shape_gradients, in, out);
1691 }
1692
1698 template <int direction, bool contract_over_rows, bool add>
1699 void
1700 hessians(const Number in[], Number out[]) const
1701 {
1702 constexpr EvaluatorQuantity hessian_type =
1703 (((variant == evaluate_general) |
1704 (variant == evaluate_symmetric_hierarchical)) ?
1705 EvaluatorQuantity::value :
1706 EvaluatorQuantity::hessian);
1707 apply<direction, contract_over_rows, add, false, hessian_type>(
1708 shape_hessians, in, out);
1709 }
1710
1718 template <int direction, bool contract_over_rows, bool add>
1719 void
1720 values_one_line(const Number in[], Number out[]) const
1721 {
1722 Assert(shape_values != nullptr, ExcNotInitialized());
1723 apply<direction, contract_over_rows, add, true, EvaluatorQuantity::value>(
1724 shape_values, in, out);
1725 }
1726
1734 template <int direction, bool contract_over_rows, bool add>
1735 void
1736 gradients_one_line(const Number in[], Number out[]) const
1737 {
1738 Assert(shape_gradients != nullptr, ExcNotInitialized());
1739 constexpr EvaluatorQuantity gradient_type =
1740 (variant == evaluate_general ? EvaluatorQuantity::value :
1741 EvaluatorQuantity::gradient);
1742 apply<direction, contract_over_rows, add, true, gradient_type>(
1743 shape_gradients, in, out);
1744 }
1745
1753 template <int direction, bool contract_over_rows, bool add>
1754 void
1755 hessians_one_line(const Number in[], Number out[]) const
1756 {
1757 Assert(shape_hessians != nullptr, ExcNotInitialized());
1758 constexpr EvaluatorQuantity hessian_type =
1759 (((variant == evaluate_general) |
1760 (variant == evaluate_symmetric_hierarchical)) ?
1761 EvaluatorQuantity::value :
1762 EvaluatorQuantity::hessian);
1763 apply<direction, contract_over_rows, add, true, hessian_type>(
1764 shape_hessians, in, out);
1765 }
1766
1803 template <int direction,
1804 bool contract_over_rows,
1805 bool add,
1806 bool one_line = false,
1807 EvaluatorQuantity quantity = EvaluatorQuantity::value,
1808 int stride = 1>
1809 static void
1810 apply(const Number2 *DEAL_II_RESTRICT shape_data,
1811 const Number *in,
1812 Number *out);
1813
1814 private:
1815 const Number2 *shape_values;
1816 const Number2 *shape_gradients;
1817 const Number2 *shape_hessians;
1818 };
1819
1820
1821
1822 template <EvaluatorVariant variant,
1823 int dim,
1824 int n_rows,
1825 int n_columns,
1826 typename Number,
1827 typename Number2>
1828 template <int direction,
1829 bool contract_over_rows,
1830 bool add,
1831 bool one_line,
1832 EvaluatorQuantity quantity,
1833 int stride>
1834 inline void
1835 EvaluatorTensorProduct<variant, dim, n_rows, n_columns, Number, Number2>::
1836 apply(const Number2 *DEAL_II_RESTRICT shape_data,
1837 const Number *in,
1838 Number *out)
1839 {
1840 static_assert(one_line == false || direction == dim - 1,
1841 "Single-line evaluation only works for direction=dim-1.");
1842 Assert(shape_data != nullptr,
1843 ExcMessage(
1844 "The given array shape_data must not be the null pointer!"));
1845 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
1846 in != out,
1847 ExcMessage("In-place operation only supported for "
1848 "n_rows==n_columns or single-line interpolation"));
1849 AssertIndexRange(direction, dim);
1850 constexpr int mm = contract_over_rows ? n_rows : n_columns,
1851 nn = contract_over_rows ? n_columns : n_rows;
1852
1853 constexpr int stride_operation = Utilities::pow(n_columns, direction);
1854 constexpr int n_blocks1 = one_line ? 1 : stride_operation;
1855 constexpr int n_blocks2 =
1856 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1857
1858 constexpr int stride_in = !contract_over_rows ? stride : 1;
1859 constexpr int stride_out = contract_over_rows ? stride : 1;
1860 for (int i2 = 0; i2 < n_blocks2; ++i2)
1861 {
1862 for (int i1 = 0; i1 < n_blocks1; ++i1)
1863 {
1864 apply_matrix_vector_product<variant,
1865 quantity,
1866 n_rows,
1867 n_columns,
1868 stride_operation * stride_in,
1869 stride_operation * stride_out,
1870 contract_over_rows,
1871 add>(shape_data, in, out);
1872
1873 if (one_line == false)
1874 {
1875 in += stride_in;
1876 out += stride_out;
1877 }
1878 }
1879 if (one_line == false)
1880 {
1881 in += stride_operation * (mm - 1) * stride_in;
1882 out += stride_operation * (nn - 1) * stride_out;
1883 }
1884 }
1885 }
1886
1887
1888
1902 template <EvaluatorVariant variant,
1903 int dim,
1904 typename Number,
1905 typename Number2>
1906 struct EvaluatorTensorProduct<variant, dim, 0, 0, Number, Number2>
1907 {
1908 static constexpr unsigned int n_rows_of_product =
1909 numbers::invalid_unsigned_int;
1910 static constexpr unsigned int n_columns_of_product =
1911 numbers::invalid_unsigned_int;
1912
1917 EvaluatorTensorProduct()
1918 : shape_values(nullptr)
1919 , shape_gradients(nullptr)
1920 , shape_hessians(nullptr)
1921 , n_rows(numbers::invalid_unsigned_int)
1922 , n_columns(numbers::invalid_unsigned_int)
1923 {}
1924
1928 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
1929 const AlignedVector<Number2> &shape_gradients,
1930 const AlignedVector<Number2> &shape_hessians,
1931 const unsigned int n_rows = 0,
1932 const unsigned int n_columns = 0)
1933 : shape_values(shape_values.begin())
1934 , shape_gradients(shape_gradients.begin())
1935 , shape_hessians(shape_hessians.begin())
1936 , n_rows(n_rows)
1937 , n_columns(n_columns)
1938 {
1939 if (variant == evaluate_evenodd)
1940 {
1941 if (!shape_values.empty())
1942 AssertDimension(shape_values.size(),
1943 n_rows * ((n_columns + 1) / 2));
1944 if (!shape_gradients.empty())
1945 AssertDimension(shape_gradients.size(),
1946 n_rows * ((n_columns + 1) / 2));
1947 if (!shape_hessians.empty())
1948 AssertDimension(shape_hessians.size(),
1949 n_rows * ((n_columns + 1) / 2));
1950 }
1951 else
1952 {
1953 Assert(shape_values.empty() ||
1954 shape_values.size() == n_rows * n_columns,
1955 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1956 Assert(shape_gradients.empty() ||
1957 shape_gradients.size() == n_rows * n_columns,
1958 ExcDimensionMismatch(shape_gradients.size(),
1959 n_rows * n_columns));
1960 Assert(shape_hessians.empty() ||
1961 shape_hessians.size() == n_rows * n_columns,
1962 ExcDimensionMismatch(shape_hessians.size(),
1963 n_rows * n_columns));
1964 }
1965 }
1966
1970 EvaluatorTensorProduct(const Number2 *shape_values,
1971 const Number2 *shape_gradients,
1972 const Number2 *shape_hessians,
1973 const unsigned int n_rows = 0,
1974 const unsigned int n_columns = 0)
1975 : shape_values(shape_values)
1976 , shape_gradients(shape_gradients)
1977 , shape_hessians(shape_hessians)
1978 , n_rows(n_rows)
1979 , n_columns(n_columns)
1980 {}
1981
1982 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1983 void
1984 values(const Number *in, Number *out) const
1985 {
1986 constexpr EvaluatorQuantity value_type = EvaluatorQuantity::value;
1987 apply<direction, contract_over_rows, add, false, value_type, stride>(
1988 shape_values, in, out);
1989 }
1990
1991 template <int direction, bool contract_over_rows, bool add, int stride = 1>
1992 void
1993 gradients(const Number *in, Number *out) const
1994 {
1995 constexpr EvaluatorQuantity gradient_type =
1996 (variant != evaluate_evenodd ? EvaluatorQuantity::value :
1997 EvaluatorQuantity::gradient);
1998 apply<direction, contract_over_rows, add, false, gradient_type, stride>(
1999 shape_gradients, in, out);
2000 }
2001
2002 template <int direction, bool contract_over_rows, bool add>
2003 void
2004 hessians(const Number *in, Number *out) const
2005 {
2006 constexpr EvaluatorQuantity hessian_type =
2007 (variant != evaluate_evenodd ? EvaluatorQuantity::value :
2008 EvaluatorQuantity::hessian);
2009 apply<direction, contract_over_rows, add, false, hessian_type>(
2010 shape_hessians, in, out);
2011 }
2012
2013 template <int direction, bool contract_over_rows, bool add>
2014 void
2015 values_one_line(const Number in[], Number out[]) const
2016 {
2017 Assert(shape_values != nullptr, ExcNotInitialized());
2018 apply<direction, contract_over_rows, add, true, EvaluatorQuantity::value>(
2019 shape_values, in, out);
2020 }
2021
2022 template <int direction, bool contract_over_rows, bool add>
2023 void
2024 gradients_one_line(const Number in[], Number out[]) const
2025 {
2026 Assert(shape_gradients != nullptr, ExcNotInitialized());
2027 constexpr EvaluatorQuantity gradient_type =
2028 (variant != evaluate_evenodd ? EvaluatorQuantity::value :
2029 EvaluatorQuantity::gradient);
2030 apply<direction, contract_over_rows, add, true, gradient_type>(
2031 shape_gradients, in, out);
2032 }
2033
2034 template <int direction, bool contract_over_rows, bool add>
2035 void
2036 hessians_one_line(const Number in[], Number out[]) const
2037 {
2038 Assert(shape_hessians != nullptr, ExcNotInitialized());
2039 constexpr EvaluatorQuantity hessian_type =
2040 (variant != evaluate_evenodd ? EvaluatorQuantity::value :
2041 EvaluatorQuantity::hessian);
2042 apply<direction, contract_over_rows, add, true, hessian_type>(
2043 shape_hessians, in, out);
2044 }
2045
2046 template <int direction,
2047 bool contract_over_rows,
2048 bool add,
2049 bool one_line = false,
2050 EvaluatorQuantity quantity = EvaluatorQuantity::value,
2051 int stride = 1>
2052 void
2053 apply(const Number2 *DEAL_II_RESTRICT shape_data,
2054 const Number *in,
2055 Number *out) const;
2056
2057 const Number2 *shape_values;
2058 const Number2 *shape_gradients;
2059 const Number2 *shape_hessians;
2060 const unsigned int n_rows;
2061 const unsigned int n_columns;
2062 };
2063
2064
2065
2066 template <EvaluatorVariant variant,
2067 int dim,
2068 typename Number,
2069 typename Number2>
2070 template <int direction,
2071 bool contract_over_rows,
2072 bool add,
2073 bool one_line,
2074 EvaluatorQuantity quantity,
2075 int stride>
2076 inline void
2077 EvaluatorTensorProduct<variant, dim, 0, 0, Number, Number2>::apply(
2078 const Number2 *DEAL_II_RESTRICT shape_data,
2079 const Number *in,
2080 Number *out) const
2081 {
2082 static_assert(one_line == false || direction == dim - 1,
2083 "Single-line evaluation only works for direction=dim-1.");
2084 Assert(shape_data != nullptr,
2085 ExcMessage(
2086 "The given array shape_data must not be the null pointer!"));
2087 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
2088 in != out,
2089 ExcMessage("In-place operation only supported for "
2090 "n_rows==n_columns or single-line interpolation"));
2091 AssertIndexRange(direction, dim);
2092 const int mm = contract_over_rows ? n_rows : n_columns,
2093 nn = contract_over_rows ? n_columns : n_rows;
2094
2095 const int stride_operation =
2096 direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
2097 const int n_blocks1 = one_line ? 1 : stride_operation;
2098 const int n_blocks2 = direction >= dim - 1 ?
2099 1 :
2100 Utilities::fixed_power<dim - direction - 1>(n_rows);
2101 Assert(n_rows <= 128, ExcNotImplemented());
2102
2103 constexpr int stride_in = !contract_over_rows ? stride : 1;
2104 constexpr int stride_out = contract_over_rows ? stride : 1;
2105 for (int i2 = 0; i2 < n_blocks2; ++i2)
2106 {
2107 for (int i1 = 0; i1 < n_blocks1; ++i1)
2108 {
2109 // the empty template case can only run the general evaluator or
2110 // evenodd
2111 constexpr EvaluatorVariant restricted_variant =
2112 variant == evaluate_evenodd ? evaluate_evenodd : evaluate_general;
2113 apply_matrix_vector_product<restricted_variant,
2114 quantity,
2115 contract_over_rows,
2116 add,
2117 (direction != 0 || stride != 1)>(
2118 shape_data,
2119 in,
2120 out,
2121 n_rows,
2122 n_columns,
2123 stride_operation * stride_in,
2124 stride_operation * stride_out);
2125
2126 if (one_line == false)
2127 {
2128 in += stride_in;
2129 out += stride_out;
2130 }
2131 }
2132 if (one_line == false)
2133 {
2134 in += stride_operation * (mm - 1) * stride_in;
2135 out += stride_operation * (nn - 1) * stride_out;
2136 }
2137 }
2138 }
2139
2140
2141
2142 template <int dim,
2143 int fe_degree,
2144 int n_q_points_1d,
2145 bool contract_over_rows,
2146 MatrixFreeFunctions::ElementType element_type =
2147 MatrixFreeFunctions::tensor_raviart_thomas,
2148 bool symmetric_evaluate = true>
2149 struct EvaluatorTensorProductAnisotropic
2150 {
2151 template <int direction,
2152 int stride = 1,
2153 typename Number = double,
2154 typename Number2 = double>
2155 static void
2156 normal(const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
2157 const Number *in,
2158 Number *out,
2159 const bool add_into_result = false,
2160 const int subface_index_1d = 0)
2161 {
2162 AssertIndexRange(direction, dim);
2163 AssertDimension(fe_degree, data.fe_degree);
2164 AssertDimension(n_q_points_1d, data.n_q_points_1d);
2165 constexpr int n_rows = fe_degree + 1;
2166 constexpr int n_columns = n_q_points_1d;
2167 constexpr int mm = contract_over_rows ? n_rows : n_columns;
2168 constexpr int nn = contract_over_rows ? n_columns : n_rows;
2169 const Number2 *shape_data =
2170 symmetric_evaluate ?
2171 data.shape_values_eo.data() :
2172 data.values_within_subface[subface_index_1d].data();
2173 Assert(shape_data != nullptr, ExcNotInitialized());
2174 Assert(contract_over_rows == false || !add_into_result,
2175 ExcMessage("Cannot add into result if contract_over_rows = true"));
2176
2177 constexpr int n_blocks1 =
2178 Utilities::pow((element_type ==
2179 MatrixFreeFunctions::tensor_raviart_thomas) ?
2180 fe_degree :
2181 fe_degree + 2,
2182 direction);
2183 constexpr int n_blocks2 =
2184 Utilities::pow((element_type ==
2185 MatrixFreeFunctions::tensor_raviart_thomas) ?
2186 fe_degree :
2187 fe_degree + 2,
2188 dim - direction - 1);
2189 constexpr int stride_in = contract_over_rows ? 1 : stride;
2190 constexpr int stride_out = contract_over_rows ? stride : 1;
2191 constexpr EvaluatorVariant variant =
2192 symmetric_evaluate ? evaluate_evenodd : evaluate_general;
2193
2194 for (int i2 = 0; i2 < n_blocks2; ++i2)
2195 {
2196 for (int i1 = 0; i1 < n_blocks1; ++i1)
2197 {
2198 if (contract_over_rows == false && add_into_result)
2199 apply_matrix_vector_product<variant,
2200 EvaluatorQuantity::value,
2201 n_rows,
2202 n_columns,
2203 n_blocks1 * stride_in,
2204 n_blocks1 * stride_out,
2205 contract_over_rows,
2206 true>(shape_data, in, out);
2207 else
2208 apply_matrix_vector_product<variant,
2209 EvaluatorQuantity::value,
2210 n_rows,
2211 n_columns,
2212 n_blocks1 * stride_in,
2213 n_blocks1 * stride_out,
2214 contract_over_rows,
2215 false>(shape_data, in, out);
2216
2217 in += stride_in;
2218 out += stride_out;
2219 }
2220 in += n_blocks1 * (mm - 1) * stride_in;
2221 out += n_blocks1 * (nn - 1) * stride_out;
2222 }
2223 }
2224
2225 template <int direction,
2226 int normal_direction,
2227 int stride = 1,
2228 typename Number = double,
2229 typename Number2 = double>
2230 static void
2231 tangential(const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
2232 const Number *in,
2233 Number *out,
2234 const int subface_index_1d = 0)
2235 {
2236 AssertIndexRange(direction, dim);
2237 AssertDimension((element_type ==
2238 MatrixFreeFunctions::tensor_raviart_thomas) ?
2239 fe_degree - 1 :
2240 fe_degree + 1,
2241 data.fe_degree);
2242 AssertDimension(n_q_points_1d, data.n_q_points_1d);
2243 static_assert(direction != normal_direction,
2244 "Cannot interpolate tangentially in normal direction");
2245
2246 constexpr int internal_dof =
2247 (element_type == MatrixFreeFunctions::tensor_raviart_thomas) ?
2248 fe_degree :
2249 fe_degree + 2;
2250 constexpr int n_rows = std::max(internal_dof, 0);
2251 constexpr int n_columns = n_q_points_1d;
2252 const Number2 *shape_data =
2253 symmetric_evaluate ?
2254 data.shape_values_eo.data() :
2255 data.values_within_subface[subface_index_1d].data();
2256 Assert(shape_data != nullptr, ExcNotInitialized());
2257
2258 constexpr int n_blocks1 =
2259 (direction > normal_direction) ?
2260 Utilities::pow(n_q_points_1d, direction) :
2261 (direction > 0 ?
2262 (Utilities::pow(internal_dof, direction - 1) * n_q_points_1d) :
2263 1);
2264 constexpr int n_blocks2 =
2265 (direction > normal_direction) ?
2266 Utilities::pow(internal_dof, dim - 1 - direction) :
2267 ((direction + 1 < dim) ?
2268 (Utilities::pow(internal_dof, dim - 2 - direction) *
2269 n_q_points_1d) :
2270 1);
2271
2272 constexpr EvaluatorVariant variant =
2273 symmetric_evaluate ? evaluate_evenodd : evaluate_general;
2274
2275 // Since we may perform an in-place interpolation, we must run the step
2276 // expanding the size of the basis backward ('contract_over_rows' aka
2277 // 'evaluate' case), so shift the pointers and decrement during the loop
2278 if (contract_over_rows)
2279 {
2280 in += (n_blocks2 - 1) * n_blocks1 * n_rows + n_blocks1 - 1;
2281 out +=
2282 stride * ((n_blocks2 - 1) * n_blocks1 * n_columns + n_blocks1 - 1);
2283 for (int i2 = 0; i2 < n_blocks2; ++i2)
2284 {
2285 for (int i1 = 0; i1 < n_blocks1; ++i1)
2286 {
2287 apply_matrix_vector_product<variant,
2288 EvaluatorQuantity::value,
2289 n_rows,
2290 n_columns,
2291 n_blocks1,
2292 n_blocks1 * stride,
2293 true,
2294 false>(shape_data, in, out);
2295
2296 --in;
2297 out -= stride;
2298 }
2299 in -= n_blocks1 * (n_rows - 1);
2300 out -= n_blocks1 * (n_columns - 1) * stride;
2301 }
2302 }
2303 else
2304 {
2305 for (int i2 = 0; i2 < n_blocks2; ++i2)
2306 {
2307 for (int i1 = 0; i1 < n_blocks1; ++i1)
2308 {
2309 apply_matrix_vector_product<variant,
2310 EvaluatorQuantity::value,
2311 n_rows,
2312 n_columns,
2313 n_blocks1 * stride,
2314 n_blocks1,
2315 false,
2316 false>(shape_data, in, out);
2317
2318 in += stride;
2319 ++out;
2320 }
2321 in += n_blocks1 * (n_columns - 1) * stride;
2322 out += n_blocks1 * (n_rows - 1);
2323 }
2324 }
2325 }
2326 };
2327
2328
2329
2375 template <int n_rows_template,
2376 int stride_template,
2377 bool contract_onto_face,
2378 bool add,
2379 int max_derivative,
2380 typename Number,
2381 typename Number2>
2382 inline std::enable_if_t<contract_onto_face, void>
2383 interpolate_to_face(const Number2 *shape_values,
2384 const std::array<int, 2> &n_blocks,
2385 const std::array<int, 2> &steps,
2386 const Number *input,
2387 Number *DEAL_II_RESTRICT output,
2388 const int n_rows_runtime = 0,
2389 const int stride_runtime = 1)
2390 {
2391 const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
2392 const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
2393
2394 Number *output1 = output + n_blocks[0] * n_blocks[1];
2395 Number *output2 = output1 + n_blocks[0] * n_blocks[1];
2396 for (int i2 = 0; i2 < n_blocks[1]; ++i2)
2397 {
2398 for (int i1 = 0; i1 < n_blocks[0]; ++i1)
2399 {
2400 Number res0 = shape_values[0] * input[0];
2401 Number res1, res2;
2402 if (max_derivative > 0)
2403 res1 = shape_values[n_rows] * input[0];
2404 if (max_derivative > 1)
2405 res2 = shape_values[2 * n_rows] * input[0];
2406 for (int ind = 1; ind < n_rows; ++ind)
2407 {
2408 res0 += shape_values[ind] * input[stride * ind];
2409 if (max_derivative > 0)
2410 res1 += shape_values[ind + n_rows] * input[stride * ind];
2411 if (max_derivative > 1)
2412 res2 += shape_values[ind + 2 * n_rows] * input[stride * ind];
2413 }
2414 if (add)
2415 {
2416 output[i1] += res0;
2417 if (max_derivative > 0)
2418 output1[i1] += res1;
2419 if (max_derivative > 1)
2420 output2[i2] += res2;
2421 }
2422 else
2423 {
2424 output[i1] = res0;
2425 if (max_derivative > 0)
2426 output1[i1] = res1;
2427 if (max_derivative > 1)
2428 output2[i1] = res2;
2429 }
2430 input += steps[0];
2431 }
2432 output += n_blocks[0];
2433 if (max_derivative > 0)
2434 output1 += n_blocks[0];
2435 if (max_derivative > 1)
2436 output2 += n_blocks[0];
2437 input += steps[1];
2438 }
2439 }
2440
2441
2442
2450 constexpr bool
2451 use_collocation_evaluation(const unsigned int fe_degree,
2452 const unsigned int n_q_points_1d)
2453 {
2454 return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
2455 (n_q_points_1d <= 3 * fe_degree / 2 + 1);
2456 }
2457
2458
2459
2465 template <int n_rows_template,
2466 int stride_template,
2467 bool contract_onto_face,
2468 bool add,
2469 int max_derivative,
2470 typename Number,
2471 typename Number2>
2472 inline std::enable_if_t<!contract_onto_face, void>
2473 interpolate_to_face(const Number2 *shape_values,
2474 const std::array<int, 2> &n_blocks,
2475 const std::array<int, 2> &steps,
2476 const Number *input,
2477 Number *DEAL_II_RESTRICT output,
2478 const int n_rows_runtime = 0,
2479 const int stride_runtime = 1)
2480 {
2481 const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
2482 const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
2483
2484 const Number *input1 = input + n_blocks[0] * n_blocks[1];
2485 const Number *input2 = input1 + n_blocks[0] * n_blocks[1];
2486 for (int i2 = 0; i2 < n_blocks[1]; ++i2)
2487 {
2488 for (int i1 = 0; i1 < n_blocks[0]; ++i1)
2489 {
2490 const Number in = input[i1];
2491 Number in1, in2;
2492 if (max_derivative > 0)
2493 in1 = input1[i1];
2494 if (max_derivative > 1)
2495 in2 = input2[i1];
2496 for (int col = 0; col < n_rows; ++col)
2497 {
2498 Number result =
2499 add ? (output[col * stride] + shape_values[col] * in) :
2500 (shape_values[col] * in);
2501 if (max_derivative > 0)
2502 result += shape_values[col + n_rows] * in1;
2503 if (max_derivative > 1)
2504 result += shape_values[col + 2 * n_rows] * in2;
2505
2506 output[col * stride] = result;
2507 }
2508 output += steps[0];
2509 }
2510 input += n_blocks[0];
2511 if (max_derivative > 0)
2512 input1 += n_blocks[0];
2513 if (max_derivative > 1)
2514 input2 += n_blocks[0];
2515 output += steps[1];
2516 }
2517 }
2518
2519 template <int dim, int n_points_1d_template, typename Number>
2520 inline void
2521 weight_fe_q_dofs_by_entity(const Number *weights,
2522 const unsigned int n_components,
2523 const int n_points_1d_non_template,
2524 Number *data)
2525 {
2526 const int n_points_1d = n_points_1d_template != -1 ?
2527 n_points_1d_template :
2528 n_points_1d_non_template;
2529
2530 Assert(n_points_1d > 0, ExcNotImplemented());
2531 Assert(n_points_1d < 100, ExcNotImplemented());
2532
2533 unsigned int compressed_index[100];
2534 compressed_index[0] = 0;
2535 for (int i = 1; i < n_points_1d - 1; ++i)
2536 compressed_index[i] = 1;
2537 compressed_index[n_points_1d - 1] = 2;
2538
2539 for (unsigned int c = 0; c < n_components; ++c)
2540 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2541 for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
2542 {
2543 const unsigned int shift =
2544 9 * compressed_index[k] + 3 * compressed_index[j];
2545 data[0] *= weights[shift];
2546 // loop bound as int avoids compiler warnings in case n_points_1d
2547 // == 1 (polynomial degree 0)
2548 const Number weight = weights[shift + 1];
2549 for (int i = 1; i < n_points_1d - 1; ++i)
2550 data[i] *= weight;
2551 data[n_points_1d - 1] *= weights[shift + 2];
2552 data += n_points_1d;
2553 }
2554 }
2555
2556
2557 template <int dim, int n_points_1d_template, typename Number>
2558 inline void
2559 weight_fe_q_dofs_by_entity_shifted(const Number *weights,
2560 const unsigned int n_components,
2561 const int n_points_1d_non_template,
2562 Number *data)
2563 {
2564 const int n_points_1d = n_points_1d_template != -1 ?
2565 n_points_1d_template :
2566 n_points_1d_non_template;
2567
2568 Assert((n_points_1d % 2) == 1,
2569 ExcMessage("The function can only with add number of points"));
2570 Assert(n_points_1d > 0, ExcNotImplemented());
2571 Assert(n_points_1d < 100, ExcNotImplemented());
2572
2573 const unsigned int n_inside_1d = n_points_1d / 2;
2574
2575 unsigned int compressed_index[100];
2576
2577 unsigned int c = 0;
2578 for (int i = 0; i < n_inside_1d; ++i)
2579 compressed_index[c++] = 0;
2580 compressed_index[c++] = 1;
2581 for (int i = 0; i < n_inside_1d; ++i)
2582 compressed_index[c++] = 2;
2583
2584 for (unsigned int c = 0; c < n_components; ++c)
2585 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2586 for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
2587 {
2588 const unsigned int shift =
2589 9 * compressed_index[k] + 3 * compressed_index[j];
2590
2591 unsigned int c = 0;
2592 const Number weight1 = weights[shift];
2593 for (int i = 0; i < n_inside_1d; ++i)
2594 data[c++] *= weight1;
2595 data[c++] *= weights[shift + 1];
2596 const Number weight2 = weights[shift + 2];
2597 for (int i = 0; i < n_inside_1d; ++i)
2598 data[c++] *= weight2;
2599 data += n_points_1d;
2600 }
2601 }
2602
2603
2604 template <int dim, int n_points_1d_template, typename Number>
2605 inline bool
2606 compute_weights_fe_q_dofs_by_entity(const Number *data,
2607 const unsigned int n_components,
2608 const int n_points_1d_non_template,
2609 Number *weights)
2610 {
2611 const int n_points_1d = n_points_1d_template != -1 ?
2612 n_points_1d_template :
2613 n_points_1d_non_template;
2614
2615 Assert(n_points_1d > 0, ExcNotImplemented());
2616 Assert(n_points_1d < 100, ExcNotImplemented());
2617
2618 unsigned int compressed_index[100];
2619 compressed_index[0] = 0;
2620 for (int i = 1; i < n_points_1d - 1; ++i)
2621 compressed_index[i] = 1;
2622 compressed_index[n_points_1d - 1] = 2;
2623
2624 // Insert the number data into a storage position for weight,
2625 // ensuring that the array has either not been touched before
2626 // or the previous content is the same. In case the previous
2627 // content has a different value, we exit this function and
2628 // signal to outer functions that the compression was not possible.
2629 const auto check_and_set = [](Number &weight, const Number &data) {
2630 if (weight == Number(-1.0) || weight == data)
2631 {
2632 weight = data;
2633 return true; // success for the entry
2634 }
2635
2636 return false; // failure for the entry
2637 };
2638
2639 for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
2640 weights[c] = Number(-1.0);
2641
2642 for (unsigned int c = 0; c < n_components; ++c)
2643 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2644 for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
2645 ++j, data += n_points_1d)
2646 {
2647 const unsigned int shift =
2648 9 * compressed_index[k] + 3 * compressed_index[j];
2649
2650 if (!check_and_set(weights[shift], data[0]))
2651 return false; // failure
2652
2653 for (int i = 1; i < n_points_1d - 1; ++i)
2654 if (!check_and_set(weights[shift + 1], data[i]))
2655 return false; // failure
2656
2657 if (!check_and_set(weights[shift + 2], data[n_points_1d - 1]))
2658 return false; // failure
2659 }
2660
2661 return true; // success
2662 }
2663
2664
2665 template <int dim, int n_points_1d_template, typename Number>
2666 inline bool
2667 compute_weights_fe_q_dofs_by_entity_shifted(
2668 const Number *data,
2669 const unsigned int n_components,
2670 const int n_points_1d_non_template,
2671 Number *weights)
2672 {
2673 const int n_points_1d = n_points_1d_template != -1 ?
2674 n_points_1d_template :
2675 n_points_1d_non_template;
2676
2677 Assert((n_points_1d % 2) == 1,
2678 ExcMessage("The function can only with add number of points"));
2679 Assert(n_points_1d > 0, ExcNotImplemented());
2680 Assert(n_points_1d < 100, ExcNotImplemented());
2681
2682 const unsigned int n_inside_1d = n_points_1d / 2;
2683
2684 unsigned int compressed_index[100];
2685
2686 unsigned int c = 0;
2687 for (int i = 0; i < n_inside_1d; ++i)
2688 compressed_index[c++] = 0;
2689 compressed_index[c++] = 1;
2690 for (int i = 0; i < n_inside_1d; ++i)
2691 compressed_index[c++] = 2;
2692
2693 // Insert the number data into a storage position for weight,
2694 // ensuring that the array has either not been touched before
2695 // or the previous content is the same. In case the previous
2696 // content has a different value, we exit this function and
2697 // signal to outer functions that the compression was not possible.
2698 const auto check_and_set = [](Number &weight, const Number &data) {
2699 if (weight == Number(-1.0) || weight == data)
2700 {
2701 weight = data;
2702 return true; // success for the entry
2703 }
2704
2705 return false; // failure for the entry
2706 };
2707
2708 for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
2709 weights[c] = Number(-1.0);
2710
2711 for (unsigned int comp = 0; comp < n_components; ++comp)
2712 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
2713 for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
2714 ++j, data += n_points_1d)
2715 {
2716 const unsigned int shift =
2717 9 * compressed_index[k] + 3 * compressed_index[j];
2718
2719 unsigned int c = 0;
2720
2721 for (int i = 0; i < n_inside_1d; ++i)
2722 if (!check_and_set(weights[shift], data[c++]))
2723 return false; // failure
2724
2725 if (!check_and_set(weights[shift + 1], data[c++]))
2726 return false; // failure
2727
2728 for (int i = 0; i < n_inside_1d; ++i)
2729 if (!check_and_set(weights[shift + 2], data[c++]))
2730 return false; // failure
2731 }
2732
2733 return true; // success
2734 }
2735
2736
2737} // end of namespace internal
2738
2739
2740DEAL_II_NAMESPACE_CLOSE
2741
2742#endif
begin
*  *  iterator begin()
Definition ALE_Finite_Strain_Plasticity.h:10397
value_type
Number value_type
Definition ALE_Finite_Strain_Plasticity.h:10345
aligned_vector.h
DEAL_II_ALWAYS_INLINE
#define DEAL_II_ALWAYS_INLINE
Definition config.h:166
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
DEAL_II_RESTRICT
#define DEAL_II_RESTRICT
Definition config.h:167
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exception_macros.h:829
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exception_macros.h:901
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)
internal::MatrixFreeFunctions::ElementType
ElementType
Definition shape_info.h:56
internal::MatrixFreeFunctions::tensor_raviart_thomas
@ tensor_raviart_thomas
Definition shape_info.h:110
data
std::vector< index_type > data
Definition mpi.cc:734
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:2451
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:2521
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:116
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:2559
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:2383
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:2606
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:2667
numbers::invalid_unsigned_int
constexpr unsigned int invalid_unsigned_int
Definition types.h:228
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6596
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6579
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6563
shape_info.h
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:2231
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:2156
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:1928
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:1970
internal::EvaluatorTensorProduct::values
void values(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1670
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:1630
internal::EvaluatorTensorProduct::values_one_line
void values_one_line(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1720
internal::EvaluatorTensorProduct::gradients
void gradients(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1684
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:1590
internal::EvaluatorTensorProduct::n_rows_of_product
static constexpr unsigned int n_rows_of_product
Definition tensor_product_kernels.h:1572
internal::EvaluatorTensorProduct::hessians
void hessians(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1700
internal::EvaluatorTensorProduct::apply
static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out)
Definition tensor_product_kernels.h:1836
internal::EvaluatorTensorProduct::n_columns_of_product
static constexpr unsigned int n_columns_of_product
Definition tensor_product_kernels.h:1574
internal::EvaluatorTensorProduct::gradients_one_line
void gradients_one_line(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1736
internal::EvaluatorTensorProduct::hessians_one_line
void hessians_one_line(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:1755
Generated by doxygen 1.9.8