Reference documentation for deal.II version 9.5.0
\(\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 - 2023 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 int normal_dir,
153 typename Number,
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
237 EvaluatorTensorProduct(const Number2 * shape_values,
238 const Number2 * shape_gradients,
239 const Number2 * shape_hessians,
240 const unsigned int dummy1 = 0,
241 const unsigned int dummy2 = 0)
242 : shape_values(shape_values)
243 , shape_gradients(shape_gradients)
244 , shape_hessians(shape_hessians)
245 {
246 (void)dummy1;
247 (void)dummy2;
248 }
249
250 template <int direction, bool contract_over_rows, bool add>
251 void
252 values(const Number in[], Number out[]) const
253 {
254 apply<direction, contract_over_rows, add>(shape_values, in, out);
255 }
256
257 template <int direction, bool contract_over_rows, bool add>
258 void
259 gradients(const Number in[], Number out[]) const
260 {
261 apply<direction, contract_over_rows, add>(shape_gradients, in, out);
262 }
263
264 template <int direction, bool contract_over_rows, bool add>
265 void
266 hessians(const Number in[], Number out[]) const
267 {
268 apply<direction, contract_over_rows, add>(shape_hessians, in, out);
269 }
270
271 template <int direction, bool contract_over_rows, bool add>
272 void
273 values_one_line(const Number in[], Number out[]) const
274 {
275 Assert(shape_values != nullptr, ExcNotInitialized());
276 apply<direction, contract_over_rows, add, true>(shape_values, in, out);
277 }
278
279 template <int direction, bool contract_over_rows, bool add>
280 void
281 gradients_one_line(const Number in[], Number out[]) const
282 {
283 Assert(shape_gradients != nullptr, ExcNotInitialized());
284 apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
285 }
286
287 template <int direction, bool contract_over_rows, bool add>
288 void
289 hessians_one_line(const Number in[], Number out[]) const
290 {
291 Assert(shape_hessians != nullptr, ExcNotInitialized());
292 apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
293 }
294
319 template <int direction,
320 bool contract_over_rows,
321 bool add,
322 bool one_line = false>
323 static void
324 apply(const Number2 *DEAL_II_RESTRICT shape_data,
325 const Number * in,
326 Number * out);
327
357 template <int face_direction,
358 bool contract_onto_face,
359 bool add,
360 int max_derivative>
361 void
362 apply_face(const Number *DEAL_II_RESTRICT in,
363 Number *DEAL_II_RESTRICT out) const;
364
365 private:
366 const Number2 *shape_values;
367 const Number2 *shape_gradients;
368 const Number2 *shape_hessians;
369 };
370
371
372
373 template <int dim,
374 int n_rows,
375 int n_columns,
376 typename Number,
377 typename Number2>
378 template <int direction, bool contract_over_rows, bool add, bool one_line>
379 inline void
380 EvaluatorTensorProduct<evaluate_general,
381 dim,
382 n_rows,
383 n_columns,
384 Number,
385 Number2>::apply(const Number2 *DEAL_II_RESTRICT
386 shape_data,
387 const Number *in,
388 Number * out)
389 {
390 static_assert(one_line == false || direction == dim - 1,
391 "Single-line evaluation only works for direction=dim-1.");
392 Assert(shape_data != nullptr,
393 ExcMessage(
394 "The given array shape_data must not be the null pointer!"));
395 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
396 in != out,
397 ExcMessage("In-place operation only supported for "
398 "n_rows==n_columns or single-line interpolation"));
399 AssertIndexRange(direction, dim);
400 constexpr int mm = contract_over_rows ? n_rows : n_columns,
401 nn = contract_over_rows ? n_columns : n_rows;
402
403 constexpr int stride = Utilities::pow(n_columns, direction);
404 constexpr int n_blocks1 = one_line ? 1 : stride;
405 constexpr int n_blocks2 =
406 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
407
408 for (int i2 = 0; i2 < n_blocks2; ++i2)
409 {
410 for (int i1 = 0; i1 < n_blocks1; ++i1)
411 {
412 Number x[mm];
413 for (int i = 0; i < mm; ++i)
414 x[i] = in[stride * i];
415 for (int col = 0; col < nn; ++col)
416 {
417 Number2 val0;
418 if (contract_over_rows == true)
419 val0 = shape_data[col];
420 else
421 val0 = shape_data[col * n_columns];
422 Number res0 = val0 * x[0];
423 for (int i = 1; i < mm; ++i)
424 {
425 if (contract_over_rows == true)
426 val0 = shape_data[i * n_columns + col];
427 else
428 val0 = shape_data[col * n_columns + i];
429 res0 += val0 * x[i];
430 }
431 if (add)
432 out[stride * col] += res0;
433 else
434 out[stride * col] = res0;
435 }
436
437 if (one_line == false)
438 {
439 ++in;
440 ++out;
441 }
442 }
443 if (one_line == false)
444 {
445 in += stride * (mm - 1);
446 out += stride * (nn - 1);
447 }
448 }
449 }
450
451
452
453 template <int dim,
454 int n_rows,
455 int n_columns,
456 typename Number,
457 typename Number2>
458 template <int face_direction,
459 bool contract_onto_face,
460 bool add,
461 int max_derivative>
462 inline void
463 EvaluatorTensorProduct<evaluate_general,
464 dim,
465 n_rows,
466 n_columns,
467 Number,
468 Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
469 Number *DEAL_II_RESTRICT
470 out) const
471 {
472 Assert(dim > 0, ExcMessage("Only dim=1,2,3 supported"));
473 static_assert(max_derivative >= 0 && max_derivative < 3,
474 "Only derivative orders 0-2 implemented");
475 Assert(shape_values != nullptr,
476 ExcMessage(
477 "The given array shape_values must not be the null pointer."));
478
479 constexpr int n_blocks1 = (dim > 1 ? n_rows : 1);
480 constexpr int n_blocks2 = (dim > 2 ? n_rows : 1);
481
482 AssertIndexRange(face_direction, dim);
483 constexpr int in_stride = Utilities::pow(n_rows, face_direction);
484 constexpr int out_stride = Utilities::pow(n_rows, dim - 1);
485 const Number2 *DEAL_II_RESTRICT shape_values = this->shape_values;
486
487 for (int i2 = 0; i2 < n_blocks2; ++i2)
488 {
489 for (int i1 = 0; i1 < n_blocks1; ++i1)
490 {
491 if (contract_onto_face == true)
492 {
493 Number res0 = shape_values[0] * in[0];
494 Number res1, res2;
495 if (max_derivative > 0)
496 res1 = shape_values[n_rows] * in[0];
497 if (max_derivative > 1)
498 res2 = shape_values[2 * n_rows] * in[0];
499 for (int ind = 1; ind < n_rows; ++ind)
500 {
501 res0 += shape_values[ind] * in[in_stride * ind];
502 if (max_derivative > 0)
503 res1 += shape_values[ind + n_rows] * in[in_stride * ind];
504 if (max_derivative > 1)
505 res2 +=
506 shape_values[ind + 2 * n_rows] * in[in_stride * ind];
507 }
508 if (add)
509 {
510 out[0] += res0;
511 if (max_derivative > 0)
512 out[out_stride] += res1;
513 if (max_derivative > 1)
514 out[2 * out_stride] += res2;
515 }
516 else
517 {
518 out[0] = res0;
519 if (max_derivative > 0)
520 out[out_stride] = res1;
521 if (max_derivative > 1)
522 out[2 * out_stride] = res2;
523 }
524 }
525 else
526 {
527 for (int col = 0; col < n_rows; ++col)
528 {
529 if (add)
530 out[col * in_stride] += shape_values[col] * in[0];
531 else
532 out[col * in_stride] = shape_values[col] * in[0];
533 if (max_derivative > 0)
534 out[col * in_stride] +=
535 shape_values[col + n_rows] * in[out_stride];
536 if (max_derivative > 1)
537 out[col * in_stride] +=
538 shape_values[col + 2 * n_rows] * in[2 * out_stride];
539 }
540 }
541
542 // increment: in regular case, just go to the next point in
543 // x-direction. If we are at the end of one chunk in x-dir, need
544 // to jump over to the next layer in z-direction
545 switch (face_direction)
546 {
547 case 0:
548 in += contract_onto_face ? n_rows : 1;
549 out += contract_onto_face ? 1 : n_rows;
550 break;
551 case 1:
552 ++in;
553 ++out;
554 // faces 2 and 3 in 3d use local coordinate system zx, which
555 // is the other way around compared to the tensor
556 // product. Need to take that into account.
557 if (dim == 3)
558 {
559 if (contract_onto_face)
560 out += n_rows - 1;
561 else
562 in += n_rows - 1;
563 }
564 break;
565 case 2:
566 ++in;
567 ++out;
568 break;
569 default:
570 Assert(false, ExcNotImplemented());
571 }
572 }
573
574 // adjust for local coordinate system zx
575 if (face_direction == 1 && dim == 3)
576 {
577 if (contract_onto_face)
578 {
579 in += n_rows * (n_rows - 1);
580 out -= n_rows * n_rows - 1;
581 }
582 else
583 {
584 out += n_rows * (n_rows - 1);
585 in -= n_rows * n_rows - 1;
586 }
587 }
588 }
589 }
590
591
592
606 template <int dim, typename Number, typename Number2>
607 struct EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>
608 {
609 static constexpr unsigned int n_rows_of_product =
610 numbers::invalid_unsigned_int;
611 static constexpr unsigned int n_columns_of_product =
612 numbers::invalid_unsigned_int;
613
618 EvaluatorTensorProduct()
619 : shape_values(nullptr)
620 , shape_gradients(nullptr)
621 , shape_hessians(nullptr)
622 , n_rows(numbers::invalid_unsigned_int)
623 , n_columns(numbers::invalid_unsigned_int)
624 {}
625
629 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
630 const AlignedVector<Number2> &shape_gradients,
631 const AlignedVector<Number2> &shape_hessians,
632 const unsigned int n_rows,
633 const unsigned int n_columns)
634 : shape_values(shape_values.begin())
635 , shape_gradients(shape_gradients.begin())
636 , shape_hessians(shape_hessians.begin())
637 , n_rows(n_rows)
638 , n_columns(n_columns)
639 {
640 // We can enter this function either for the apply() path that has
641 // n_rows * n_columns entries or for the apply_face() path that only has
642 // n_rows * 3 entries in the array. Since we cannot decide about the use
643 // we must allow for both here.
644 Assert(shape_values.size() == 0 ||
645 shape_values.size() == n_rows * n_columns ||
646 shape_values.size() == n_rows * 3,
647 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
648 Assert(shape_gradients.size() == 0 ||
649 shape_gradients.size() == n_rows * n_columns,
650 ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
651 Assert(shape_hessians.size() == 0 ||
652 shape_hessians.size() == n_rows * n_columns,
653 ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
654 }
655
659 EvaluatorTensorProduct(const Number2 * shape_values,
660 const Number2 * shape_gradients,
661 const Number2 * shape_hessians,
662 const unsigned int n_rows,
663 const unsigned int n_columns)
664 : shape_values(shape_values)
665 , shape_gradients(shape_gradients)
666 , shape_hessians(shape_hessians)
667 , n_rows(n_rows)
668 , n_columns(n_columns)
669 {}
670
671 template <int direction, bool contract_over_rows, bool add>
672 void
673 values(const Number *in, Number *out) const
674 {
675 apply<direction, contract_over_rows, add>(shape_values, in, out);
676 }
677
678 template <int direction, bool contract_over_rows, bool add>
679 void
680 gradients(const Number *in, Number *out) const
681 {
682 apply<direction, contract_over_rows, add>(shape_gradients, in, out);
683 }
684
685 template <int direction, bool contract_over_rows, bool add>
686 void
687 hessians(const Number *in, Number *out) const
688 {
689 apply<direction, contract_over_rows, add>(shape_hessians, in, out);
690 }
691
692 template <int direction, bool contract_over_rows, bool add>
693 void
694 values_one_line(const Number in[], Number out[]) const
695 {
696 Assert(shape_values != nullptr, ExcNotInitialized());
697 apply<direction, contract_over_rows, add, true>(shape_values, in, out);
698 }
699
700 template <int direction, bool contract_over_rows, bool add>
701 void
702 gradients_one_line(const Number in[], Number out[]) const
703 {
704 Assert(shape_gradients != nullptr, ExcNotInitialized());
705 apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
706 }
707
708 template <int direction, bool contract_over_rows, bool add>
709 void
710 hessians_one_line(const Number in[], Number out[]) const
711 {
712 Assert(shape_hessians != nullptr, ExcNotInitialized());
713 apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
714 }
715
716 template <int direction,
717 bool contract_over_rows,
718 bool add,
719 bool one_line = false>
720 void
721 apply(const Number2 *DEAL_II_RESTRICT shape_data,
722 const Number * in,
723 Number * out) const;
724
725 template <int face_direction,
726 bool contract_onto_face,
727 bool add,
728 int max_derivative>
729 void
730 apply_face(const Number *DEAL_II_RESTRICT in,
731 Number *DEAL_II_RESTRICT out) const;
732
733 const Number2 * shape_values;
734 const Number2 * shape_gradients;
735 const Number2 * shape_hessians;
736 const unsigned int n_rows;
737 const unsigned int n_columns;
738 };
739
740
741
742 template <int dim, typename Number, typename Number2>
743 template <int direction, bool contract_over_rows, bool add, bool one_line>
744 inline void
745 EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>::apply(
746 const Number2 *DEAL_II_RESTRICT shape_data,
747 const Number * in,
748 Number * out) const
749 {
750 static_assert(one_line == false || direction == dim - 1,
751 "Single-line evaluation only works for direction=dim-1.");
752 Assert(shape_data != nullptr,
753 ExcMessage(
754 "The given array shape_data must not be the null pointer!"));
755 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
756 in != out,
757 ExcMessage("In-place operation only supported for "
758 "n_rows==n_columns or single-line interpolation"));
759 AssertIndexRange(direction, dim);
760 const int mm = contract_over_rows ? n_rows : n_columns,
761 nn = contract_over_rows ? n_columns : n_rows;
762
763 const int stride =
764 direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
765 const int n_blocks1 = one_line ? 1 : stride;
766 const int n_blocks2 = direction >= dim - 1 ?
767 1 :
768 Utilities::fixed_power<dim - direction - 1>(n_rows);
769 Assert(n_rows <= 128, ExcNotImplemented());
770
771 // specialization for n_rows = 2 that manually unrolls the innermost loop
772 // to make the operation perform better (not completely as good as the
773 // templated one, but much better than the generic version down below,
774 // because the loop over col can be more effectively unrolled by the
775 // compiler)
776 if (contract_over_rows && n_rows == 2)
777 {
778 const Number2 *shape_data_1 = shape_data + n_columns;
779 for (int i2 = 0; i2 < n_blocks2; ++i2)
780 {
781 for (int i1 = 0; i1 < n_blocks1; ++i1)
782 {
783 const Number x0 = in[0], x1 = in[stride];
784 for (int col = 0; col < nn; ++col)
785 {
786 const Number result =
787 shape_data[col] * x0 + shape_data_1[col] * x1;
788 if (add)
789 out[stride * col] += result;
790 else
791 out[stride * col] = result;
792 }
793
794 if (one_line == false)
795 {
796 ++in;
797 ++out;
798 }
799 }
800 if (one_line == false)
801 {
802 in += stride * (mm - 1);
803 out += stride * (nn - 1);
804 }
805 }
806 }
807 // specialization for n = 3
808 else if (contract_over_rows && n_rows == 3)
809 {
810 const Number2 *shape_data_1 = shape_data + n_columns;
811 const Number2 *shape_data_2 = shape_data + 2 * n_columns;
812 for (int i2 = 0; i2 < n_blocks2; ++i2)
813 {
814 for (int i1 = 0; i1 < n_blocks1; ++i1)
815 {
816 const Number x0 = in[0], x1 = in[stride], x2 = in[2 * stride];
817 for (int col = 0; col < nn; ++col)
818 {
819 const Number result = shape_data[col] * x0 +
820 shape_data_1[col] * x1 +
821 shape_data_2[col] * x2;
822 if (add)
823 out[stride * col] += result;
824 else
825 out[stride * col] = result;
826 }
827
828 if (one_line == false)
829 {
830 ++in;
831 ++out;
832 }
833 }
834 if (one_line == false)
835 {
836 in += stride * (mm - 1);
837 out += stride * (nn - 1);
838 }
839 }
840 }
841 // general loop for all other cases
842 else
843 for (int i2 = 0; i2 < n_blocks2; ++i2)
844 {
845 for (int i1 = 0; i1 < n_blocks1; ++i1)
846 {
847 Number x[129];
848 for (int i = 0; i < mm; ++i)
849 x[i] = in[stride * i];
850 for (int col = 0; col < nn; ++col)
851 {
852 Number2 val0;
853 if (contract_over_rows == true)
854 val0 = shape_data[col];
855 else
856 val0 = shape_data[col * n_columns];
857 Number res0 = val0 * x[0];
858 for (int i = 1; i < mm; ++i)
859 {
860 if (contract_over_rows == true)
861 val0 = shape_data[i * n_columns + col];
862 else
863 val0 = shape_data[col * n_columns + i];
864 res0 += val0 * x[i];
865 }
866 if (add)
867 out[stride * col] += res0;
868 else
869 out[stride * col] = res0;
870 }
871
872 if (one_line == false)
873 {
874 ++in;
875 ++out;
876 }
877 }
878 if (one_line == false)
879 {
880 in += stride * (mm - 1);
881 out += stride * (nn - 1);
882 }
883 }
884 }
885
886
887
888 template <int dim, typename Number, typename Number2>
889 template <int face_direction,
890 bool contract_onto_face,
891 bool add,
892 int max_derivative>
893 inline void
894 EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>::
895 apply_face(const Number *DEAL_II_RESTRICT in,
896 Number *DEAL_II_RESTRICT out) const
897 {
898 Assert(shape_values != nullptr,
899 ExcMessage(
900 "The given array shape_data must not be the null pointer!"));
901 static_assert(dim > 0 && dim < 4, "Only dim=1,2,3 supported");
902 const int n_blocks1 = dim > 1 ? n_rows : 1;
903 const int n_blocks2 = dim > 2 ? n_rows : 1;
904
905 AssertIndexRange(face_direction, dim);
906 const int in_stride =
907 face_direction > 0 ? Utilities::fixed_power<face_direction>(n_rows) : 1;
908 const int out_stride =
909 dim > 1 ? Utilities::fixed_power<dim - 1>(n_rows) : 1;
910
911 for (int i2 = 0; i2 < n_blocks2; ++i2)
912 {
913 for (int i1 = 0; i1 < n_blocks1; ++i1)
914 {
915 if (contract_onto_face == true)
916 {
917 Number res0 = shape_values[0] * in[0];
918 Number res1, res2;
919 if (max_derivative > 0)
920 res1 = shape_values[n_rows] * in[0];
921 if (max_derivative > 1)
922 res2 = shape_values[2 * n_rows] * in[0];
923 for (unsigned int ind = 1; ind < n_rows; ++ind)
924 {
925 res0 += shape_values[ind] * in[in_stride * ind];
926 if (max_derivative > 0)
927 res1 += shape_values[ind + n_rows] * in[in_stride * ind];
928 if (max_derivative > 1)
929 res2 +=
930 shape_values[ind + 2 * n_rows] * in[in_stride * ind];
931 }
932 if (add)
933 {
934 out[0] += res0;
935 if (max_derivative > 0)
936 out[out_stride] += res1;
937 if (max_derivative > 1)
938 out[2 * out_stride] += res2;
939 }
940 else
941 {
942 out[0] = res0;
943 if (max_derivative > 0)
944 out[out_stride] = res1;
945 if (max_derivative > 1)
946 out[2 * out_stride] = res2;
947 }
948 }
949 else
950 {
951 for (unsigned int col = 0; col < n_rows; ++col)
952 {
953 if (add)
954 out[col * in_stride] += shape_values[col] * in[0];
955 else
956 out[col * in_stride] = shape_values[col] * in[0];
957 if (max_derivative > 0)
958 out[col * in_stride] +=
959 shape_values[col + n_rows] * in[out_stride];
960 if (max_derivative > 1)
961 out[col * in_stride] +=
962 shape_values[col + 2 * n_rows] * in[2 * out_stride];
963 }
964 }
965
966 // increment: in regular case, just go to the next point in
967 // x-direction. If we are at the end of one chunk in x-dir, need
968 // to jump over to the next layer in z-direction
969 switch (face_direction)
970 {
971 case 0:
972 in += contract_onto_face ? n_rows : 1;
973 out += contract_onto_face ? 1 : n_rows;
974 break;
975 case 1:
976 ++in;
977 ++out;
978 // faces 2 and 3 in 3d use local coordinate system zx, which
979 // is the other way around compared to the tensor
980 // product. Need to take that into account.
981 if (dim == 3)
982 {
983 if (contract_onto_face)
984 out += n_rows - 1;
985 else
986 in += n_rows - 1;
987 }
988 break;
989 case 2:
990 ++in;
991 ++out;
992 break;
993 default:
994 Assert(false, ExcNotImplemented());
995 }
996 }
997 if (face_direction == 1 && dim == 3)
998 {
999 // adjust for local coordinate system zx
1000 if (contract_onto_face)
1001 {
1002 in += n_rows * (n_rows - 1);
1003 out -= n_rows * n_rows - 1;
1004 }
1005 else
1006 {
1007 out += n_rows * (n_rows - 1);
1008 in -= n_rows * n_rows - 1;
1009 }
1010 }
1011 }
1012 }
1013
1014
1015
1036 template <int dim,
1037 int n_rows,
1038 int n_columns,
1039 typename Number,
1040 typename Number2>
1041 struct EvaluatorTensorProduct<evaluate_symmetric,
1042 dim,
1043 n_rows,
1044 n_columns,
1045 Number,
1046 Number2>
1047 {
1048 static constexpr unsigned int n_rows_of_product =
1049 Utilities::pow(n_rows, dim);
1050 static constexpr unsigned int n_columns_of_product =
1051 Utilities::pow(n_columns, dim);
1052
1056 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
1057 const AlignedVector<Number2> &shape_gradients,
1058 const AlignedVector<Number2> &shape_hessians,
1059 const unsigned int dummy1 = 0,
1060 const unsigned int dummy2 = 0)
1061 : shape_values(shape_values.begin())
1062 , shape_gradients(shape_gradients.begin())
1063 , shape_hessians(shape_hessians.begin())
1064 {
1065 Assert(shape_values.size() == 0 ||
1066 shape_values.size() == n_rows * n_columns,
1067 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1068 Assert(shape_gradients.size() == 0 ||
1069 shape_gradients.size() == n_rows * n_columns,
1070 ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
1071 Assert(shape_hessians.size() == 0 ||
1072 shape_hessians.size() == n_rows * n_columns,
1073 ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
1074 (void)dummy1;
1075 (void)dummy2;
1076 }
1077
1078 template <int direction, bool contract_over_rows, bool add>
1079 void
1080 values(const Number in[], Number out[]) const;
1081
1082 template <int direction, bool contract_over_rows, bool add>
1083 void
1084 gradients(const Number in[], Number out[]) const;
1085
1086 template <int direction, bool contract_over_rows, bool add>
1087 void
1088 hessians(const Number in[], Number out[]) const;
1089
1090 private:
1091 const Number2 *shape_values;
1092 const Number2 *shape_gradients;
1093 const Number2 *shape_hessians;
1094 };
1095
1096
1097
1098 // In this case, the 1d shape values read (sorted lexicographically, rows
1099 // run over 1d dofs, columns over quadrature points):
1100 // Q2 --> [ 0.687 0 -0.087 ]
1101 // [ 0.4 1 0.4 ]
1102 // [-0.087 0 0.687 ]
1103 // Q3 --> [ 0.66 0.003 0.002 0.049 ]
1104 // [ 0.521 1.005 -0.01 -0.230 ]
1105 // [-0.230 -0.01 1.005 0.521 ]
1106 // [ 0.049 0.002 0.003 0.66 ]
1107 // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ]
1108 // [ 0.608 1.059 0 0.039 0.176 ]
1109 // [-0.409 -0.113 1 -0.113 -0.409 ]
1110 // [ 0.176 0.039 0 1.059 0.608 ]
1111 // [-0.032 -0.007 0 0.022 0.658 ]
1112 //
1113 // In these matrices, we want to use avoid computations involving zeros and
1114 // ones and in addition use the symmetry in entries to reduce the number of
1115 // read operations.
1116 template <int dim,
1117 int n_rows,
1118 int n_columns,
1119 typename Number,
1120 typename Number2>
1121 template <int direction, bool contract_over_rows, bool add>
1122 inline void
1123 EvaluatorTensorProduct<evaluate_symmetric,
1124 dim,
1125 n_rows,
1126 n_columns,
1127 Number,
1128 Number2>::values(const Number in[], Number out[]) const
1129 {
1130 Assert(shape_values != nullptr, ExcNotInitialized());
1131 AssertIndexRange(direction, dim);
1132 constexpr int mm = contract_over_rows ? n_rows : n_columns,
1133 nn = contract_over_rows ? n_columns : n_rows;
1134 constexpr int n_cols = nn / 2;
1135 constexpr int mid = mm / 2;
1136
1137 constexpr int stride = Utilities::pow(n_columns, direction);
1138 constexpr int n_blocks1 = stride;
1139 constexpr int n_blocks2 =
1140 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1141
1142 for (int i2 = 0; i2 < n_blocks2; ++i2)
1143 {
1144 for (int i1 = 0; i1 < n_blocks1; ++i1)
1145 {
1146 for (int col = 0; col < n_cols; ++col)
1147 {
1148 Number2 val0, val1;
1149 Number in0, in1, res0, res1;
1150 if (contract_over_rows == true)
1151 {
1152 val0 = shape_values[col];
1153 val1 = shape_values[nn - 1 - col];
1154 }
1155 else
1156 {
1157 val0 = shape_values[col * n_columns];
1158 val1 = shape_values[(col + 1) * n_columns - 1];
1159 }
1160 if (mid > 0)
1161 {
1162 in0 = in[0];
1163 in1 = in[stride * (mm - 1)];
1164 res0 = val0 * in0;
1165 res1 = val1 * in0;
1166 res0 += val1 * in1;
1167 res1 += val0 * in1;
1168 for (int ind = 1; ind < mid; ++ind)
1169 {
1170 if (contract_over_rows == true)
1171 {
1172 val0 = shape_values[ind * n_columns + col];
1173 val1 = shape_values[ind * n_columns + nn - 1 - col];
1174 }
1175 else
1176 {
1177 val0 = shape_values[col * n_columns + ind];
1178 val1 =
1179 shape_values[(col + 1) * n_columns - 1 - ind];
1180 }
1181 in0 = in[stride * ind];
1182 in1 = in[stride * (mm - 1 - ind)];
1183 res0 += val0 * in0;
1184 res1 += val1 * in0;
1185 res0 += val1 * in1;
1186 res1 += val0 * in1;
1187 }
1188 }
1189 else
1190 res0 = res1 = Number();
1191 if (contract_over_rows == true)
1192 {
1193 if (mm % 2 == 1)
1194 {
1195 val0 = shape_values[mid * n_columns + col];
1196 in1 = val0 * in[stride * mid];
1197 res0 += in1;
1198 res1 += in1;
1199 }
1200 }
1201 else
1202 {
1203 if (mm % 2 == 1 && nn % 2 == 0)
1204 {
1205 val0 = shape_values[col * n_columns + mid];
1206 in1 = val0 * in[stride * mid];
1207 res0 += in1;
1208 res1 += in1;
1209 }
1210 }
1211 if (add)
1212 {
1213 out[stride * col] += res0;
1214 out[stride * (nn - 1 - col)] += res1;
1215 }
1216 else
1217 {
1218 out[stride * col] = res0;
1219 out[stride * (nn - 1 - col)] = res1;
1220 }
1221 }
1222 if (contract_over_rows == true && nn % 2 == 1 && mm % 2 == 1)
1223 {
1224 if (add)
1225 out[stride * n_cols] += in[stride * mid];
1226 else
1227 out[stride * n_cols] = in[stride * mid];
1228 }
1229 else if (contract_over_rows == true && nn % 2 == 1)
1230 {
1231 Number res0;
1232 Number2 val0 = shape_values[n_cols];
1233 if (mid > 0)
1234 {
1235 res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1236 for (int ind = 1; ind < mid; ++ind)
1237 {
1238 val0 = shape_values[ind * n_columns + n_cols];
1239 res0 += val0 * (in[stride * ind] +
1240 in[stride * (mm - 1 - ind)]);
1241 }
1242 }
1243 else
1244 res0 = Number();
1245 if (add)
1246 out[stride * n_cols] += res0;
1247 else
1248 out[stride * n_cols] = res0;
1249 }
1250 else if (contract_over_rows == false && nn % 2 == 1)
1251 {
1252 Number res0;
1253 if (mid > 0)
1254 {
1255 Number2 val0 = shape_values[n_cols * n_columns];
1256 res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1257 for (int ind = 1; ind < mid; ++ind)
1258 {
1259 val0 = shape_values[n_cols * n_columns + ind];
1260 Number in1 = val0 * (in[stride * ind] +
1261 in[stride * (mm - 1 - ind)]);
1262 res0 += in1;
1263 }
1264 if (mm % 2)
1265 res0 += in[stride * mid];
1266 }
1267 else
1268 res0 = in[0];
1269 if (add)
1270 out[stride * n_cols] += res0;
1271 else
1272 out[stride * n_cols] = res0;
1273 }
1274
1275 ++in;
1276 ++out;
1277 }
1278 in += stride * (mm - 1);
1279 out += stride * (nn - 1);
1280 }
1281 }
1282
1283
1284
1285 // For the specialized loop used for the gradient computation in
1286 // here, the 1d shape values read (sorted lexicographically, rows
1287 // run over 1d dofs, columns over quadrature points):
1288 // Q2 --> [-2.549 -1 0.549 ]
1289 // [ 3.098 0 -3.098 ]
1290 // [-0.549 1 2.549 ]
1291 // Q3 --> [-4.315 -1.03 0.5 -0.44 ]
1292 // [ 6.07 -1.44 -2.97 2.196 ]
1293 // [-2.196 2.97 1.44 -6.07 ]
1294 // [ 0.44 -0.5 1.03 4.315 ]
1295 // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ]
1296 // [10.111 -2.76 -2.667 2.066 -2.306 ]
1297 // [-5.688 5.773 0 -5.773 5.688 ]
1298 // [ 2.306 -2.066 2.667 2.76 -10.111 ]
1299 // [-0.413 0.353 -0.333 -0.353 0.413 ]
1300 //
1301 // In these matrices, we want to use avoid computations involving
1302 // zeros and ones and in addition use the symmetry in entries to
1303 // reduce the number of read operations.
1304 template <int dim,
1305 int n_rows,
1306 int n_columns,
1307 typename Number,
1308 typename Number2>
1309 template <int direction, bool contract_over_rows, bool add>
1310 inline void
1311 EvaluatorTensorProduct<evaluate_symmetric,
1312 dim,
1313 n_rows,
1314 n_columns,
1315 Number,
1316 Number2>::gradients(const Number in[],
1317 Number out[]) const
1318 {
1319 Assert(shape_gradients != nullptr, ExcNotInitialized());
1320 AssertIndexRange(direction, dim);
1321 constexpr int mm = contract_over_rows ? n_rows : n_columns,
1322 nn = contract_over_rows ? n_columns : n_rows;
1323 constexpr int n_cols = nn / 2;
1324 constexpr int mid = mm / 2;
1325
1326 constexpr int stride = Utilities::pow(n_columns, direction);
1327 constexpr int n_blocks1 = stride;
1328 constexpr int n_blocks2 =
1329 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1330
1331 for (int i2 = 0; i2 < n_blocks2; ++i2)
1332 {
1333 for (int i1 = 0; i1 < n_blocks1; ++i1)
1334 {
1335 for (int col = 0; col < n_cols; ++col)
1336 {
1337 Number2 val0, val1;
1338 Number in0, in1, res0, res1;
1339 if (contract_over_rows == true)
1340 {
1341 val0 = shape_gradients[col];
1342 val1 = shape_gradients[nn - 1 - col];
1343 }
1344 else
1345 {
1346 val0 = shape_gradients[col * n_columns];
1347 val1 = shape_gradients[(nn - col - 1) * n_columns];
1348 }
1349 if (mid > 0)
1350 {
1351 in0 = in[0];
1352 in1 = in[stride * (mm - 1)];
1353 res0 = val0 * in0;
1354 res1 = val1 * in0;
1355 res0 -= val1 * in1;
1356 res1 -= val0 * in1;
1357 for (int ind = 1; ind < mid; ++ind)
1358 {
1359 if (contract_over_rows == true)
1360 {
1361 val0 = shape_gradients[ind * n_columns + col];
1362 val1 =
1363 shape_gradients[ind * n_columns + nn - 1 - col];
1364 }
1365 else
1366 {
1367 val0 = shape_gradients[col * n_columns + ind];
1368 val1 =
1369 shape_gradients[(nn - col - 1) * n_columns + ind];
1370 }
1371 in0 = in[stride * ind];
1372 in1 = in[stride * (mm - 1 - ind)];
1373 res0 += val0 * in0;
1374 res1 += val1 * in0;
1375 res0 -= val1 * in1;
1376 res1 -= val0 * in1;
1377 }
1378 }
1379 else
1380 res0 = res1 = Number();
1381 if (mm % 2 == 1)
1382 {
1383 if (contract_over_rows == true)
1384 val0 = shape_gradients[mid * n_columns + col];
1385 else
1386 val0 = shape_gradients[col * n_columns + mid];
1387 in1 = val0 * in[stride * mid];
1388 res0 += in1;
1389 res1 -= in1;
1390 }
1391 if (add)
1392 {
1393 out[stride * col] += res0;
1394 out[stride * (nn - 1 - col)] += res1;
1395 }
1396 else
1397 {
1398 out[stride * col] = res0;
1399 out[stride * (nn - 1 - col)] = res1;
1400 }
1401 }
1402 if (nn % 2 == 1)
1403 {
1404 Number2 val0;
1405 Number res0;
1406 if (contract_over_rows == true)
1407 val0 = shape_gradients[n_cols];
1408 else
1409 val0 = shape_gradients[n_cols * n_columns];
1410 res0 = val0 * (in[0] - in[stride * (mm - 1)]);
1411 for (int ind = 1; ind < mid; ++ind)
1412 {
1413 if (contract_over_rows == true)
1414 val0 = shape_gradients[ind * n_columns + n_cols];
1415 else
1416 val0 = shape_gradients[n_cols * n_columns + ind];
1417 Number in1 =
1418 val0 * (in[stride * ind] - in[stride * (mm - 1 - ind)]);
1419 res0 += in1;
1420 }
1421 if (add)
1422 out[stride * n_cols] += res0;
1423 else
1424 out[stride * n_cols] = res0;
1425 }
1426
1427 ++in;
1428 ++out;
1429 }
1430 in += stride * (mm - 1);
1431 out += stride * (nn - 1);
1432 }
1433 }
1434
1435
1436
1437 // evaluates the given shape data in 1d-3d using the tensor product
1438 // form assuming the symmetries of unit cell shape hessians for
1439 // finite elements in FEEvaluation
1440 template <int dim,
1441 int n_rows,
1442 int n_columns,
1443 typename Number,
1444 typename Number2>
1445 template <int direction, bool contract_over_rows, bool add>
1446 inline void
1447 EvaluatorTensorProduct<evaluate_symmetric,
1448 dim,
1449 n_rows,
1450 n_columns,
1451 Number,
1452 Number2>::hessians(const Number in[],
1453 Number out[]) const
1454 {
1455 Assert(shape_hessians != nullptr, ExcNotInitialized());
1456 AssertIndexRange(direction, dim);
1457 constexpr int mm = contract_over_rows ? n_rows : n_columns;
1458 constexpr int nn = contract_over_rows ? n_columns : n_rows;
1459 constexpr int n_cols = nn / 2;
1460 constexpr int mid = mm / 2;
1461
1462 constexpr int stride = Utilities::pow(n_columns, direction);
1463 constexpr int n_blocks1 = stride;
1464 constexpr int n_blocks2 =
1465 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1466
1467 for (int i2 = 0; i2 < n_blocks2; ++i2)
1468 {
1469 for (int i1 = 0; i1 < n_blocks1; ++i1)
1470 {
1471 for (int col = 0; col < n_cols; ++col)
1472 {
1473 Number2 val0, val1;
1474 Number in0, in1, res0, res1;
1475 if (contract_over_rows == true)
1476 {
1477 val0 = shape_hessians[col];
1478 val1 = shape_hessians[nn - 1 - col];
1479 }
1480 else
1481 {
1482 val0 = shape_hessians[col * n_columns];
1483 val1 = shape_hessians[(col + 1) * n_columns - 1];
1484 }
1485 if (mid > 0)
1486 {
1487 in0 = in[0];
1488 in1 = in[stride * (mm - 1)];
1489 res0 = val0 * in0;
1490 res1 = val1 * in0;
1491 res0 += val1 * in1;
1492 res1 += val0 * in1;
1493 for (int ind = 1; ind < mid; ++ind)
1494 {
1495 if (contract_over_rows == true)
1496 {
1497 val0 = shape_hessians[ind * n_columns + col];
1498 val1 =
1499 shape_hessians[ind * n_columns + nn - 1 - col];
1500 }
1501 else
1502 {
1503 val0 = shape_hessians[col * n_columns + ind];
1504 val1 =
1505 shape_hessians[(col + 1) * n_columns - 1 - ind];
1506 }
1507 in0 = in[stride * ind];
1508 in1 = in[stride * (mm - 1 - ind)];
1509 res0 += val0 * in0;
1510 res1 += val1 * in0;
1511 res0 += val1 * in1;
1512 res1 += val0 * in1;
1513 }
1514 }
1515 else
1516 res0 = res1 = Number();
1517 if (mm % 2 == 1)
1518 {
1519 if (contract_over_rows == true)
1520 val0 = shape_hessians[mid * n_columns + col];
1521 else
1522 val0 = shape_hessians[col * n_columns + mid];
1523 in1 = val0 * in[stride * mid];
1524 res0 += in1;
1525 res1 += in1;
1526 }
1527 if (add)
1528 {
1529 out[stride * col] += res0;
1530 out[stride * (nn - 1 - col)] += res1;
1531 }
1532 else
1533 {
1534 out[stride * col] = res0;
1535 out[stride * (nn - 1 - col)] = res1;
1536 }
1537 }
1538 if (nn % 2 == 1)
1539 {
1540 Number2 val0;
1541 Number res0;
1542 if (contract_over_rows == true)
1543 val0 = shape_hessians[n_cols];
1544 else
1545 val0 = shape_hessians[n_cols * n_columns];
1546 if (mid > 0)
1547 {
1548 res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1549 for (int ind = 1; ind < mid; ++ind)
1550 {
1551 if (contract_over_rows == true)
1552 val0 = shape_hessians[ind * n_columns + n_cols];
1553 else
1554 val0 = shape_hessians[n_cols * n_columns + ind];
1555 Number in1 = val0 * (in[stride * ind] +
1556 in[stride * (mm - 1 - ind)]);
1557 res0 += in1;
1558 }
1559 }
1560 else
1561 res0 = Number();
1562 if (mm % 2 == 1)
1563 {
1564 if (contract_over_rows == true)
1565 val0 = shape_hessians[mid * n_columns + n_cols];
1566 else
1567 val0 = shape_hessians[n_cols * n_columns + mid];
1568 res0 += val0 * in[stride * mid];
1569 }
1570 if (add)
1571 out[stride * n_cols] += res0;
1572 else
1573 out[stride * n_cols] = res0;
1574 }
1575
1576 ++in;
1577 ++out;
1578 }
1579 in += stride * (mm - 1);
1580 out += stride * (nn - 1);
1581 }
1582 }
1583
1584
1585
1586 template <int dim,
1587 int n_rows_static,
1588 int n_columns_static,
1589 typename Number,
1590 typename Number2,
1591 int direction,
1592 bool contract_over_rows,
1593 bool add,
1594 int type,
1595 bool one_line>
1596 inline void
1597 even_odd_apply(const int n_rows_in,
1598 const int n_columns_in,
1599 const Number2 *DEAL_II_RESTRICT shapes,
1600 const Number * in,
1601 Number * out)
1602 {
1603 static_assert(type < 3, "Only three variants type=0,1,2 implemented");
1604 static_assert(one_line == false || direction == dim - 1,
1605 "Single-line evaluation only works for direction=dim-1.");
1606
1607 const int n_rows = n_rows_static == -1 ? n_rows_in : n_rows_static;
1608 const int n_columns =
1609 n_columns_static == -1 ? n_columns_in : n_columns_static;
1610
1611 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
1612 in != out,
1613 ExcMessage("In-place operation only supported for "
1614 "n_rows==n_columns or single-line interpolation"));
1615
1616 // We cannot statically assert that direction is less than dim, so must do
1617 // an additional dynamic check
1618 AssertIndexRange(direction, dim);
1619
1620 const int nn = contract_over_rows ? n_columns : n_rows;
1621 const int mm = contract_over_rows ? n_rows : n_columns;
1622 constexpr int mm_static =
1623 contract_over_rows ? n_rows_static : n_columns_static;
1624 const int n_cols = nn / 2;
1625 const int mid = mm / 2;
1626 constexpr int mid_static = mm_static / 2;
1627 constexpr int max_mid = 15; // for non-templated execution
1628
1629 Assert((n_rows_static != -1 && n_columns_static != -1) || mid <= max_mid,
1630 ExcNotImplemented());
1631
1632 const int stride = Utilities::pow(n_columns, direction);
1633 const int n_blocks1 = one_line ? 1 : stride;
1634 const int n_blocks2 =
1635 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1636
1637 const int offset = (n_columns + 1) / 2;
1638
1639 // this code may look very inefficient at first sight due to the many
1640 // different cases with if's at the innermost loop part, but all of the
1641 // conditionals can be evaluated at compile time because they are
1642 // templates, so the compiler should optimize everything away
1643 for (int i2 = 0; i2 < n_blocks2; ++i2)
1644 {
1645 for (int i1 = 0; i1 < n_blocks1; ++i1)
1646 {
1647 constexpr unsigned int mid_size =
1648 (n_rows_static == -1 || n_columns_static == -1) ?
1649 max_mid :
1650 (mid_static > 0 ? mid_static : 1);
1651 Number xp[mid_size], xm[mid_size];
1652 for (int i = 0; i < mid; ++i)
1653 {
1654 if (contract_over_rows == true && type == 1)
1655 {
1656 xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
1657 xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
1658 }
1659 else
1660 {
1661 xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
1662 xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
1663 }
1664 }
1665 Number xmid = in[stride * mid];
1666 for (int col = 0; col < n_cols; ++col)
1667 {
1668 Number r0, r1;
1669 if (mid > 0)
1670 {
1671 if (contract_over_rows == true)
1672 {
1673 r0 = shapes[col] * xp[0];
1674 r1 = shapes[(n_rows - 1) * offset + col] * xm[0];
1675 }
1676 else
1677 {
1678 r0 = shapes[col * offset] * xp[0];
1679 r1 = shapes[(n_rows - 1 - col) * offset] * xm[0];
1680 }
1681 for (int ind = 1; ind < mid; ++ind)
1682 {
1683 if (contract_over_rows == true)
1684 {
1685 r0 += shapes[ind * offset + col] * xp[ind];
1686 r1 += shapes[(n_rows - 1 - ind) * offset + col] *
1687 xm[ind];
1688 }
1689 else
1690 {
1691 r0 += shapes[col * offset + ind] * xp[ind];
1692 r1 += shapes[(n_rows - 1 - col) * offset + ind] *
1693 xm[ind];
1694 }
1695 }
1696 }
1697 else
1698 r0 = r1 = Number();
1699 if (mm % 2 == 1 && contract_over_rows == true)
1700 {
1701 if (type == 1)
1702 r1 += shapes[mid * offset + col] * xmid;
1703 else
1704 r0 += shapes[mid * offset + col] * xmid;
1705 }
1706 else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0 || mm == 3))
1707 r0 += shapes[col * offset + mid] * xmid;
1708
1709 if (add)
1710 {
1711 out[stride * col] += r0 + r1;
1712 if (type == 1 && contract_over_rows == false)
1713 out[stride * (nn - 1 - col)] += r1 - r0;
1714 else
1715 out[stride * (nn - 1 - col)] += r0 - r1;
1716 }
1717 else
1718 {
1719 out[stride * col] = r0 + r1;
1720 if (type == 1 && contract_over_rows == false)
1721 out[stride * (nn - 1 - col)] = r1 - r0;
1722 else
1723 out[stride * (nn - 1 - col)] = r0 - r1;
1724 }
1725 }
1726 if (type == 0 && contract_over_rows == true && nn % 2 == 1 &&
1727 mm % 2 == 1 && mm > 3)
1728 {
1729 if (add)
1730 out[stride * n_cols] += shapes[mid * offset + n_cols] * xmid;
1731 else
1732 out[stride * n_cols] = shapes[mid * offset + n_cols] * xmid;
1733 }
1734 else if (contract_over_rows == true && nn % 2 == 1)
1735 {
1736 Number r0;
1737 if (mid > 0)
1738 {
1739 r0 = shapes[n_cols] * xp[0];
1740 for (int ind = 1; ind < mid; ++ind)
1741 r0 += shapes[ind * offset + n_cols] * xp[ind];
1742 }
1743 else
1744 r0 = Number();
1745 if (type != 1 && mm % 2 == 1)
1746 r0 += shapes[mid * offset + n_cols] * xmid;
1747
1748 if (add)
1749 out[stride * n_cols] += r0;
1750 else
1751 out[stride * n_cols] = r0;
1752 }
1753 else if (contract_over_rows == false && nn % 2 == 1)
1754 {
1755 Number r0;
1756 if (mid > 0)
1757 {
1758 if (type == 1)
1759 {
1760 r0 = shapes[n_cols * offset] * xm[0];
1761 for (int ind = 1; ind < mid; ++ind)
1762 r0 += shapes[n_cols * offset + ind] * xm[ind];
1763 }
1764 else
1765 {
1766 r0 = shapes[n_cols * offset] * xp[0];
1767 for (int ind = 1; ind < mid; ++ind)
1768 r0 += shapes[n_cols * offset + ind] * xp[ind];
1769 }
1770 }
1771 else
1772 r0 = Number();
1773
1774 if ((type == 0 || type == 2) && mm % 2 == 1)
1775 r0 += shapes[n_cols * offset + mid] * xmid;
1776
1777 if (add)
1778 out[stride * n_cols] += r0;
1779 else
1780 out[stride * n_cols] = r0;
1781 }
1782 if (one_line == false)
1783 {
1784 in += 1;
1785 out += 1;
1786 }
1787 }
1788 if (one_line == false)
1789 {
1790 in += stride * (mm - 1);
1791 out += stride * (nn - 1);
1792 }
1793 }
1794 }
1795
1796
1797
1829 template <int dim,
1830 int n_rows,
1831 int n_columns,
1832 typename Number,
1833 typename Number2>
1834 struct EvaluatorTensorProduct<evaluate_evenodd,
1835 dim,
1836 n_rows,
1837 n_columns,
1838 Number,
1839 Number2>
1840 {
1841 static constexpr unsigned int n_rows_of_product =
1842 Utilities::pow(n_rows, dim);
1843 static constexpr unsigned int n_columns_of_product =
1844 Utilities::pow(n_columns, dim);
1845
1851 EvaluatorTensorProduct()
1852 : shape_values(nullptr)
1853 , shape_gradients(nullptr)
1854 , shape_hessians(nullptr)
1855 {}
1856
1861 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values)
1862 : shape_values(shape_values.begin())
1863 , shape_gradients(nullptr)
1864 , shape_hessians(nullptr)
1865 {
1866 AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
1867 }
1868
1873 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
1874 const AlignedVector<Number2> &shape_gradients,
1875 const AlignedVector<Number2> &shape_hessians,
1876 const unsigned int dummy1 = 0,
1877 const unsigned int dummy2 = 0)
1878 : shape_values(shape_values.begin())
1879 , shape_gradients(shape_gradients.begin())
1880 , shape_hessians(shape_hessians.begin())
1881 {
1882 // In this function, we allow for dummy pointers if some of values,
1883 // gradients or hessians should not be computed
1884 if (!shape_values.empty())
1885 AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
1886 if (!shape_gradients.empty())
1887 AssertDimension(shape_gradients.size(), n_rows * ((n_columns + 1) / 2));
1888 if (!shape_hessians.empty())
1889 AssertDimension(shape_hessians.size(), n_rows * ((n_columns + 1) / 2));
1890 (void)dummy1;
1891 (void)dummy2;
1892 }
1893
1894 template <int direction, bool contract_over_rows, bool add>
1895 void
1896 values(const Number in[], Number out[]) const
1897 {
1898 Assert(shape_values != nullptr, ExcNotInitialized());
1899 apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
1900 }
1901
1902 template <int direction, bool contract_over_rows, bool add>
1903 void
1904 gradients(const Number in[], Number out[]) const
1905 {
1906 Assert(shape_gradients != nullptr, ExcNotInitialized());
1907 apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
1908 }
1909
1910 template <int direction, bool contract_over_rows, bool add>
1911 void
1912 hessians(const Number in[], Number out[]) const
1913 {
1914 Assert(shape_hessians != nullptr, ExcNotInitialized());
1915 apply<direction, contract_over_rows, add, 2>(shape_hessians, in, out);
1916 }
1917
1918 template <int direction, bool contract_over_rows, bool add>
1919 void
1920 values_one_line(const Number in[], Number out[]) const
1921 {
1922 Assert(shape_values != nullptr, ExcNotInitialized());
1923 apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
1924 }
1925
1926 template <int direction, bool contract_over_rows, bool add>
1927 void
1928 gradients_one_line(const Number in[], Number out[]) const
1929 {
1930 Assert(shape_gradients != nullptr, ExcNotInitialized());
1931 apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
1932 in,
1933 out);
1934 }
1935
1936 template <int direction, bool contract_over_rows, bool add>
1937 void
1938 hessians_one_line(const Number in[], Number out[]) const
1939 {
1940 Assert(shape_hessians != nullptr, ExcNotInitialized());
1941 apply<direction, contract_over_rows, add, 2, true>(shape_hessians,
1942 in,
1943 out);
1944 }
1945
1974 template <int direction,
1975 bool contract_over_rows,
1976 bool add,
1977 int type,
1978 bool one_line = false>
1979 static void
1980 apply(const Number2 *DEAL_II_RESTRICT shape_data,
1981 const Number * in,
1982 Number * out)
1983 {
1984 even_odd_apply<dim,
1985 n_rows,
1986 n_columns,
1987 Number,
1988 Number2,
1989 direction,
1990 contract_over_rows,
1991 add,
1992 type,
1993 one_line>(n_rows, n_columns, shape_data, in, out);
1994 }
1995
1996 private:
1997 const Number2 *shape_values;
1998 const Number2 *shape_gradients;
1999 const Number2 *shape_hessians;
2000 };
2001
2002
2009 template <int dim, typename Number, typename Number2>
2010 struct EvaluatorTensorProduct<evaluate_evenodd, dim, 0, 0, Number, Number2>
2011 {
2012 EvaluatorTensorProduct()
2013 : shape_values(nullptr)
2014 , shape_gradients(nullptr)
2015 , shape_hessians(nullptr)
2016 , n_rows(numbers::invalid_unsigned_int)
2017 , n_columns(numbers::invalid_unsigned_int)
2018 {}
2019
2020 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
2021 const unsigned int n_rows = 0,
2022 const unsigned int n_columns = 0)
2023 : shape_values(shape_values.begin())
2024 , shape_gradients(nullptr)
2025 , shape_hessians(nullptr)
2026 , n_rows(n_rows)
2027 , n_columns(n_columns)
2028 {
2029 AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
2030 }
2031
2032 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
2033 const AlignedVector<Number2> &shape_gradients,
2034 const AlignedVector<Number2> &shape_hessians,
2035 const unsigned int n_rows = 0,
2036 const unsigned int n_columns = 0)
2037 : shape_values(shape_values.begin())
2038 , shape_gradients(shape_gradients.begin())
2039 , shape_hessians(shape_hessians.begin())
2040 , n_rows(n_rows)
2041 , n_columns(n_columns)
2042 {
2043 if (!shape_values.empty())
2044 AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
2045 if (!shape_gradients.empty())
2046 AssertDimension(shape_gradients.size(), n_rows * ((n_columns + 1) / 2));
2047 if (!shape_hessians.empty())
2048 AssertDimension(shape_hessians.size(), n_rows * ((n_columns + 1) / 2));
2049 }
2050
2051 template <int direction, bool contract_over_rows, bool add>
2052 void
2053 values(const Number in[], Number out[]) const
2054 {
2055 Assert(shape_values != nullptr, ExcNotInitialized());
2056 apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
2057 }
2058
2059 template <int direction, bool contract_over_rows, bool add>
2060 void
2061 gradients(const Number in[], Number out[]) const
2062 {
2063 Assert(shape_gradients != nullptr, ExcNotInitialized());
2064 apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
2065 }
2066
2067 template <int direction, bool contract_over_rows, bool add>
2068 void
2069 hessians(const Number in[], Number out[]) const
2070 {
2071 Assert(shape_hessians != nullptr, ExcNotInitialized());
2072 apply<direction, contract_over_rows, add, 2>(shape_hessians, in, out);
2073 }
2074
2075 template <int direction, bool contract_over_rows, bool add>
2076 void
2077 values_one_line(const Number in[], Number out[]) const
2078 {
2079 Assert(shape_values != nullptr, ExcNotInitialized());
2080 apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
2081 }
2082
2083 template <int direction, bool contract_over_rows, bool add>
2084 void
2085 gradients_one_line(const Number in[], Number out[]) const
2086 {
2087 Assert(shape_gradients != nullptr, ExcNotInitialized());
2088 apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
2089 in,
2090 out);
2091 }
2092
2093 template <int direction, bool contract_over_rows, bool add>
2094 void
2095 hessians_one_line(const Number in[], Number out[]) const
2096 {
2097 Assert(shape_hessians != nullptr, ExcNotInitialized());
2098 apply<direction, contract_over_rows, add, 2, true>(shape_hessians,
2099 in,
2100 out);
2101 }
2102
2103 template <int direction,
2104 bool contract_over_rows,
2105 bool add,
2106 int type,
2107 bool one_line = false>
2108 void
2109 apply(const Number2 *DEAL_II_RESTRICT shape_data,
2110 const Number * in,
2111 Number * out) const
2112 {
2113 even_odd_apply<dim,
2114 -1,
2115 -1,
2116 Number,
2117 Number2,
2118 direction,
2119 contract_over_rows,
2120 add,
2121 type,
2122 one_line>(n_rows, n_columns, shape_data, in, out);
2123 }
2124
2125 private:
2126 const Number2 * shape_values;
2127 const Number2 * shape_gradients;
2128 const Number2 * shape_hessians;
2129 const unsigned int n_rows;
2130 const unsigned int n_columns;
2131 };
2132
2133
2134
2163 template <int dim,
2164 int n_rows,
2165 int n_columns,
2166 typename Number,
2167 typename Number2>
2168 struct EvaluatorTensorProduct<evaluate_symmetric_hierarchical,
2169 dim,
2170 n_rows,
2171 n_columns,
2172 Number,
2173 Number2>
2174 {
2175 static constexpr unsigned int n_rows_of_product =
2176 Utilities::pow(n_rows, dim);
2177 static constexpr unsigned int n_columns_of_product =
2178 Utilities::pow(n_columns, dim);
2179
2185 EvaluatorTensorProduct()
2186 : shape_values(nullptr)
2187 , shape_gradients(nullptr)
2188 , shape_hessians(nullptr)
2189 {}
2190
2195 EvaluatorTensorProduct(const AlignedVector<Number> &shape_values)
2196 : shape_values(shape_values.begin())
2197 , shape_gradients(nullptr)
2198 , shape_hessians(nullptr)
2199 {}
2200
2205 EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
2206 const AlignedVector<Number2> &shape_gradients,
2207 const AlignedVector<Number2> &shape_hessians,
2208 const unsigned int dummy1 = 0,
2209 const unsigned int dummy2 = 0)
2210 : shape_values(shape_values.begin())
2211 , shape_gradients(shape_gradients.begin())
2212 , shape_hessians(shape_hessians.begin())
2213 {
2214 (void)dummy1;
2215 (void)dummy2;
2216 }
2217
2218 template <int direction, bool contract_over_rows, bool add>
2219 void
2220 values(const Number in[], Number out[]) const
2221 {
2222 Assert(shape_values != nullptr, ExcNotInitialized());
2223 apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
2224 }
2225
2226 template <int direction, bool contract_over_rows, bool add>
2227 void
2228 gradients(const Number in[], Number out[]) const
2229 {
2230 Assert(shape_gradients != nullptr, ExcNotInitialized());
2231 apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
2232 }
2233
2234 template <int direction, bool contract_over_rows, bool add>
2235 void
2236 hessians(const Number in[], Number out[]) const
2237 {
2238 Assert(shape_hessians != nullptr, ExcNotInitialized());
2239 apply<direction, contract_over_rows, add, 0>(shape_hessians, in, out);
2240 }
2241
2242 template <int direction, bool contract_over_rows, bool add>
2243 void
2244 values_one_line(const Number in[], Number out[]) const
2245 {
2246 Assert(shape_values != nullptr, ExcNotInitialized());
2247 apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
2248 }
2249
2250 template <int direction, bool contract_over_rows, bool add>
2251 void
2252 gradients_one_line(const Number in[], Number out[]) const
2253 {
2254 Assert(shape_gradients != nullptr, ExcNotInitialized());
2255 apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
2256 in,
2257 out);
2258 }
2259
2260 template <int direction, bool contract_over_rows, bool add>
2261 void
2262 hessians_one_line(const Number in[], Number out[]) const
2263 {
2264 Assert(shape_hessians != nullptr, ExcNotInitialized());
2265 apply<direction, contract_over_rows, add, 0, true>(shape_hessians,
2266 in,
2267 out);
2268 }
2269
2297 template <int direction,
2298 bool contract_over_rows,
2299 bool add,
2300 int type,
2301 bool one_line = false>
2302 static void
2303 apply(const Number2 *DEAL_II_RESTRICT shape_data,
2304 const Number * in,
2305 Number * out);
2306
2307 private:
2308 const Number2 *shape_values;
2309 const Number2 *shape_gradients;
2310 const Number2 *shape_hessians;
2311 };
2312
2313
2314
2315 template <int dim,
2316 int n_rows,
2317 int n_columns,
2318 typename Number,
2319 typename Number2>
2320 template <int direction,
2321 bool contract_over_rows,
2322 bool add,
2323 int type,
2324 bool one_line>
2325 inline void
2326 EvaluatorTensorProduct<evaluate_symmetric_hierarchical,
2327 dim,
2328 n_rows,
2329 n_columns,
2330 Number,
2331 Number2>::apply(const Number2 *DEAL_II_RESTRICT shapes,
2332 const Number * in,
2333 Number * out)
2334 {
2335 static_assert(one_line == false || direction == dim - 1,
2336 "Single-line evaluation only works for direction=dim-1.");
2337 static_assert(
2338 type == 0 || type == 1,
2339 "Only types 0 and 1 implemented for evaluate_symmetric_hierarchical.");
2340 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
2341 in != out,
2342 ExcMessage("In-place operation only supported for "
2343 "n_rows==n_columns or single-line interpolation"));
2344
2345 // We cannot statically assert that direction is less than dim, so must do
2346 // an additional dynamic check
2347 AssertIndexRange(direction, dim);
2348
2349 constexpr int nn = contract_over_rows ? n_columns : n_rows;
2350 constexpr int mm = contract_over_rows ? n_rows : n_columns;
2351 constexpr int n_cols = nn / 2;
2352 constexpr int mid = mm / 2;
2353
2354 constexpr int stride = Utilities::pow(n_columns, direction);
2355 constexpr int n_blocks1 = one_line ? 1 : stride;
2356 constexpr int n_blocks2 =
2357 Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
2358
2359 // this code may look very inefficient at first sight due to the many
2360 // different cases with if's at the innermost loop part, but all of the
2361 // conditionals can be evaluated at compile time because they are
2362 // templates, so the compiler should optimize everything away
2363 for (int i2 = 0; i2 < n_blocks2; ++i2)
2364 {
2365 for (int i1 = 0; i1 < n_blocks1; ++i1)
2366 {
2367 if (contract_over_rows)
2368 {
2369 Number x[mm];
2370 for (unsigned int i = 0; i < mm; ++i)
2371 x[i] = in[stride * i];
2372 for (unsigned int col = 0; col < n_cols; ++col)
2373 {
2374 Number r0, r1;
2375 if (mid > 0)
2376 {
2377 r0 = shapes[col] * x[0];
2378 r1 = shapes[col + n_columns] * x[1];
2379 for (unsigned int ind = 1; ind < mid; ++ind)
2380 {
2381 r0 +=
2382 shapes[col + 2 * ind * n_columns] * x[2 * ind];
2383 r1 += shapes[col + (2 * ind + 1) * n_columns] *
2384 x[2 * ind + 1];
2385 }
2386 }
2387 else
2388 r0 = r1 = Number();
2389 if (mm % 2 == 1)
2390 r0 += shapes[col + (mm - 1) * n_columns] * x[mm - 1];
2391 if (add)
2392 {
2393 out[stride * col] += r0 + r1;
2394 if (type == 1)
2395 out[stride * (nn - 1 - col)] += r1 - r0;
2396 else
2397 out[stride * (nn - 1 - col)] += r0 - r1;
2398 }
2399 else
2400 {
2401 out[stride * col] = r0 + r1;
2402 if (type == 1)
2403 out[stride * (nn - 1 - col)] = r1 - r0;
2404 else
2405 out[stride * (nn - 1 - col)] = r0 - r1;
2406 }
2407 }
2408 if (nn % 2 == 1)
2409 {
2410 Number r0;
2411 const unsigned int shift = type == 1 ? 1 : 0;
2412 if (mid > 0)
2413 {
2414 r0 = shapes[n_cols + shift * n_columns] * x[shift];
2415 for (unsigned int ind = 1; ind < mid; ++ind)
2416 r0 += shapes[n_cols + (2 * ind + shift) * n_columns] *
2417 x[2 * ind + shift];
2418 }
2419 else
2420 r0 = 0;
2421 if (type != 1 && mm % 2 == 1)
2422 r0 += shapes[n_cols + (mm - 1) * n_columns] * x[mm - 1];
2423 if (add)
2424 out[stride * n_cols] += r0;
2425 else
2426 out[stride * n_cols] = r0;
2427 }
2428 }
2429 else
2430 {
2431 Number xp[mid + 1], xm[mid > 0 ? mid : 1];
2432 for (int i = 0; i < mid; ++i)
2433 if (type == 0)
2434 {
2435 xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
2436 xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
2437 }
2438 else
2439 {
2440 xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
2441 xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
2442 }
2443 if (mm % 2 == 1)
2444 xp[mid] = in[stride * mid];
2445 for (unsigned int col = 0; col < n_cols; ++col)
2446 {
2447 Number r0, r1;
2448 if (mid > 0)
2449 {
2450 r0 = shapes[2 * col * n_columns] * xp[0];
2451 r1 = shapes[(2 * col + 1) * n_columns] * xm[0];
2452 for (unsigned int ind = 1; ind < mid; ++ind)
2453 {
2454 r0 += shapes[2 * col * n_columns + ind] * xp[ind];
2455 r1 +=
2456 shapes[(2 * col + 1) * n_columns + ind] * xm[ind];
2457 }
2458 }
2459 else
2460 r0 = r1 = Number();
2461 if (mm % 2 == 1)
2462 {
2463 if (type == 1)
2464 r1 +=
2465 shapes[(2 * col + 1) * n_columns + mid] * xp[mid];
2466 else
2467 r0 += shapes[2 * col * n_columns + mid] * xp[mid];
2468 }
2469 if (add)
2470 {
2471 out[stride * (2 * col)] += r0;
2472 out[stride * (2 * col + 1)] += r1;
2473 }
2474 else
2475 {
2476 out[stride * (2 * col)] = r0;
2477 out[stride * (2 * col + 1)] = r1;
2478 }
2479 }
2480 if (nn % 2 == 1)
2481 {
2482 Number r0;
2483 if (mid > 0)
2484 {
2485 r0 = shapes[(nn - 1) * n_columns] * xp[0];
2486 for (unsigned int ind = 1; ind < mid; ++ind)
2487 r0 += shapes[(nn - 1) * n_columns + ind] * xp[ind];
2488 }
2489 else
2490 r0 = Number();
2491 if (mm % 2 == 1 && type == 0)
2492 r0 += shapes[(nn - 1) * n_columns + mid] * xp[mid];
2493 if (add)
2494 out[stride * (nn - 1)] += r0;
2495 else
2496 out[stride * (nn - 1)] = r0;
2497 }
2498 }
2499 if (one_line == false)
2500 {
2501 in += 1;
2502 out += 1;
2503 }
2504 }
2505 if (one_line == false)
2506 {
2507 in += stride * (mm - 1);
2508 out += stride * (nn - 1);
2509 }
2510 }
2511 }
2512
2513
2514
2534 template <int dim,
2535 int n_rows,
2536 int n_columns,
2537 int normal_dir,
2538 typename Number,
2539 typename Number2>
2540 struct EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
2541 dim,
2542 n_rows,
2543 n_columns,
2544 normal_dir,
2545 Number,
2546 Number2>
2547 {
2548 static constexpr unsigned int n_rows_of_product =
2549 numbers::invalid_unsigned_int;
2550 static constexpr unsigned int n_columns_of_product =
2551 numbers::invalid_unsigned_int;
2552
2557 EvaluatorTensorProductAnisotropic()
2558 : shape_values(nullptr)
2559 , shape_gradients(nullptr)
2560 , shape_hessians(nullptr)
2561 {}
2562
2566 EvaluatorTensorProductAnisotropic(
2567 const AlignedVector<Number2> &shape_values,
2568 const AlignedVector<Number2> &shape_gradients,
2569 const AlignedVector<Number2> &shape_hessians,
2570 const unsigned int dummy1 = 0,
2571 const unsigned int dummy2 = 0)
2572 : shape_values(shape_values.begin())
2573 , shape_gradients(shape_gradients.begin())
2574 , shape_hessians(shape_hessians.begin())
2575 {
2576 // We can enter this function either for the apply() path that has
2577 // n_rows * n_columns entries or for the apply_face() path that only has
2578 // n_rows * 3 entries in the array. Since we cannot decide about the use
2579 // we must allow for both here.
2580 Assert(shape_values.size() == 0 ||
2581 shape_values.size() == n_rows * n_columns ||
2582 shape_values.size() == 3 * n_rows,
2583 ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
2584 Assert(shape_gradients.size() == 0 ||
2585 shape_gradients.size() == n_rows * n_columns,
2586 ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
2587 Assert(shape_hessians.size() == 0 ||
2588 shape_hessians.size() == n_rows * n_columns,
2589 ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
2590 (void)dummy1;
2591 (void)dummy2;
2592 }
2593
2594 template <int direction, bool contract_over_rows, bool add>
2595 void
2596 values(const Number in[], Number out[]) const
2597 {
2598 apply<direction, contract_over_rows, add>(shape_values, in, out);
2599 }
2600
2601 template <int direction, bool contract_over_rows, bool add>
2602 void
2603 gradients(const Number in[], Number out[]) const
2604 {
2605 apply<direction, contract_over_rows, add>(shape_gradients, in, out);
2606 }
2607
2608 template <int direction, bool contract_over_rows, bool add>
2609 void
2610 hessians(const Number in[], Number out[]) const
2611 {
2612 apply<direction, contract_over_rows, add>(shape_hessians, in, out);
2613 }
2614
2643 template <int direction,
2644 bool contract_over_rows,
2645 bool add,
2646 bool one_line = false>
2647 static void
2648 apply(const Number2 *DEAL_II_RESTRICT shape_data,
2649 const Number * in,
2650 Number * out);
2651
2652 template <int face_direction,
2653 bool contract_onto_face,
2654 bool add,
2655 int max_derivative>
2656 void
2657 apply_face(const Number *DEAL_II_RESTRICT in,
2658 Number *DEAL_II_RESTRICT out) const;
2659
2660 private:
2661 const Number2 *shape_values;
2662 const Number2 *shape_gradients;
2663 const Number2 *shape_hessians;
2664 };
2665
2666 template <int dim,
2667 int n_rows,
2668 int n_columns,
2669 int normal_dir,
2670 typename Number,
2671 typename Number2>
2672 template <int direction, bool contract_over_rows, bool add, bool one_line>
2673 inline void
2674 EvaluatorTensorProductAnisotropic<
2675 evaluate_raviart_thomas,
2676 dim,
2677 n_rows,
2678 n_columns,
2679 normal_dir,
2680 Number,
2681 Number2>::apply(const Number2 *DEAL_II_RESTRICT shape_data,
2682 const Number * in,
2683 Number * out)
2684 {
2685 static_assert(one_line == false || direction == dim - 1,
2686 "Single-line evaluation only works for direction=dim-1.");
2687 Assert(shape_data != nullptr,
2688 ExcMessage(
2689 "The given array shape_data must not be the null pointer!"));
2690 Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
2691 in != out,
2692 ExcMessage("In-place operation only supported for "
2693 "n_rows==n_columns or single-line interpolation"));
2694 AssertIndexRange(direction, dim);
2695 constexpr int mm = contract_over_rows ? n_rows : n_columns,
2696 nn = contract_over_rows ? n_columns : n_rows;
2697
2698 constexpr int stride = Utilities::pow(n_columns, direction);
2699 constexpr int n_blocks1 = one_line ? 1 : stride;
2700
2701 // The number of blocks depend on both direction and dimension.
2702 constexpr int n_blocks2 =
2703 (dim - direction - 1 == 0) ?
2704 1 :
2705 ((direction == normal_dir) ?
2706 Utilities::pow((n_rows - 1),
2707 (direction >= dim) ? 0 : dim - direction - 1) :
2708 (((direction < normal_dir) ? (n_rows + 1) : n_rows) *
2709 ((dim - direction == 3) ? n_rows : 1)));
2710
2711 for (int i2 = 0; i2 < n_blocks2; ++i2)
2712 {
2713 for (int i1 = 0; i1 < n_blocks1; ++i1)
2714 {
2715 Number x[mm];
2716 for (int i = 0; i < mm; ++i)
2717 x[i] = in[stride * i];
2718
2719 for (int col = 0; col < nn; ++col)
2720 {
2721 Number2 val0;
2722
2723 if (contract_over_rows)
2724 val0 = shape_data[col];
2725 else
2726 val0 = shape_data[col * n_columns];
2727
2728 Number res0 = val0 * x[0];
2729 for (int i = 1; i < mm; ++i)
2730 {
2731 if (contract_over_rows)
2732 val0 = shape_data[i * n_columns + col];
2733 else
2734 val0 = shape_data[col * n_columns + i];
2735
2736 res0 += val0 * x[i];
2737 }
2738 if (add)
2739 out[stride * col] += res0;
2740
2741 else
2742 out[stride * col] = res0;
2743 }
2744
2745 if (one_line == false)
2746 {
2747 ++in;
2748 ++out;
2749 }
2750 }
2751 if (one_line == false)
2752 {
2753 in += stride * (mm - 1);
2754 out += stride * (nn - 1);
2755 }
2756 }
2757 }
2758
2759 template <int dim,
2760 int n_rows,
2761 int n_columns,
2762 int normal_dir,
2763 typename Number,
2764 typename Number2>
2765 template <int face_direction,
2766 bool contract_onto_face,
2767 bool add,
2768 int max_derivative>
2769 inline void
2770 EvaluatorTensorProductAnisotropic<
2771 evaluate_raviart_thomas,
2772 dim,
2773 n_rows,
2774 n_columns,
2775 normal_dir,
2776 Number,
2777 Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
2778 Number *DEAL_II_RESTRICT out) const
2779 {
2780 Assert(dim > 1 && dim < 4, ExcMessage("Only dim=2,3 supported"));
2781 static_assert(max_derivative >= 0 && max_derivative < 3,
2782 "Only derivative orders 0-2 implemented");
2783 Assert(shape_values != nullptr,
2784 ExcMessage(
2785 "The given array shape_values must not be the null pointer."));
2786
2787 // Determine the number of blocks depending on the face and normaldirection,
2788 // as well as dimension.
2789 constexpr int n_blocks1 = (face_direction == normal_dir) ? (n_rows - 1) :
2790 ((face_direction == 0 && normal_dir == 2) ||
2791 (face_direction == 1 && normal_dir == 2) ||
2792 (face_direction == 2 && normal_dir == 1)) ?
2793 n_rows :
2794 (n_rows + 1);
2795 constexpr int n_blocks2 = (dim == 2) ?
2796 1 :
2797 ((face_direction == normal_dir) ?
2798 (n_rows - 1) :
2799 (((face_direction == 0 && normal_dir == 1) ||
2800 (face_direction == 1 && normal_dir == 0) ||
2801 (face_direction == 2 && normal_dir == 0)) ?
2802 n_rows :
2803 (n_rows + 1)));
2804
2805 AssertIndexRange(face_direction, dim);
2806
2807 constexpr int in_stride =
2808 (face_direction == normal_dir) ?
2809 Utilities::pow(n_rows - 1, face_direction) :
2810 ((face_direction == 0) ?
2811 1 :
2812 ((face_direction == 2) ?
2813 n_rows * (n_rows + 1) :
2814 ((face_direction == 1 && normal_dir == 0) ? (n_rows + 1) :
2815 n_rows)));
2816 constexpr int out_stride = n_blocks1 * n_blocks2;
2817
2818 const Number2 *DEAL_II_RESTRICT shape_values = this->shape_values;
2819
2820 for (int i2 = 0; i2 < n_blocks2; ++i2)
2821 {
2822 for (int i1 = 0; i1 < n_blocks1; ++i1)
2823 {
2824 if (contract_onto_face == true)
2825 {
2826 Number res0 = shape_values[0] * in[0];
2827 Number res1, res2;
2828
2829 if (max_derivative > 0)
2830 res1 = shape_values[n_rows] * in[0];
2831
2832 if (max_derivative > 1)
2833 res2 = shape_values[2 * n_rows] * in[0];
2834
2835 for (int ind = 1; ind < n_rows; ++ind)
2836 {
2837 res0 += shape_values[ind] * in[in_stride * ind];
2838 if (max_derivative > 0)
2839 res1 += shape_values[ind + n_rows] * in[in_stride * ind];
2840
2841 if (max_derivative > 1)
2842 res2 +=
2843 shape_values[ind + 2 * n_rows] * in[in_stride * ind];
2844 }
2845 if (add)
2846 {
2847 out[0] += res0;
2848
2849 if (max_derivative > 0)
2850 out[out_stride] += res1;
2851
2852 if (max_derivative > 1)
2853 out[2 * out_stride] += res2;
2854 }
2855 else
2856 {
2857 out[0] = res0;
2858
2859 if (max_derivative > 0)
2860 out[out_stride] = res1;
2861
2862 if (max_derivative > 1)
2863 out[2 * out_stride] = res2;
2864 }
2865 }
2866 else
2867 {
2868 for (int col = 0; col < n_rows; ++col)
2869 {
2870 if (add)
2871 out[col * in_stride] += shape_values[col] * in[0];
2872 else
2873 out[col * in_stride] = shape_values[col] * in[0];
2874
2875 if (max_derivative > 0)
2876 out[col * in_stride] +=
2877 shape_values[col + n_rows] * in[out_stride];
2878
2879 if (max_derivative > 1)
2880 out[col * in_stride] +=
2881 shape_values[col + 2 * n_rows] * in[2 * out_stride];
2882 }
2883 }
2884
2885 // increment: in regular case, just go to the next point in
2886 // x-direction. If we are at the end of one chunk in x-dir, need
2887 // to jump over to the next layer in z-direction
2888 switch (face_direction)
2889 {
2890 case 0:
2891 in += contract_onto_face ? n_rows : 1;
2892 out += contract_onto_face ? 1 : n_rows;
2893 break;
2894
2895 case 1:
2896 ++in;
2897 ++out;
2898 // faces 2 and 3 in 3d use local coordinate system zx, which
2899 // is the other way around compared to the tensor
2900 // product. Need to take that into account.
2901 if (dim == 3)
2902 {
2903 if (normal_dir == 0)
2904 {
2905 if (contract_onto_face)
2906 out += n_rows - 1;
2907 else
2908 in += n_rows - 1;
2909 }
2910 if (normal_dir == 1)
2911 {
2912 if (contract_onto_face)
2913 out += n_rows - 2;
2914 else
2915 in += n_rows - 2;
2916 }
2917 if (normal_dir == 2)
2918 {
2919 if (contract_onto_face)
2920 out += n_rows;
2921 else
2922 in += n_rows;
2923 }
2924 }
2925 break;
2926
2927 case 2:
2928 ++in;
2929 ++out;
2930 break;
2931
2932 default:
2933 Assert(false, ExcNotImplemented());
2934 }
2935 }
2936 if (face_direction == 1 && dim == 3)
2937 {
2938 // adjust for local coordinate system zx
2939 if (contract_onto_face)
2940 {
2941 if (normal_dir == 0)
2942 {
2943 in += (n_rows + 1) * (n_rows - 1);
2944 out -= n_rows * (n_rows + 1) - 1;
2945 }
2946 if (normal_dir == 1)
2947 {
2948 in += (n_rows - 1) * (n_rows - 1);
2949 out -= (n_rows - 1) * (n_rows - 1) - 1;
2950 }
2951 if (normal_dir == 2)
2952 {
2953 in += (n_rows - 1) * (n_rows);
2954 out -= (n_rows) * (n_rows + 1) - 1;
2955 }
2956 }
2957 else
2958 {
2959 if (normal_dir == 0)
2960 {
2961 out += (n_rows + 1) * (n_rows - 1);
2962 in -= n_rows * (n_rows + 1) - 1;
2963 }
2964 if (normal_dir == 1)
2965 {
2966 out += (n_rows - 1) * (n_rows - 1);
2967 in -= (n_rows - 1) * (n_rows - 1) - 1;
2968 }
2969 if (normal_dir == 2)
2970 {
2971 out += (n_rows - 1) * (n_rows);
2972 in -= (n_rows) * (n_rows + 1) - 1;
2973 }
2974 }
2975 }
2976 }
2977 }
2978
2979
2980
2987 template <typename Number, typename Number2>
2988 struct ProductTypeNoPoint
2989 {
2990 using type = typename ProductType<Number, Number2>::type;
2991 };
2992
2993 template <int dim, typename Number, typename Number2>
2994 struct ProductTypeNoPoint<Point<dim, Number>, Number2>
2995 {
2996 using type = typename ProductType<Tensor<1, dim, Number>, Number2>::type;
2997 };
2998
2999
3000
3005 template <int dim, typename Number>
3006 inline void
3007 compute_values_of_array(
3008 ::ndarray<Number, 2, dim> * shapes,
3009 const std::vector<Polynomials::Polynomial<double>> &poly,
3010 const Point<dim, Number> & p,
3011 const unsigned int derivative = 1)
3012 {
3013 const int n_shapes = poly.size();
3014
3015 // Evaluate 1d polynomials and their derivatives
3016 std::array<Number, dim> point;
3017 for (unsigned int d = 0; d < dim; ++d)
3018 point[d] = p[d];
3019 for (int i = 0; i < n_shapes; ++i)
3020 poly[i].values_of_array(point, derivative, shapes[i].data());
3021 }
3022
3023
3024
3028 template <typename Number>
3029 inline void
3030 compute_values_of_array(::ndarray<Number, 2, 0> *,
3031 const std::vector<Polynomials::Polynomial<double>> &,
3032 const Point<0, Number> &,
3033 const unsigned int)
3034 {
3035 Assert(false, ExcInternalError());
3036 }
3037
3038
3039
3043 template <int dim,
3044 int length,
3045 typename Number2,
3046 typename Number,
3047 int n_values = 1,
3048 bool do_renumber = true>
3049 inline
3050#ifndef DEBUG
3051 DEAL_II_ALWAYS_INLINE
3052#endif
3053 std::array<typename ProductTypeNoPoint<Number, Number2>::type,
3054 2 + n_values>
3055 do_interpolate_xy(const Number * values,
3056 const std::vector<unsigned int> & renumber,
3057 const ::ndarray<Number2, 2, dim> *shapes,
3058 const int n_shapes_runtime,
3059 int & i)
3060 {
3061 static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
3062 static_assert(1 <= n_values && n_values <= 2,
3063 "Only n_values=1,2 implemented");
3064
3065 const int n_shapes = length > 0 ? length : n_shapes_runtime;
3066
3067 // If n_values > 1, we want to interpolate from a second array,
3068 // placed in the same array immediately after the main data. This
3069 // is used to interpolate normal derivatives onto faces.
3070 const Number *values_2 =
3071 n_values > 1 ?
3072 values + (length > 0 ? Utilities::pow(length, dim) :
3073 Utilities::fixed_power<dim>(n_shapes_runtime)) :
3074 nullptr;
3075 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
3076 std::array<Number3, 2 + n_values> result = {};
3077 for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
3078 {
3079 // Interpolation + derivative x direction
3080 std::array<Number3, 1 + n_values> inner_result = {};
3081
3082 // Distinguish the inner loop based on whether we have a
3083 // renumbering or not
3084 if (do_renumber && !renumber.empty())
3085 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3086 {
3087 // gradient
3088 inner_result[0] += shapes[i0][1][0] * values[renumber[i]];
3089 // values
3090 inner_result[1] += shapes[i0][0][0] * values[renumber[i]];
3091 if (n_values > 1)
3092 inner_result[2] += shapes[i0][0][0] * values_2[renumber[i]];
3093 }
3094 else
3095 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3096 {
3097 // gradient
3098 inner_result[0] += shapes[i0][1][0] * values[i];
3099 // values
3100 inner_result[1] += shapes[i0][0][0] * values[i];
3101 if (n_values > 1)
3102 inner_result[2] += shapes[i0][0][0] * values_2[i];
3103 }
3104
3105 if (dim > 1)
3106 {
3107 // Interpolation + derivative in y direction
3108 // gradient
3109 result[0] += inner_result[0] * shapes[i1][0][1];
3110 result[1] += inner_result[1] * shapes[i1][1][1];
3111 // values
3112 result[2] += inner_result[1] * shapes[i1][0][1];
3113 if (n_values > 1)
3114 result[3] += inner_result[2] * shapes[i1][0][1];
3115 }
3116 else
3117 {
3118 // gradient
3119 result[0] = inner_result[0];
3120 // values
3121 result[1] = inner_result[1];
3122 if (n_values > 1)
3123 result[2] = inner_result[2];
3124 }
3125 }
3126 return result;
3127 }
3128
3129
3130
3135 template <int dim,
3136 typename Number,
3137 typename Number2,
3138 int n_values = 1,
3139 bool do_renumber = true>
3140 inline std::array<typename ProductTypeNoPoint<Number, Number2>::type,
3141 dim + n_values>
3142 evaluate_tensor_product_value_and_gradient_shapes(
3143 const ::ndarray<Number2, 2, dim> *shapes,
3144 const int n_shapes,
3145 const Number * values,
3146 const std::vector<unsigned int> & renumber = {})
3147 {
3148 static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
3149 static_assert(1 <= n_values && n_values <= 2,
3150 "Only n_values=1,2 implemented");
3151
3152 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
3153
3154 std::array<Number3, dim + n_values> result = {};
3155 if (dim == 0)
3156 {
3157 // We only need the interpolation of the value and normal derivatives on
3158 // faces of a 1d element. As the interpolation is the value at the
3159 // point, simply set the result vector accordingly.
3160 result[0] = values[0];
3161 if (n_values > 1)
3162 result[1] = values[1];
3163 return result;
3164 }
3165
3166 // Go through the tensor product of shape functions and interpolate
3167 // with optimal algorithm
3168 for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
3169 {
3170 std::array<Number3, 2 + n_values> inner_result;
3171 // Generate separate code with known loop bounds for the most common
3172 // cases
3173 if (n_shapes == 2)
3174 inner_result =
3175 do_interpolate_xy<dim, 2, Number2, Number, n_values, do_renumber>(
3176 values, renumber, shapes, n_shapes, i);
3177 else if (n_shapes == 3)
3178 inner_result =
3179 do_interpolate_xy<dim, 3, Number2, Number, n_values, do_renumber>(
3180 values, renumber, shapes, n_shapes, i);
3181 else if (n_shapes == 4)
3182 inner_result =
3183 do_interpolate_xy<dim, 4, Number2, Number, n_values, do_renumber>(
3184 values, renumber, shapes, n_shapes, i);
3185 else if (n_shapes == 5)
3186 inner_result =
3187 do_interpolate_xy<dim, 5, Number2, Number, n_values, do_renumber>(
3188 values, renumber, shapes, n_shapes, i);
3189 else if (n_shapes == 6)
3190 inner_result =
3191 do_interpolate_xy<dim, 6, Number2, Number, n_values, do_renumber>(
3192 values, renumber, shapes, n_shapes, i);
3193 else
3194 inner_result =
3195 do_interpolate_xy<dim, -1, Number2, Number, n_values, do_renumber>(
3196 values, renumber, shapes, n_shapes, i);
3197 if (dim == 3)
3198 {
3199 // derivative + interpolation in z direction
3200 // gradient
3201 result[0] += inner_result[0] * shapes[i2][0][2];
3202 result[1] += inner_result[1] * shapes[i2][0][2];
3203 result[2] += inner_result[2] * shapes[i2][1][2];
3204 // values
3205 result[3] += inner_result[2] * shapes[i2][0][2];
3206 if (n_values > 1)
3207 result[4] += inner_result[3] * shapes[i2][0][2];
3208 }
3209 else if (dim == 2)
3210 {
3211 // gradient
3212 result[0] = inner_result[0];
3213 result[1] = inner_result[1];
3214 // values
3215 result[2] = inner_result[2];
3216 if (n_values > 1)
3217 result[3] = inner_result[3];
3218 }
3219 else
3220 {
3221 // gradient
3222 result[0] = inner_result[0];
3223 // values
3224 result[1] = inner_result[1];
3225 if (n_values > 1)
3226 result[2] = inner_result[2];
3227 }
3228 }
3229
3230 return result;
3231 }
3232
3233
3234
3239 template <int dim, typename Number, typename Number2, int n_values = 1>
3240 inline std::array<typename ProductTypeNoPoint<Number, Number2>::type,
3241 dim + n_values>
3242 evaluate_tensor_product_value_and_gradient_linear(
3243 const unsigned int n_shapes,
3244 const Number * values,
3245 const Point<dim, Number2> & p,
3246 const std::vector<unsigned int> &renumber = {})
3247 {
3248 static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
3249 static_assert(1 <= n_values && n_values <= 2,
3250 "Only n_values=1,2 implemented");
3251
3252 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
3253
3254 // If n_values > 1, we want to interpolate from a second array,
3255 // placed in the same array immediately after the main data. This
3256 // is used to interpolate normal derivatives onto faces.
3257 const Number *values_2 =
3258 n_values > 1 ? values + Utilities::fixed_power<dim>(n_shapes) : nullptr;
3259
3260 AssertDimension(n_shapes, 2);
3261 for (unsigned int i = 0; i < renumber.size(); ++i)
3262 AssertDimension(renumber[i], i);
3263
3264 std::array<Number3, dim + n_values> result;
3265 if (dim == 0)
3266 {
3267 // we only need the value on faces of a 1d element
3268 result[0] = values[0];
3269 if (n_values > 1)
3270 result[1] = values_2[0];
3271 }
3272 else if (dim == 1)
3273 {
3274 // gradient
3275 result[0] = Number3(values[1] - values[0]);
3276 // values
3277 result[1] = Number3(values[0]) + p[0] * result[0];
3278 if (n_values > 1)
3279 result[2] = Number3(values_2[0]) + p[0] * (values_2[1] - values_2[0]);
3280 }
3281 else if (dim == 2)
3282 {
3283 const Number3 val10 = Number3(values[1] - values[0]);
3284 const Number3 val32 = Number3(values[3] - values[2]);
3285 const Number3 tmp0 = Number3(values[0]) + p[0] * val10;
3286 const Number3 tmp1 = Number3(values[2]) + p[0] * val32;
3287
3288 // gradient
3289 result[0] = val10 + p[1] * (val32 - val10);
3290 result[1] = tmp1 - tmp0;
3291
3292 // values
3293 result[2] = tmp0 + p[1] * result[1];
3294
3295 if (n_values > 1)
3296 {
3297 const Number3 tmp0_2 =
3298 Number3(values_2[0]) + p[0] * (values_2[1] - values_2[0]);
3299 const Number3 tmp1_2 =
3300 Number3(values_2[2]) + p[0] * (values_2[3] - values_2[0]);
3301 result[3] = tmp0_2 + p[1] * (tmp1_2 - tmp0_2);
3302 }
3303 }
3304 else if (dim == 3)
3305 {
3306 const Number3 val10 = Number3(values[1] - values[0]);
3307 const Number3 val32 = Number3(values[3] - values[2]);
3308 const Number3 tmp0 = Number3(values[0]) + p[0] * val10;
3309 const Number3 tmp1 = Number3(values[2]) + p[0] * val32;
3310 const Number3 tmp10 = tmp1 - tmp0;
3311 const Number3 tmpy0 = tmp0 + p[1] * tmp10;
3312
3313 const Number3 val54 = Number3(values[5] - values[4]);
3314 const Number3 val76 = Number3(values[7] - values[6]);
3315 const Number3 tmp2 = Number3(values[4]) + p[0] * val54;
3316 const Number3 tmp3 = Number3(values[6]) + p[0] * val76;
3317 const Number3 tmp32 = tmp3 - tmp2;
3318 const Number3 tmpy1 = tmp2 + p[1] * tmp32;
3319
3320 // gradient
3321 result[2] = tmpy1 - tmpy0;
3322 result[1] = tmp10 + p[2] * (tmp32 - tmp10);
3323 const Number3 tmpz0 = val10 + p[1] * (val32 - val10);
3324 result[0] = tmpz0 + p[2] * (val54 + p[1] * (val76 - val54) - tmpz0);
3325
3326 // value
3327 result[3] = tmpy0 + p[2] * result[2];
3328 Assert(n_values == 1, ExcNotImplemented());
3329 }
3330
3331 return result;
3332 }
3333
3334
3335
3370 template <int dim, typename Number, typename Number2>
3371 inline std::pair<
3372 typename ProductTypeNoPoint<Number, Number2>::type,
3373 Tensor<1, dim, typename ProductTypeNoPoint<Number, Number2>::type>>
3374 evaluate_tensor_product_value_and_gradient(
3375 const std::vector<Polynomials::Polynomial<double>> &poly,
3376 const std::vector<Number> & values,
3377 const Point<dim, Number2> & p,
3378 const bool d_linear = false,
3379 const std::vector<unsigned int> & renumber = {})
3380 {
3381 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
3382
3383 std::array<Number3, dim + 1> result;
3384 if (d_linear)
3385 {
3386 result = evaluate_tensor_product_value_and_gradient_linear(
3387 poly.size(), values.data(), p, renumber);
3388 }
3389 else
3390 {
3391 AssertIndexRange(poly.size(), 200);
3392 std::array<::ndarray<Number2, 2, dim>, 200> shapes;
3393 compute_values_of_array(shapes.data(), poly, p);
3394 result = evaluate_tensor_product_value_and_gradient_shapes<dim,
3395 Number,
3396 Number2>(
3397 shapes.data(), poly.size(), values.data(), renumber);
3398 }
3399 return std::make_pair(result[dim],
3400 Tensor<1, dim, Number3>(
3401 ArrayView<Number3>(result.data(), dim)));
3402 }
3403
3404
3405
3406 template <int dim,
3407 int length,
3408 typename Number2,
3409 typename Number,
3410 bool do_renumber = true>
3411 inline
3412#ifndef DEBUG
3413 DEAL_II_ALWAYS_INLINE
3414#endif
3415 typename ProductTypeNoPoint<Number, Number2>::type
3416 do_interpolate_xy_value(const Number * values,
3417 const std::vector<unsigned int> & renumber,
3418 const ::ndarray<Number2, 2, dim> *shapes,
3419 const int n_shapes_runtime,
3420 int & i)
3421 {
3422 const int n_shapes = length > 0 ? length : n_shapes_runtime;
3423 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
3424 Number3 result = {};
3425 for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
3426 {
3427 // Interpolation x direction
3428 Number3 value = {};
3429
3430 // Distinguish the inner loop based on whether we have a
3431 // renumbering or not
3432 if (do_renumber && !renumber.empty())
3433 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3434 value += shapes[i0][0][0] * values[renumber[i]];
3435 else
3436 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3437 value += shapes[i0][0][0] * values[i];
3438
3439 if (dim > 1)
3440 result += value * shapes[i1][0][1];
3441 else
3442 result = value;
3443 }
3444 return result;
3445 }
3446
3447
3448
3449 template <int dim, typename Number, typename Number2, bool do_renumber = true>
3450 inline typename ProductTypeNoPoint<Number, Number2>::type
3451 evaluate_tensor_product_value_shapes(
3452 const ::ndarray<Number2, 2, dim> *shapes,
3453 const int n_shapes,
3454 const Number * values,
3455 const std::vector<unsigned int> & renumber = {})
3456 {
3457 static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
3458
3459 // we only need the value on faces of a 1d element
3460 if (dim == 0)
3461 {
3462 return values[0];
3463 }
3464
3465 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
3466
3467 // Go through the tensor product of shape functions and interpolate
3468 // with optimal algorithm
3469 Number3 result = {};
3470 for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
3471 {
3472 Number3 inner_result;
3473 // Generate separate code with known loop bounds for the most common
3474 // cases
3475 if (n_shapes == 2)
3476 inner_result =
3477 do_interpolate_xy_value<dim, 2, Number2, Number, do_renumber>(
3478 values, renumber, shapes, n_shapes, i);
3479 else if (n_shapes == 3)
3480 inner_result =
3481 do_interpolate_xy_value<dim, 3, Number2, Number, do_renumber>(
3482 values, renumber, shapes, n_shapes, i);
3483 else if (n_shapes == 4)
3484 inner_result =
3485 do_interpolate_xy_value<dim, 4, Number2, Number, do_renumber>(
3486 values, renumber, shapes, n_shapes, i);
3487 else if (n_shapes == 5)
3488 inner_result =
3489 do_interpolate_xy_value<dim, 5, Number2, Number, do_renumber>(
3490 values, renumber, shapes, n_shapes, i);
3491 else if (n_shapes == 6)
3492 inner_result =
3493 do_interpolate_xy_value<dim, 6, Number2, Number, do_renumber>(
3494 values, renumber, shapes, n_shapes, i);
3495 else
3496 inner_result =
3497 do_interpolate_xy_value<dim, -1, Number2, Number, do_renumber>(
3498 values, renumber, shapes, n_shapes, i);
3499 if (dim == 3)
3500 {
3501 // Interpolation + derivative in z direction
3502 result += inner_result * shapes[i2][0][2];
3503 }
3504 else
3505 {
3506 result = inner_result;
3507 }
3508 }
3509
3510 return result;
3511 }
3512
3513
3514
3515 template <int dim, typename Number, typename Number2>
3516 inline typename ProductTypeNoPoint<Number, Number2>::type
3517 evaluate_tensor_product_value_linear(
3518 const unsigned int n_shapes,
3519 const Number * values,
3520 const Point<dim, Number2> & p,
3521 const std::vector<unsigned int> &renumber = {})
3522 {
3523 (void)n_shapes;
3524 static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
3525
3526 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
3527
3528 AssertDimension(n_shapes, 2);
3529 for (unsigned int i = 0; i < renumber.size(); ++i)
3530 AssertDimension(renumber[i], i);
3531
3532 if (dim == 0)
3533 {
3534 // we only need the value on faces of a 1d element
3535 return values[0];
3536 }
3537 else if (dim == 1)
3538 {
3539 return Number3(values[0]) + p[0] * Number3(values[1] - values[0]);
3540 }
3541 else if (dim == 2)
3542 {
3543 const Number3 val10 = Number3(values[1] - values[0]);
3544 const Number3 val32 = Number3(values[3] - values[2]);
3545 const Number3 tmp0 = Number3(values[0]) + p[0] * val10;
3546 const Number3 tmp1 = Number3(values[2]) + p[0] * val32;
3547 return tmp0 + p[1] * (tmp1 - tmp0);
3548 }
3549 else if (dim == 3)
3550 {
3551 const Number3 val10 = Number3(values[1] - values[0]);
3552 const Number3 val32 = Number3(values[3] - values[2]);
3553 const Number3 tmp0 = Number3(values[0]) + p[0] * val10;
3554 const Number3 tmp1 = Number3(values[2]) + p[0] * val32;
3555 const Number3 tmpy0 = tmp0 + p[1] * (tmp1 - tmp0);
3556
3557 const Number3 val54 = Number3(values[5] - values[4]);
3558 const Number3 val76 = Number3(values[7] - values[6]);
3559 const Number3 tmp2 = Number3(values[4]) + p[0] * val54;
3560 const Number3 tmp3 = Number3(values[6]) + p[0] * val76;
3561 const Number3 tmpy1 = tmp2 + p[1] * (tmp3 - tmp2);
3562
3563 return tmpy0 + p[2] * (tmpy1 - tmpy0);
3564 }
3565
3566 // work around a compile error: missing return statement
3567 return Number3();
3568 }
3569
3570
3571
3572 template <int dim, typename Number, typename Number2>
3573 inline typename ProductTypeNoPoint<Number, Number2>::type
3574 evaluate_tensor_product_value(
3575 const std::vector<Polynomials::Polynomial<double>> &poly,
3576 const std::vector<Number> & values,
3577 const Point<dim, Number2> & p,
3578 const bool d_linear = false,
3579 const std::vector<unsigned int> & renumber = {})
3580 {
3581 typename ProductTypeNoPoint<Number, Number2>::type result;
3582 if (d_linear)
3583 {
3584 result = evaluate_tensor_product_value_linear(poly.size(),
3585 values.data(),
3586 p,
3587 renumber);
3588 }
3589 else
3590 {
3591 AssertIndexRange(poly.size(), 200);
3592 std::array<::ndarray<Number2, 2, dim>, 200> shapes;
3593 const int n_shapes = poly.size();
3594 std::array<Number2, dim> point;
3595 for (unsigned int d = 0; d < dim; ++d)
3596 point[d] = p[d];
3597 for (int i = 0; i < n_shapes; ++i)
3598 poly[i].values_of_array(point, 0, shapes[i].data());
3599 result = evaluate_tensor_product_value_shapes<dim, Number, Number2>(
3600 shapes.data(), n_shapes, values.data(), renumber);
3601 }
3602 return result;
3603 }
3604
3605
3606
3611 template <int derivative_order, typename Number, typename Number2>
3612 inline Tensor<1, 1, typename ProductTypeNoPoint<Number, Number2>::type>
3613 evaluate_tensor_product_higher_derivatives(
3614 const std::vector<Polynomials::Polynomial<double>> &poly,
3615 const std::vector<Number> & values,
3616 const Point<1, Number2> & p,
3617 const std::vector<unsigned int> & renumber = {})
3618 {
3619 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
3620
3621 const int n_shapes = poly.size();
3622 AssertDimension(n_shapes, values.size());
3623 Assert(renumber.empty() || renumber.size() == values.size(),
3624 ExcDimensionMismatch(renumber.size(), values.size()));
3625
3626 std::array<Number2, derivative_order + 1> shapes;
3627 Tensor<1, 1, Number3> result;
3628 if (renumber.empty())
3629 for (int i = 0; i < n_shapes; ++i)
3630 {
3631 poly[i].value(p[0], derivative_order, shapes.data());
3632 result[0] += shapes[derivative_order] * values[i];
3633 }
3634 else
3635 for (int i = 0; i < n_shapes; ++i)
3636 {
3637 poly[i].value(p[0], derivative_order, shapes.data());
3638 result[0] += shapes[derivative_order] * values[renumber[i]];
3639 }
3640 return result;
3641 }
3642
3643
3644
3649 template <int derivative_order, typename Number, typename Number2>
3650 inline Tensor<1,
3651 derivative_order + 1,
3652 typename ProductTypeNoPoint<Number, Number2>::type>
3653 evaluate_tensor_product_higher_derivatives(
3654 const std::vector<Polynomials::Polynomial<double>> &poly,
3655 const std::vector<Number> & values,
3656 const Point<2, Number2> & p,
3657 const std::vector<unsigned int> & renumber = {})
3658 {
3659 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
3660 constexpr int dim = 2;
3661
3662 const int n_shapes = poly.size();
3663 AssertDimension(Utilities::pow(n_shapes, 2), values.size());
3664 Assert(renumber.empty() || renumber.size() == values.size(),
3665 ExcDimensionMismatch(renumber.size(), values.size()));
3666
3667 AssertIndexRange(n_shapes, 100);
3668 ::ndarray<Number2, 100, derivative_order + 1, dim> shapes;
3669 // Evaluate 1d polynomials and their derivatives
3670 std::array<Number2, dim> point;
3671 for (unsigned int d = 0; d < dim; ++d)
3672 point[d] = p[d];
3673 for (int i = 0; i < n_shapes; ++i)
3674 poly[i].values_of_array(point, derivative_order, &shapes[i][0]);
3675
3676 Tensor<1, derivative_order + 1, Number3> result;
3677 for (int i1 = 0, i = 0; i1 < n_shapes; ++i1)
3678 {
3679 Tensor<1, derivative_order + 1, Number3> result_x;
3680 if (renumber.empty())
3681 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3682 for (unsigned int d = 0; d <= derivative_order; ++d)
3683 result_x[d] += shapes[i0][d][0] * values[i];
3684 else
3685 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3686 for (unsigned int d = 0; d <= derivative_order; ++d)
3687 result_x[d] += shapes[i0][d][0] * values[renumber[i]];
3688
3689 for (unsigned int d = 0; d <= derivative_order; ++d)
3690 result[d] += shapes[i1][d][1] * result_x[derivative_order - d];
3691 }
3692 return result;
3693 }
3694
3695
3696
3701 template <int derivative_order, typename Number, typename Number2>
3702 inline Tensor<1,
3703 ((derivative_order + 1) * (derivative_order + 2)) / 2,
3704 typename ProductTypeNoPoint<Number, Number2>::type>
3705 evaluate_tensor_product_higher_derivatives(
3706 const std::vector<Polynomials::Polynomial<double>> &poly,
3707 const std::vector<Number> & values,
3708 const Point<3, Number2> & p,
3709 const std::vector<unsigned int> & renumber = {})
3710 {
3711 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
3712 constexpr int dim = 3;
3713 constexpr int n_derivatives =
3714 ((derivative_order + 1) * (derivative_order + 2)) / 2;
3715
3716 const int n_shapes = poly.size();
3717 AssertDimension(Utilities::pow(n_shapes, 3), values.size());
3718 Assert(renumber.empty() || renumber.size() == values.size(),
3719 ExcDimensionMismatch(renumber.size(), values.size()));
3720
3721 AssertIndexRange(n_shapes, 100);
3722 ::ndarray<Number2, 100, derivative_order + 1, dim> shapes;
3723 // Evaluate 1d polynomials and their derivatives
3724 std::array<Number2, dim> point;
3725 for (unsigned int d = 0; d < dim; ++d)
3726 point[d] = p[d];
3727 for (int i = 0; i < n_shapes; ++i)
3728 poly[i].values_of_array(point, derivative_order, &shapes[i][0]);
3729
3730 Tensor<1, n_derivatives, Number3> result;
3731 for (int i2 = 0, i = 0; i2 < n_shapes; ++i2)
3732 {
3733 Tensor<1, n_derivatives, Number3> result_xy;
3734 for (int i1 = 0; i1 < n_shapes; ++i1)
3735 {
3736 // apply x derivatives
3737 Tensor<1, derivative_order + 1, Number3> result_x;
3738 if (renumber.empty())
3739 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3740 for (unsigned int d = 0; d <= derivative_order; ++d)
3741 result_x[d] += shapes[i0][d][0] * values[i];
3742 else
3743 for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3744 for (unsigned int d = 0; d <= derivative_order; ++d)
3745 result_x[d] += shapes[i0][d][0] * values[renumber[i]];
3746
3747 // multiply by y derivatives, sorting them in upper triangular
3748 // matrix, starting with highest global derivative order,
3749 // decreasing the combined order of xy derivatives by one in each
3750 // row, to be combined with z derivatives in the next step
3751 for (unsigned int d = 0, c = 0; d <= derivative_order; ++d)
3752 for (unsigned int e = d; e <= derivative_order; ++e, ++c)
3753 result_xy[c] +=
3754 shapes[i1][e - d][1] * result_x[derivative_order - e];
3755 }
3756
3757 // multiply by z derivatives, starting with highest x derivative
3758 for (unsigned int d = 0, c = 0; d <= derivative_order; ++d)
3759 for (unsigned int e = d; e <= derivative_order; ++e, ++c)
3760 result[c] += shapes[i2][d][2] * result_xy[c];
3761 }
3762 return result;
3763 }
3764
3765
3766
3767 template <int dim, typename Number, typename Number2>
3768 SymmetricTensor<2, dim, typename ProductTypeNoPoint<Number, Number2>::type>
3769 evaluate_tensor_product_hessian(
3770 const std::vector<Polynomials::Polynomial<double>> &poly,
3771 const std::vector<Number> & values,
3772 const Point<dim, Number2> & p,
3773 const std::vector<unsigned int> & renumber = {})
3774 {
3775 static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
3776
3777 const auto hessian =
3778 evaluate_tensor_product_higher_derivatives<2>(poly, values, p, renumber);
3779
3780 using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
3781 SymmetricTensor<2, dim, Number3> result;
3782 if (dim == 1)
3783 result[0][0] = hessian[0];
3784 else if (dim >= 2)
3785 {
3786 // derivatives in Hessian are xx, xy, yy, xz, yz, zz, so must re-order
3787 // them for 3D
3788 for (unsigned int d = 0, c = 0; d < 2; ++d)
3789 for (unsigned int e = d; e < 2; ++e, ++c)
3790 result[d][e] = hessian[c];
3791 if (dim == 3)
3792 {
3793 for (unsigned int d = 0; d < 2; ++d)
3794 result[d][2] = hessian[3 + d];
3795 result[2][2] = hessian[5];
3796 }
3797 }
3798
3799 return result;
3800 }
3801
3802
3803
3807 template <int dim,
3808 int length,
3809 typename Number2,
3810 typename Number,
3811 bool add,
3812 int n_values = 1>
3813 inline
3814#ifndef DEBUG
3815 DEAL_II_ALWAYS_INLINE
3816#endif
3817 void
3818 do_apply_test_functions_xy(
3819 Number2 * values,
3820 const ::ndarray<Number, 2, dim> * shapes,
3821 const std::array<Number2, 2 + n_values> &test_grads_value,
3822 const int n_shapes_runtime,
3823 int & i)
3824 {
3825 static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
3826 static_assert(1 <= n_values && n_values <= 2,
3827 "Only n_values=1,2 implemented");
3828
3829 // Note that 'add' is a template argument, so the compiler will remove
3830 // these checks
3831 if (length > 0)
3832 {
3833 constexpr unsigned int array_size = length > 0 ? length : 1;
3834 std::array<Number, array_size> shape_values_x;
3835 std::array<Number, array_size> shape_derivs_x;
3836 for (unsigned int j = 0; j < array_size; ++j)
3837 {
3838 shape_values_x[j] = shapes[j][0][0];
3839 shape_derivs_x[j] = shapes[j][1][0];
3840 }
3841 for (int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1)
3842 {
3843 const Number2 test_value_y =
3844 dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] +
3845 test_grads_value[1] * shapes[i1][1][1]) :
3846 test_grads_value[2];
3847 const Number2 test_grad_xy =
3848 dim > 1 ? test_grads_value[0] * shapes[i1][0][1] :
3849 test_grads_value[0];
3850 Number2 test_value_y_2;
3851 if (n_values > 1)
3852 test_value_y_2 = dim > 1 ?
3853 test_grads_value[3] * shapes[i1][0][1] :
3854 test_grads_value[3];
3855
3856 Number2 *values_ptr = values + i + i1 * length;
3857 Number2 *values_ptr_2 =
3858 n_values > 1 ? values_ptr + Utilities::pow(length, dim) : nullptr;
3859 for (int i0 = 0; i0 < length; ++i0)
3860 {
3861 if (add)
3862 values_ptr[i0] += shape_values_x[i0] * test_value_y;
3863 else
3864 values_ptr[i0] = shape_values_x[i0] * test_value_y;
3865 values_ptr[i0] += shape_derivs_x[i0] * test_grad_xy;
3866 if (n_values > 1)
3867 {
3868 if (add)
3869 values_ptr_2[i0] += shape_values_x[i0] * test_value_y_2;
3870 else
3871 values_ptr_2[i0] = shape_values_x[i0] * test_value_y_2;
3872 }
3873 }
3874 }
3875 i += (dim > 1 ? length * length : length);
3876 }
3877 else
3878 {
3879 for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1)
3880 {
3881 const Number2 test_value_y =
3882 dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] +
3883 test_grads_value[1] * shapes[i1][1][1]) :
3884 test_grads_value[2];
3885 const Number2 test_grad_xy =
3886 dim > 1 ? test_grads_value[0] * shapes[i1][0][1] :
3887 test_grads_value[0];
3888 Number2 test_value_y_2;
3889 if (n_values > 1)
3890 test_value_y_2 = dim > 1 ?
3891 test_grads_value[3] * shapes[i1][0][1] :
3892 test_grads_value[3];
3893
3894 Number2 *values_ptr = values + i + i1 * n_shapes_runtime;
3895 Number2 *values_ptr_2 =
3896 n_values > 1 ?
3897 values_ptr + Utilities::fixed_power<dim>(n_shapes_runtime) :
3898 nullptr;
3899 for (int i0 = 0; i0 < n_shapes_runtime; ++i0)
3900 {
3901 if (add)
3902 values_ptr[i0] += shapes[i0][0][0] * test_value_y;
3903 else
3904 values_ptr[i0] = shapes[i0][0][0] * test_value_y;
3905 values_ptr[i0] += shapes[i0][1][0] * test_grad_xy;
3906 if (n_values > 1)
3907 {
3908 if (add)
3909 values_ptr_2[i0] += shapes[i0][0][0] * test_value_y_2;
3910 else
3911 values_ptr_2[i0] = shapes[i0][0][0] * test_value_y_2;
3912 }
3913 }
3914 }
3915 i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime);
3916 }
3917 }
3918
3919
3920
3925 template <int dim,
3926 typename Number,
3927 typename Number2,
3928 bool add,
3929 int n_values = 1>
3930 inline void
3931 integrate_add_tensor_product_value_and_gradient_shapes(
3932 const ::ndarray<Number, 2, dim> *shapes,
3933 const int n_shapes,
3934 const Number2 * value,
3935 const Tensor<1, dim, Number2> & gradient,
3936 Number2 * values)
3937 {
3938 static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
3939 static_assert(1 <= n_values && n_values <= 2,
3940 "Only n_values=1,2 implemented");
3941
3942 // Note that 'add' is a template argument, so the compiler will remove
3943 // these checks
3944 if (dim == 0)
3945 {
3946 if (add)
3947 values[0] += value[0];
3948 else
3949 values[0] = value[0];
3950 if (n_values > 1)
3951 {
3952 if (add)
3953 values[1] += value[1];
3954 else
3955 values[1] = value[1];
3956 }
3957 return;
3958 }
3959
3960 // Implement the transpose of the function above
3961 // as in evaluate, use `int` type to produce better code in this context
3962 std::array<Number2, 2 + n_values> test_grads_value;
3963 for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
3964 {
3965 // test grad x
3966 test_grads_value[0] =
3967 dim > 2 ? gradient[0] * shapes[i2][0][2] : gradient[0];
3968 // test grad y
3969 test_grads_value[1] = dim > 2 ? gradient[1] * shapes[i2][0][2] :
3970 (dim > 1 ? gradient[1] : Number2());
3971 // test value z
3972 test_grads_value[2] =
3973 dim > 2 ?
3974 (value[0] * shapes[i2][0][2] + gradient[2] * shapes[i2][1][2]) :
3975 value[0];
3976
3977 if (n_values > 1)
3978 test_grads_value[3] =
3979 dim > 2 ? value[1] * shapes[i2][0][2] : value[1];
3980 // Generate separate code with known loop bounds for the most common
3981 // cases
3982 if (n_shapes == 2)
3983 do_apply_test_functions_xy<dim, 2, Number2, Number, add, n_values>(
3984 values, shapes, test_grads_value, n_shapes, i);
3985 else if (n_shapes == 3)
3986 do_apply_test_functions_xy<dim, 3, Number2, Number, add, n_values>(
3987 values, shapes, test_grads_value, n_shapes, i);
3988 else if (n_shapes == 4)
3989 do_apply_test_functions_xy<dim, 4, Number2, Number, add, n_values>(
3990 values, shapes, test_grads_value, n_shapes, i);
3991 else if (n_shapes == 5)
3992 do_apply_test_functions_xy<dim, 5, Number2, Number, add, n_values>(
3993 values, shapes, test_grads_value, n_shapes, i);
3994 else if (n_shapes == 6)
3995 do_apply_test_functions_xy<dim, 6, Number2, Number, add, n_values>(
3996 values, shapes, test_grads_value, n_shapes, i);
3997 else
3998 do_apply_test_functions_xy<dim, -1, Number2, Number, add, n_values>(
3999 values, shapes, test_grads_value, n_shapes, i);
4000 }
4001 }
4002
4003
4004
4009 template <int dim,
4010 typename Number,
4011 typename Number2,
4012 bool add,
4013 int n_values = 1>
4014 inline void
4015 integrate_add_tensor_product_value_and_gradient_linear(
4016 const unsigned int n_shapes,
4017 const Number2 * value,
4018 const Tensor<1, dim, Number2> &gradient,
4019 Number2 * values,
4020 const Point<dim, Number> & p)
4021 {
4022 (void)n_shapes;
4023 static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
4024 static_assert(1 <= n_values && n_values <= 2,
4025 "Only n_values=1,2 implemented");
4026
4027 AssertDimension(n_shapes, 2);
4028
4029 // Note that 'add' is a template argument, so the compiler will remove
4030 // these checks
4031 if (dim == 0)
4032 {
4033 if (add)
4034 values[0] += value[0];
4035 else
4036 values[0] = value[0];
4037 if (n_values > 1)
4038 {
4039 if (add)
4040 values[1] += value[1];
4041 else
4042 values[1] = value[1];
4043 }
4044 }
4045 else if (dim == 1)
4046 {
4047 const Number2 difference = value[0] * p[0] + gradient[0];
4048 if (add)
4049 {
4050 values[0] += value[0] - difference;
4051 values[1] += difference;
4052 }
4053 else
4054 {
4055 values[0] = value[0] - difference;
4056 values[1] = difference;
4057 }
4058 if (n_values > 1)
4059 {
4060 const Number2 product = value[1] * p[0];
4061 if (add)
4062 {
4063 values[2] += value[1] - product;
4064 values[3] += product;
4065 }
4066 else
4067 {
4068 values[2] = value[1] - product;
4069 values[3] = product;
4070 }
4071 }
4072 }
4073 else if (dim == 2)
4074 {
4075 const Number2 test_value_y1 = value[0] * p[1] + gradient[1];
4076 const Number2 test_value_y0 = value[0] - test_value_y1;
4077 const Number2 test_grad_xy1 = gradient[0] * p[1];
4078 const Number2 test_grad_xy0 = gradient[0] - test_grad_xy1;
4079 const Number2 value0 = p[0] * test_value_y0 + test_grad_xy0;
4080 const Number2 value1 = p[0] * test_value_y1 + test_grad_xy1;
4081
4082 if (add)
4083 {
4084 values[0] += test_value_y0 - value0;
4085 values[1] += value0;
4086 values[2] += test_value_y1 - value1;
4087 values[3] += value1;
4088 }
4089 else
4090 {
4091 values[0] = test_value_y0 - value0;
4092 values[1] = value0;
4093 values[2] = test_value_y1 - value1;
4094 values[3] = value1;
4095 }
4096
4097 if (n_values > 1)
4098 {
4099 const Number2 test_value_y1_2 = value[1] * p[1];
4100 const Number2 test_value_y0_2 = value[1] - test_value_y1_2;
4101 const Number2 value0_2 = p[0] * test_value_y1_2;
4102 const Number2 value1_2 = p[0] * test_value_y1_2;
4103
4104 if (add)
4105 {
4106 values[4] += test_value_y0_2 - value0_2;
4107 values[5] += value0_2;
4108 values[6] += test_value_y1_2 - value1_2;
4109 values[7] += value1_2;
4110 }
4111 else
4112 {
4113 values[4] = test_value_y0_2 - value0_2;
4114 values[5] = value0_2;
4115 values[6] = test_value_y1_2 - value1_2;
4116 values[7] = value1_2;
4117 }
4118 }
4119 }
4120 else if (dim == 3)
4121 {
4122 Assert(n_values == 1, ExcNotImplemented());
4123
4124 const Number2 test_value_z1 = value[0] * p[2] + gradient[2];
4125 const Number2 test_value_z0 = value[0] - test_value_z1;
4126 const Number2 test_grad_x1 = gradient[0] * p[2];
4127 const Number2 test_grad_x0 = gradient[0] - test_grad_x1;
4128 const Number2 test_grad_y1 = gradient[1] * p[2];
4129 const Number2 test_grad_y0 = gradient[1] - test_grad_y1;
4130
4131 const Number2 test_value_y01 = test_value_z0 * p[1] + test_grad_y0;
4132 const Number2 test_value_y00 = test_value_z0 - test_value_y01;
4133 const Number2 test_grad_xy01 = test_grad_x0 * p[1];
4134 const Number2 test_grad_xy00 = test_grad_x0 - test_grad_xy01;
4135 const Number2 test_value_y11 = test_value_z1 * p[1] + test_grad_y1;
4136 const Number2 test_value_y10 = test_value_z1 - test_value_y11;
4137 const Number2 test_grad_xy11 = test_grad_x1 * p[1];
4138 const Number2 test_grad_xy10 = test_grad_x1 - test_grad_xy11;
4139
4140 const Number2 value00 = p[0] * test_value_y00 + test_grad_xy00;
4141 const Number2 value01 = p[0] * test_value_y01 + test_grad_xy01;
4142 const Number2 value10 = p[0] * test_value_y10 + test_grad_xy10;
4143 const Number2 value11 = p[0] * test_value_y11 + test_grad_xy11;
4144
4145 if (add)
4146 {
4147 values[0] += test_value_y00 - value00;
4148 values[1] += value00;
4149 values[2] += test_value_y01 - value01;
4150 values[3] += value01;
4151 values[4] += test_value_y10 - value10;
4152 values[5] += value10;
4153 values[6] += test_value_y11 - value11;
4154 values[7] += value11;
4155 }
4156 else
4157 {
4158 values[0] = test_value_y00 - value00;
4159 values[1] = value00;
4160 values[2] = test_value_y01 - value01;
4161 values[3] = value01;
4162 values[4] = test_value_y10 - value10;
4163 values[5] = value10;
4164 values[6] = test_value_y11 - value11;
4165 values[7] = value11;
4166 }
4167 }
4168 }
4169
4170
4171
4177 template <int dim, typename Number, typename Number2, int n_values = 1>
4178 inline void
4179 integrate_tensor_product_value_and_gradient(
4180 const ::ndarray<Number, 2, dim> *shapes,
4181 const unsigned int n_shapes,
4182 const Number2 * value,
4183 const Tensor<1, dim, Number2> & gradient,
4184 Number2 * values,
4185 const Point<dim, Number> & p,
4186 const bool is_linear,
4187 const bool do_add)
4188 {
4189 if (do_add)
4190 {
4191 if (is_linear)
4192 internal::integrate_add_tensor_product_value_and_gradient_linear<
4193 dim,
4194 Number,
4195 Number2,
4196 true,
4197 n_values>(n_shapes, value, gradient, values, p);
4198 else
4199 internal::integrate_add_tensor_product_value_and_gradient_shapes<
4200 dim,
4201 Number,
4202 Number2,
4203 true,
4204 n_values>(shapes, n_shapes, value, gradient, values);
4205 }
4206 else
4207 {
4208 if (is_linear)
4209 internal::integrate_add_tensor_product_value_and_gradient_linear<
4210 dim,
4211 Number,
4212 Number2,
4213 false,
4214 n_values>(n_shapes, value, gradient, values, p);
4215 else
4216 internal::integrate_add_tensor_product_value_and_gradient_shapes<
4217 dim,
4218 Number,
4219 Number2,
4220 false,
4221 n_values>(shapes, n_shapes, value, gradient, values);
4222 }
4223 }
4224
4225
4226
4230 template <int dim, int length, typename Number2, typename Number, bool add>
4231 inline
4232#ifndef DEBUG
4233 DEAL_II_ALWAYS_INLINE
4234#endif
4235 void
4236 do_apply_test_functions_xy_value(
4237 Number2 * values,
4238 const ::ndarray<Number, 2, dim> *shapes,
4239 const Number2 & test_value,
4240 const int n_shapes_runtime,
4241 int & i)
4242 {
4243 if (length > 0)
4244 {
4245 constexpr unsigned int array_size = length > 0 ? length : 1;
4246 std::array<Number, array_size> shape_values_x;
4247 for (unsigned int i1 = 0; i1 < array_size; ++i1)
4248 shape_values_x[i1] = shapes[i1][0][0];
4249 for (unsigned int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1)
4250 {
4251 const Number2 test_value_y =
4252 dim > 1 ? test_value * shapes[i1][0][1] : test_value;
4253
4254 Number2 *values_ptr = values + i + i1 * length;
4255 for (unsigned int i0 = 0; i0 < length; ++i0)
4256 {
4257 if (add)
4258 values_ptr[i0] += shape_values_x[i0] * test_value_y;
4259 else
4260 values_ptr[i0] = shape_values_x[i0] * test_value_y;
4261 }
4262 }
4263 i += (dim > 1 ? length * length : length);
4264 }
4265 else
4266 {
4267 for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1)
4268 {
4269 const Number2 test_value_y =
4270 dim > 1 ? test_value * shapes[i1][0][1] : test_value;
4271
4272 Number2 *values_ptr = values + i + i1 * n_shapes_runtime;
4273 for (int i0 = 0; i0 < n_shapes_runtime; ++i0)
4274 {
4275 if (add)
4276 values_ptr[i0] += shapes[i0][0][0] * test_value_y;
4277 else
4278 values_ptr[i0] = shapes[i0][0][0] * test_value_y;
4279 }
4280 }
4281 i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime);
4282 }
4283 }
4284
4285
4286
4290 template <int dim, typename Number, typename Number2, bool add>
4291 inline void
4292 integrate_add_tensor_product_value_shapes(
4293 const ::ndarray<Number, 2, dim> *shapes,
4294 const int n_shapes,
4295 const Number2 & value,
4296 Number2 * values)
4297 {
4298 static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
4299
4300 // as in evaluate, use `int` type to produce better code in this context
4301
4302 if (dim == 0)
4303 {
4304 if (add)
4305 values[0] += value;
4306 else
4307 values[0] = value;
4308 return;
4309 }
4310
4311 // Implement the transpose of the function above
4312 Number2 test_value;
4313 for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
4314 {
4315 // test value z
4316 test_value = dim > 2 ? value * shapes[i2][0][2] : value;
4317
4318 // Generate separate code with known loop bounds for the most common
4319 // cases
4320 if (n_shapes == 2)
4321 do_apply_test_functions_xy_value<dim, 2, Number2, Number, add>(
4322 values, shapes, test_value, n_shapes, i);
4323 else if (n_shapes == 3)
4324 do_apply_test_functions_xy_value<dim, 3, Number2, Number, add>(
4325 values, shapes, test_value, n_shapes, i);
4326 else if (n_shapes == 4)
4327 do_apply_test_functions_xy_value<dim, 4, Number2, Number, add>(
4328 values, shapes, test_value, n_shapes, i);
4329 else if (n_shapes == 5)
4330 do_apply_test_functions_xy_value<dim, 5, Number2, Number, add>(
4331 values, shapes, test_value, n_shapes, i);
4332 else if (n_shapes == 6)
4333 do_apply_test_functions_xy_value<dim, 6, Number2, Number, add>(
4334 values, shapes, test_value, n_shapes, i);
4335 else
4336 do_apply_test_functions_xy_value<dim, -1, Number2, Number, add>(
4337 values, shapes, test_value, n_shapes, i);
4338 }
4339 }
4340
4341
4342
4347 template <int dim, typename Number, typename Number2, bool add>
4348 inline void
4349 integrate_add_tensor_product_value_linear(const unsigned int n_shapes,
4350 const Number2 & value,
4351 Number2 * values,
4352 const Point<dim, Number> &p)
4353 {
4354 (void)n_shapes;
4355 static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
4356
4357 AssertDimension(n_shapes, 2);
4358
4359 if (dim == 0)
4360 {
4361 if (add)
4362 values[0] += value;
4363 else
4364 values[0] = value;
4365 }
4366 else if (dim == 1)
4367 {
4368 const auto x0 = 1. - p[0], x1 = p[0];
4369
4370 if (add)
4371 {
4372 values[0] += value * x0;
4373 values[1] += value * x1;
4374 }
4375 else
4376 {
4377 values[0] = value * x0;
4378 values[1] = value * x1;
4379 }
4380 }
4381 else if (dim == 2)
4382 {
4383 const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1];
4384
4385 const auto test_value_y0 = value * y0;
4386 const auto test_value_y1 = value * y1;
4387
4388 if (add)
4389 {
4390 values[0] += x0 * test_value_y0;
4391 values[1] += x1 * test_value_y0;
4392 values[2] += x0 * test_value_y1;
4393 values[3] += x1 * test_value_y1;
4394 }
4395 else
4396 {
4397 values[0] = x0 * test_value_y0;
4398 values[1] = x1 * test_value_y0;
4399 values[2] = x0 * test_value_y1;
4400 values[3] = x1 * test_value_y1;
4401 }
4402 }
4403 else if (dim == 3)
4404 {
4405 const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
4406 z0 = 1. - p[2], z1 = p[2];
4407
4408 const auto test_value_z0 = value * z0;
4409 const auto test_value_z1 = value * z1;
4410
4411 const auto test_value_y00 = test_value_z0 * y0;
4412 const auto test_value_y01 = test_value_z0 * y1;
4413 const auto test_value_y10 = test_value_z1 * y0;
4414 const auto test_value_y11 = test_value_z1 * y1;
4415
4416 if (add)
4417 {
4418 values[0] += x0 * test_value_y00;
4419 values[1] += x1 * test_value_y00;
4420 values[2] += x0 * test_value_y01;
4421 values[3] += x1 * test_value_y01;
4422 values[4] += x0 * test_value_y10;
4423 values[5] += x1 * test_value_y10;
4424 values[6] += x0 * test_value_y11;
4425 values[7] += x1 * test_value_y11;
4426 }
4427 else
4428 {
4429 values[0] = x0 * test_value_y00;
4430 values[1] = x1 * test_value_y00;
4431 values[2] = x0 * test_value_y01;
4432 values[3] = x1 * test_value_y01;
4433 values[4] = x0 * test_value_y10;
4434 values[5] = x1 * test_value_y10;
4435 values[6] = x0 * test_value_y11;
4436 values[7] = x1 * test_value_y11;
4437 }
4438 }
4439 }
4440
4441
4442
4448 template <int dim, typename Number, typename Number2>
4449 inline void
4450 integrate_tensor_product_value(const ::ndarray<Number, 2, dim> *shapes,
4451 const unsigned int n_shapes,
4452 const Number2 & value,
4453 Number2 * values,
4454 const Point<dim, Number> &p,
4455 const bool is_linear,
4456 const bool do_add)
4457 {
4458 if (do_add)
4459 {
4460 if (is_linear)
4461 internal::integrate_add_tensor_product_value_linear<dim,
4462 Number,
4463 Number2,
4464 true>(n_shapes,
4465 value,
4466 values,
4467 p);
4468 else
4469 internal::integrate_add_tensor_product_value_shapes<dim,
4470 Number,
4471 Number2,
4472 true>(shapes,
4473 n_shapes,
4474 value,
4475 values);
4476 }
4477 else
4478 {
4479 if (is_linear)
4480 internal::integrate_add_tensor_product_value_linear<dim,
4481 Number,
4482 Number2,
4483 false>(n_shapes,
4484 value,
4485 values,
4486 p);
4487 else
4488 internal::integrate_add_tensor_product_value_shapes<dim,
4489 Number,
4490 Number2,
4491 false>(shapes,
4492 n_shapes,
4493 value,
4494 values);
4495 }
4496 }
4497
4498
4499
4500 template <int dim, int n_points_1d_template, typename Number>
4501 inline void
4502 weight_fe_q_dofs_by_entity(const Number * weights,
4503 const unsigned int n_components,
4504 const int n_points_1d_non_template,
4505 Number * data)
4506 {
4507 const int n_points_1d = n_points_1d_template != -1 ?
4508 n_points_1d_template :
4509 n_points_1d_non_template;
4510
4511 Assert(n_points_1d > 0, ExcNotImplemented());
4512 Assert(n_points_1d < 100, ExcNotImplemented());
4513
4514 unsigned int compressed_index[100];
4515 compressed_index[0] = 0;
4516 for (int i = 1; i < n_points_1d - 1; ++i)
4517 compressed_index[i] = 1;
4518 compressed_index[n_points_1d - 1] = 2;
4519
4520 for (unsigned int c = 0; c < n_components; ++c)
4521 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
4522 for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
4523 {
4524 const unsigned int shift =
4525 9 * compressed_index[k] + 3 * compressed_index[j];
4526 data[0] *= weights[shift];
4527 // loop bound as int avoids compiler warnings in case n_points_1d
4528 // == 1 (polynomial degree 0)
4529 const Number weight = weights[shift + 1];
4530 for (int i = 1; i < n_points_1d - 1; ++i)
4531 data[i] *= weight;
4532 data[n_points_1d - 1] *= weights[shift + 2];
4533 data += n_points_1d;
4534 }
4535 }
4536
4537
4538 template <int dim, int n_points_1d_template, typename Number>
4539 inline void
4540 weight_fe_q_dofs_by_entity_shifted(const Number * weights,
4541 const unsigned int n_components,
4542 const int n_points_1d_non_template,
4543 Number * data)
4544 {
4545 const int n_points_1d = n_points_1d_template != -1 ?
4546 n_points_1d_template :
4547 n_points_1d_non_template;
4548
4549 Assert((n_points_1d % 2) == 1,
4550 ExcMessage("The function can only with add number of points"));
4551 Assert(n_points_1d > 0, ExcNotImplemented());
4552 Assert(n_points_1d < 100, ExcNotImplemented());
4553
4554 const unsigned int n_inside_1d = n_points_1d / 2;
4555
4556 unsigned int compressed_index[100];
4557
4558 unsigned int c = 0;
4559 for (int i = 0; i < n_inside_1d; ++i)
4560 compressed_index[c++] = 0;
4561 compressed_index[c++] = 1;
4562 for (int i = 0; i < n_inside_1d; ++i)
4563 compressed_index[c++] = 2;
4564
4565 for (unsigned int c = 0; c < n_components; ++c)
4566 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
4567 for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
4568 {
4569 const unsigned int shift =
4570 9 * compressed_index[k] + 3 * compressed_index[j];
4571
4572 unsigned int c = 0;
4573 const Number weight1 = weights[shift];
4574 for (int i = 0; i < n_inside_1d; ++i)
4575 data[c++] *= weight1;
4576 data[c++] *= weights[shift + 1];
4577 const Number weight2 = weights[shift + 2];
4578 for (int i = 0; i < n_inside_1d; ++i)
4579 data[c++] *= weight2;
4580 data += n_points_1d;
4581 }
4582 }
4583
4584
4585 template <int dim, int n_points_1d_template, typename Number>
4586 inline bool
4587 compute_weights_fe_q_dofs_by_entity(const Number * data,
4588 const unsigned int n_components,
4589 const int n_points_1d_non_template,
4590 Number * weights)
4591 {
4592 const int n_points_1d = n_points_1d_template != -1 ?
4593 n_points_1d_template :
4594 n_points_1d_non_template;
4595
4596 Assert(n_points_1d > 0, ExcNotImplemented());
4597 Assert(n_points_1d < 100, ExcNotImplemented());
4598
4599 unsigned int compressed_index[100];
4600 compressed_index[0] = 0;
4601 for (int i = 1; i < n_points_1d - 1; ++i)
4602 compressed_index[i] = 1;
4603 compressed_index[n_points_1d - 1] = 2;
4604
4605 // Insert the number data into a storage position for weight,
4606 // ensuring that the array has either not been touched before
4607 // or the previous content is the same. In case the previous
4608 // content has a different value, we exit this function and
4609 // signal to outer functions that the compression was not possible.
4610 const auto check_and_set = [](Number &weight, const Number &data) {
4611 if (weight == Number(-1.0) || weight == data)
4612 {
4613 weight = data;
4614 return true; // success for the entry
4615 }
4616
4617 return false; // failure for the entry
4618 };
4619
4620 for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
4621 weights[c] = Number(-1.0);
4622
4623 for (unsigned int c = 0; c < n_components; ++c)
4624 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
4625 for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
4626 ++j, data += n_points_1d)
4627 {
4628 const unsigned int shift =
4629 9 * compressed_index[k] + 3 * compressed_index[j];
4630
4631 if (!check_and_set(weights[shift], data[0]))
4632 return false; // failure
4633
4634 for (int i = 1; i < n_points_1d - 1; ++i)
4635 if (!check_and_set(weights[shift + 1], data[i]))
4636 return false; // failure
4637
4638 if (!check_and_set(weights[shift + 2], data[n_points_1d - 1]))
4639 return false; // failure
4640 }
4641
4642 return true; // success
4643 }
4644
4645
4646 template <int dim, int n_points_1d_template, typename Number>
4647 inline bool
4648 compute_weights_fe_q_dofs_by_entity_shifted(
4649 const Number * data,
4650 const unsigned int n_components,
4651 const int n_points_1d_non_template,
4652 Number * weights)
4653 {
4654 const int n_points_1d = n_points_1d_template != -1 ?
4655 n_points_1d_template :
4656 n_points_1d_non_template;
4657
4658 Assert((n_points_1d % 2) == 1,
4659 ExcMessage("The function can only with add number of points"));
4660 Assert(n_points_1d > 0, ExcNotImplemented());
4661 Assert(n_points_1d < 100, ExcNotImplemented());
4662
4663 const unsigned int n_inside_1d = n_points_1d / 2;
4664
4665 unsigned int compressed_index[100];
4666
4667 unsigned int c = 0;
4668 for (int i = 0; i < n_inside_1d; ++i)
4669 compressed_index[c++] = 0;
4670 compressed_index[c++] = 1;
4671 for (int i = 0; i < n_inside_1d; ++i)
4672 compressed_index[c++] = 2;
4673
4674 // Insert the number data into a storage position for weight,
4675 // ensuring that the array has either not been touched before
4676 // or the previous content is the same. In case the previous
4677 // content has a different value, we exit this function and
4678 // signal to outer functions that the compression was not possible.
4679 const auto check_and_set = [](Number &weight, const Number &data) {
4680 if (weight == Number(-1.0) || weight == data)
4681 {
4682 weight = data;
4683 return true; // success for the entry
4684 }
4685
4686 return false; // failure for the entry
4687 };
4688
4689 for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
4690 weights[c] = Number(-1.0);
4691
4692 for (unsigned int comp = 0; comp < n_components; ++comp)
4693 for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
4694 for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
4695 ++j, data += n_points_1d)
4696 {
4697 const unsigned int shift =
4698 9 * compressed_index[k] + 3 * compressed_index[j];
4699
4700 unsigned int c = 0;
4701
4702 for (int i = 0; i < n_inside_1d; ++i)
4703 if (!check_and_set(weights[shift], data[c++]))
4704 return false; // failure
4705
4706 if (!check_and_set(weights[shift + 1], data[c++]))
4707 return false; // failure
4708
4709 for (int i = 0; i < n_inside_1d; ++i)
4710 if (!check_and_set(weights[shift + 2], data[c++]))
4711 return false; // failure
4712 }
4713
4714 return true; // success
4715 }
4716
4717
4718} // end of namespace internal
4719
4720
4721DEAL_II_NAMESPACE_CLOSE
4722
4723#endif
aligned_vector.h
AlignedVector::empty
bool empty() const
AlignedVector::size
size_type size() const
Point
Definition point.h:112
DEAL_II_ALWAYS_INLINE
#define DEAL_II_ALWAYS_INLINE
Definition config.h:106
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
DEAL_II_RESTRICT
#define DEAL_II_RESTRICT
Definition config.h:107
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1614
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1787
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:1855
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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:447
Utilities::fixed_power
T fixed_power(const T t)
Definition utilities.h:983
internal::evaluate_tensor_product_value
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value(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:3574
internal::evaluate_tensor_product_value_and_gradient_linear
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, dim+n_values > evaluate_tensor_product_value_and_gradient_linear(const unsigned int n_shapes, const Number *values, const Point< dim, Number2 > &p, const std::vector< unsigned int > &renumber={})
Definition tensor_product_kernels.h:3242
internal::weight_fe_q_dofs_by_entity
void weight_fe_q_dofs_by_entity(const Number *weights, const unsigned int n_components, const int n_points_1d_non_template, Number *data)
Definition tensor_product_kernels.h:4502
internal::do_apply_test_functions_xy_value
void do_apply_test_functions_xy_value(Number2 *values, const ::ndarray< Number, 2, dim > *shapes, const Number2 &test_value, const int n_shapes_runtime, int &i)
Definition tensor_product_kernels.h:4236
internal::integrate_add_tensor_product_value_shapes
void integrate_add_tensor_product_value_shapes(const ::ndarray< Number, 2, dim > *shapes, const int n_shapes, const Number2 &value, Number2 *values)
Definition tensor_product_kernels.h:4292
internal::even_odd_apply
void even_odd_apply(const int n_rows_in, const int n_columns_in, const Number2 *DEAL_II_RESTRICT shapes, const Number *in, Number *out)
Definition tensor_product_kernels.h:1597
internal::compute_values_of_array
void compute_values_of_array(::ndarray< Number, 2, dim > *shapes, const std::vector< Polynomials::Polynomial< double > > &poly, const Point< dim, Number > &p, const unsigned int derivative=1)
Definition tensor_product_kernels.h:3007
internal::integrate_add_tensor_product_value_and_gradient_linear
void integrate_add_tensor_product_value_and_gradient_linear(const unsigned int n_shapes, const Number2 *value, const Tensor< 1, dim, Number2 > &gradient, Number2 *values, const Point< dim, Number > &p)
Definition tensor_product_kernels.h:4015
internal::integrate_add_tensor_product_value_and_gradient_shapes
void integrate_add_tensor_product_value_and_gradient_shapes(const ::ndarray< Number, 2, dim > *shapes, const int n_shapes, const Number2 *value, const Tensor< 1, dim, Number2 > &gradient, Number2 *values)
Definition tensor_product_kernels.h:3931
internal::integrate_tensor_product_value
void integrate_tensor_product_value(const ::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 &value, Number2 *values, const Point< dim, Number > &p, const bool is_linear, const bool do_add)
Definition tensor_product_kernels.h:4450
internal::weight_fe_q_dofs_by_entity_shifted
void weight_fe_q_dofs_by_entity_shifted(const Number *weights, const unsigned int n_components, const int n_points_1d_non_template, Number *data)
Definition tensor_product_kernels.h:4540
internal::evaluate_tensor_product_value_linear
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_linear(const unsigned int n_shapes, const Number *values, const Point< dim, Number2 > &p, const std::vector< unsigned int > &renumber={})
Definition tensor_product_kernels.h:3517
internal::do_apply_test_functions_xy
void do_apply_test_functions_xy(Number2 *values, const ::ndarray< Number, 2, dim > *shapes, const std::array< Number2, 2+n_values > &test_grads_value, const int n_shapes_runtime, int &i)
Definition tensor_product_kernels.h:3818
internal::do_interpolate_xy_value
ProductTypeNoPoint< Number, Number2 >::type do_interpolate_xy_value(const Number *values, const std::vector< unsigned int > &renumber, const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes_runtime, int &i)
Definition tensor_product_kernels.h:3416
internal::evaluate_tensor_product_value_shapes
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_shapes(const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes, const Number *values, const std::vector< unsigned int > &renumber={})
Definition tensor_product_kernels.h:3451
internal::evaluate_tensor_product_higher_derivatives
Tensor< 1, 1, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_higher_derivatives(const std::vector< Polynomials::Polynomial< double > > &poly, const std::vector< Number > &values, const Point< 1, Number2 > &p, const std::vector< unsigned int > &renumber={})
Definition tensor_product_kernels.h:3613
internal::compute_weights_fe_q_dofs_by_entity
bool compute_weights_fe_q_dofs_by_entity(const Number *data, const unsigned int n_components, const int n_points_1d_non_template, Number *weights)
Definition tensor_product_kernels.h:4587
internal::compute_weights_fe_q_dofs_by_entity_shifted
bool compute_weights_fe_q_dofs_by_entity_shifted(const Number *data, const unsigned int n_components, const int n_points_1d_non_template, Number *weights)
Definition tensor_product_kernels.h:4648
internal::do_interpolate_xy
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, 2+n_values > do_interpolate_xy(const Number *values, const std::vector< unsigned int > &renumber, const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes_runtime, int &i)
Definition tensor_product_kernels.h:3055
internal::integrate_add_tensor_product_value_linear
void integrate_add_tensor_product_value_linear(const unsigned int n_shapes, const Number2 &value, Number2 *values, const Point< dim, Number > &p)
Definition tensor_product_kernels.h:4349
internal::evaluate_tensor_product_value_and_gradient_shapes
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, dim+n_values > evaluate_tensor_product_value_and_gradient_shapes(const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes, const Number *values, const std::vector< unsigned int > &renumber={})
Definition tensor_product_kernels.h:3142
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:3374
internal::integrate_tensor_product_value_and_gradient
void integrate_tensor_product_value_and_gradient(const ::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 *value, const Tensor< 1, dim, Number2 > &gradient, Number2 *values, const Point< dim, Number > &p, const bool is_linear, const bool do_add)
Definition tensor_product_kernels.h:4179
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:3769
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:213
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:485
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, normal_dir, Number, Number2 >::values
void values(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:2596
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, normal_dir, Number, Number2 >::shape_values
const Number2 * shape_values
Definition tensor_product_kernels.h:2661
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, normal_dir, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition tensor_product_kernels.h:2663
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, normal_dir, Number, 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:2566
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, normal_dir, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition tensor_product_kernels.h:2662
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, normal_dir, Number, Number2 >::EvaluatorTensorProductAnisotropic
EvaluatorTensorProductAnisotropic()
Definition tensor_product_kernels.h:2557
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, normal_dir, Number, Number2 >::gradients
void gradients(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:2603
internal::EvaluatorTensorProductAnisotropic< evaluate_raviart_thomas, dim, n_rows, n_columns, normal_dir, Number, Number2 >::hessians
void hessians(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:2610
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::hessians
void hessians(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:2069
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::gradients_one_line
void gradients_one_line(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:2085
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int n_rows=0, const unsigned int n_columns=0)
Definition tensor_product_kernels.h:2032
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::n_columns
const unsigned int n_columns
Definition tensor_product_kernels.h:2130
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const unsigned int n_rows=0, const unsigned int n_columns=0)
Definition tensor_product_kernels.h:2020
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::gradients
void gradients(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:2061
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::hessians_one_line
void hessians_one_line(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:2095
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition tensor_product_kernels.h:2012
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::shape_values
const Number2 * shape_values
Definition tensor_product_kernels.h:2126
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::values
void values(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:2053
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition tensor_product_kernels.h:2127
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::n_rows
const unsigned int n_rows
Definition tensor_product_kernels.h:2129
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::values_one_line
void values_one_line(const Number in[], Number out[]) const
Definition tensor_product_kernels.h:2077
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::apply
void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out) const
Definition tensor_product_kernels.h:2109
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, 0, 0, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition tensor_product_kernels.h:2128
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:1873
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:1938
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values)
Definition tensor_product_kernels.h:1861
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:1928
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition tensor_product_kernels.h:1997
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition tensor_product_kernels.h:1999
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition tensor_product_kernels.h:1998
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:1920
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::apply
static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out)
Definition tensor_product_kernels.h:1980
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:1896
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:1904
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition tensor_product_kernels.h:1851
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:1912
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:659
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::n_columns
const unsigned int n_columns
Definition tensor_product_kernels.h:737
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition tensor_product_kernels.h:734
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:710
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition tensor_product_kernels.h:618
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::gradients
void gradients(const Number *in, Number *out) const
Definition tensor_product_kernels.h:680
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:694
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:702
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::shape_values
const Number2 * shape_values
Definition tensor_product_kernels.h:733
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::values
void values(const Number *in, Number *out) const
Definition tensor_product_kernels.h:673
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition tensor_product_kernels.h:735
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:629
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::n_rows
const unsigned int n_rows
Definition tensor_product_kernels.h:736
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::hessians
void hessians(const Number *in, Number *out) const
Definition tensor_product_kernels.h:687
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:259
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition tensor_product_kernels.h:366
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:367
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:252
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:266
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:289
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:273
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition tensor_product_kernels.h:368
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const Number2 *shape_values, const Number2 *shape_gradients, const Number2 *shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition tensor_product_kernels.h:237
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:281
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition tensor_product_kernels.h:1092
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition tensor_product_kernels.h:1093
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:1056
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition tensor_product_kernels.h:1091
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition tensor_product_kernels.h:2310
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:2228
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:2220
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:2262
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number > &shape_values)
Definition tensor_product_kernels.h:2195
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:2244
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:2205
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:2252
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:2236
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition tensor_product_kernels.h:2309
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition tensor_product_kernels.h:2308
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition tensor_product_kernels.h:2185
internal::ProductTypeNoPoint< Point< dim, Number >, Number2 >::type
typename ProductType< Tensor< 1, dim, Number >, Number2 >::type type
Definition tensor_product_kernels.h:2996
internal::ProductTypeNoPoint::type
typename ProductType< Number, Number2 >::type type
Definition tensor_product_kernels.h:2990
Generated by doxygen 1.9.7