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