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