Reference documentation for deal.II version 9.4.1
\(\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// Copyright (C) 2017 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_matrix_free_tensor_product_kernels_h
18#define dealii_matrix_free_tensor_product_kernels_h
19
20#include <deal.II/base/config.h>
21
22#include <deal.II/base/aligned_vector.h>
23#include <deal.II/base/ndarray.h>
24#include <deal.II/base/polynomial.h>
25#include <deal.II/base/utilities.h>
26
27
28DEAL_II_NAMESPACE_OPEN
29
30
31
32namespace internal
33{
38 enum EvaluatorVariant
39 {
44 evaluate_general,
50 evaluate_symmetric,
56 evaluate_evenodd,
67 evaluate_symmetric_hierarchical,
71 evaluate_raviart_thomas
72 };
73
74
75
79 enum class EvaluatorQuantity
80 {
84 value,
88 gradient,
92 hessian
93 };
94
95
96
117 template <EvaluatorVariant variant,
118 int dim,
119 int n_rows,
120 int n_columns,
121 typename Number,
122 typename Number2 = Number>
123 struct EvaluatorTensorProduct
124 {};
125
148 template <EvaluatorVariant variant,
149 int dim,
150 int n_rows,
151 int n_columns,
152 typename Number,
153 int normal_dir,
154 typename Number2 = Number>
155 struct EvaluatorTensorProductAnisotropic
156 {};
157
158
159
177 template <int dim,
178 int n_rows,
179 int n_columns,
180 typename Number,
181 typename Number2>
182 struct EvaluatorTensorProduct<evaluate_general,
183 dim,
184 n_rows,
185 n_columns,
186 Number,
187 Number2>
188 {
189 static constexpr unsigned int n_rows_of_product =
190 Utilities::pow(n_rows, dim);
191 static constexpr unsigned int n_columns_of_product =
192 Utilities::pow(n_columns, dim);
193
198 EvaluatorTensorProduct()
199 : shape_values(nullptr)
200 , shape_gradients(nullptr)
201 , shape_hessians(nullptr)
202 {}
203
207 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
208 const AlignedVector<Number2> &shape_gradients,
209 const AlignedVector<Number2> &shape_hessians,
210 const unsigned int dummy1 = 0,
211 const unsigned int dummy2 = 0)
212 : shape_values(shape_values.begin())
213 , shape_gradients(shape_gradients.begin())
214 , shape_hessians(shape_hessians.begin())
215 {
216 // We can enter this function either for the apply() path that has
217 // n_rows * n_columns entries or for the apply_face() path that only has
218 // n_rows * 3 entries in the array. Since we cannot decide about the use
219 // we must allow for both here.
220 Assert(shape_values.size() == 0 ||
221 shape_values.size() == n_rows * n_columns ||
222 shape_values.size() == 3 * n_rows,
223 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
224 Assert(shape_gradients.size() == 0 ||
225 shape_gradients.size() == n_rows * n_columns,
226 ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
227 Assert(shape_hessians.size() == 0 ||
228 shape_hessians.size() == n_rows * n_columns,
229 ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
230 (void)dummy1;
231 (void)dummy2;
232 }
233
234 template <int direction, bool contract_over_rows, bool add>
235 void
236 values(const Number in[], Number out[]) const
237 {
238 apply<direction, contract_over_rows, add>(shape_values, in, out);
239 }
240
241 template <int direction, bool contract_over_rows, bool add>
242 void
243 gradients(const Number in[], Number out[]) const
244 {
245 apply<direction, contract_over_rows, add>(shape_gradients, in, out);
246 }
247
248 template <int direction, bool contract_over_rows, bool add>
249 void
250 hessians(const Number in[], Number out[]) const
251 {
252 apply<direction, contract_over_rows, add>(shape_hessians, in, out);
253 }
254
255 template <int direction, bool contract_over_rows, bool add>
256 void
257 values_one_line(const Number in[], Number out[]) const
258 {
259 Assert(shape_values != nullptr, ExcNotInitialized());
260 apply<direction, contract_over_rows, add, true>(shape_values, in, out);
261 }
262
263 template <int direction, bool contract_over_rows, bool add>
264 void
265 gradients_one_line(const Number in[], Number out[]) const
266 {
267 Assert(shape_gradients != nullptr, ExcNotInitialized());
268 apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
269 }
270
271 template <int direction, bool contract_over_rows, bool add>
272 void
273 hessians_one_line(const Number in[], Number out[]) const
274 {
275 Assert(shape_hessians != nullptr, ExcNotInitialized());
276 apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
277 }
278
303 template <int direction,
304 bool contract_over_rows,
305 bool add,
306 bool one_line = false>
307 static void
308 apply(const Number2 *DEAL_II_RESTRICT shape_data,
309 const Number * in,
310 Number * out);
311
341 template <int face_direction,
342 bool contract_onto_face,
343 bool add,
344 int max_derivative>
345 void
346 apply_face(const Number *DEAL_II_RESTRICT in,
347 Number *DEAL_II_RESTRICT out) const;
348
349 private:
350 const Number2 *shape_values;
351 const Number2 *shape_gradients;
352 const Number2 *shape_hessians;
353 };
354
355
356
357 template <int dim,
358 int n_rows,
359 int n_columns,
360 typename Number,
361 typename Number2>
362 template <int direction, bool contract_over_rows, bool add, bool one_line>
363 inline void
364 EvaluatorTensorProduct<evaluate_general,
365 dim,
366 n_rows,
367 n_columns,
368 Number,
369 Number2>::apply(const Number2 *DEAL_II_RESTRICT
370 shape_data,
371 const Number *in,
372 Number * out)
373 {
374 static_assert(one_line == false || direction == dim - 1,
375 "Single-line evaluation only works for direction=dim-1.");
376 Assert(shape_data != nullptr,
377 ExcMessage(
378 "The given array shape_data must not be the null pointer!"));
379 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
380 in != out,
381 ExcMessage("In-place operation only supported for "
382 "n_rows==n_columns or single-line interpolation"));
383 AssertIndexRange(direction, dim);
384 constexpr int mm = contract_over_rows ? n_rows : n_columns,
385 nn = contract_over_rows ? n_columns : n_rows;
386
387 constexpr int stride = Utilities::pow(n_columns, direction);
388 constexpr int n_blocks1 = one_line ? 1 : stride;
389 constexpr int n_blocks2 =
390 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
391
392 for (int i2 = 0; i2 < n_blocks2; ++i2)
393 {
394 for (int i1 = 0; i1 < n_blocks1; ++i1)
395 {
396 Number x[mm];
397 for (int i = 0; i < mm; ++i)
398 x[i] = in[stride * i];
399 for (int col = 0; col < nn; ++col)
400 {
401 Number2 val0;
402 if (contract_over_rows == true)
403 val0 = shape_data[col];
404 else
405 val0 = shape_data[col * n_columns];
406 Number res0 = val0 * x[0];
407 for (int i = 1; i < mm; ++i)
408 {
409 if (contract_over_rows == true)
410 val0 = shape_data[i * n_columns + col];
411 else
412 val0 = shape_data[col * n_columns + i];
413 res0 += val0 * x[i];
414 }
415 if (add)
416 out[stride * col] += res0;
417 else
418 out[stride * col] = res0;
419 }
420
421 if (one_line == false)
422 {
423 ++in;
424 ++out;
425 }
426 }
427 if (one_line == false)
428 {
429 in += stride * (mm - 1);
430 out += stride * (nn - 1);
431 }
432 }
433 }
434
435
436
437 template <int dim,
438 int n_rows,
439 int n_columns,
440 typename Number,
441 typename Number2>
442 template <int face_direction,
443 bool contract_onto_face,
444 bool add,
445 int max_derivative>
446 inline void
447 EvaluatorTensorProduct<evaluate_general,
448 dim,
449 n_rows,
450 n_columns,
451 Number,
452 Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
453 Number *DEAL_II_RESTRICT
454 out) const
455 {
456 Assert(dim > 0, ExcMessage("Only dim=1,2,3 supported"));
457 static_assert(max_derivative >= 0 && max_derivative < 3,
458 "Only derivative orders 0-2 implemented");
459 Assert(shape_values != nullptr,
460 ExcMessage(
461 "The given array shape_values must not be the null pointer."));
462
463 constexpr int n_blocks1 = (dim > 1 ? n_rows : 1);
464 constexpr int n_blocks2 = (dim > 2 ? n_rows : 1);
465
466 AssertIndexRange(face_direction, dim);
467 constexpr int in_stride = Utilities::pow(n_rows, face_direction);
468 constexpr int out_stride = Utilities::pow(n_rows, dim - 1);
469 const Number *DEAL_II_RESTRICT shape_values = this->shape_values;
470
471 for (int i2 = 0; i2 < n_blocks2; ++i2)
472 {
473 for (int i1 = 0; i1 < n_blocks1; ++i1)
474 {
475 if (contract_onto_face == true)
476 {
477 Number res0 = shape_values[0] * in[0];
478 Number res1, res2;
479 if (max_derivative > 0)
480 res1 = shape_values[n_rows] * in[0];
481 if (max_derivative > 1)
482 res2 = shape_values[2 * n_rows] * in[0];
483 for (int ind = 1; ind < n_rows; ++ind)
484 {
485 res0 += shape_values[ind] * in[in_stride * ind];
486 if (max_derivative > 0)
487 res1 += shape_values[ind + n_rows] * in[in_stride * ind];
488 if (max_derivative > 1)
489 res2 +=
490 shape_values[ind + 2 * n_rows] * in[in_stride * ind];
491 }
492 if (add)
493 {
494 out[0] += res0;
495 if (max_derivative > 0)
496 out[out_stride] += res1;
497 if (max_derivative > 1)
498 out[2 * out_stride] += res2;
499 }
500 else
501 {
502 out[0] = res0;
503 if (max_derivative > 0)
504 out[out_stride] = res1;
505 if (max_derivative > 1)
506 out[2 * out_stride] = res2;
507 }
508 }
509 else
510 {
511 for (int col = 0; col < n_rows; ++col)
512 {
513 if (add)
514 out[col * in_stride] += shape_values[col] * in[0];
515 else
516 out[col * in_stride] = shape_values[col] * in[0];
517 if (max_derivative > 0)
518 out[col * in_stride] +=
519 shape_values[col + n_rows] * in[out_stride];
520 if (max_derivative > 1)
521 out[col * in_stride] +=
522 shape_values[col + 2 * n_rows] * in[2 * out_stride];
523 }
524 }
525
526 // increment: in regular case, just go to the next point in
527 // x-direction. If we are at the end of one chunk in x-dir, need
528 // to jump over to the next layer in z-direction
529 switch (face_direction)
530 {
531 case 0:
532 in += contract_onto_face ? n_rows : 1;
533 out += contract_onto_face ? 1 : n_rows;
534 break;
535 case 1:
536 ++in;
537 ++out;
538 // faces 2 and 3 in 3D use local coordinate system zx, which
539 // is the other way around compared to the tensor
540 // product. Need to take that into account.
541 if (dim == 3)
542 {
543 if (contract_onto_face)
544 out += n_rows - 1;
545 else
546 in += n_rows - 1;
547 }
548 break;
549 case 2:
550 ++in;
551 ++out;
552 break;
553 default:
554 Assert(false, ExcNotImplemented());
555 }
556 }
557
558 // adjust for local coordinate system zx
559 if (face_direction == 1 && dim == 3)
560 {
561 if (contract_onto_face)
562 {
563 in += n_rows * (n_rows - 1);
564 out -= n_rows * n_rows - 1;
565 }
566 else
567 {
568 out += n_rows * (n_rows - 1);
569 in -= n_rows * n_rows - 1;
570 }
571 }
572 }
573 }
574
575
576
590 template <int dim, typename Number, typename Number2>
591 struct EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>
592 {
593 static constexpr unsigned int n_rows_of_product =
594 numbers::invalid_unsigned_int;
595 static constexpr unsigned int n_columns_of_product =
596 numbers::invalid_unsigned_int;
597
602 EvaluatorTensorProduct()
603 : shape_values(nullptr)
604 , shape_gradients(nullptr)
605 , shape_hessians(nullptr)
606 , n_rows(numbers::invalid_unsigned_int)
607 , n_columns(numbers::invalid_unsigned_int)
608 {}
609
613 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
614 const AlignedVector<Number2> &shape_gradients,
615 const AlignedVector<Number2> &shape_hessians,
616 const unsigned int n_rows,
617 const unsigned int n_columns)
618 : shape_values(shape_values.begin())
619 , shape_gradients(shape_gradients.begin())
620 , shape_hessians(shape_hessians.begin())
621 , n_rows(n_rows)
622 , n_columns(n_columns)
623 {
624 // We can enter this function either for the apply() path that has
625 // n_rows * n_columns entries or for the apply_face() path that only has
626 // n_rows * 3 entries in the array. Since we cannot decide about the use
627 // we must allow for both here.
628 Assert(shape_values.size() == 0 ||
629 shape_values.size() == n_rows * n_columns ||
630 shape_values.size() == n_rows * 3,
631 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
632 Assert(shape_gradients.size() == 0 ||
633 shape_gradients.size() == n_rows * n_columns,
634 ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
635 Assert(shape_hessians.size() == 0 ||
636 shape_hessians.size() == n_rows * n_columns,
637 ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
638 }
639
643 EvaluatorTensorProduct(const Number2 * shape_values,
644 const Number2 * shape_gradients,
645 const Number2 * shape_hessians,
646 const unsigned int n_rows,
647 const unsigned int n_columns)
648 : shape_values(shape_values)
649 , shape_gradients(shape_gradients)
650 , shape_hessians(shape_hessians)
651 , n_rows(n_rows)
652 , n_columns(n_columns)
653 {}
654
655 template <int direction, bool contract_over_rows, bool add>
656 void
657 values(const Number *in, Number *out) const
658 {
659 apply<direction, contract_over_rows, add>(shape_values, in, out);
660 }
661
662 template <int direction, bool contract_over_rows, bool add>
663 void
664 gradients(const Number *in, Number *out) const
665 {
666 apply<direction, contract_over_rows, add>(shape_gradients, in, out);
667 }
668
669 template <int direction, bool contract_over_rows, bool add>
670 void
671 hessians(const Number *in, Number *out) const
672 {
673 apply<direction, contract_over_rows, add>(shape_hessians, in, out);
674 }
675
676 template <int direction, bool contract_over_rows, bool add>
677 void
678 values_one_line(const Number in[], Number out[]) const
679 {
680 Assert(shape_values != nullptr, ExcNotInitialized());
681 apply<direction, contract_over_rows, add, true>(shape_values, in, out);
682 }
683
684 template <int direction, bool contract_over_rows, bool add>
685 void
686 gradients_one_line(const Number in[], Number out[]) const
687 {
688 Assert(shape_gradients != nullptr, ExcNotInitialized());
689 apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
690 }
691
692 template <int direction, bool contract_over_rows, bool add>
693 void
694 hessians_one_line(const Number in[], Number out[]) const
695 {
696 Assert(shape_hessians != nullptr, ExcNotInitialized());
697 apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
698 }
699
700 template <int direction,
701 bool contract_over_rows,
702 bool add,
703 bool one_line = false>
704 void
705 apply(const Number2 *DEAL_II_RESTRICT shape_data,
706 const Number * in,
707 Number * out) const;
708
709 template <int face_direction,
710 bool contract_onto_face,
711 bool add,
712 int max_derivative>
713 void
714 apply_face(const Number *DEAL_II_RESTRICT in,
715 Number *DEAL_II_RESTRICT out) const;
716
717 const Number2 * shape_values;
718 const Number2 * shape_gradients;
719 const Number2 * shape_hessians;
720 const unsigned int n_rows;
721 const unsigned int n_columns;
722 };
723
724
725
726 template <int dim, typename Number, typename Number2>
727 template <int direction, bool contract_over_rows, bool add, bool one_line>
728 inline void
729 EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>::apply(
730 const Number2 *DEAL_II_RESTRICT shape_data,
731 const Number * in,
732 Number * out) const
733 {
734 static_assert(one_line == false || direction == dim - 1,
735 "Single-line evaluation only works for direction=dim-1.");
736 Assert(shape_data != nullptr,
737 ExcMessage(
738 "The given array shape_data must not be the null pointer!"));
739 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
740 in != out,
741 ExcMessage("In-place operation only supported for "
742 "n_rows==n_columns or single-line interpolation"));
743 AssertIndexRange(direction, dim);
744 const int mm = contract_over_rows ? n_rows : n_columns,
745 nn = contract_over_rows ? n_columns : n_rows;
746
747 const int stride =
748 direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
749 const int n_blocks1 = one_line ? 1 : stride;
750 const int n_blocks2 = direction >= dim - 1 ?
751 1 :
752 Utilities::fixed_power<dim - direction - 1>(n_rows);
753 Assert(n_rows <= 128, ExcNotImplemented());
754
755 // specialization for n_rows = 2 that manually unrolls the innermost loop
756 // to make the operation perform better (not completely as good as the
757 // templated one, but much better than the generic version down below,
758 // because the loop over col can be more effectively unrolled by the
759 // compiler)
760 if (contract_over_rows && n_rows == 2)
761 {
762 const Number2 *shape_data_1 = shape_data + n_columns;
763 for (int i2 = 0; i2 < n_blocks2; ++i2)
764 {
765 for (int i1 = 0; i1 < n_blocks1; ++i1)
766 {
767 const Number x0 = in[0], x1 = in[stride];
768 for (int col = 0; col < nn; ++col)
769 {
770 const Number result =
771 shape_data[col] * x0 + shape_data_1[col] * x1;
772 if (add)
773 out[stride * col] += result;
774 else
775 out[stride * col] = result;
776 }
777
778 if (one_line == false)
779 {
780 ++in;
781 ++out;
782 }
783 }
784 if (one_line == false)
785 {
786 in += stride * (mm - 1);
787 out += stride * (nn - 1);
788 }
789 }
790 }
791 // specialization for n = 3
792 else if (contract_over_rows && n_rows == 3)
793 {
794 const Number2 *shape_data_1 = shape_data + n_columns;
795 const Number2 *shape_data_2 = shape_data + 2 * n_columns;
796 for (int i2 = 0; i2 < n_blocks2; ++i2)
797 {
798 for (int i1 = 0; i1 < n_blocks1; ++i1)
799 {
800 const Number x0 = in[0], x1 = in[stride], x2 = in[2 * stride];
801 for (int col = 0; col < nn; ++col)
802 {
803 const Number result = shape_data[col] * x0 +
804 shape_data_1[col] * x1 +
805 shape_data_2[col] * x2;
806 if (add)
807 out[stride * col] += result;
808 else
809 out[stride * col] = result;
810 }
811
812 if (one_line == false)
813 {
814 ++in;
815 ++out;
816 }
817 }
818 if (one_line == false)
819 {
820 in += stride * (mm - 1);
821 out += stride * (nn - 1);
822 }
823 }
824 }
825 // general loop for all other cases
826 else
827 for (int i2 = 0; i2 < n_blocks2; ++i2)
828 {
829 for (int i1 = 0; i1 < n_blocks1; ++i1)
830 {
831 Number x[129];
832 for (int i = 0; i < mm; ++i)
833 x[i] = in[stride * i];
834 for (int col = 0; col < nn; ++col)
835 {
836 Number2 val0;
837 if (contract_over_rows == true)
838 val0 = shape_data[col];
839 else
840 val0 = shape_data[col * n_columns];
841 Number res0 = val0 * x[0];
842 for (int i = 1; i < mm; ++i)
843 {
844 if (contract_over_rows == true)
845 val0 = shape_data[i * n_columns + col];
846 else
847 val0 = shape_data[col * n_columns + i];
848 res0 += val0 * x[i];
849 }
850 if (add)
851 out[stride * col] += res0;
852 else
853 out[stride * col] = res0;
854 }
855
856 if (one_line == false)
857 {
858 ++in;
859 ++out;
860 }
861 }
862 if (one_line == false)
863 {
864 in += stride * (mm - 1);
865 out += stride * (nn - 1);
866 }
867 }
868 }
869
870
871
872 template <int dim, typename Number, typename Number2>
873 template <int face_direction,
874 bool contract_onto_face,
875 bool add,
876 int max_derivative>
877 inline void
878 EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>::
879 apply_face(const Number *DEAL_II_RESTRICT in,
880 Number *DEAL_II_RESTRICT out) const
881 {
882 Assert(shape_values != nullptr,
883 ExcMessage(
884 "The given array shape_data must not be the null pointer!"));
885 static_assert(dim > 0 && dim < 4, "Only dim=1,2,3 supported");
886 const int n_blocks1 = dim > 1 ? n_rows : 1;
887 const int n_blocks2 = dim > 2 ? n_rows : 1;
888
889 AssertIndexRange(face_direction, dim);
890 const int in_stride =
891 face_direction > 0 ? Utilities::fixed_power<face_direction>(n_rows) : 1;
892 const int out_stride =
893 dim > 1 ? Utilities::fixed_power<dim - 1>(n_rows) : 1;
894
895 for (int i2 = 0; i2 < n_blocks2; ++i2)
896 {
897 for (int i1 = 0; i1 < n_blocks1; ++i1)
898 {
899 if (contract_onto_face == true)
900 {
901 Number res0 = shape_values[0] * in[0];
902 Number res1, res2;
903 if (max_derivative > 0)
904 res1 = shape_values[n_rows] * in[0];
905 if (max_derivative > 1)
906 res2 = shape_values[2 * n_rows] * in[0];
907 for (unsigned int ind = 1; ind < n_rows; ++ind)
908 {
909 res0 += shape_values[ind] * in[in_stride * ind];
910 if (max_derivative > 0)
911 res1 += shape_values[ind + n_rows] * in[in_stride * ind];
912 if (max_derivative > 1)
913 res2 +=
914 shape_values[ind + 2 * n_rows] * in[in_stride * ind];
915 }
916 if (add)
917 {
918 out[0] += res0;
919 if (max_derivative > 0)
920 out[out_stride] += res1;
921 if (max_derivative > 1)
922 out[2 * out_stride] += res2;
923 }
924 else
925 {
926 out[0] = res0;
927 if (max_derivative > 0)
928 out[out_stride] = res1;
929 if (max_derivative > 1)
930 out[2 * out_stride] = res2;
931 }
932 }
933 else
934 {
935 for (unsigned int col = 0; col < n_rows; ++col)
936 {
937 if (add)
938 out[col * in_stride] += shape_values[col] * in[0];
939 else
940 out[col * in_stride] = shape_values[col] * in[0];
941 if (max_derivative > 0)
942 out[col * in_stride] +=
943 shape_values[col + n_rows] * in[out_stride];
944 if (max_derivative > 1)
945 out[col * in_stride] +=
946 shape_values[col + 2 * n_rows] * in[2 * out_stride];
947 }
948 }
949
950 // increment: in regular case, just go to the next point in
951 // x-direction. If we are at the end of one chunk in x-dir, need
952 // to jump over to the next layer in z-direction
953 switch (face_direction)
954 {
955 case 0:
956 in += contract_onto_face ? n_rows : 1;
957 out += contract_onto_face ? 1 : n_rows;
958 break;
959 case 1:
960 ++in;
961 ++out;
962 // faces 2 and 3 in 3D use local coordinate system zx, which
963 // is the other way around compared to the tensor
964 // product. Need to take that into account.
965 if (dim == 3)
966 {
967 if (contract_onto_face)
968 out += n_rows - 1;
969 else
970 in += n_rows - 1;
971 }
972 break;
973 case 2:
974 ++in;
975 ++out;
976 break;
977 default:
978 Assert(false, ExcNotImplemented());
979 }
980 }
981 if (face_direction == 1 && dim == 3)
982 {
983 // adjust for local coordinate system zx
984 if (contract_onto_face)
985 {
986 in += n_rows * (n_rows - 1);
987 out -= n_rows * n_rows - 1;
988 }
989 else
990 {
991 out += n_rows * (n_rows - 1);
992 in -= n_rows * n_rows - 1;
993 }
994 }
995 }
996 }
997
998
999
1020 template <int dim,
1021 int n_rows,
1022 int n_columns,
1023 typename Number,
1024 typename Number2>
1025 struct EvaluatorTensorProduct<evaluate_symmetric,
1026 dim,
1027 n_rows,
1028 n_columns,
1029 Number,
1030 Number2>
1031 {
1032 static constexpr unsigned int n_rows_of_product =
1033 Utilities::pow(n_rows, dim);
1034 static constexpr unsigned int n_columns_of_product =
1035 Utilities::pow(n_columns, dim);
1036
1040 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
1041 const AlignedVector<Number2> &shape_gradients,
1042 const AlignedVector<Number2> &shape_hessians,
1043 const unsigned int dummy1 = 0,
1044 const unsigned int dummy2 = 0)
1045 : shape_values(shape_values.begin())
1046 , shape_gradients(shape_gradients.begin())
1047 , shape_hessians(shape_hessians.begin())
1048 {
1049 Assert(shape_values.size() == 0 ||
1050 shape_values.size() == n_rows * n_columns,
1051 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1052 Assert(shape_gradients.size() == 0 ||
1053 shape_gradients.size() == n_rows * n_columns,
1054 ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
1055 Assert(shape_hessians.size() == 0 ||
1056 shape_hessians.size() == n_rows * n_columns,
1057 ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
1058 (void)dummy1;
1059 (void)dummy2;
1060 }
1061
1062 template <int direction, bool contract_over_rows, bool add>
1063 void
1064 values(const Number in[], Number out[]) const;
1065
1066 template <int direction, bool contract_over_rows, bool add>
1067 void
1068 gradients(const Number in[], Number out[]) const;
1069
1070 template <int direction, bool contract_over_rows, bool add>
1071 void
1072 hessians(const Number in[], Number out[]) const;
1073
1074 private:
1075 const Number2 *shape_values;
1076 const Number2 *shape_gradients;
1077 const Number2 *shape_hessians;
1078 };
1079
1080
1081
1082 // In this case, the 1D shape values read (sorted lexicographically, rows
1083 // run over 1D dofs, columns over quadrature points):
1084 // Q2 --> [ 0.687 0 -0.087 ]
1085 // [ 0.4 1 0.4 ]
1086 // [-0.087 0 0.687 ]
1087 // Q3 --> [ 0.66 0.003 0.002 0.049 ]
1088 // [ 0.521 1.005 -0.01 -0.230 ]
1089 // [-0.230 -0.01 1.005 0.521 ]
1090 // [ 0.049 0.002 0.003 0.66 ]
1091 // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ]
1092 // [ 0.608 1.059 0 0.039 0.176 ]
1093 // [-0.409 -0.113 1 -0.113 -0.409 ]
1094 // [ 0.176 0.039 0 1.059 0.608 ]
1095 // [-0.032 -0.007 0 0.022 0.658 ]
1096 //
1097 // In these matrices, we want to use avoid computations involving zeros and
1098 // ones and in addition use the symmetry in entries to reduce the number of
1099 // read operations.
1100 template <int dim,
1101 int n_rows,
1102 int n_columns,
1103 typename Number,
1104 typename Number2>
1105 template <int direction, bool contract_over_rows, bool add>
1106 inline void
1107 EvaluatorTensorProduct<evaluate_symmetric,
1108 dim,
1109 n_rows,
1110 n_columns,
1111 Number,
1112 Number2>::values(const Number in[], Number out[]) const
1113 {
1114 Assert(shape_values != nullptr, ExcNotInitialized());
1115 AssertIndexRange(direction, dim);
1116 constexpr int mm = contract_over_rows ? n_rows : n_columns,
1117 nn = contract_over_rows ? n_columns : n_rows;
1118 constexpr int n_cols = nn / 2;
1119 constexpr int mid = mm / 2;
1120
1121 constexpr int stride = Utilities::pow(n_columns, direction);
1122 constexpr int n_blocks1 = stride;
1123 constexpr int n_blocks2 =
1124 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1125
1126 for (int i2 = 0; i2 < n_blocks2; ++i2)
1127 {
1128 for (int i1 = 0; i1 < n_blocks1; ++i1)
1129 {
1130 for (int col = 0; col < n_cols; ++col)
1131 {
1132 Number2 val0, val1;
1133 Number in0, in1, res0, res1;
1134 if (contract_over_rows == true)
1135 {
1136 val0 = shape_values[col];
1137 val1 = shape_values[nn - 1 - col];
1138 }
1139 else
1140 {
1141 val0 = shape_values[col * n_columns];
1142 val1 = shape_values[(col + 1) * n_columns - 1];
1143 }
1144 if (mid > 0)
1145 {
1146 in0 = in[0];
1147 in1 = in[stride * (mm - 1)];
1148 res0 = val0 * in0;
1149 res1 = val1 * in0;
1150 res0 += val1 * in1;
1151 res1 += val0 * in1;
1152 for (int ind = 1; ind < mid; ++ind)
1153 {
1154 if (contract_over_rows == true)
1155 {
1156 val0 = shape_values[ind * n_columns + col];
1157 val1 = shape_values[ind * n_columns + nn - 1 - col];
1158 }
1159 else
1160 {
1161 val0 = shape_values[col * n_columns + ind];
1162 val1 =
1163 shape_values[(col + 1) * n_columns - 1 - ind];
1164 }
1165 in0 = in[stride * ind];
1166 in1 = in[stride * (mm - 1 - ind)];
1167 res0 += val0 * in0;
1168 res1 += val1 * in0;
1169 res0 += val1 * in1;
1170 res1 += val0 * in1;
1171 }
1172 }
1173 else
1174 res0 = res1 = Number();
1175 if (contract_over_rows == true)
1176 {
1177 if (mm % 2 == 1)
1178 {
1179 val0 = shape_values[mid * n_columns + col];
1180 in1 = val0 * in[stride * mid];
1181 res0 += in1;
1182 res1 += in1;
1183 }
1184 }
1185 else
1186 {
1187 if (mm % 2 == 1 && nn % 2 == 0)
1188 {
1189 val0 = shape_values[col * n_columns + mid];
1190 in1 = val0 * in[stride * mid];
1191 res0 += in1;
1192 res1 += in1;
1193 }
1194 }
1195 if (add)
1196 {
1197 out[stride * col] += res0;
1198 out[stride * (nn - 1 - col)] += res1;
1199 }
1200 else
1201 {
1202 out[stride * col] = res0;
1203 out[stride * (nn - 1 - col)] = res1;
1204 }
1205 }
1206 if (contract_over_rows == true && nn % 2 == 1 && mm % 2 == 1)
1207 {
1208 if (add)
1209 out[stride * n_cols] += in[stride * mid];
1210 else
1211 out[stride * n_cols] = in[stride * mid];
1212 }
1213 else if (contract_over_rows == true && nn % 2 == 1)
1214 {
1215 Number res0;
1216 Number2 val0 = shape_values[n_cols];
1217 if (mid > 0)
1218 {
1219 res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1220 for (int ind = 1; ind < mid; ++ind)
1221 {
1222 val0 = shape_values[ind * n_columns + n_cols];
1223 res0 += val0 * (in[stride * ind] +
1224 in[stride * (mm - 1 - ind)]);
1225 }
1226 }
1227 else
1228 res0 = Number();
1229 if (add)
1230 out[stride * n_cols] += res0;
1231 else
1232 out[stride * n_cols] = res0;
1233 }
1234 else if (contract_over_rows == false && nn % 2 == 1)
1235 {
1236 Number res0;
1237 if (mid > 0)
1238 {
1239 Number2 val0 = shape_values[n_cols * n_columns];
1240 res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1241 for (int ind = 1; ind < mid; ++ind)
1242 {
1243 val0 = shape_values[n_cols * n_columns + ind];
1244 Number in1 = val0 * (in[stride * ind] +
1245 in[stride * (mm - 1 - ind)]);
1246 res0 += in1;
1247 }
1248 if (mm % 2)
1249 res0 += in[stride * mid];
1250 }
1251 else
1252 res0 = in[0];
1253 if (add)
1254 out[stride * n_cols] += res0;
1255 else
1256 out[stride * n_cols] = res0;
1257 }
1258
1259 ++in;
1260 ++out;
1261 }
1262 in += stride * (mm - 1);
1263 out += stride * (nn - 1);
1264 }
1265 }
1266
1267
1268
1269 // For the specialized loop used for the gradient computation in
1270 // here, the 1D shape values read (sorted lexicographically, rows
1271 // run over 1D dofs, columns over quadrature points):
1272 // Q2 --> [-2.549 -1 0.549 ]
1273 // [ 3.098 0 -3.098 ]
1274 // [-0.549 1 2.549 ]
1275 // Q3 --> [-4.315 -1.03 0.5 -0.44 ]
1276 // [ 6.07 -1.44 -2.97 2.196 ]
1277 // [-2.196 2.97 1.44 -6.07 ]
1278 // [ 0.44 -0.5 1.03 4.315 ]
1279 // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ]
1280 // [10.111 -2.76 -2.667 2.066 -2.306 ]
1281 // [-5.688 5.773 0 -5.773 5.688 ]
1282 // [ 2.306 -2.066 2.667 2.76 -10.111 ]
1283 // [-0.413 0.353 -0.333 -0.353 0.413 ]
1284 //
1285 // In these matrices, we want to use avoid computations involving
1286 // zeros and ones and in addition use the symmetry in entries to
1287 // reduce the number of read operations.
1288 template <int dim,
1289 int n_rows,
1290 int n_columns,
1291 typename Number,
1292 typename Number2>
1293 template <int direction, bool contract_over_rows, bool add>
1294 inline void
1295 EvaluatorTensorProduct<evaluate_symmetric,
1296 dim,
1297 n_rows,
1298 n_columns,
1299 Number,
1300 Number2>::gradients(const Number in[],
1301 Number out[]) const
1302 {
1303 Assert(shape_gradients != nullptr, ExcNotInitialized());
1304 AssertIndexRange(direction, dim);
1305 constexpr int mm = contract_over_rows ? n_rows : n_columns,
1306 nn = contract_over_rows ? n_columns : n_rows;
1307 constexpr int n_cols = nn / 2;
1308 constexpr int mid = mm / 2;
1309
1310 constexpr int stride = Utilities::pow(n_columns, direction);
1311 constexpr int n_blocks1 = stride;
1312 constexpr int n_blocks2 =
1313 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1314
1315 for (int i2 = 0; i2 < n_blocks2; ++i2)
1316 {
1317 for (int i1 = 0; i1 < n_blocks1; ++i1)
1318 {
1319 for (int col = 0; col < n_cols; ++col)
1320 {
1321 Number2 val0, val1;
1322 Number in0, in1, res0, res1;
1323 if (contract_over_rows == true)
1324 {
1325 val0 = shape_gradients[col];
1326 val1 = shape_gradients[nn - 1 - col];
1327 }
1328 else
1329 {
1330 val0 = shape_gradients[col * n_columns];
1331 val1 = shape_gradients[(nn - col - 1) * n_columns];
1332 }
1333 if (mid > 0)
1334 {
1335 in0 = in[0];
1336 in1 = in[stride * (mm - 1)];
1337 res0 = val0 * in0;
1338 res1 = val1 * in0;
1339 res0 -= val1 * in1;
1340 res1 -= val0 * in1;
1341 for (int ind = 1; ind < mid; ++ind)
1342 {
1343 if (contract_over_rows == true)
1344 {
1345 val0 = shape_gradients[ind * n_columns + col];
1346 val1 =
1347 shape_gradients[ind * n_columns + nn - 1 - col];
1348 }
1349 else
1350 {
1351 val0 = shape_gradients[col * n_columns + ind];
1352 val1 =
1353 shape_gradients[(nn - col - 1) * n_columns + ind];
1354 }
1355 in0 = in[stride * ind];
1356 in1 = in[stride * (mm - 1 - ind)];
1357 res0 += val0 * in0;
1358 res1 += val1 * in0;
1359 res0 -= val1 * in1;
1360 res1 -= val0 * in1;
1361 }
1362 }
1363 else
1364 res0 = res1 = Number();
1365 if (mm % 2 == 1)
1366 {
1367 if (contract_over_rows == true)
1368 val0 = shape_gradients[mid * n_columns + col];
1369 else
1370 val0 = shape_gradients[col * n_columns + mid];
1371 in1 = val0 * in[stride * mid];
1372 res0 += in1;
1373 res1 -= in1;
1374 }
1375 if (add)
1376 {
1377 out[stride * col] += res0;
1378 out[stride * (nn - 1 - col)] += res1;
1379 }
1380 else
1381 {
1382 out[stride * col] = res0;
1383 out[stride * (nn - 1 - col)] = res1;
1384 }
1385 }
1386 if (nn % 2 == 1)
1387 {
1388 Number2 val0;
1389 Number res0;
1390 if (contract_over_rows == true)
1391 val0 = shape_gradients[n_cols];
1392 else
1393 val0 = shape_gradients[n_cols * n_columns];
1394 res0 = val0 * (in[0] - in[stride * (mm - 1)]);
1395 for (int ind = 1; ind < mid; ++ind)
1396 {
1397 if (contract_over_rows == true)
1398 val0 = shape_gradients[ind * n_columns + n_cols];
1399 else
1400 val0 = shape_gradients[n_cols * n_columns + ind];
1401 Number in1 =
1402 val0 * (in[stride * ind] - in[stride * (mm - 1 - ind)]);
1403 res0 += in1;
1404 }
1405 if (add)
1406 out[stride * n_cols] += res0;
1407 else
1408 out[stride * n_cols] = res0;
1409 }
1410
1411 ++in;
1412 ++out;
1413 }
1414 in += stride * (mm - 1);
1415 out += stride * (nn - 1);
1416 }
1417 }
1418
1419
1420
1421 // evaluates the given shape data in 1d-3d using the tensor product
1422 // form assuming the symmetries of unit cell shape hessians for
1423 // finite elements in FEEvaluation
1424 template <int dim,
1425 int n_rows,
1426 int n_columns,
1427 typename Number,
1428 typename Number2>
1429 template <int direction, bool contract_over_rows, bool add>
1430 inline void
1431 EvaluatorTensorProduct<evaluate_symmetric,
1432 dim,
1433 n_rows,
1434 n_columns,
1435 Number,
1436 Number2>::hessians(const Number in[],
1437 Number out[]) const
1438 {
1439 Assert(shape_hessians != nullptr, ExcNotInitialized());
1440 AssertIndexRange(direction, dim);
1441 constexpr int mm = contract_over_rows ? n_rows : n_columns;
1442 constexpr int nn = contract_over_rows ? n_columns : n_rows;
1443 constexpr int n_cols = nn / 2;
1444 constexpr int mid = mm / 2;
1445
1446 constexpr int stride = Utilities::pow(n_columns, direction);
1447 constexpr int n_blocks1 = stride;
1448 constexpr int n_blocks2 =
1449 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1450
1451 for (int i2 = 0; i2 < n_blocks2; ++i2)
1452 {
1453 for (int i1 = 0; i1 < n_blocks1; ++i1)
1454 {
1455 for (int col = 0; col < n_cols; ++col)
1456 {
1457 Number2 val0, val1;
1458 Number in0, in1, res0, res1;
1459 if (contract_over_rows == true)
1460 {
1461 val0 = shape_hessians[col];
1462 val1 = shape_hessians[nn - 1 - col];
1463 }
1464 else
1465 {
1466 val0 = shape_hessians[col * n_columns];
1467 val1 = shape_hessians[(col + 1) * n_columns - 1];
1468 }
1469 if (mid > 0)
1470 {
1471 in0 = in[0];
1472 in1 = in[stride * (mm - 1)];
1473 res0 = val0 * in0;
1474 res1 = val1 * in0;
1475 res0 += val1 * in1;
1476 res1 += val0 * in1;
1477 for (int ind = 1; ind < mid; ++ind)
1478 {
1479 if (contract_over_rows == true)
1480 {
1481 val0 = shape_hessians[ind * n_columns + col];
1482 val1 =
1483 shape_hessians[ind * n_columns + nn - 1 - col];
1484 }
1485 else
1486 {
1487 val0 = shape_hessians[col * n_columns + ind];
1488 val1 =
1489 shape_hessians[(col + 1) * n_columns - 1 - ind];
1490 }
1491 in0 = in[stride * ind];
1492 in1 = in[stride * (mm - 1 - ind)];
1493 res0 += val0 * in0;
1494 res1 += val1 * in0;
1495 res0 += val1 * in1;
1496 res1 += val0 * in1;
1497 }
1498 }
1499 else
1500 res0 = res1 = Number();
1501 if (mm % 2 == 1)
1502 {
1503 if (contract_over_rows == true)
1504 val0 = shape_hessians[mid * n_columns + col];
1505 else
1506 val0 = shape_hessians[col * n_columns + mid];
1507 in1 = val0 * in[stride * mid];
1508 res0 += in1;
1509 res1 += in1;
1510 }
1511 if (add)
1512 {
1513 out[stride * col] += res0;
1514 out[stride * (nn - 1 - col)] += res1;
1515 }
1516 else
1517 {
1518 out[stride * col] = res0;
1519 out[stride * (nn - 1 - col)] = res1;
1520 }
1521 }
1522 if (nn % 2 == 1)
1523 {
1524 Number2 val0;
1525 Number res0;
1526 if (contract_over_rows == true)
1527 val0 = shape_hessians[n_cols];
1528 else
1529 val0 = shape_hessians[n_cols * n_columns];
1530 if (mid > 0)
1531 {
1532 res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1533 for (int ind = 1; ind < mid; ++ind)
1534 {
1535 if (contract_over_rows == true)
1536 val0 = shape_hessians[ind * n_columns + n_cols];
1537 else
1538 val0 = shape_hessians[n_cols * n_columns + ind];
1539 Number in1 = val0 * (in[stride * ind] +
1540 in[stride * (mm - 1 - ind)]);
1541 res0 += in1;
1542 }
1543 }
1544 else
1545 res0 = Number();
1546 if (mm % 2 == 1)
1547 {
1548 if (contract_over_rows == true)
1549 val0 = shape_hessians[mid * n_columns + n_cols];
1550 else
1551 val0 = shape_hessians[n_cols * n_columns + mid];
1552 res0 += val0 * in[stride * mid];
1553 }
1554 if (add)
1555 out[stride * n_cols] += res0;
1556 else
1557 out[stride * n_cols] = res0;
1558 }
1559
1560 ++in;
1561 ++out;
1562 }
1563 in += stride * (mm - 1);
1564 out += stride * (nn - 1);
1565 }
1566 }
1567
1568
1569
1601 template <int dim,
1602 int n_rows,
1603 int n_columns,
1604 typename Number,
1605 typename Number2>
1606 struct EvaluatorTensorProduct<evaluate_evenodd,
1607 dim,
1608 n_rows,
1609 n_columns,
1610 Number,
1611 Number2>
1612 {
1613 static constexpr unsigned int n_rows_of_product =
1614 Utilities::pow(n_rows, dim);
1615 static constexpr unsigned int n_columns_of_product =
1616 Utilities::pow(n_columns, dim);
1617
1623 EvaluatorTensorProduct()
1624 : shape_values(nullptr)
1625 , shape_gradients(nullptr)
1626 , shape_hessians(nullptr)
1627 {}
1628
1633 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values)
1634 : shape_values(shape_values.begin())
1635 , shape_gradients(nullptr)
1636 , shape_hessians(nullptr)
1637 {
1638 AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
1639 }
1640
1645 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
1646 const AlignedVector<Number2> &shape_gradients,
1647 const AlignedVector<Number2> &shape_hessians,
1648 const unsigned int dummy1 = 0,
1649 const unsigned int dummy2 = 0)
1650 : shape_values(shape_values.begin())
1651 , shape_gradients(shape_gradients.begin())
1652 , shape_hessians(shape_hessians.begin())
1653 {
1654 // In this function, we allow for dummy pointers if some of values,
1655 // gradients or hessians should not be computed
1656 if (!shape_values.empty())
1657 AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
1658 if (!shape_gradients.empty())
1659 AssertDimension(shape_gradients.size(), n_rows * ((n_columns + 1) / 2));
1660 if (!shape_hessians.empty())
1661 AssertDimension(shape_hessians.size(), n_rows * ((n_columns + 1) / 2));
1662 (void)dummy1;
1663 (void)dummy2;
1664 }
1665
1666 template <int direction, bool contract_over_rows, bool add>
1667 void
1668 values(const Number in[], Number out[]) const
1669 {
1670 Assert(shape_values != nullptr, ExcNotInitialized());
1671 apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
1672 }
1673
1674 template <int direction, bool contract_over_rows, bool add>
1675 void
1676 gradients(const Number in[], Number out[]) const
1677 {
1678 Assert(shape_gradients != nullptr, ExcNotInitialized());
1679 apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
1680 }
1681
1682 template <int direction, bool contract_over_rows, bool add>
1683 void
1684 hessians(const Number in[], Number out[]) const
1685 {
1686 Assert(shape_hessians != nullptr, ExcNotInitialized());
1687 apply<direction, contract_over_rows, add, 2>(shape_hessians, in, out);
1688 }
1689
1690 template <int direction, bool contract_over_rows, bool add>
1691 void
1692 values_one_line(const Number in[], Number out[]) const
1693 {
1694 Assert(shape_values != nullptr, ExcNotInitialized());
1695 apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
1696 }
1697
1698 template <int direction, bool contract_over_rows, bool add>
1699 void
1700 gradients_one_line(const Number in[], Number out[]) const
1701 {
1702 Assert(shape_gradients != nullptr, ExcNotInitialized());
1703 apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
1704 in,
1705 out);
1706 }
1707
1708 template <int direction, bool contract_over_rows, bool add>
1709 void
1710 hessians_one_line(const Number in[], Number out[]) const
1711 {
1712 Assert(shape_hessians != nullptr, ExcNotInitialized());
1713 apply<direction, contract_over_rows, add, 2, true>(shape_hessians,
1714 in,
1715 out);
1716 }
1717
1746 template <int direction,
1747 bool contract_over_rows,
1748 bool add,
1749 int type,
1750 bool one_line = false>
1751 static void
1752 apply(const Number2 *DEAL_II_RESTRICT shape_data,
1753 const Number * in,
1754 Number * out);
1755
1756 private:
1757 const Number2 *shape_values;
1758 const Number2 *shape_gradients;
1759 const Number2 *shape_hessians;
1760 };
1761
1762
1763
1764 template <int dim,
1765 int n_rows,
1766 int n_columns,
1767 typename Number,
1768 typename Number2>
1769 template <int direction,
1770 bool contract_over_rows,
1771 bool add,
1772 int type,
1773 bool one_line>
1774 inline void
1775 EvaluatorTensorProduct<evaluate_evenodd,
1776 dim,
1777 n_rows,
1778 n_columns,
1779 Number,
1780 Number2>::apply(const Number2 *DEAL_II_RESTRICT shapes,
1781 const Number * in,
1782 Number * out)
1783 {
1784 static_assert(type < 3, "Only three variants type=0,1,2 implemented");
1785 static_assert(one_line == false || direction == dim - 1,
1786 "Single-line evaluation only works for direction=dim-1.");
1787 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
1788 in != out,
1789 ExcMessage("In-place operation only supported for "
1790 "n_rows==n_columns or single-line interpolation"));
1791
1792 // We cannot statically assert that direction is less than dim, so must do
1793 // an additional dynamic check
1794 AssertIndexRange(direction, dim);
1795
1796 constexpr int nn = contract_over_rows ? n_columns : n_rows;
1797 constexpr int mm = contract_over_rows ? n_rows : n_columns;
1798 constexpr int n_cols = nn / 2;
1799 constexpr int mid = mm / 2;
1800
1801 constexpr int stride = Utilities::pow(n_columns, direction);
1802 constexpr int n_blocks1 = one_line ? 1 : stride;
1803 constexpr int n_blocks2 =
1804 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1805
1806 constexpr int offset = (n_columns + 1) / 2;
1807
1808 // this code may look very inefficient at first sight due to the many
1809 // different cases with if's at the innermost loop part, but all of the
1810 // conditionals can be evaluated at compile time because they are
1811 // templates, so the compiler should optimize everything away
1812 for (int i2 = 0; i2 < n_blocks2; ++i2)
1813 {
1814 for (int i1 = 0; i1 < n_blocks1; ++i1)
1815 {
1816 Number xp[mid > 0 ? mid : 1], xm[mid > 0 ? mid : 1];
1817 for (int i = 0; i < mid; ++i)
1818 {
1819 if (contract_over_rows == true && type == 1)
1820 {
1821 xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
1822 xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
1823 }
1824 else
1825 {
1826 xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
1827 xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
1828 }
1829 }
1830 Number xmid = in[stride * mid];
1831 for (int col = 0; col < n_cols; ++col)
1832 {
1833 Number r0, r1;
1834 if (mid > 0)
1835 {
1836 if (contract_over_rows == true)
1837 {
1838 r0 = shapes[col] * xp[0];
1839 r1 = shapes[(n_rows - 1) * offset + col] * xm[0];
1840 }
1841 else
1842 {
1843 r0 = shapes[col * offset] * xp[0];
1844 r1 = shapes[(n_rows - 1 - col) * offset] * xm[0];
1845 }
1846 for (int ind = 1; ind < mid; ++ind)
1847 {
1848 if (contract_over_rows == true)
1849 {
1850 r0 += shapes[ind * offset + col] * xp[ind];
1851 r1 += shapes[(n_rows - 1 - ind) * offset + col] *
1852 xm[ind];
1853 }
1854 else
1855 {
1856 r0 += shapes[col * offset + ind] * xp[ind];
1857 r1 += shapes[(n_rows - 1 - col) * offset + ind] *
1858 xm[ind];
1859 }
1860 }
1861 }
1862 else
1863 r0 = r1 = Number();
1864 if (mm % 2 == 1 && contract_over_rows == true)
1865 {
1866 if (type == 1)
1867 r1 += shapes[mid * offset + col] * xmid;
1868 else
1869 r0 += shapes[mid * offset + col] * xmid;
1870 }
1871 else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0 || mm == 3))
1872 r0 += shapes[col * offset + mid] * xmid;
1873
1874 if (add)
1875 {
1876 out[stride * col] += r0 + r1;
1877 if (type == 1 && contract_over_rows == false)
1878 out[stride * (nn - 1 - col)] += r1 - r0;
1879 else
1880 out[stride * (nn - 1 - col)] += r0 - r1;
1881 }
1882 else
1883 {
1884 out[stride * col] = r0 + r1;
1885 if (type == 1 && contract_over_rows == false)
1886 out[stride * (nn - 1 - col)] = r1 - r0;
1887 else
1888 out[stride * (nn - 1 - col)] = r0 - r1;
1889 }
1890 }
1891 if (type == 0 && contract_over_rows == true && nn % 2 == 1 &&
1892 mm % 2 == 1 && mm > 3)
1893 {
1894 if (add)
1895 out[stride * n_cols] += shapes[mid * offset + n_cols] * xmid;
1896 else
1897 out[stride * n_cols] = shapes[mid * offset + n_cols] * xmid;
1898 }
1899 else if (contract_over_rows == true && nn % 2 == 1)
1900 {
1901 Number r0;
1902 if (mid > 0)
1903 {
1904 r0 = shapes[n_cols] * xp[0];
1905 for (int ind = 1; ind < mid; ++ind)
1906 r0 += shapes[ind * offset + n_cols] * xp[ind];
1907 }
1908 else
1909 r0 = Number();
1910 if (type != 1 && mm % 2 == 1)
1911 r0 += shapes[mid * offset + n_cols] * xmid;
1912
1913 if (add)
1914 out[stride * n_cols] += r0;
1915 else
1916 out[stride * n_cols] = r0;
1917 }
1918 else if (contract_over_rows == false && nn % 2 == 1)
1919 {
1920 Number r0;
1921 if (mid > 0)
1922 {
1923 if (type == 1)
1924 {
1925 r0 = shapes[n_cols * offset] * xm[0];
1926 for (int ind = 1; ind < mid; ++ind)
1927 r0 += shapes[n_cols * offset + ind] * xm[ind];
1928 }
1929 else
1930 {
1931 r0 = shapes[n_cols * offset] * xp[0];
1932 for (int ind = 1; ind < mid; ++ind)
1933 r0 += shapes[n_cols * offset + ind] * xp[ind];
1934 }
1935 }
1936 else
1937 r0 = Number();
1938
1939 if ((type == 0 || type == 2) && mm % 2 == 1)
1940 r0 += shapes[n_cols * offset + mid] * xmid;
1941
1942 if (add)
1943 out[stride * n_cols] += r0;
1944 else
1945 out[stride * n_cols] = r0;
1946 }
1947 if (one_line == false)
1948 {
1949 in += 1;
1950 out += 1;
1951 }
1952 }
1953 if (one_line == false)
1954 {
1955 in += stride * (mm - 1);
1956 out += stride * (nn - 1);
1957 }
1958 }
1959 }
1960
1961
1962
1991 template <int dim,
1992 int n_rows,
1993 int n_columns,
1994 typename Number,
1995 typename Number2>
1996 struct EvaluatorTensorProduct<evaluate_symmetric_hierarchical,
1997 dim,
1998 n_rows,
1999 n_columns,
2000 Number,
2001 Number2>
2002 {
2003 static constexpr unsigned int n_rows_of_product =
2004 Utilities::pow(n_rows, dim);
2005 static constexpr unsigned int n_columns_of_product =
2006 Utilities::pow(n_columns, dim);
2007
2013 EvaluatorTensorProduct()
2014 : shape_values(nullptr)
2015 , shape_gradients(nullptr)
2016 , shape_hessians(nullptr)
2017 {}
2018
2023 EvaluatorTensorProduct(const AlignedVector<Number> &shape_values)
2024 : shape_values(shape_values.begin())
2025 , shape_gradients(nullptr)
2026 , shape_hessians(nullptr)
2027 {}
2028
2033 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
2034 const AlignedVector<Number2> &shape_gradients,
2035 const AlignedVector<Number2> &shape_hessians,
2036 const unsigned int dummy1 = 0,
2037 const unsigned int dummy2 = 0)
2038 : shape_values(shape_values.begin())
2039 , shape_gradients(shape_gradients.begin())
2040 , shape_hessians(shape_hessians.begin())
2041 {
2042 (void)dummy1;
2043 (void)dummy2;
2044 }
2045
2046 template <int direction, bool contract_over_rows, bool add>
2047 void
2048 values(const Number in[], Number out[]) const
2049 {
2050 Assert(shape_values != nullptr, ExcNotInitialized());
2051 apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
2052 }
2053
2054 template <int direction, bool contract_over_rows, bool add>
2055 void
2056 gradients(const Number in[], Number out[]) const
2057 {
2058 Assert(shape_gradients != nullptr, ExcNotInitialized());
2059 apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
2060 }
2061
2062 template <int direction, bool contract_over_rows, bool add>
2063 void
2064 hessians(const Number in[], Number out[]) const
2065 {
2066 Assert(shape_hessians != nullptr, ExcNotInitialized());
2067 apply<direction, contract_over_rows, add, 0>(shape_hessians, in, out);
2068 }
2069
2070 template <int direction, bool contract_over_rows, bool add>
2071 void
2072 values_one_line(const Number in[], Number out[]) const
2073 {
2074 Assert(shape_values != nullptr, ExcNotInitialized());
2075 apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
2076 }
2077
2078 template <int direction, bool contract_over_rows, bool add>
2079 void
2080 gradients_one_line(const Number in[], Number out[]) const
2081 {
2082 Assert(shape_gradients != nullptr, ExcNotInitialized());
2083 apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
2084 in,
2085 out);
2086 }
2087
2088 template <int direction, bool contract_over_rows, bool add>
2089 void
2090 hessians_one_line(const Number in[], Number out[]) const
2091 {
2092 Assert(shape_hessians != nullptr, ExcNotInitialized());
2093 apply<direction, contract_over_rows, add, 0, true>(shape_hessians,
2094 in,
2095 out);
2096 }
2097
2125 template <int direction,
2126 bool contract_over_rows,
2127 bool add,
2128 int type,
2129 bool one_line = false>
2130 static void
2131 apply(const Number2 *DEAL_II_RESTRICT shape_data,
2132 const Number * in,
2133 Number * out);
2134
2135 private:
2136 const Number2 *shape_values;
2137 const Number2 *shape_gradients;
2138 const Number2 *shape_hessians;
2139 };
2140
2141
2142
2143 template <int dim,
2144 int n_rows,
2145 int n_columns,
2146 typename Number,
2147 typename Number2>
2148 template <int direction,
2149 bool contract_over_rows,
2150 bool add,
2151 int type,
2152 bool one_line>
2153 inline void
2154 EvaluatorTensorProduct<evaluate_symmetric_hierarchical,
2155 dim,
2156 n_rows,
2157 n_columns,
2158 Number,
2159 Number2>::apply(const Number2 *DEAL_II_RESTRICT shapes,
2160 const Number * in,
2161 Number * out)
2162 {
2163 static_assert(one_line == false || direction == dim - 1,
2164 "Single-line evaluation only works for direction=dim-1.");
2165 static_assert(
2166 type == 0 || type == 1,
2167 "Only types 0 and 1 implemented for evaluate_symmetric_hierarchical.");
2168 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
2169 in != out,
2170 ExcMessage("In-place operation only supported for "
2171 "n_rows==n_columns or single-line interpolation"));
2172
2173 // We cannot statically assert that direction is less than dim, so must do
2174 // an additional dynamic check
2175 AssertIndexRange(direction, dim);
2176
2177 constexpr int nn = contract_over_rows ? n_columns : n_rows;
2178 constexpr int mm = contract_over_rows ? n_rows : n_columns;
2179 constexpr int n_cols = nn / 2;
2180 constexpr int mid = mm / 2;
2181
2182 constexpr int stride = Utilities::pow(n_columns, direction);
2183 constexpr int n_blocks1 = one_line ? 1 : stride;
2184 constexpr int n_blocks2 =
2185 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
2186
2187 // this code may look very inefficient at first sight due to the many
2188 // different cases with if's at the innermost loop part, but all of the
2189 // conditionals can be evaluated at compile time because they are
2190 // templates, so the compiler should optimize everything away
2191 for (int i2 = 0; i2 < n_blocks2; ++i2)
2192 {
2193 for (int i1 = 0; i1 < n_blocks1; ++i1)
2194 {
2195 if (contract_over_rows)
2196 {
2197 Number x[mm];
2198 for (unsigned int i = 0; i < mm; ++i)
2199 x[i] = in[stride * i];
2200 for (unsigned int col = 0; col < n_cols; ++col)
2201 {
2202 Number r0, r1;
2203 if (mid > 0)
2204 {
2205 r0 = shapes[col] * x[0];
2206 r1 = shapes[col + n_columns] * x[1];
2207 for (unsigned int ind = 1; ind < mid; ++ind)
2208 {
2209 r0 +=
2210 shapes[col + 2 * ind * n_columns] * x[2 * ind];
2211 r1 += shapes[col + (2 * ind + 1) * n_columns] *
2212 x[2 * ind + 1];
2213 }
2214 }
2215 else
2216 r0 = r1 = Number();
2217 if (mm % 2 == 1)
2218 r0 += shapes[col + (mm - 1) * n_columns] * x[mm - 1];
2219 if (add)
2220 {
2221 out[stride * col] += r0 + r1;
2222 if (type == 1)
2223 out[stride * (nn - 1 - col)] += r1 - r0;
2224 else
2225 out[stride * (nn - 1 - col)] += r0 - r1;
2226 }
2227 else
2228 {
2229 out[stride * col] = r0 + r1;
2230 if (type == 1)
2231 out[stride * (nn - 1 - col)] = r1 - r0;
2232 else
2233 out[stride * (nn - 1 - col)] = r0 - r1;
2234 }
2235 }
2236 if (nn % 2 == 1)
2237 {
2238 Number r0;
2239 const unsigned int shift = type == 1 ? 1 : 0;
2240 if (mid > 0)
2241 {
2242 r0 = shapes[n_cols + shift * n_columns] * x[shift];
2243 for (unsigned int ind = 1; ind < mid; ++ind)
2244 r0 += shapes[n_cols + (2 * ind + shift) * n_columns] *
2245 x[2 * ind + shift];
2246 }
2247 else
2248 r0 = 0;
2249 if (type != 1 && mm % 2 == 1)
2250 r0 += shapes[n_cols + (mm - 1) * n_columns] * x[mm - 1];
2251 if (add)
2252 out[stride * n_cols] += r0;
2253 else
2254 out[stride * n_cols] = r0;
2255 }
2256 }
2257 else
2258 {
2259 Number xp[mid + 1], xm[mid > 0 ? mid : 1];
2260 for (int i = 0; i < mid; ++i)
2261 if (type == 0)
2262 {
2263 xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
2264 xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
2265 }
2266 else
2267 {
2268 xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
2269 xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
2270 }
2271 if (mm % 2 == 1)
2272 xp[mid] = in[stride * mid];
2273 for (unsigned int col = 0; col < n_cols; ++col)
2274 {
2275 Number r0, r1;
2276 if (mid > 0)
2277 {
2278 r0 = shapes[2 * col * n_columns] * xp[0];
2279 r1 = shapes[(2 * col + 1) * n_columns] * xm[0];
2280 for (unsigned int ind = 1; ind < mid; ++ind)
2281 {
2282 r0 += shapes[2 * col * n_columns + ind] * xp[ind];
2283 r1 +=
2284 shapes[(2 * col + 1) * n_columns + ind] * xm[ind];
2285 }
2286 }
2287 else
2288 r0 = r1 = Number();
2289 if (mm % 2 == 1)
2290 {
2291 if (type == 1)
2292 r1 +=
2293 shapes[(2 * col + 1) * n_columns + mid] * xp[mid];
2294 else
2295 r0 += shapes[2 * col * n_columns + mid] * xp[mid];
2296 }
2297 if (add)
2298 {
2299 out[stride * (2 * col)] += r0;
2300 out[stride * (2 * col + 1)] += r1;
2301 }
2302 else
2303 {
2304 out[stride * (2 * col)] = r0;
2305 out[stride * (2 * col + 1)] = r1;
2306 }
2307 }
2308 if (nn % 2 == 1)
2309 {
2310 Number r0;
2311 if (mid > 0)
2312 {
2313 r0 = shapes[(nn - 1) * n_columns] * xp[0];
2314 for (unsigned int ind = 1; ind < mid; ++ind)
2315 r0 += shapes[(nn - 1) * n_columns + ind] * xp[ind];
2316 }
2317 else
2318 r0 = Number();
2319 if (mm % 2 == 1 && type == 0)
2320 r0 += shapes[(nn - 1) * n_columns + mid] * xp[mid];
2321 if (add)
2322 out[stride * (nn - 1)] += r0;
2323 else
2324 out[stride * (nn - 1)] = r0;
2325 }
2326 }
2327 if (one_line == false)
2328 {
2329 in += 1;
2330 out += 1;
2331 }
2332 }
2333 if (one_line == false)
2334 {
2335 in += stride * (mm - 1);
2336 out += stride * (nn - 1);
2337 }
2338 }
2339 }
2340
2341
2342
2362 template <int dim,
2363 int n_rows,
2364 int n_columns,
2365 typename Number,
2366 int normal_dir,
2367 typename Number2>
2368 struct EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
2369 dim,
2370 n_rows,
2371 n_columns,
2372 Number,
2373 normal_dir,
2374 Number2>
2375 {
2376 static constexpr unsigned int n_rows_of_product =
2377 numbers::invalid_unsigned_int;
2378 static constexpr unsigned int n_columns_of_product =
2379 numbers::invalid_unsigned_int;
2380
2385 EvaluatorTensorProductAnisotropic()
2386 : shape_values(nullptr)
2387 , shape_gradients(nullptr)
2388 , shape_hessians(nullptr)
2389 {}
2390
2394 EvaluatorTensorProductAnisotropic(
2395 const AlignedVector<Number2> &shape_values,
2396 const AlignedVector<Number2> &shape_gradients,
2397 const AlignedVector<Number2> &shape_hessians,
2398 const unsigned int dummy1 = 0,
2399 const unsigned int dummy2 = 0)
2400 : shape_values(shape_values.begin())
2401 , shape_gradients(shape_gradients.begin())
2402 , shape_hessians(shape_hessians.begin())
2403 {
2404 // We can enter this function either for the apply() path that has
2405 // n_rows * n_columns entries or for the apply_face() path that only has
2406 // n_rows * 3 entries in the array. Since we cannot decide about the use
2407 // we must allow for both here.
2408 Assert(shape_values.size() == 0 ||
2409 shape_values.size() == n_rows * n_columns ||
2410 shape_values.size() == 3 * n_rows,
2411 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
2412 Assert(shape_gradients.size() == 0 ||
2413 shape_gradients.size() == n_rows * n_columns,
2414 ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
2415 Assert(shape_hessians.size() == 0 ||
2416 shape_hessians.size() == n_rows * n_columns,
2417 ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
2418 (void)dummy1;
2419 (void)dummy2;
2420 }
2421
2422 template <int direction, bool contract_over_rows, bool add>
2423 void
2424 values(const Number in[], Number out[]) const
2425 {
2426 apply<direction, contract_over_rows, add>(shape_values, in, out);
2427 }
2428
2429 template <int direction, bool contract_over_rows, bool add>
2430 void
2431 gradients(const Number in[], Number out[]) const
2432 {
2433 apply<direction, contract_over_rows, add>(shape_gradients, in, out);
2434 }
2435
2436 template <int direction, bool contract_over_rows, bool add>
2437 void
2438 hessians(const Number in[], Number out[]) const
2439 {
2440 apply<direction, contract_over_rows, add>(shape_hessians, in, out);
2441 }
2442
2471 template <int direction,
2472 bool contract_over_rows,
2473 bool add,
2474 bool one_line = false>
2475 static void
2476 apply(const Number2 *DEAL_II_RESTRICT shape_data,
2477 const Number * in,
2478 Number * out);
2479
2480 template <int face_direction,
2481 bool contract_onto_face,
2482 bool add,
2483 int max_derivative>
2484 void
2485 apply_face(const Number *DEAL_II_RESTRICT in,
2486 Number *DEAL_II_RESTRICT out) const;
2487
2488 private:
2489 const Number2 *shape_values;
2490 const Number2 *shape_gradients;
2491 const Number2 *shape_hessians;
2492 };
2493
2494 template <int dim,
2495 int n_rows,
2496 int n_columns,
2497 typename Number,
2498 int normal_dir,
2499 typename Number2>
2500 template <int direction, bool contract_over_rows, bool add, bool one_line>
2501 inline void
2502 EvaluatorTensorProductAnisotropic<
2503 evaluate_raviart_thomas,
2504 dim,
2505 n_rows,
2506 n_columns,
2507 Number,
2508 normal_dir,
2509 Number2>::apply(const Number2 *DEAL_II_RESTRICT shape_data,
2510 const Number * in,
2511 Number * out)
2512 {
2513 static_assert(one_line == false || direction == dim - 1,
2514 "Single-line evaluation only works for direction=dim-1.");
2515 Assert(shape_data != nullptr,
2516 ExcMessage(
2517 "The given array shape_data must not be the null pointer!"));
2518 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
2519 in != out,
2520 ExcMessage("In-place operation only supported for "
2521 "n_rows==n_columns or single-line interpolation"));
2522 AssertIndexRange(direction, dim);
2523 constexpr int mm = contract_over_rows ? n_rows : n_columns,
2524 nn = contract_over_rows ? n_columns : n_rows;
2525
2526 constexpr int stride = Utilities::pow(n_columns, direction);
2527 constexpr int n_blocks1 = one_line ? 1 : stride;
2528
2529 // The number of blocks depend on both direction and dimension.
2530 constexpr int n_blocks2 =
2531 (dim - direction - 1 == 0) ?
2532 1 :
2533 ((direction == normal_dir) ?
2534 Utilities::pow((n_rows - 1),
2535 (direction >= dim) ? 0 : dim - direction - 1) :
2536 (((direction < normal_dir) ? (n_rows + 1) : n_rows) *
2537 ((dim - direction == 3) ? n_rows : 1)));
2538
2539 for (int i2 = 0; i2 < n_blocks2; ++i2)
2540 {
2541 for (int i1 = 0; i1 < n_blocks1; ++i1)
2542 {
2543 Number x[mm];
2544 for (int i = 0; i < mm; ++i)
2545 x[i] = in[stride * i];
2546
2547 for (int col = 0; col < nn; ++col)
2548 {
2549 Number2 val0;
2550
2551 if (contract_over_rows)
2552 val0 = shape_data[col];
2553 else
2554 val0 = shape_data[col * n_columns];
2555
2556 Number res0 = val0 * x[0];
2557 for (int i = 1; i < mm; ++i)
2558 {
2559 if (contract_over_rows)
2560 val0 = shape_data[i * n_columns + col];
2561 else
2562 val0 = shape_data[col * n_columns + i];
2563
2564 res0 += val0 * x[i];
2565 }
2566 if (add)
2567 out[stride * col] += res0;
2568
2569 else
2570 out[stride * col] = res0;
2571 }
2572
2573 if (one_line == false)
2574 {
2575 ++in;
2576 ++out;
2577 }
2578 }
2579 if (one_line == false)
2580 {
2581 in += stride * (mm - 1);
2582 out += stride * (nn - 1);
2583 }
2584 }
2585 }
2586
2587 template <int dim,
2588 int n_rows,
2589 int n_columns,
2590 typename Number,
2591 int normal_dir,
2592 typename Number2>
2593 template <int face_direction,
2594 bool contract_onto_face,
2595 bool add,
2596 int max_derivative>
2597 inline void
2598 EvaluatorTensorProductAnisotropic<
2599 evaluate_raviart_thomas,
2600 dim,
2601 n_rows,
2602 n_columns,
2603 Number,
2604 normal_dir,
2605 Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
2606 Number *DEAL_II_RESTRICT out) const
2607 {
2608 Assert(dim > 1 && dim < 4, ExcMessage("Only dim=2,3 supported"));
2609 static_assert(max_derivative >= 0 && max_derivative < 3,
2610 "Only derivative orders 0-2 implemented");
2611 Assert(shape_values != nullptr,
2612 ExcMessage(
2613 "The given array shape_values must not be the null pointer."));
2614
2615 // Determine the number of blocks depending on the face and normaldirection,
2616 // as well as dimension.
2617 constexpr int n_blocks1 = (face_direction == normal_dir) ? (n_rows - 1) :
2618 ((face_direction == 0 && normal_dir == 2) ||
2619 (face_direction == 1 && normal_dir == 2) ||
2620 (face_direction == 2 && normal_dir == 1)) ?
2621 n_rows :
2622 (n_rows + 1);
2623 constexpr int n_blocks2 = (dim == 2) ?
2624 1 :
2625 ((face_direction == normal_dir) ?
2626 (n_rows - 1) :
2627 (((face_direction == 0 && normal_dir == 1) ||
2628 (face_direction == 1 && normal_dir == 0) ||
2629 (face_direction == 2 && normal_dir == 0)) ?
2630 n_rows :
2631 (n_rows + 1)));
2632
2633 AssertIndexRange(face_direction, dim);
2634
2635 constexpr int in_stride =
2636 (face_direction == normal_dir) ?
2637 Utilities::pow(n_rows - 1, face_direction) :
2638 ((face_direction == 0) ?
2639 1 :
2640 ((face_direction == 2) ?
2641 n_rows * (n_rows + 1) :
2642 ((face_direction == 1 && normal_dir == 0) ? (n_rows + 1) :
2643 n_rows)));
2644 constexpr int out_stride = n_blocks1 * n_blocks2;
2645
2646 const Number *DEAL_II_RESTRICT shape_values = this->shape_values;
2647
2648 for (int i2 = 0; i2 < n_blocks2; ++i2)
2649 {
2650 for (int i1 = 0; i1 < n_blocks1; ++i1)
2651 {
2652 if (contract_onto_face == true)
2653 {
2654 Number res0 = shape_values[0] * in[0];
2655 Number res1, res2;
2656
2657 if (max_derivative > 0)
2658 res1 = shape_values[n_rows] * in[0];
2659
2660 if (max_derivative > 1)
2661 res2 = shape_values[2 * n_rows] * in[0];
2662
2663 for (int ind = 1; ind < n_rows; ++ind)
2664 {
2665 res0 += shape_values[ind] * in[in_stride * ind];
2666 if (max_derivative > 0)
2667 res1 += shape_values[ind + n_rows] * in[in_stride * ind];
2668
2669 if (max_derivative > 1)
2670 res2 +=
2671 shape_values[ind + 2 * n_rows] * in[in_stride * ind];
2672 }
2673 if (add)
2674 {
2675 out[0] += res0;
2676
2677 if (max_derivative > 0)
2678 out[out_stride] += res1;
2679
2680 if (max_derivative > 1)
2681 out[2 * out_stride] += res2;
2682 }
2683 else
2684 {
2685 out[0] = res0;
2686
2687 if (max_derivative > 0)
2688 out[out_stride] = res1;
2689
2690 if (max_derivative > 1)
2691 out[2 * out_stride] = res2;
2692 }
2693 }
2694 else
2695 {
2696 for (int col = 0; col < n_rows; ++col)
2697 {
2698 if (add)
2699 out[col * in_stride] += shape_values[col] * in[0];
2700 else
2701 out[col * in_stride] = shape_values[col] * in[0];
2702
2703 if (max_derivative > 0)
2704 out[col * in_stride] +=
2705 shape_values[col + n_rows] * in[out_stride];
2706
2707 if (max_derivative > 1)
2708 out[col * in_stride] +=
2709 shape_values[col + 2 * n_rows] * in[2 * out_stride];
2710 }
2711 }
2712
2713 // increment: in regular case, just go to the next point in
2714 // x-direction. If we are at the end of one chunk in x-dir, need
2715 // to jump over to the next layer in z-direction
2716 switch (face_direction)
2717 {
2718 case 0:
2719 in += contract_onto_face ? n_rows : 1;
2720 out += contract_onto_face ? 1 : n_rows;
2721 break;
2722
2723 case 1:
2724 ++in;
2725 ++out;
2726 // faces 2 and 3 in 3D use local coordinate system zx, which
2727 // is the other way around compared to the tensor
2728 // product. Need to take that into account.
2729 if (dim == 3)
2730 {
2731 if (normal_dir == 0)
2732 {
2733 if (contract_onto_face)
2734 out += n_rows - 1;
2735 else
2736 in += n_rows - 1;
2737 }
2738 if (normal_dir == 1)
2739 {
2740 if (contract_onto_face)
2741 out += n_rows - 2;
2742 else
2743 in += n_rows - 2;
2744 }
2745 if (normal_dir == 2)
2746 {
2747 if (contract_onto_face)
2748 out += n_rows;
2749 else
2750 in += n_rows;
2751 }
2752 }
2753 break;
2754
2755 case 2:
2756 ++in;
2757 ++out;
2758 break;
2759
2760 default:
2761 Assert(false, ExcNotImplemented());
2762 }
2763 }
2764 if (face_direction == 1 && dim == 3)
2765 {
2766 // adjust for local coordinate system zx
2767 if (contract_onto_face)
2768 {
2769 if (normal_dir == 0)
2770 {
2771 in += (n_rows + 1) * (n_rows - 1);
2772 out -= n_rows * (n_rows + 1) - 1;
2773 }
2774 if (normal_dir == 1)
2775 {
2776 in += (n_rows - 1) * (n_rows - 1);
2777 out -= (n_rows - 1) * (n_rows - 1) - 1;
2778 }
2779 if (normal_dir == 2)
2780 {
2781 in += (n_rows - 1) * (n_rows);
2782 out -= (n_rows) * (n_rows + 1) - 1;
2783 }
2784 }
2785 else
2786 {
2787 if (normal_dir == 0)
2788 {
2789 out += (n_rows + 1) * (n_rows - 1);
2790 in -= n_rows * (n_rows + 1) - 1;
2791 }
2792 if (normal_dir == 1)
2793 {
2794 out += (n_rows - 1) * (n_rows - 1);
2795 in -= (n_rows - 1) * (n_rows - 1) - 1;
2796 }
2797 if (normal_dir == 2)
2798 {
2799 out += (n_rows - 1) * (n_rows);
2800 in -= (n_rows) * (n_rows + 1) - 1;
2801 }
2802 }
2803 }
2804 }
2805 }
2806
2807
2808
2815 template <typename Number, typename Number2>
2816 struct ProductTypeNoPoint
2817 {
2818 using type = typename ProductType<Number, Number2>::type;
2819 };
2820
2821 template <int dim, typename Number, typename Number2>
2822 struct ProductTypeNoPoint<Point<dim, Number>, Number2>
2823 {
2824 using type = typename ProductType<Tensor<1, dim, Number>, Number2>::type;
2825 };
2826
2827
2828
2863 template <int dim, typename Number, typename Number2>
2864 inline std::pair<
2865 typename ProductTypeNoPoint<Number, Number2>::type,
2866 Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>>
2867 evaluate_tensor_product_value_and_gradient(
2868 const std::vector<Polynomials::Polynomial<double>> &poly,
2869 const std::vector<Number> & values,
2870 const Point<dim, Number2> & p,
2871 const bool d_linear = false,
2872 const std::vector<unsigned int> & renumber = {})
2873 {
2874 static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
2875
2876 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
2877
2878 // use `int` type for this variable and the loops below to inform the
2879 // compiler that the loops below will never overflow, which allows it to
2880 // generate more optimized code for the variable loop bounds in the
2881 // present context
2882 const int n_shapes = poly.size();
2883 AssertDimension(Utilities::pow(n_shapes, dim), values.size());
2884 Assert(renumber.empty() || renumber.size() == values.size(),
2885 ExcDimensionMismatch(renumber.size(), values.size()));
2886
2887 // shortcut for linear interpolation to speed up evaluation
2888 if (d_linear)
2889 {
2890 AssertDimension(poly.size(), 2);
2891 for (unsigned int i = 0; i < renumber.size(); ++i)
2892 AssertDimension(renumber[i], i);
2893
2894 if (dim == 1)
2895 {
2896 Tensor<1, dim, Number3> derivative;
2897 derivative[0] = values[1] - values[0];
2898 return std::make_pair((1. - p[0]) * values[0] + p[0] * values[1],
2899 derivative);
2900 }
2901 else if (dim == 2)
2902 {
2903 const Number2 x0 = 1. - p[0], x1 = p[0];
2904 const Number3 tmp0 = x0 * values[0] + x1 * values[1];
2905 const Number3 tmp1 = x0 * values[2] + x1 * values[3];
2906 const Number3 mapped = (1. - p[1]) * tmp0 + p[1] * tmp1;
2907 Tensor<1, dim, Number3> derivative;
2908 derivative[0] = (1. - p[1]) * (values[1] - values[0]) +
2909 p[1] * (values[3] - values[2]);
2910 derivative[1] = tmp1 - tmp0;
2911 return std::make_pair(mapped, derivative);
2912 }
2913 else if (dim == 3)
2914 {
2915 const Number2 x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
2916 z0 = 1. - p[2], z1 = p[2];
2917 const Number3 tmp0 = x0 * values[0] + x1 * values[1];
2918 const Number3 tmp1 = x0 * values[2] + x1 * values[3];
2919 const Number3 tmpy0 = y0 * tmp0 + y1 * tmp1;
2920 const Number3 tmp2 = x0 * values[4] + x1 * values[5];
2921 const Number3 tmp3 = x0 * values[6] + x1 * values[7];
2922 const Number3 tmpy1 = y0 * tmp2 + y1 * tmp3;
2923 const Number3 mapped = z0 * tmpy0 + z1 * tmpy1;
2924 Tensor<1, dim, Number3> derivative;
2925 derivative[2] = tmpy1 - tmpy0;
2926 derivative[1] = z0 * (tmp1 - tmp0) + z1 * (tmp3 - tmp2);
2927 derivative[0] =
2928 z0 *
2929 (y0 * (values[1] - values[0]) + y1 * (values[3] - values[2])) +
2930 z1 *
2931 (y0 * (values[5] - values[4]) + y1 * (values[7] - values[6]));
2932 return std::make_pair(mapped, derivative);
2933 }
2934 }
2935
2936 AssertIndexRange(n_shapes, 200);
2937 ::ndarray<Number2, 200, 2, dim> shapes;
2938
2939 // Evaluate 1D polynomials and their derivatives
2940 std::array<Number2, dim> point;
2941 for (unsigned int d = 0; d < dim; ++d)
2942 point[d] = p[d];
2943 for (int i = 0; i < n_shapes; ++i)
2944 poly[i].values_of_array(point, 1, &shapes[i][0]);
2945
2946 // Go through the tensor product of shape functions and interpolate
2947 // with optimal algorithm
2948 std::pair<Number3, Tensor<1, dim, Number3>> result = {};
2949 for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
2950 {
2951 Number3 value_y = {}, deriv_x = {}, deriv_y = {};
2952 for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
2953 {
2954 // Interpolation + derivative x direction
2955 Number3 value = {}, deriv = {};
2956
2957 // Distinguish the inner loop based on whether we have a
2958 // renumbering or not
2959 if (renumber.empty())
2960 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
2961 {
2962 value += shapes[i0][0][0] * values[i];
2963 deriv += shapes[i0][1][0] * values[i];
2964 }
2965 else
2966 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
2967 {
2968 value += shapes[i0][0][0] * values[renumber[i]];
2969 deriv += shapes[i0][1][0] * values[renumber[i]];
2970 }
2971
2972 // Interpolation + derivative in y direction
2973 if (dim > 1)
2974 {
2975 value_y += value * shapes[i1][0][1];
2976 deriv_x += deriv * shapes[i1][0][1];
2977 deriv_y += value * shapes[i1][1][1];
2978 }
2979 else
2980 {
2981 result.first = value;
2982 result.second[0] = deriv;
2983 }
2984 }
2985 if (dim == 3)
2986 {
2987 // Interpolation + derivative in z direction
2988 result.first += value_y * shapes[i2][0][2];
2989 result.second[0] += deriv_x * shapes[i2][0][2];
2990 result.second[1] += deriv_y * shapes[i2][0][2];
2991 result.second[2] += value_y * shapes[i2][1][2];
2992 }
2993 else if (dim == 2)
2994 {
2995 result.first = value_y;
2996 result.second[0] = deriv_x;
2997 result.second[1] = deriv_y;
2998 }
2999 }
3000
3001 return result;
3002 }
3003
3004
3005
3006 template <int dim, typename Number, typename Number2>
3007 SymmetricTensor<2, dim, typename ProductTypeNoPoint<Number, Number2>::type>
3008 evaluate_tensor_product_hessian(
3009 const std::vector<Polynomials::Polynomial<double>> &poly,
3010 const std::vector<Number> & values,
3011 const Point<dim, Number2> & p,
3012 const std::vector<unsigned int> & renumber = {})
3013 {
3014 static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
3015
3016 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
3017
3018 // use `int` type for this variable and the loops below to inform the
3019 // compiler that the loops below will never overflow, which allows it to
3020 // generate more optimized code for the variable loop bounds in the
3021 // present context
3022 const int n_shapes = poly.size();
3023 AssertDimension(Utilities::pow(n_shapes, dim), values.size());
3024 Assert(renumber.empty() || renumber.size() == values.size(),
3025 ExcDimensionMismatch(renumber.size(), values.size()));
3026
3027 AssertIndexRange(n_shapes, 200);
3028 ::ndarray<Number2, 200, 3, dim> shapes;
3029
3030 // Evaluate 1D polynomials and their derivatives
3031 std::array<Number2, dim> point;
3032 for (unsigned int d = 0; d < dim; ++d)
3033 point[d] = p[d];
3034 for (int i = 0; i < n_shapes; ++i)
3035 poly[i].values_of_array(point, 2, &shapes[i][0]);
3036
3037 // Go through the tensor product of shape functions and interpolate
3038 // with optimal algorithm
3039 SymmetricTensor<2, dim, Number3> result;
3040 for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
3041 {
3042 Number3 value_y = {}, deriv_x = {}, deriv_y = {}, deriv_xx = {},
3043 deriv_xy = {}, deriv_yy = {};
3044 for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
3045 {
3046 // Interpolation + derivative x direction
3047 Number3 value = {}, deriv_1 = {}, deriv_2 = {};
3048
3049 // Distinguish the inner loop based on whether we have a
3050 // renumbering or not
3051 if (renumber.empty())
3052 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3053 {
3054 value += shapes[i0][0][0] * values[i];
3055 deriv_1 += shapes[i0][1][0] * values[i];
3056 deriv_2 += shapes[i0][2][0] * values[i];
3057 }
3058 else
3059 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3060 {
3061 value += shapes[i0][0][0] * values[renumber[i]];
3062 deriv_1 += shapes[i0][1][0] * values[renumber[i]];
3063 deriv_2 += shapes[i0][2][0] * values[renumber[i]];
3064 }
3065
3066 // Interpolation + derivative in y direction
3067 if (dim > 1)
3068 {
3069 if (dim > 2)
3070 {
3071 value_y += value * shapes[i1][0][1];
3072 deriv_x += deriv_1 * shapes[i1][0][1];
3073 deriv_y += value * shapes[i1][1][1];
3074 }
3075 deriv_xx += deriv_2 * shapes[i1][0][1];
3076 deriv_xy += deriv_1 * shapes[i1][1][1];
3077 deriv_yy += value * shapes[i1][2][1];
3078 }
3079 else
3080 {
3081 result[0][0] = deriv_2;
3082 }
3083 }
3084 if (dim == 3)
3085 {
3086 // Interpolation + derivative in z direction
3087 result[0][0] += deriv_xx * shapes[i2][0][2];
3088 result[0][1] += deriv_xy * shapes[i2][0][2];
3089 result[0][2] += deriv_x * shapes[i2][1][2];
3090 result[1][1] += deriv_yy * shapes[i2][0][2];
3091 result[1][2] += deriv_y * shapes[i2][1][2];
3092 result[2][2] += value_y * shapes[i2][2][2];
3093 }
3094 else if (dim == 2)
3095 {
3096 result[0][0] = deriv_xx;
3097 result[1][0] = deriv_xy;
3098 result[1][1] = deriv_yy;
3099 }
3100 }
3101
3102 return result;
3103 }
3104
3105
3106
3110 template <int dim, typename Number, typename Number2>
3111 inline void
3112 integrate_add_tensor_product_value_and_gradient(
3113 const std::vector<Polynomials::Polynomial<double>> &poly,
3114 const Number2 & value,
3115 const Tensor<1, dim, Number2> & gradient,
3116 const Point<dim, Number> & p,
3117 AlignedVector<Number2> & values,
3118 const std::vector<unsigned int> & renumber = {})
3119 {
3120 static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
3121
3122 // as in evaluate, use `int` type to produce better code in this context
3123 const int n_shapes = poly.size();
3124 AssertDimension(Utilities::pow(n_shapes, dim), values.size());
3125 Assert(renumber.empty() || renumber.size() == values.size(),
3126 ExcDimensionMismatch(renumber.size(), values.size()));
3127
3128 AssertIndexRange(n_shapes, 200);
3129 ::ndarray<Number, 200, 2, dim> shapes;
3130
3131 // Evaluate 1D polynomials and their derivatives
3132 std::array<Number, dim> point;
3133 for (unsigned int d = 0; d < dim; ++d)
3134 point[d] = p[d];
3135 for (int i = 0; i < n_shapes; ++i)
3136 poly[i].values_of_array(point, 1, &shapes[i][0]);
3137
3138 // Implement the transpose of the function above
3139 for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
3140 {
3141 const Number2 test_value_z =
3142 dim > 2 ?
3143 (value * shapes[i2][0][2] + gradient[2] * shapes[i2][1][2]) :
3144 value;
3145 const Number2 test_grad_x =
3146 dim > 2 ? gradient[0] * shapes[i2][0][2] : gradient[0];
3147 const Number2 test_grad_y = dim > 2 ?
3148 gradient[1] * shapes[i2][0][2] :
3149 (dim > 1 ? gradient[1] : Number2());
3150 for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
3151 {
3152 const Number2 test_value_y = dim > 1 ?
3153 (test_value_z * shapes[i1][0][1] +
3154 test_grad_y * shapes[i1][1][1]) :
3155 test_value_z;
3156 const Number2 test_grad_xy =
3157 dim > 1 ? test_grad_x * shapes[i1][0][1] : test_grad_x;
3158 if (renumber.empty())
3159 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3160 values[i] += shapes[i0][0][0] * test_value_y +
3161 shapes[i0][1][0] * test_grad_xy;
3162 else
3163 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3164 values[renumber[i]] += shapes[i0][0][0] * test_value_y +
3165 shapes[i0][1][0] * test_grad_xy;
3166 }
3167 }
3168 }
3169
3170
3171 template <int dim, int loop_length_template, typename Number>
3172 inline void
3173 weight_fe_q_dofs_by_entity(const VectorizedArray<Number> *weights,
3174 const unsigned int n_components,
3175 const int loop_length_non_template,
3176 VectorizedArray<Number> *data)
3177 {
3178 const int loop_length = loop_length_template != -1 ?
3179 loop_length_template :
3180 loop_length_non_template;
3181
3182 Assert(loop_length > 0, ExcNotImplemented());
3183 Assert(loop_length < 100, ExcNotImplemented());
3184 unsigned int degree_to_3[100];
3185 degree_to_3[0] = 0;
3186 for (int i = 1; i < loop_length - 1; ++i)
3187 degree_to_3[i] = 1;
3188 degree_to_3[loop_length - 1] = 2;
3189 for (unsigned int c = 0; c < n_components; ++c)
3190 for (int k = 0; k < (dim > 2 ? loop_length : 1); ++k)
3191 for (int j = 0; j < (dim > 1 ? loop_length : 1); ++j)
3192 {
3193 const unsigned int shift = 9 * degree_to_3[k] + 3 * degree_to_3[j];
3194 data[0] *= weights[shift];
3195 // loop bound as int avoids compiler warnings in case loop_length
3196 // == 1 (polynomial degree 0)
3197 for (int i = 1; i < loop_length - 1; ++i)
3198 data[i] *= weights[shift + 1];
3199 data[loop_length - 1] *= weights[shift + 2];
3200 data += loop_length;
3201 }
3202 }
3203
3204
3205} // end of namespace internal
3206
3207
3208DEAL_II_NAMESPACE_CLOSE
3209
3210#endif
aligned_vector.h
AlignedVector::empty
bool empty() const
AlignedVector::size
size_type size() const
Point
Definition: point.h:111
Tensor
Definition: tensor.h:503
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
DEAL_II_RESTRICT
#define DEAL_II_RESTRICT
Definition: config.h:103
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1473
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
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)
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:190
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:462
Utilities::fixed_power
T fixed_power(const T t)
Definition: utilities.h:1123
internal::integrate_add_tensor_product_value_and_gradient
void integrate_add_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double > > &poly, const Number2 &value, const Tensor< 1, dim, Number2 > &gradient, const Point< dim, Number > &p, AlignedVector< Number2 > &values, const std::vector< unsigned int > &renumber={})
Definition: tensor_product_kernels.h:3112
internal::weight_fe_q_dofs_by_entity
void weight_fe_q_dofs_by_entity(const VectorizedArray< Number > *weights, const unsigned int n_components, const int loop_length_non_template, VectorizedArray< Number > *data)
Definition: tensor_product_kernels.h:3173
internal::evaluate_tensor_product_value_and_gradient
std::pair< typename ProductTypeNoPoint< Number, Number2 >::type, Tensor< 1, dim, typename ProductTypeNoPoint< Number, Number2 >::type > > evaluate_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double > > &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
Definition: tensor_product_kernels.h:2867
internal::evaluate_tensor_product_hessian
SymmetricTensor< 2, dim, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_hessian(const std::vector< Polynomials::Polynomial< double > > &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const std::vector< unsigned int > &renumber={})
Definition: tensor_product_kernels.h:3008
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:201
ndarray
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108
ProductType::type
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
Definition: template_constraints.h:488
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, Number, normal_dir, Number2 >::gradients
void gradients(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:2431
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, Number, normal_dir, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:2491
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, Number, normal_dir, Number2 >::hessians
void hessians(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:2438
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, Number, normal_dir, Number2 >::EvaluatorTensorProductAnisotropic
EvaluatorTensorProductAnisotropic()
Definition: tensor_product_kernels.h:2385
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, Number, normal_dir, Number2 >::values
void values(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:2424
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, Number, normal_dir, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:2490
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, Number, normal_dir, Number2 >::EvaluatorTensorProductAnisotropic
EvaluatorTensorProductAnisotropic(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition: tensor_product_kernels.h:2394
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, Number, normal_dir, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:2489
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition: tensor_product_kernels.h:1645
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::hessians_one_line
void hessians_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1710
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values)
Definition: tensor_product_kernels.h:1633
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::gradients_one_line
void gradients_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1700
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:1757
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:1759
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:1758
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::values_one_line
void values_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1692
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::values
void values(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1668
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::gradients
void gradients(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1676
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:1623
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::hessians
void hessians(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1684
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const Number2 *shape_values, const Number2 *shape_gradients, const Number2 *shape_hessians, const unsigned int n_rows, const unsigned int n_columns)
Definition: tensor_product_kernels.h:643
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::n_columns
const unsigned int n_columns
Definition: tensor_product_kernels.h:721
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:718
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::hessians_one_line
void hessians_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:694
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:602
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::gradients
void gradients(const Number *in, Number *out) const
Definition: tensor_product_kernels.h:664
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::values_one_line
void values_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:678
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::gradients_one_line
void gradients_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:686
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:717
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::values
void values(const Number *in, Number *out) const
Definition: tensor_product_kernels.h:657
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:719
internal::EvaluatorTensorProduct< evaluate_general, 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, const unsigned int n_columns)
Definition: tensor_product_kernels.h:613
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::n_rows
const unsigned int n_rows
Definition: tensor_product_kernels.h:720
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::hessians
void hessians(const Number *in, Number *out) const
Definition: tensor_product_kernels.h:671
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::gradients
void gradients(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:243
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:350
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition: tensor_product_kernels.h:207
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:351
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:198
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::values
void values(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:236
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::hessians
void hessians(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:250
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::hessians_one_line
void hessians_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:273
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::values_one_line
void values_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:257
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:352
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::gradients_one_line
void gradients_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:265
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:1076
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:1077
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition: tensor_product_kernels.h:1040
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:1075
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:2138
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::gradients
void gradients(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:2056
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::values
void values(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:2048
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::hessians_one_line
void hessians_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:2090
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number > &shape_values)
Definition: tensor_product_kernels.h:2023
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::values_one_line
void values_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:2072
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition: tensor_product_kernels.h:2033
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::gradients_one_line
void gradients_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:2080
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::hessians
void hessians(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:2064
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:2137
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:2136
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:2013
internal::ProductTypeNoPoint< Point< dim, Number >, Number2 >::type
typename ProductType< Tensor< 1, dim, Number >, Number2 >::type type
Definition: tensor_product_kernels.h:2824
internal::ProductTypeNoPoint::type
typename ProductType< Number, Number2 >::type type
Definition: tensor_product_kernels.h:2818
Generated by doxygen 1.9.5