Reference documentation for deal.II version 9.3.3
\(\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\}}\)
tensor_product_kernels.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2017 - 2021 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/polynomial.h>
24#include <deal.II/base/utilities.h>
25
26
27DEAL_II_NAMESPACE_OPEN
28
29
30
31namespace internal
32{
37 enum EvaluatorVariant
38 {
43 evaluate_general,
49 evaluate_symmetric,
55 evaluate_evenodd,
66 evaluate_symmetric_hierarchical
67 };
68
69
70
74 enum class EvaluatorQuantity
75 {
79 value,
83 gradient,
87 hessian
88 };
89
90
91
112 template <EvaluatorVariant variant,
113 int dim,
114 int n_rows,
115 int n_columns,
116 typename Number,
117 typename Number2 = Number>
118 struct EvaluatorTensorProduct
119 {};
120
121
122
140 template <int dim,
141 int n_rows,
142 int n_columns,
143 typename Number,
144 typename Number2>
145 struct EvaluatorTensorProduct<evaluate_general,
146 dim,
147 n_rows,
148 n_columns,
149 Number,
150 Number2>
151 {
152 static constexpr unsigned int n_rows_of_product =
153 Utilities::pow(n_rows, dim);
154 static constexpr unsigned int n_columns_of_product =
155 Utilities::pow(n_columns, dim);
156
161 EvaluatorTensorProduct()
162 : shape_values(nullptr)
163 , shape_gradients(nullptr)
164 , shape_hessians(nullptr)
165 {}
166
170 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
171 const AlignedVector<Number2> &shape_gradients,
172 const AlignedVector<Number2> &shape_hessians,
173 const unsigned int dummy1 = 0,
174 const unsigned int dummy2 = 0)
175 : shape_values(shape_values.begin())
176 , shape_gradients(shape_gradients.begin())
177 , shape_hessians(shape_hessians.begin())
178 {
179 // We can enter this function either for the apply() path that has
180 // n_rows * n_columns entries or for the apply_face() path that only has
181 // n_rows * 3 entries in the array. Since we cannot decide about the use
182 // we must allow for both here.
183 Assert(shape_values.size() == 0 ||
184 shape_values.size() == n_rows * n_columns ||
185 shape_values.size() == 3 * n_rows,
186 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
187 Assert(shape_gradients.size() == 0 ||
188 shape_gradients.size() == n_rows * n_columns,
189 ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
190 Assert(shape_hessians.size() == 0 ||
191 shape_hessians.size() == n_rows * n_columns,
192 ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
193 (void)dummy1;
194 (void)dummy2;
195 }
196
197 template <int direction, bool contract_over_rows, bool add>
198 void
199 values(const Number in[], Number out[]) const
200 {
201 apply<direction, contract_over_rows, add>(shape_values, in, out);
202 }
203
204 template <int direction, bool contract_over_rows, bool add>
205 void
206 gradients(const Number in[], Number out[]) const
207 {
208 apply<direction, contract_over_rows, add>(shape_gradients, in, out);
209 }
210
211 template <int direction, bool contract_over_rows, bool add>
212 void
213 hessians(const Number in[], Number out[]) const
214 {
215 apply<direction, contract_over_rows, add>(shape_hessians, in, out);
216 }
217
218 template <int direction, bool contract_over_rows, bool add>
219 void
220 values_one_line(const Number in[], Number out[]) const
221 {
222 Assert(shape_values != nullptr, ExcNotInitialized());
223 apply<direction, contract_over_rows, add, true>(shape_values, in, out);
224 }
225
226 template <int direction, bool contract_over_rows, bool add>
227 void
228 gradients_one_line(const Number in[], Number out[]) const
229 {
230 Assert(shape_gradients != nullptr, ExcNotInitialized());
231 apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
232 }
233
234 template <int direction, bool contract_over_rows, bool add>
235 void
236 hessians_one_line(const Number in[], Number out[]) const
237 {
238 Assert(shape_hessians != nullptr, ExcNotInitialized());
239 apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
240 }
241
266 template <int direction,
267 bool contract_over_rows,
268 bool add,
269 bool one_line = false>
270 static void
271 apply(const Number2 *DEAL_II_RESTRICT shape_data,
272 const Number * in,
273 Number * out);
274
309 template <int face_direction,
310 bool contract_onto_face,
311 bool add,
312 int max_derivative,
313 bool lex_faces = false>
314 void
315 apply_face(const Number *DEAL_II_RESTRICT in,
316 Number *DEAL_II_RESTRICT out) const;
317
318 const Number2 *shape_values;
319 const Number2 *shape_gradients;
320 const Number2 *shape_hessians;
321 };
322
323
324
325 template <int dim,
326 int n_rows,
327 int n_columns,
328 typename Number,
329 typename Number2>
330 template <int direction, bool contract_over_rows, bool add, bool one_line>
331 inline void
332 EvaluatorTensorProduct<evaluate_general,
333 dim,
334 n_rows,
335 n_columns,
336 Number,
337 Number2>::apply(const Number2 *DEAL_II_RESTRICT
338 shape_data,
339 const Number * in,
340 Number * out)
341 {
342 static_assert(one_line == false || direction == dim - 1,
343 "Single-line evaluation only works for direction=dim-1.");
344 Assert(shape_data != nullptr,
345 ExcMessage(
346 "The given array shape_data must not be the null pointer!"));
347 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
348 in != out,
349 ExcMessage("In-place operation only supported for "
350 "n_rows==n_columns or single-line interpolation"));
351 AssertIndexRange(direction, dim);
352 constexpr int mm = contract_over_rows ? n_rows : n_columns,
353 nn = contract_over_rows ? n_columns : n_rows;
354
355 constexpr int stride = Utilities::pow(n_columns, direction);
356 constexpr int n_blocks1 = one_line ? 1 : stride;
357 constexpr int n_blocks2 =
358 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
359
360 for (int i2 = 0; i2 < n_blocks2; ++i2)
361 {
362 for (int i1 = 0; i1 < n_blocks1; ++i1)
363 {
364 Number x[mm];
365 for (int i = 0; i < mm; ++i)
366 x[i] = in[stride * i];
367 for (int col = 0; col < nn; ++col)
368 {
369 Number2 val0;
370 if (contract_over_rows == true)
371 val0 = shape_data[col];
372 else
373 val0 = shape_data[col * n_columns];
374 Number res0 = val0 * x[0];
375 for (int i = 1; i < mm; ++i)
376 {
377 if (contract_over_rows == true)
378 val0 = shape_data[i * n_columns + col];
379 else
380 val0 = shape_data[col * n_columns + i];
381 res0 += val0 * x[i];
382 }
383 if (add == false)
384 out[stride * col] = res0;
385 else
386 out[stride * col] += res0;
387 }
388
389 if (one_line == false)
390 {
391 ++in;
392 ++out;
393 }
394 }
395 if (one_line == false)
396 {
397 in += stride * (mm - 1);
398 out += stride * (nn - 1);
399 }
400 }
401 }
402
403
404
405 template <int dim,
406 int n_rows,
407 int n_columns,
408 typename Number,
409 typename Number2>
410 template <int face_direction,
411 bool contract_onto_face,
412 bool add,
413 int max_derivative,
414 bool lex_faces>
415 inline void
416 EvaluatorTensorProduct<evaluate_general,
417 dim,
418 n_rows,
419 n_columns,
420 Number,
421 Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
422 Number *DEAL_II_RESTRICT
423 out) const
424 {
425 Assert(dim > 0 && (lex_faces || dim < 4),
426 ExcMessage("Only dim=1,2,3 supported"));
427 static_assert(max_derivative >= 0 && max_derivative < 3,
428 "Only derivative orders 0-2 implemented");
429 Assert(shape_values != nullptr,
430 ExcMessage(
431 "The given array shape_values must not be the null pointer."));
432
433 constexpr int n_blocks1 =
434 lex_faces ? ::Utilities::pow<unsigned int>(n_rows, face_direction) :
435 (dim > 1 ? n_rows : 1);
436 constexpr int n_blocks2 =
437 lex_faces ? ::Utilities::pow<unsigned int>(
438 n_rows, std::max(dim - face_direction - 1, 0)) :
439 (dim > 2 ? n_rows : 1);
440
441 AssertIndexRange(face_direction, dim);
442 constexpr int stride = Utilities::pow(n_rows, face_direction);
443 constexpr int out_stride = Utilities::pow(n_rows, dim - 1);
444 const Number *DEAL_II_RESTRICT shape_values = this->shape_values;
445
446 for (int i2 = 0; i2 < n_blocks2; ++i2)
447 {
448 for (int i1 = 0; i1 < n_blocks1; ++i1)
449 {
450 if (contract_onto_face == true)
451 {
452 Number res0 = shape_values[0] * in[0];
453 Number res1, res2;
454 if (max_derivative > 0)
455 res1 = shape_values[n_rows] * in[0];
456 if (max_derivative > 1)
457 res2 = shape_values[2 * n_rows] * in[0];
458 for (int ind = 1; ind < n_rows; ++ind)
459 {
460 res0 += shape_values[ind] * in[stride * ind];
461 if (max_derivative > 0)
462 res1 += shape_values[ind + n_rows] * in[stride * ind];
463 if (max_derivative > 1)
464 res2 += shape_values[ind + 2 * n_rows] * in[stride * ind];
465 }
466 if (add == false)
467 {
468 out[0] = res0;
469 if (max_derivative > 0)
470 out[out_stride] = res1;
471 if (max_derivative > 1)
472 out[2 * out_stride] = res2;
473 }
474 else
475 {
476 out[0] += res0;
477 if (max_derivative > 0)
478 out[out_stride] += res1;
479 if (max_derivative > 1)
480 out[2 * out_stride] += res2;
481 }
482 }
483 else
484 {
485 for (int col = 0; col < n_rows; ++col)
486 {
487 if (add == false)
488 out[col * stride] = shape_values[col] * in[0];
489 else
490 out[col * stride] += shape_values[col] * in[0];
491 if (max_derivative > 0)
492 out[col * stride] +=
493 shape_values[col + n_rows] * in[out_stride];
494 if (max_derivative > 1)
495 out[col * stride] +=
496 shape_values[col + 2 * n_rows] * in[2 * out_stride];
497 }
498 }
499
500 if (lex_faces)
501 {
502 ++out;
503 ++in;
504 }
505 else
506 // increment: in regular case, just go to the next point in
507 // x-direction. If we are at the end of one chunk in x-dir, need
508 // to jump over to the next layer in z-direction
509 switch (face_direction)
510 {
511 case 0:
512 in += contract_onto_face ? n_rows : 1;
513 out += contract_onto_face ? 1 : n_rows;
514 break;
515 case 1:
516 ++in;
517 ++out;
518 // faces 2 and 3 in 3D use local coordinate system zx, which
519 // is the other way around compared to the tensor
520 // product. Need to take that into account.
521 if (dim == 3)
522 {
523 if (contract_onto_face)
524 out += n_rows - 1;
525 else
526 in += n_rows - 1;
527 }
528 break;
529 case 2:
530 ++in;
531 ++out;
532 break;
533 default:
534 Assert(false, ExcNotImplemented());
535 }
536 }
537 if (lex_faces)
538 {
539 if (contract_onto_face)
540 in += (::Utilities::pow(n_rows, face_direction + 1) -
541 n_blocks1);
542 else
543 out += (::Utilities::pow(n_rows, face_direction + 1) -
544 n_blocks1);
545 }
546 else if (face_direction == 1 && dim == 3)
547 {
548 // adjust for local coordinate system zx
549 if (contract_onto_face)
550 {
551 in += n_rows * (n_rows - 1);
552 out -= n_rows * n_rows - 1;
553 }
554 else
555 {
556 out += n_rows * (n_rows - 1);
557 in -= n_rows * n_rows - 1;
558 }
559 }
560 }
561 }
562
563
564
578 template <int dim, typename Number, typename Number2>
579 struct EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>
580 {
581 static constexpr unsigned int n_rows_of_product =
582 numbers::invalid_unsigned_int;
583 static constexpr unsigned int n_columns_of_product =
584 numbers::invalid_unsigned_int;
585
590 EvaluatorTensorProduct()
591 : shape_values(nullptr)
592 , shape_gradients(nullptr)
593 , shape_hessians(nullptr)
594 , n_rows(numbers::invalid_unsigned_int)
595 , n_columns(numbers::invalid_unsigned_int)
596 {}
597
601 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
602 const AlignedVector<Number2> &shape_gradients,
603 const AlignedVector<Number2> &shape_hessians,
604 const unsigned int n_rows,
605 const unsigned int n_columns)
606 : shape_values(shape_values.begin())
607 , shape_gradients(shape_gradients.begin())
608 , shape_hessians(shape_hessians.begin())
609 , n_rows(n_rows)
610 , n_columns(n_columns)
611 {
612 // We can enter this function either for the apply() path that has
613 // n_rows * n_columns entries or for the apply_face() path that only has
614 // n_rows * 3 entries in the array. Since we cannot decide about the use
615 // we must allow for both here.
616 Assert(shape_values.size() == 0 ||
617 shape_values.size() == n_rows * n_columns ||
618 shape_values.size() == n_rows * 3,
619 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
620 Assert(shape_gradients.size() == 0 ||
621 shape_gradients.size() == n_rows * n_columns,
622 ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
623 Assert(shape_hessians.size() == 0 ||
624 shape_hessians.size() == n_rows * n_columns,
625 ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
626 }
627
631 EvaluatorTensorProduct(const Number2 * shape_values,
632 const Number2 * shape_gradients,
633 const Number2 * shape_hessians,
634 const unsigned int n_rows,
635 const unsigned int n_columns)
636 : shape_values(shape_values)
637 , shape_gradients(shape_gradients)
638 , shape_hessians(shape_hessians)
639 , n_rows(n_rows)
640 , n_columns(n_columns)
641 {}
642
643 template <int direction, bool contract_over_rows, bool add>
644 void
645 values(const Number *in, Number *out) const
646 {
647 apply<direction, contract_over_rows, add>(shape_values, in, out);
648 }
649
650 template <int direction, bool contract_over_rows, bool add>
651 void
652 gradients(const Number *in, Number *out) const
653 {
654 apply<direction, contract_over_rows, add>(shape_gradients, in, out);
655 }
656
657 template <int direction, bool contract_over_rows, bool add>
658 void
659 hessians(const Number *in, Number *out) const
660 {
661 apply<direction, contract_over_rows, add>(shape_hessians, in, out);
662 }
663
664 template <int direction, bool contract_over_rows, bool add>
665 void
666 values_one_line(const Number in[], Number out[]) const
667 {
668 Assert(shape_values != nullptr, ExcNotInitialized());
669 apply<direction, contract_over_rows, add, true>(shape_values, in, out);
670 }
671
672 template <int direction, bool contract_over_rows, bool add>
673 void
674 gradients_one_line(const Number in[], Number out[]) const
675 {
676 Assert(shape_gradients != nullptr, ExcNotInitialized());
677 apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
678 }
679
680 template <int direction, bool contract_over_rows, bool add>
681 void
682 hessians_one_line(const Number in[], Number out[]) const
683 {
684 Assert(shape_hessians != nullptr, ExcNotInitialized());
685 apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
686 }
687
688 template <int direction,
689 bool contract_over_rows,
690 bool add,
691 bool one_line = false>
692 void
693 apply(const Number2 *DEAL_II_RESTRICT shape_data,
694 const Number * in,
695 Number * out) const;
696
697 template <int face_direction,
698 bool contract_onto_face,
699 bool add,
700 int max_derivative,
701 bool lex_faces = false>
702 void
703 apply_face(const Number *DEAL_II_RESTRICT in,
704 Number *DEAL_II_RESTRICT out) const;
705
706 const Number2 * shape_values;
707 const Number2 * shape_gradients;
708 const Number2 * shape_hessians;
709 const unsigned int n_rows;
710 const unsigned int n_columns;
711 };
712
713
714
715 template <int dim, typename Number, typename Number2>
716 template <int direction, bool contract_over_rows, bool add, bool one_line>
717 inline void
718 EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>::apply(
719 const Number2 *DEAL_II_RESTRICT shape_data,
720 const Number * in,
721 Number * out) const
722 {
723 static_assert(one_line == false || direction == dim - 1,
724 "Single-line evaluation only works for direction=dim-1.");
725 Assert(shape_data != nullptr,
726 ExcMessage(
727 "The given array shape_data must not be the null pointer!"));
728 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
729 in != out,
730 ExcMessage("In-place operation only supported for "
731 "n_rows==n_columns or single-line interpolation"));
732 AssertIndexRange(direction, dim);
733 const int mm = contract_over_rows ? n_rows : n_columns,
734 nn = contract_over_rows ? n_columns : n_rows;
735
736 const int stride =
737 direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
738 const int n_blocks1 = one_line ? 1 : stride;
739 const int n_blocks2 = direction >= dim - 1 ?
740 1 :
741 Utilities::fixed_power<dim - direction - 1>(n_rows);
742 Assert(n_rows <= 128, ExcNotImplemented());
743
744 // specialization for n_rows = 2 that manually unrolls the innermost loop
745 // to make the operation perform better (not completely as good as the
746 // templated one, but much better than the generic version down below,
747 // because the loop over col can be more effectively unrolled by the
748 // compiler)
749 if (contract_over_rows && n_rows == 2)
750 {
751 const Number2 *shape_data_1 = shape_data + n_columns;
752 for (int i2 = 0; i2 < n_blocks2; ++i2)
753 {
754 for (int i1 = 0; i1 < n_blocks1; ++i1)
755 {
756 const Number x0 = in[0], x1 = in[stride];
757 for (int col = 0; col < nn; ++col)
758 {
759 const Number result =
760 shape_data[col] * x0 + shape_data_1[col] * x1;
761 if (add == false)
762 out[stride * col] = result;
763 else
764 out[stride * col] += result;
765 }
766
767 if (one_line == false)
768 {
769 ++in;
770 ++out;
771 }
772 }
773 if (one_line == false)
774 {
775 in += stride * (mm - 1);
776 out += stride * (nn - 1);
777 }
778 }
779 }
780 // specialization for n = 3
781 else if (contract_over_rows && n_rows == 3)
782 {
783 const Number2 *shape_data_1 = shape_data + n_columns;
784 const Number2 *shape_data_2 = shape_data + 2 * n_columns;
785 for (int i2 = 0; i2 < n_blocks2; ++i2)
786 {
787 for (int i1 = 0; i1 < n_blocks1; ++i1)
788 {
789 const Number x0 = in[0], x1 = in[stride], x2 = in[2 * stride];
790 for (int col = 0; col < nn; ++col)
791 {
792 const Number result = shape_data[col] * x0 +
793 shape_data_1[col] * x1 +
794 shape_data_2[col] * x2;
795 if (add == false)
796 out[stride * col] = result;
797 else
798 out[stride * col] += result;
799 }
800
801 if (one_line == false)
802 {
803 ++in;
804 ++out;
805 }
806 }
807 if (one_line == false)
808 {
809 in += stride * (mm - 1);
810 out += stride * (nn - 1);
811 }
812 }
813 }
814 // general loop for all other cases
815 else
816 for (int i2 = 0; i2 < n_blocks2; ++i2)
817 {
818 for (int i1 = 0; i1 < n_blocks1; ++i1)
819 {
820 Number x[129];
821 for (int i = 0; i < mm; ++i)
822 x[i] = in[stride * i];
823 for (int col = 0; col < nn; ++col)
824 {
825 Number2 val0;
826 if (contract_over_rows == true)
827 val0 = shape_data[col];
828 else
829 val0 = shape_data[col * n_columns];
830 Number res0 = val0 * x[0];
831 for (int i = 1; i < mm; ++i)
832 {
833 if (contract_over_rows == true)
834 val0 = shape_data[i * n_columns + col];
835 else
836 val0 = shape_data[col * n_columns + i];
837 res0 += val0 * x[i];
838 }
839 if (add == false)
840 out[stride * col] = res0;
841 else
842 out[stride * col] += res0;
843 }
844
845 if (one_line == false)
846 {
847 ++in;
848 ++out;
849 }
850 }
851 if (one_line == false)
852 {
853 in += stride * (mm - 1);
854 out += stride * (nn - 1);
855 }
856 }
857 }
858
859
860
861 template <int dim, typename Number, typename Number2>
862 template <int face_direction,
863 bool contract_onto_face,
864 bool add,
865 int max_derivative,
866 bool lex_faces>
867 inline void
868 EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>::
869 apply_face(const Number *DEAL_II_RESTRICT in,
870 Number *DEAL_II_RESTRICT out) const
871 {
872 static_assert(lex_faces == false, "Not implemented yet.");
873
874 Assert(shape_values != nullptr,
875 ExcMessage(
876 "The given array shape_data must not be the null pointer!"));
877 static_assert(dim > 0 && dim < 4, "Only dim=1,2,3 supported");
878 const int n_blocks1 = dim > 1 ? n_rows : 1;
879 const int n_blocks2 = dim > 2 ? n_rows : 1;
880
881 AssertIndexRange(face_direction, dim);
882 const int stride =
883 face_direction > 0 ? Utilities::fixed_power<face_direction>(n_rows) : 1;
884 const int out_stride =
885 dim > 1 ? Utilities::fixed_power<dim - 1>(n_rows) : 1;
886
887 for (int i2 = 0; i2 < n_blocks2; ++i2)
888 {
889 for (int i1 = 0; i1 < n_blocks1; ++i1)
890 {
891 if (contract_onto_face == true)
892 {
893 Number res0 = shape_values[0] * in[0];
894 Number res1, res2;
895 if (max_derivative > 0)
896 res1 = shape_values[n_rows] * in[0];
897 if (max_derivative > 1)
898 res2 = shape_values[2 * n_rows] * in[0];
899 for (unsigned int ind = 1; ind < n_rows; ++ind)
900 {
901 res0 += shape_values[ind] * in[stride * ind];
902 if (max_derivative > 0)
903 res1 += shape_values[ind + n_rows] * in[stride * ind];
904 if (max_derivative > 1)
905 res2 += shape_values[ind + 2 * n_rows] * in[stride * ind];
906 }
907 if (add == false)
908 {
909 out[0] = res0;
910 if (max_derivative > 0)
911 out[out_stride] = res1;
912 if (max_derivative > 1)
913 out[2 * out_stride] = res2;
914 }
915 else
916 {
917 out[0] += res0;
918 if (max_derivative > 0)
919 out[out_stride] += res1;
920 if (max_derivative > 1)
921 out[2 * out_stride] += res2;
922 }
923 }
924 else
925 {
926 for (unsigned int col = 0; col < n_rows; ++col)
927 {
928 if (add == false)
929 out[col * stride] = shape_values[col] * in[0];
930 else
931 out[col * stride] += shape_values[col] * in[0];
932 if (max_derivative > 0)
933 out[col * stride] +=
934 shape_values[col + n_rows] * in[out_stride];
935 if (max_derivative > 1)
936 out[col * stride] +=
937 shape_values[col + 2 * n_rows] * in[2 * out_stride];
938 }
939 }
940
941 // increment: in regular case, just go to the next point in
942 // x-direction. If we are at the end of one chunk in x-dir, need
943 // to jump over to the next layer in z-direction
944 switch (face_direction)
945 {
946 case 0:
947 in += contract_onto_face ? n_rows : 1;
948 out += contract_onto_face ? 1 : n_rows;
949 break;
950 case 1:
951 ++in;
952 ++out;
953 // faces 2 and 3 in 3D use local coordinate system zx, which
954 // is the other way around compared to the tensor
955 // product. Need to take that into account.
956 if (dim == 3)
957 {
958 if (contract_onto_face)
959 out += n_rows - 1;
960 else
961 in += n_rows - 1;
962 }
963 break;
964 case 2:
965 ++in;
966 ++out;
967 break;
968 default:
969 Assert(false, ExcNotImplemented());
970 }
971 }
972 if (face_direction == 1 && dim == 3)
973 {
974 // adjust for local coordinate system zx
975 if (contract_onto_face)
976 {
977 in += n_rows * (n_rows - 1);
978 out -= n_rows * n_rows - 1;
979 }
980 else
981 {
982 out += n_rows * (n_rows - 1);
983 in -= n_rows * n_rows - 1;
984 }
985 }
986 }
987 }
988
989
990
1011 template <int dim,
1012 int n_rows,
1013 int n_columns,
1014 typename Number,
1015 typename Number2>
1016 struct EvaluatorTensorProduct<evaluate_symmetric,
1017 dim,
1018 n_rows,
1019 n_columns,
1020 Number,
1021 Number2>
1022 {
1023 static constexpr unsigned int n_rows_of_product =
1024 Utilities::pow(n_rows, dim);
1025 static constexpr unsigned int n_columns_of_product =
1026 Utilities::pow(n_columns, dim);
1027
1031 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
1032 const AlignedVector<Number2> &shape_gradients,
1033 const AlignedVector<Number2> &shape_hessians,
1034 const unsigned int dummy1 = 0,
1035 const unsigned int dummy2 = 0)
1036 : shape_values(shape_values.begin())
1037 , shape_gradients(shape_gradients.begin())
1038 , shape_hessians(shape_hessians.begin())
1039 {
1040 Assert(shape_values.size() == 0 ||
1041 shape_values.size() == n_rows * n_columns,
1042 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1043 Assert(shape_gradients.size() == 0 ||
1044 shape_gradients.size() == n_rows * n_columns,
1045 ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
1046 Assert(shape_hessians.size() == 0 ||
1047 shape_hessians.size() == n_rows * n_columns,
1048 ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
1049 (void)dummy1;
1050 (void)dummy2;
1051 }
1052
1053 template <int direction, bool contract_over_rows, bool add>
1054 void
1055 values(const Number in[], Number out[]) const;
1056
1057 template <int direction, bool contract_over_rows, bool add>
1058 void
1059 gradients(const Number in[], Number out[]) const;
1060
1061 template <int direction, bool contract_over_rows, bool add>
1062 void
1063 hessians(const Number in[], Number out[]) const;
1064
1065 const Number2 *shape_values;
1066 const Number2 *shape_gradients;
1067 const Number2 *shape_hessians;
1068 };
1069
1070
1071
1072 // In this case, the 1D shape values read (sorted lexicographically, rows
1073 // run over 1D dofs, columns over quadrature points):
1074 // Q2 --> [ 0.687 0 -0.087 ]
1075 // [ 0.4 1 0.4 ]
1076 // [-0.087 0 0.687 ]
1077 // Q3 --> [ 0.66 0.003 0.002 0.049 ]
1078 // [ 0.521 1.005 -0.01 -0.230 ]
1079 // [-0.230 -0.01 1.005 0.521 ]
1080 // [ 0.049 0.002 0.003 0.66 ]
1081 // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ]
1082 // [ 0.608 1.059 0 0.039 0.176 ]
1083 // [-0.409 -0.113 1 -0.113 -0.409 ]
1084 // [ 0.176 0.039 0 1.059 0.608 ]
1085 // [-0.032 -0.007 0 0.022 0.658 ]
1086 //
1087 // In these matrices, we want to use avoid computations involving zeros and
1088 // ones and in addition use the symmetry in entries to reduce the number of
1089 // read operations.
1090 template <int dim,
1091 int n_rows,
1092 int n_columns,
1093 typename Number,
1094 typename Number2>
1095 template <int direction, bool contract_over_rows, bool add>
1096 inline void
1097 EvaluatorTensorProduct<evaluate_symmetric,
1098 dim,
1099 n_rows,
1100 n_columns,
1101 Number,
1102 Number2>::values(const Number in[], Number out[]) const
1103 {
1104 Assert(shape_values != nullptr, ExcNotInitialized());
1105 AssertIndexRange(direction, dim);
1106 constexpr int mm = contract_over_rows ? n_rows : n_columns,
1107 nn = contract_over_rows ? n_columns : n_rows;
1108 constexpr int n_cols = nn / 2;
1109 constexpr int mid = mm / 2;
1110
1111 constexpr int stride = Utilities::pow(n_columns, direction);
1112 constexpr int n_blocks1 = stride;
1113 constexpr int n_blocks2 =
1114 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1115
1116 for (int i2 = 0; i2 < n_blocks2; ++i2)
1117 {
1118 for (int i1 = 0; i1 < n_blocks1; ++i1)
1119 {
1120 for (int col = 0; col < n_cols; ++col)
1121 {
1122 Number2 val0, val1;
1123 Number in0, in1, res0, res1;
1124 if (contract_over_rows == true)
1125 {
1126 val0 = shape_values[col];
1127 val1 = shape_values[nn - 1 - col];
1128 }
1129 else
1130 {
1131 val0 = shape_values[col * n_columns];
1132 val1 = shape_values[(col + 1) * n_columns - 1];
1133 }
1134 if (mid > 0)
1135 {
1136 in0 = in[0];
1137 in1 = in[stride * (mm - 1)];
1138 res0 = val0 * in0;
1139 res1 = val1 * in0;
1140 res0 += val1 * in1;
1141 res1 += val0 * in1;
1142 for (int ind = 1; ind < mid; ++ind)
1143 {
1144 if (contract_over_rows == true)
1145 {
1146 val0 = shape_values[ind * n_columns + col];
1147 val1 = shape_values[ind * n_columns + nn - 1 - col];
1148 }
1149 else
1150 {
1151 val0 = shape_values[col * n_columns + ind];
1152 val1 =
1153 shape_values[(col + 1) * n_columns - 1 - ind];
1154 }
1155 in0 = in[stride * ind];
1156 in1 = in[stride * (mm - 1 - ind)];
1157 res0 += val0 * in0;
1158 res1 += val1 * in0;
1159 res0 += val1 * in1;
1160 res1 += val0 * in1;
1161 }
1162 }
1163 else
1164 res0 = res1 = Number();
1165 if (contract_over_rows == true)
1166 {
1167 if (mm % 2 == 1)
1168 {
1169 val0 = shape_values[mid * n_columns + col];
1170 in1 = val0 * in[stride * mid];
1171 res0 += in1;
1172 res1 += in1;
1173 }
1174 }
1175 else
1176 {
1177 if (mm % 2 == 1 && nn % 2 == 0)
1178 {
1179 val0 = shape_values[col * n_columns + mid];
1180 in1 = val0 * in[stride * mid];
1181 res0 += in1;
1182 res1 += in1;
1183 }
1184 }
1185 if (add == false)
1186 {
1187 out[stride * col] = res0;
1188 out[stride * (nn - 1 - col)] = res1;
1189 }
1190 else
1191 {
1192 out[stride * col] += res0;
1193 out[stride * (nn - 1 - col)] += res1;
1194 }
1195 }
1196 if (contract_over_rows == true && nn % 2 == 1 && mm % 2 == 1)
1197 {
1198 if (add == false)
1199 out[stride * n_cols] = in[stride * mid];
1200 else
1201 out[stride * n_cols] += in[stride * mid];
1202 }
1203 else if (contract_over_rows == true && nn % 2 == 1)
1204 {
1205 Number res0;
1206 Number2 val0 = shape_values[n_cols];
1207 if (mid > 0)
1208 {
1209 res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1210 for (int ind = 1; ind < mid; ++ind)
1211 {
1212 val0 = shape_values[ind * n_columns + n_cols];
1213 res0 += val0 * (in[stride * ind] +
1214 in[stride * (mm - 1 - ind)]);
1215 }
1216 }
1217 else
1218 res0 = Number();
1219 if (add == false)
1220 out[stride * n_cols] = res0;
1221 else
1222 out[stride * n_cols] += res0;
1223 }
1224 else if (contract_over_rows == false && nn % 2 == 1)
1225 {
1226 Number res0;
1227 if (mid > 0)
1228 {
1229 Number2 val0 = shape_values[n_cols * n_columns];
1230 res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1231 for (int ind = 1; ind < mid; ++ind)
1232 {
1233 val0 = shape_values[n_cols * n_columns + ind];
1234 Number in1 = val0 * (in[stride * ind] +
1235 in[stride * (mm - 1 - ind)]);
1236 res0 += in1;
1237 }
1238 if (mm % 2)
1239 res0 += in[stride * mid];
1240 }
1241 else
1242 res0 = in[0];
1243 if (add == false)
1244 out[stride * n_cols] = res0;
1245 else
1246 out[stride * n_cols] += res0;
1247 }
1248
1249 ++in;
1250 ++out;
1251 }
1252 in += stride * (mm - 1);
1253 out += stride * (nn - 1);
1254 }
1255 }
1256
1257
1258
1259 // For the specialized loop used for the gradient computation in
1260 // here, the 1D shape values read (sorted lexicographically, rows
1261 // run over 1D dofs, columns over quadrature points):
1262 // Q2 --> [-2.549 -1 0.549 ]
1263 // [ 3.098 0 -3.098 ]
1264 // [-0.549 1 2.549 ]
1265 // Q3 --> [-4.315 -1.03 0.5 -0.44 ]
1266 // [ 6.07 -1.44 -2.97 2.196 ]
1267 // [-2.196 2.97 1.44 -6.07 ]
1268 // [ 0.44 -0.5 1.03 4.315 ]
1269 // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ]
1270 // [10.111 -2.76 -2.667 2.066 -2.306 ]
1271 // [-5.688 5.773 0 -5.773 5.688 ]
1272 // [ 2.306 -2.066 2.667 2.76 -10.111 ]
1273 // [-0.413 0.353 -0.333 -0.353 0.413 ]
1274 //
1275 // In these matrices, we want to use avoid computations involving
1276 // zeros and ones and in addition use the symmetry in entries to
1277 // reduce the number of read operations.
1278 template <int dim,
1279 int n_rows,
1280 int n_columns,
1281 typename Number,
1282 typename Number2>
1283 template <int direction, bool contract_over_rows, bool add>
1284 inline void
1285 EvaluatorTensorProduct<evaluate_symmetric,
1286 dim,
1287 n_rows,
1288 n_columns,
1289 Number,
1290 Number2>::gradients(const Number in[],
1291 Number out[]) const
1292 {
1293 Assert(shape_gradients != nullptr, ExcNotInitialized());
1294 AssertIndexRange(direction, dim);
1295 constexpr int mm = contract_over_rows ? n_rows : n_columns,
1296 nn = contract_over_rows ? n_columns : n_rows;
1297 constexpr int n_cols = nn / 2;
1298 constexpr int mid = mm / 2;
1299
1300 constexpr int stride = Utilities::pow(n_columns, direction);
1301 constexpr int n_blocks1 = stride;
1302 constexpr int n_blocks2 =
1303 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1304
1305 for (int i2 = 0; i2 < n_blocks2; ++i2)
1306 {
1307 for (int i1 = 0; i1 < n_blocks1; ++i1)
1308 {
1309 for (int col = 0; col < n_cols; ++col)
1310 {
1311 Number2 val0, val1;
1312 Number in0, in1, res0, res1;
1313 if (contract_over_rows == true)
1314 {
1315 val0 = shape_gradients[col];
1316 val1 = shape_gradients[nn - 1 - col];
1317 }
1318 else
1319 {
1320 val0 = shape_gradients[col * n_columns];
1321 val1 = shape_gradients[(nn - col - 1) * n_columns];
1322 }
1323 if (mid > 0)
1324 {
1325 in0 = in[0];
1326 in1 = in[stride * (mm - 1)];
1327 res0 = val0 * in0;
1328 res1 = val1 * in0;
1329 res0 -= val1 * in1;
1330 res1 -= val0 * in1;
1331 for (int ind = 1; ind < mid; ++ind)
1332 {
1333 if (contract_over_rows == true)
1334 {
1335 val0 = shape_gradients[ind * n_columns + col];
1336 val1 =
1337 shape_gradients[ind * n_columns + nn - 1 - col];
1338 }
1339 else
1340 {
1341 val0 = shape_gradients[col * n_columns + ind];
1342 val1 =
1343 shape_gradients[(nn - col - 1) * n_columns + ind];
1344 }
1345 in0 = in[stride * ind];
1346 in1 = in[stride * (mm - 1 - ind)];
1347 res0 += val0 * in0;
1348 res1 += val1 * in0;
1349 res0 -= val1 * in1;
1350 res1 -= val0 * in1;
1351 }
1352 }
1353 else
1354 res0 = res1 = Number();
1355 if (mm % 2 == 1)
1356 {
1357 if (contract_over_rows == true)
1358 val0 = shape_gradients[mid * n_columns + col];
1359 else
1360 val0 = shape_gradients[col * n_columns + mid];
1361 in1 = val0 * in[stride * mid];
1362 res0 += in1;
1363 res1 -= in1;
1364 }
1365 if (add == false)
1366 {
1367 out[stride * col] = res0;
1368 out[stride * (nn - 1 - col)] = res1;
1369 }
1370 else
1371 {
1372 out[stride * col] += res0;
1373 out[stride * (nn - 1 - col)] += res1;
1374 }
1375 }
1376 if (nn % 2 == 1)
1377 {
1378 Number2 val0;
1379 Number res0;
1380 if (contract_over_rows == true)
1381 val0 = shape_gradients[n_cols];
1382 else
1383 val0 = shape_gradients[n_cols * n_columns];
1384 res0 = val0 * (in[0] - in[stride * (mm - 1)]);
1385 for (int ind = 1; ind < mid; ++ind)
1386 {
1387 if (contract_over_rows == true)
1388 val0 = shape_gradients[ind * n_columns + n_cols];
1389 else
1390 val0 = shape_gradients[n_cols * n_columns + ind];
1391 Number in1 =
1392 val0 * (in[stride * ind] - in[stride * (mm - 1 - ind)]);
1393 res0 += in1;
1394 }
1395 if (add == false)
1396 out[stride * n_cols] = res0;
1397 else
1398 out[stride * n_cols] += res0;
1399 }
1400
1401 ++in;
1402 ++out;
1403 }
1404 in += stride * (mm - 1);
1405 out += stride * (nn - 1);
1406 }
1407 }
1408
1409
1410
1411 // evaluates the given shape data in 1d-3d using the tensor product
1412 // form assuming the symmetries of unit cell shape hessians for
1413 // finite elements in FEEvaluation
1414 template <int dim,
1415 int n_rows,
1416 int n_columns,
1417 typename Number,
1418 typename Number2>
1419 template <int direction, bool contract_over_rows, bool add>
1420 inline void
1421 EvaluatorTensorProduct<evaluate_symmetric,
1422 dim,
1423 n_rows,
1424 n_columns,
1425 Number,
1426 Number2>::hessians(const Number in[],
1427 Number out[]) const
1428 {
1429 Assert(shape_hessians != nullptr, ExcNotInitialized());
1430 AssertIndexRange(direction, dim);
1431 constexpr int mm = contract_over_rows ? n_rows : n_columns;
1432 constexpr int nn = contract_over_rows ? n_columns : n_rows;
1433 constexpr int n_cols = nn / 2;
1434 constexpr int mid = mm / 2;
1435
1436 constexpr int stride = Utilities::pow(n_columns, direction);
1437 constexpr int n_blocks1 = stride;
1438 constexpr int n_blocks2 =
1439 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1440
1441 for (int i2 = 0; i2 < n_blocks2; ++i2)
1442 {
1443 for (int i1 = 0; i1 < n_blocks1; ++i1)
1444 {
1445 for (int col = 0; col < n_cols; ++col)
1446 {
1447 Number2 val0, val1;
1448 Number in0, in1, res0, res1;
1449 if (contract_over_rows == true)
1450 {
1451 val0 = shape_hessians[col];
1452 val1 = shape_hessians[nn - 1 - col];
1453 }
1454 else
1455 {
1456 val0 = shape_hessians[col * n_columns];
1457 val1 = shape_hessians[(col + 1) * n_columns - 1];
1458 }
1459 if (mid > 0)
1460 {
1461 in0 = in[0];
1462 in1 = in[stride * (mm - 1)];
1463 res0 = val0 * in0;
1464 res1 = val1 * in0;
1465 res0 += val1 * in1;
1466 res1 += val0 * in1;
1467 for (int ind = 1; ind < mid; ++ind)
1468 {
1469 if (contract_over_rows == true)
1470 {
1471 val0 = shape_hessians[ind * n_columns + col];
1472 val1 =
1473 shape_hessians[ind * n_columns + nn - 1 - col];
1474 }
1475 else
1476 {
1477 val0 = shape_hessians[col * n_columns + ind];
1478 val1 =
1479 shape_hessians[(col + 1) * n_columns - 1 - ind];
1480 }
1481 in0 = in[stride * ind];
1482 in1 = in[stride * (mm - 1 - ind)];
1483 res0 += val0 * in0;
1484 res1 += val1 * in0;
1485 res0 += val1 * in1;
1486 res1 += val0 * in1;
1487 }
1488 }
1489 else
1490 res0 = res1 = Number();
1491 if (mm % 2 == 1)
1492 {
1493 if (contract_over_rows == true)
1494 val0 = shape_hessians[mid * n_columns + col];
1495 else
1496 val0 = shape_hessians[col * n_columns + mid];
1497 in1 = val0 * in[stride * mid];
1498 res0 += in1;
1499 res1 += in1;
1500 }
1501 if (add == false)
1502 {
1503 out[stride * col] = res0;
1504 out[stride * (nn - 1 - col)] = res1;
1505 }
1506 else
1507 {
1508 out[stride * col] += res0;
1509 out[stride * (nn - 1 - col)] += res1;
1510 }
1511 }
1512 if (nn % 2 == 1)
1513 {
1514 Number2 val0;
1515 Number res0;
1516 if (contract_over_rows == true)
1517 val0 = shape_hessians[n_cols];
1518 else
1519 val0 = shape_hessians[n_cols * n_columns];
1520 if (mid > 0)
1521 {
1522 res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1523 for (int ind = 1; ind < mid; ++ind)
1524 {
1525 if (contract_over_rows == true)
1526 val0 = shape_hessians[ind * n_columns + n_cols];
1527 else
1528 val0 = shape_hessians[n_cols * n_columns + ind];
1529 Number in1 = val0 * (in[stride * ind] +
1530 in[stride * (mm - 1 - ind)]);
1531 res0 += in1;
1532 }
1533 }
1534 else
1535 res0 = Number();
1536 if (mm % 2 == 1)
1537 {
1538 if (contract_over_rows == true)
1539 val0 = shape_hessians[mid * n_columns + n_cols];
1540 else
1541 val0 = shape_hessians[n_cols * n_columns + mid];
1542 res0 += val0 * in[stride * mid];
1543 }
1544 if (add == false)
1545 out[stride * n_cols] = res0;
1546 else
1547 out[stride * n_cols] += res0;
1548 }
1549
1550 ++in;
1551 ++out;
1552 }
1553 in += stride * (mm - 1);
1554 out += stride * (nn - 1);
1555 }
1556 }
1557
1558
1559
1591 template <int dim,
1592 int n_rows,
1593 int n_columns,
1594 typename Number,
1595 typename Number2>
1596 struct EvaluatorTensorProduct<evaluate_evenodd,
1597 dim,
1598 n_rows,
1599 n_columns,
1600 Number,
1601 Number2>
1602 {
1603 static constexpr unsigned int n_rows_of_product =
1604 Utilities::pow(n_rows, dim);
1605 static constexpr unsigned int n_columns_of_product =
1606 Utilities::pow(n_columns, dim);
1607
1613 EvaluatorTensorProduct()
1614 : shape_values(nullptr)
1615 , shape_gradients(nullptr)
1616 , shape_hessians(nullptr)
1617 {}
1618
1623 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values)
1624 : shape_values(shape_values.begin())
1625 , shape_gradients(nullptr)
1626 , shape_hessians(nullptr)
1627 {
1628 AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
1629 }
1630
1635 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
1636 const AlignedVector<Number2> &shape_gradients,
1637 const AlignedVector<Number2> &shape_hessians,
1638 const unsigned int dummy1 = 0,
1639 const unsigned int dummy2 = 0)
1640 : shape_values(shape_values.begin())
1641 , shape_gradients(shape_gradients.begin())
1642 , shape_hessians(shape_hessians.begin())
1643 {
1644 // In this function, we allow for dummy pointers if some of values,
1645 // gradients or hessians should not be computed
1646 if (!shape_values.empty())
1647 AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
1648 if (!shape_gradients.empty())
1649 AssertDimension(shape_gradients.size(), n_rows * ((n_columns + 1) / 2));
1650 if (!shape_hessians.empty())
1651 AssertDimension(shape_hessians.size(), n_rows * ((n_columns + 1) / 2));
1652 (void)dummy1;
1653 (void)dummy2;
1654 }
1655
1656 template <int direction, bool contract_over_rows, bool add>
1657 void
1658 values(const Number in[], Number out[]) const
1659 {
1660 Assert(shape_values != nullptr, ExcNotInitialized());
1661 apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
1662 }
1663
1664 template <int direction, bool contract_over_rows, bool add>
1665 void
1666 gradients(const Number in[], Number out[]) const
1667 {
1668 Assert(shape_gradients != nullptr, ExcNotInitialized());
1669 apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
1670 }
1671
1672 template <int direction, bool contract_over_rows, bool add>
1673 void
1674 hessians(const Number in[], Number out[]) const
1675 {
1676 Assert(shape_hessians != nullptr, ExcNotInitialized());
1677 apply<direction, contract_over_rows, add, 2>(shape_hessians, in, out);
1678 }
1679
1680 template <int direction, bool contract_over_rows, bool add>
1681 void
1682 values_one_line(const Number in[], Number out[]) const
1683 {
1684 Assert(shape_values != nullptr, ExcNotInitialized());
1685 apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
1686 }
1687
1688 template <int direction, bool contract_over_rows, bool add>
1689 void
1690 gradients_one_line(const Number in[], Number out[]) const
1691 {
1692 Assert(shape_gradients != nullptr, ExcNotInitialized());
1693 apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
1694 in,
1695 out);
1696 }
1697
1698 template <int direction, bool contract_over_rows, bool add>
1699 void
1700 hessians_one_line(const Number in[], Number out[]) const
1701 {
1702 Assert(shape_hessians != nullptr, ExcNotInitialized());
1703 apply<direction, contract_over_rows, add, 2, true>(shape_hessians,
1704 in,
1705 out);
1706 }
1707
1736 template <int direction,
1737 bool contract_over_rows,
1738 bool add,
1739 int type,
1740 bool one_line = false>
1741 static void
1742 apply(const Number2 *DEAL_II_RESTRICT shape_data,
1743 const Number * in,
1744 Number * out);
1745
1746 const Number2 *shape_values;
1747 const Number2 *shape_gradients;
1748 const Number2 *shape_hessians;
1749 };
1750
1751
1752
1753 template <int dim,
1754 int n_rows,
1755 int n_columns,
1756 typename Number,
1757 typename Number2>
1758 template <int direction,
1759 bool contract_over_rows,
1760 bool add,
1761 int type,
1762 bool one_line>
1763 inline void
1764 EvaluatorTensorProduct<evaluate_evenodd,
1765 dim,
1766 n_rows,
1767 n_columns,
1768 Number,
1769 Number2>::apply(const Number2 *DEAL_II_RESTRICT shapes,
1770 const Number * in,
1771 Number * out)
1772 {
1773 static_assert(type < 3, "Only three variants type=0,1,2 implemented");
1774 static_assert(one_line == false || direction == dim - 1,
1775 "Single-line evaluation only works for direction=dim-1.");
1776 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
1777 in != out,
1778 ExcMessage("In-place operation only supported for "
1779 "n_rows==n_columns or single-line interpolation"));
1780
1781 // We cannot statically assert that direction is less than dim, so must do
1782 // an additional dynamic check
1783 AssertIndexRange(direction, dim);
1784
1785 constexpr int nn = contract_over_rows ? n_columns : n_rows;
1786 constexpr int mm = contract_over_rows ? n_rows : n_columns;
1787 constexpr int n_cols = nn / 2;
1788 constexpr int mid = mm / 2;
1789
1790 constexpr int stride = Utilities::pow(n_columns, direction);
1791 constexpr int n_blocks1 = one_line ? 1 : stride;
1792 constexpr int n_blocks2 =
1793 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1794
1795 constexpr int offset = (n_columns + 1) / 2;
1796
1797 // this code may look very inefficient at first sight due to the many
1798 // different cases with if's at the innermost loop part, but all of the
1799 // conditionals can be evaluated at compile time because they are
1800 // templates, so the compiler should optimize everything away
1801 for (int i2 = 0; i2 < n_blocks2; ++i2)
1802 {
1803 for (int i1 = 0; i1 < n_blocks1; ++i1)
1804 {
1805 Number xp[mid > 0 ? mid : 1], xm[mid > 0 ? mid : 1];
1806 for (int i = 0; i < mid; ++i)
1807 {
1808 if (contract_over_rows == true && type == 1)
1809 {
1810 xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
1811 xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
1812 }
1813 else
1814 {
1815 xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
1816 xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
1817 }
1818 }
1819 Number xmid = in[stride * mid];
1820 for (int col = 0; col < n_cols; ++col)
1821 {
1822 Number r0, r1;
1823 if (mid > 0)
1824 {
1825 if (contract_over_rows == true)
1826 {
1827 r0 = shapes[col] * xp[0];
1828 r1 = shapes[(n_rows - 1) * offset + col] * xm[0];
1829 }
1830 else
1831 {
1832 r0 = shapes[col * offset] * xp[0];
1833 r1 = shapes[(n_rows - 1 - col) * offset] * xm[0];
1834 }
1835 for (int ind = 1; ind < mid; ++ind)
1836 {
1837 if (contract_over_rows == true)
1838 {
1839 r0 += shapes[ind * offset + col] * xp[ind];
1840 r1 += shapes[(n_rows - 1 - ind) * offset + col] *
1841 xm[ind];
1842 }
1843 else
1844 {
1845 r0 += shapes[col * offset + ind] * xp[ind];
1846 r1 += shapes[(n_rows - 1 - col) * offset + ind] *
1847 xm[ind];
1848 }
1849 }
1850 }
1851 else
1852 r0 = r1 = Number();
1853 if (mm % 2 == 1 && contract_over_rows == true)
1854 {
1855 if (type == 1)
1856 r1 += shapes[mid * offset + col] * xmid;
1857 else
1858 r0 += shapes[mid * offset + col] * xmid;
1859 }
1860 else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0 || mm == 3))
1861 r0 += shapes[col * offset + mid] * xmid;
1862
1863 if (add == false)
1864 {
1865 out[stride * col] = r0 + r1;
1866 if (type == 1 && contract_over_rows == false)
1867 out[stride * (nn - 1 - col)] = r1 - r0;
1868 else
1869 out[stride * (nn - 1 - col)] = r0 - r1;
1870 }
1871 else
1872 {
1873 out[stride * col] += r0 + r1;
1874 if (type == 1 && contract_over_rows == false)
1875 out[stride * (nn - 1 - col)] += r1 - r0;
1876 else
1877 out[stride * (nn - 1 - col)] += r0 - r1;
1878 }
1879 }
1880 if (type == 0 && contract_over_rows == true && nn % 2 == 1 &&
1881 mm % 2 == 1 && mm > 3)
1882 {
1883 if (add == false)
1884 out[stride * n_cols] = shapes[mid * offset + n_cols] * xmid;
1885 else
1886 out[stride * n_cols] += shapes[mid * offset + n_cols] * xmid;
1887 }
1888 else if (contract_over_rows == true && nn % 2 == 1)
1889 {
1890 Number r0;
1891 if (mid > 0)
1892 {
1893 r0 = shapes[n_cols] * xp[0];
1894 for (int ind = 1; ind < mid; ++ind)
1895 r0 += shapes[ind * offset + n_cols] * xp[ind];
1896 }
1897 else
1898 r0 = Number();
1899 if (type != 1 && mm % 2 == 1)
1900 r0 += shapes[mid * offset + n_cols] * xmid;
1901
1902 if (add == false)
1903 out[stride * n_cols] = r0;
1904 else
1905 out[stride * n_cols] += r0;
1906 }
1907 else if (contract_over_rows == false && nn % 2 == 1)
1908 {
1909 Number r0;
1910 if (mid > 0)
1911 {
1912 if (type == 1)
1913 {
1914 r0 = shapes[n_cols * offset] * xm[0];
1915 for (int ind = 1; ind < mid; ++ind)
1916 r0 += shapes[n_cols * offset + ind] * xm[ind];
1917 }
1918 else
1919 {
1920 r0 = shapes[n_cols * offset] * xp[0];
1921 for (int ind = 1; ind < mid; ++ind)
1922 r0 += shapes[n_cols * offset + ind] * xp[ind];
1923 }
1924 }
1925 else
1926 r0 = Number();
1927
1928 if ((type == 0 || type == 2) && mm % 2 == 1)
1929 r0 += shapes[n_cols * offset + mid] * xmid;
1930
1931 if (add == false)
1932 out[stride * n_cols] = r0;
1933 else
1934 out[stride * n_cols] += r0;
1935 }
1936 if (one_line == false)
1937 {
1938 in += 1;
1939 out += 1;
1940 }
1941 }
1942 if (one_line == false)
1943 {
1944 in += stride * (mm - 1);
1945 out += stride * (nn - 1);
1946 }
1947 }
1948 }
1949
1950
1951
1980 template <int dim,
1981 int n_rows,
1982 int n_columns,
1983 typename Number,
1984 typename Number2>
1985 struct EvaluatorTensorProduct<evaluate_symmetric_hierarchical,
1986 dim,
1987 n_rows,
1988 n_columns,
1989 Number,
1990 Number2>
1991 {
1992 static constexpr unsigned int n_rows_of_product =
1993 Utilities::pow(n_rows, dim);
1994 static constexpr unsigned int n_columns_of_product =
1995 Utilities::pow(n_columns, dim);
1996
2002 EvaluatorTensorProduct()
2003 : shape_values(nullptr)
2004 , shape_gradients(nullptr)
2005 , shape_hessians(nullptr)
2006 {}
2007
2012 EvaluatorTensorProduct(const AlignedVector<Number> &shape_values)
2013 : shape_values(shape_values.begin())
2014 , shape_gradients(nullptr)
2015 , shape_hessians(nullptr)
2016 {}
2017
2022 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
2023 const AlignedVector<Number2> &shape_gradients,
2024 const AlignedVector<Number2> &shape_hessians,
2025 const unsigned int dummy1 = 0,
2026 const unsigned int dummy2 = 0)
2027 : shape_values(shape_values.begin())
2028 , shape_gradients(shape_gradients.begin())
2029 , shape_hessians(shape_hessians.begin())
2030 {
2031 (void)dummy1;
2032 (void)dummy2;
2033 }
2034
2035 template <int direction, bool contract_over_rows, bool add>
2036 void
2037 values(const Number in[], Number out[]) const
2038 {
2039 Assert(shape_values != nullptr, ExcNotInitialized());
2040 apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
2041 }
2042
2043 template <int direction, bool contract_over_rows, bool add>
2044 void
2045 gradients(const Number in[], Number out[]) const
2046 {
2047 Assert(shape_gradients != nullptr, ExcNotInitialized());
2048 apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
2049 }
2050
2051 template <int direction, bool contract_over_rows, bool add>
2052 void
2053 hessians(const Number in[], Number out[]) const
2054 {
2055 Assert(shape_hessians != nullptr, ExcNotInitialized());
2056 apply<direction, contract_over_rows, add, 0>(shape_hessians, in, out);
2057 }
2058
2059 template <int direction, bool contract_over_rows, bool add>
2060 void
2061 values_one_line(const Number in[], Number out[]) const
2062 {
2063 Assert(shape_values != nullptr, ExcNotInitialized());
2064 apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
2065 }
2066
2067 template <int direction, bool contract_over_rows, bool add>
2068 void
2069 gradients_one_line(const Number in[], Number out[]) const
2070 {
2071 Assert(shape_gradients != nullptr, ExcNotInitialized());
2072 apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
2073 in,
2074 out);
2075 }
2076
2077 template <int direction, bool contract_over_rows, bool add>
2078 void
2079 hessians_one_line(const Number in[], Number out[]) const
2080 {
2081 Assert(shape_hessians != nullptr, ExcNotInitialized());
2082 apply<direction, contract_over_rows, add, 0, true>(shape_hessians,
2083 in,
2084 out);
2085 }
2086
2114 template <int direction,
2115 bool contract_over_rows,
2116 bool add,
2117 int type,
2118 bool one_line = false>
2119 static void
2120 apply(const Number2 *DEAL_II_RESTRICT shape_data,
2121 const Number * in,
2122 Number * out);
2123
2124 const Number2 *shape_values;
2125 const Number2 *shape_gradients;
2126 const Number2 *shape_hessians;
2127 };
2128
2129
2130
2131 template <int dim,
2132 int n_rows,
2133 int n_columns,
2134 typename Number,
2135 typename Number2>
2136 template <int direction,
2137 bool contract_over_rows,
2138 bool add,
2139 int type,
2140 bool one_line>
2141 inline void
2142 EvaluatorTensorProduct<evaluate_symmetric_hierarchical,
2143 dim,
2144 n_rows,
2145 n_columns,
2146 Number,
2147 Number2>::apply(const Number2 *DEAL_II_RESTRICT shapes,
2148 const Number * in,
2149 Number * out)
2150 {
2151 static_assert(one_line == false || direction == dim - 1,
2152 "Single-line evaluation only works for direction=dim-1.");
2153 static_assert(
2154 type == 0 || type == 1,
2155 "Only types 0 and 1 implemented for evaluate_symmetric_hierarchical.");
2156 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
2157 in != out,
2158 ExcMessage("In-place operation only supported for "
2159 "n_rows==n_columns or single-line interpolation"));
2160
2161 // We cannot statically assert that direction is less than dim, so must do
2162 // an additional dynamic check
2163 AssertIndexRange(direction, dim);
2164
2165 constexpr int nn = contract_over_rows ? n_columns : n_rows;
2166 constexpr int mm = contract_over_rows ? n_rows : n_columns;
2167 constexpr int n_cols = nn / 2;
2168 constexpr int mid = mm / 2;
2169
2170 constexpr int stride = Utilities::pow(n_columns, direction);
2171 constexpr int n_blocks1 = one_line ? 1 : stride;
2172 constexpr int n_blocks2 =
2173 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
2174
2175 // this code may look very inefficient at first sight due to the many
2176 // different cases with if's at the innermost loop part, but all of the
2177 // conditionals can be evaluated at compile time because they are
2178 // templates, so the compiler should optimize everything away
2179 for (int i2 = 0; i2 < n_blocks2; ++i2)
2180 {
2181 for (int i1 = 0; i1 < n_blocks1; ++i1)
2182 {
2183 if (contract_over_rows)
2184 {
2185 Number x[mm];
2186 for (unsigned int i = 0; i < mm; ++i)
2187 x[i] = in[stride * i];
2188 for (unsigned int col = 0; col < n_cols; ++col)
2189 {
2190 Number r0, r1;
2191 if (mid > 0)
2192 {
2193 r0 = shapes[col] * x[0];
2194 r1 = shapes[col + n_columns] * x[1];
2195 for (unsigned int ind = 1; ind < mid; ++ind)
2196 {
2197 r0 +=
2198 shapes[col + 2 * ind * n_columns] * x[2 * ind];
2199 r1 += shapes[col + (2 * ind + 1) * n_columns] *
2200 x[2 * ind + 1];
2201 }
2202 }
2203 else
2204 r0 = r1 = Number();
2205 if (mm % 2 == 1)
2206 r0 += shapes[col + (mm - 1) * n_columns] * x[mm - 1];
2207 if (add == false)
2208 {
2209 out[stride * col] = r0 + r1;
2210 if (type == 1)
2211 out[stride * (nn - 1 - col)] = r1 - r0;
2212 else
2213 out[stride * (nn - 1 - col)] = r0 - r1;
2214 }
2215 else
2216 {
2217 out[stride * col] += r0 + r1;
2218 if (type == 1)
2219 out[stride * (nn - 1 - col)] += r1 - r0;
2220 else
2221 out[stride * (nn - 1 - col)] += r0 - r1;
2222 }
2223 }
2224 if (nn % 2 == 1)
2225 {
2226 Number r0;
2227 const unsigned int shift = type == 1 ? 1 : 0;
2228 if (mid > 0)
2229 {
2230 r0 = shapes[n_cols + shift * n_columns] * x[shift];
2231 for (unsigned int ind = 1; ind < mid; ++ind)
2232 r0 += shapes[n_cols + (2 * ind + shift) * n_columns] *
2233 x[2 * ind + shift];
2234 }
2235 else
2236 r0 = 0;
2237 if (type != 1 && mm % 2 == 1)
2238 r0 += shapes[n_cols + (mm - 1) * n_columns] * x[mm - 1];
2239 if (add == false)
2240 out[stride * n_cols] = r0;
2241 else
2242 out[stride * n_cols] += r0;
2243 }
2244 }
2245 else
2246 {
2247 Number xp[mid + 1], xm[mid > 0 ? mid : 1];
2248 for (int i = 0; i < mid; ++i)
2249 if (type == 0)
2250 {
2251 xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
2252 xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
2253 }
2254 else
2255 {
2256 xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
2257 xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
2258 }
2259 if (mm % 2 == 1)
2260 xp[mid] = in[stride * mid];
2261 for (unsigned int col = 0; col < n_cols; ++col)
2262 {
2263 Number r0, r1;
2264 if (mid > 0)
2265 {
2266 r0 = shapes[2 * col * n_columns] * xp[0];
2267 r1 = shapes[(2 * col + 1) * n_columns] * xm[0];
2268 for (unsigned int ind = 1; ind < mid; ++ind)
2269 {
2270 r0 += shapes[2 * col * n_columns + ind] * xp[ind];
2271 r1 +=
2272 shapes[(2 * col + 1) * n_columns + ind] * xm[ind];
2273 }
2274 }
2275 else
2276 r0 = r1 = Number();
2277 if (mm % 2 == 1)
2278 {
2279 if (type == 1)
2280 r1 +=
2281 shapes[(2 * col + 1) * n_columns + mid] * xp[mid];
2282 else
2283 r0 += shapes[2 * col * n_columns + mid] * xp[mid];
2284 }
2285 if (add == false)
2286 {
2287 out[stride * (2 * col)] = r0;
2288 out[stride * (2 * col + 1)] = r1;
2289 }
2290 else
2291 {
2292 out[stride * (2 * col)] += r0;
2293 out[stride * (2 * col + 1)] += r1;
2294 }
2295 }
2296 if (nn % 2 == 1)
2297 {
2298 Number r0;
2299 if (mid > 0)
2300 {
2301 r0 = shapes[(nn - 1) * n_columns] * xp[0];
2302 for (unsigned int ind = 1; ind < mid; ++ind)
2303 r0 += shapes[(nn - 1) * n_columns + ind] * xp[ind];
2304 }
2305 else
2306 r0 = Number();
2307 if (mm % 2 == 1 && type == 0)
2308 r0 += shapes[(nn - 1) * n_columns + mid] * xp[mid];
2309 if (add == false)
2310 out[stride * (nn - 1)] = r0;
2311 else
2312 out[stride * (nn - 1)] += r0;
2313 }
2314 }
2315 if (one_line == false)
2316 {
2317 in += 1;
2318 out += 1;
2319 }
2320 }
2321 if (one_line == false)
2322 {
2323 in += stride * (mm - 1);
2324 out += stride * (nn - 1);
2325 }
2326 }
2327 }
2328
2329
2330
2337 template <typename Number, typename Number2>
2338 struct ProductTypeNoPoint
2339 {
2340 using type = typename ProductType<Number, Number2>::type;
2341 };
2342
2343 template <int dim, typename Number, typename Number2>
2344 struct ProductTypeNoPoint<Point<dim, Number>, Number2>
2345 {
2346 using type = typename ProductType<Tensor<1, dim, Number>, Number2>::type;
2347 };
2348
2349
2350
2385 template <int dim, typename Number, typename Number2>
2386 inline std::pair<
2387 typename ProductTypeNoPoint<Number, Number2>::type,
2388 Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>>
2389 evaluate_tensor_product_value_and_gradient(
2390 const std::vector<Polynomials::Polynomial<double>> &poly,
2391 const std::vector<Number> & values,
2392 const Point<dim, Number2> & p,
2393 const bool d_linear = false,
2394 const std::vector<unsigned int> & renumber = {})
2395 {
2396 static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
2397
2398 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
2399
2400 // use `int` type for this variable and the loops below to inform the
2401 // compiler that the loops below will never overflow, which allows it to
2402 // generate more optimized code for the variable loop bounds in the
2403 // present context
2404 const int n_shapes = poly.size();
2405 AssertDimension(Utilities::pow(n_shapes, dim), values.size());
2406 Assert(renumber.empty() || renumber.size() == values.size(),
2407 ExcDimensionMismatch(renumber.size(), values.size()));
2408
2409 // shortcut for linear interpolation to speed up evaluation
2410 if (d_linear)
2411 {
2412 AssertDimension(poly.size(), 2);
2413 for (unsigned int i = 0; i < renumber.size(); ++i)
2414 AssertDimension(renumber[i], i);
2415
2416 if (dim == 1)
2417 {
2418 Tensor<1, dim, Number3> derivative;
2419 derivative[0] = values[1] - values[0];
2420 return std::make_pair((1. - p[0]) * values[0] + p[0] * values[1],
2421 derivative);
2422 }
2423 else if (dim == 2)
2424 {
2425 const Number2 x0 = 1. - p[0], x1 = p[0];
2426 const Number3 tmp0 = x0 * values[0] + x1 * values[1];
2427 const Number3 tmp1 = x0 * values[2] + x1 * values[3];
2428 const Number3 mapped = (1. - p[1]) * tmp0 + p[1] * tmp1;
2429 Tensor<1, dim, Number3> derivative;
2430 derivative[0] = (1. - p[1]) * (values[1] - values[0]) +
2431 p[1] * (values[3] - values[2]);
2432 derivative[1] = tmp1 - tmp0;
2433 return std::make_pair(mapped, derivative);
2434 }
2435 else if (dim == 3)
2436 {
2437 const Number2 x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
2438 z0 = 1. - p[2], z1 = p[2];
2439 const Number3 tmp0 = x0 * values[0] + x1 * values[1];
2440 const Number3 tmp1 = x0 * values[2] + x1 * values[3];
2441 const Number3 tmpy0 = y0 * tmp0 + y1 * tmp1;
2442 const Number3 tmp2 = x0 * values[4] + x1 * values[5];
2443 const Number3 tmp3 = x0 * values[6] + x1 * values[7];
2444 const Number3 tmpy1 = y0 * tmp2 + y1 * tmp3;
2445 const Number3 mapped = z0 * tmpy0 + z1 * tmpy1;
2446 Tensor<1, dim, Number3> derivative;
2447 derivative[2] = tmpy1 - tmpy0;
2448 derivative[1] = z0 * (tmp1 - tmp0) + z1 * (tmp3 - tmp2);
2449 derivative[0] =
2450 z0 *
2451 (y0 * (values[1] - values[0]) + y1 * (values[3] - values[2])) +
2452 z1 *
2453 (y0 * (values[5] - values[4]) + y1 * (values[7] - values[6]));
2454 return std::make_pair(mapped, derivative);
2455 }
2456 }
2457
2458 AssertIndexRange(n_shapes, 200);
2459 std::array<Number2, 2 * dim * 200> shapes;
2460
2461 // Evaluate 1D polynomials and their derivatives
2462 for (unsigned int d = 0; d < dim; ++d)
2463 for (int i = 0; i < n_shapes; ++i)
2464 poly[i].value(p[d], 1, shapes.data() + 2 * (d * n_shapes + i));
2465
2466 // Go through the tensor product of shape functions and interpolate
2467 // with optimal algorithm
2468 std::pair<Number3, Tensor<1, dim, Number3>> result = {};
2469 for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
2470 {
2471 Number3 value_y = {}, deriv_x = {}, deriv_y = {};
2472 for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
2473 {
2474 // Interpolation + derivative x direction
2475 Number3 value = {}, deriv = {};
2476
2477 // Distinguish the inner loop based on whether we have a
2478 // renumbering or not
2479 if (renumber.empty())
2480 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
2481 {
2482 value += shapes[2 * i0] * values[i];
2483 deriv += shapes[2 * i0 + 1] * values[i];
2484 }
2485 else
2486 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
2487 {
2488 value += shapes[2 * i0] * values[renumber[i]];
2489 deriv += shapes[2 * i0 + 1] * values[renumber[i]];
2490 }
2491
2492 // Interpolation + derivative in y direction
2493 if (dim > 1)
2494 {
2495 value_y += value * shapes[2 * n_shapes + 2 * i1];
2496 deriv_x += deriv * shapes[2 * n_shapes + 2 * i1];
2497 deriv_y += value * shapes[2 * n_shapes + 2 * i1 + 1];
2498 }
2499 else
2500 {
2501 result.first = value;
2502 result.second[0] = deriv;
2503 }
2504 }
2505 if (dim == 3)
2506 {
2507 // Interpolation + derivative in z direction
2508 result.first += value_y * shapes[4 * n_shapes + 2 * i2];
2509 result.second[0] += deriv_x * shapes[4 * n_shapes + 2 * i2];
2510 result.second[1] += deriv_y * shapes[4 * n_shapes + 2 * i2];
2511 result.second[2] += value_y * shapes[4 * n_shapes + 2 * i2 + 1];
2512 }
2513 else if (dim == 2)
2514 {
2515 result.first = value_y;
2516 result.second[0] = deriv_x;
2517 result.second[1] = deriv_y;
2518 }
2519 }
2520
2521 return result;
2522 }
2523
2524
2525
2529 template <int dim, typename Number, typename Number2>
2530 inline void
2531 integrate_add_tensor_product_value_and_gradient(
2532 const std::vector<Polynomials::Polynomial<double>> &poly,
2533 const Number2 & value,
2534 const Tensor<1, dim, Number2> & gradient,
2535 const Point<dim, Number> & p,
2536 AlignedVector<Number2> & values,
2537 const std::vector<unsigned int> & renumber = {})
2538 {
2539 static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
2540
2541 // as in evaluate, use `int` type to produce better code in this context
2542 const int n_shapes = poly.size();
2543 AssertDimension(Utilities::pow(n_shapes, dim), values.size());
2544 Assert(renumber.empty() || renumber.size() == values.size(),
2545 ExcDimensionMismatch(renumber.size(), values.size()));
2546
2547 AssertIndexRange(n_shapes, 200);
2548 std::array<Number, 2 * dim * 200> shapes;
2549
2550 // Evaluate 1D polynomials and their derivatives
2551 for (unsigned int d = 0; d < dim; ++d)
2552 for (int i = 0; i < n_shapes; ++i)
2553 poly[i].value(p[d], 1, shapes.data() + 2 * (d * n_shapes + i));
2554
2555 // Implement the transpose of the function above
2556 for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
2557 {
2558 const Number2 test_value_z =
2559 dim > 2 ? (value * shapes[4 * n_shapes + 2 * i2] +
2560 gradient[2] * shapes[4 * n_shapes + 2 * i2 + 1]) :
2561 value;
2562 const Number2 test_grad_x =
2563 dim > 2 ? gradient[0] * shapes[4 * n_shapes + 2 * i2] : gradient[0];
2564 const Number2 test_grad_y =
2565 dim > 2 ? gradient[1] * shapes[4 * n_shapes + 2 * i2] :
2566 (dim > 1 ? gradient[1] : Number2());
2567 for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
2568 {
2569 const Number2 test_value_y =
2570 dim > 1 ? (test_value_z * shapes[2 * n_shapes + 2 * i1] +
2571 test_grad_y * shapes[2 * n_shapes + 2 * i1 + 1]) :
2572 test_value_z;
2573 const Number2 test_grad_xy =
2574 dim > 1 ? test_grad_x * shapes[2 * n_shapes + 2 * i1] :
2575 test_grad_x;
2576 if (renumber.empty())
2577 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
2578 values[i] += shapes[2 * i0] * test_value_y +
2579 shapes[2 * i0 + 1] * test_grad_xy;
2580 else
2581 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
2582 values[renumber[i]] += shapes[2 * i0] * test_value_y +
2583 shapes[2 * i0 + 1] * test_grad_xy;
2584 }
2585 }
2586 }
2587
2588
2589} // end of namespace internal
2590
2591
2592DEAL_II_NAMESPACE_CLOSE
2593
2594#endif
aligned_vector.h
AlignedVector::empty
bool empty() const
AlignedVector::size
size_type size() const
Point
Definition: point.h:111
Tensor
Definition: tensor.h:472
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
DEAL_II_RESTRICT
#define DEAL_II_RESTRICT
Definition: config.h:101
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1465
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
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)
GridTools::shift
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2022
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition: trilinos_sparse_matrix.cc:53
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:461
Utilities::fixed_power
T fixed_power(const T t)
Definition: utilities.h:1081
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:2531
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:2389
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:196
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5484
std::pow
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
Definition: vectorization.h:5422
ProductType::type
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
Definition: template_constraints.h:323
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:1635
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:1700
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values)
Definition: tensor_product_kernels.h:1623
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:1690
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:1746
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:1748
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:1747
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:1682
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:1658
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:1666
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:1613
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:1674
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:631
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::n_columns
const unsigned int n_columns
Definition: tensor_product_kernels.h:710
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:707
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:682
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:590
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::gradients
void gradients(const Number *in, Number *out) const
Definition: tensor_product_kernels.h:652
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:666
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:674
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:706
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::values
void values(const Number *in, Number *out) const
Definition: tensor_product_kernels.h:645
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:708
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:601
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::n_rows
const unsigned int n_rows
Definition: tensor_product_kernels.h:709
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::hessians
void hessians(const Number *in, Number *out) const
Definition: tensor_product_kernels.h:659
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:206
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:318
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:170
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:319
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:161
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:199
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:213
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:236
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:220
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:320
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:228
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:1066
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:1067
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:1031
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:1065
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:2126
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:2045
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:2037
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:2079
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number > &shape_values)
Definition: tensor_product_kernels.h:2012
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:2061
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:2022
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:2069
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:2053
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:2125
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:2124
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:2002
internal::ProductTypeNoPoint< Point< dim, Number >, Number2 >::type
typename ProductType< Tensor< 1, dim, Number >, Number2 >::type type
Definition: tensor_product_kernels.h:2346
internal::ProductTypeNoPoint::type
typename ProductType< Number, Number2 >::type type
Definition: tensor_product_kernels.h:2340
Generated by doxygen 1.9.4