Reference documentation for deal.II version Git f0919993dd 2020-09-21 18:25:06 -0600
\(\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 - 2020 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/utilities.h>
24 
25 
27 
28 
29 
30 namespace internal
31 {
37  {
66  };
67 
68 
69 
73  enum class EvaluatorQuantity
74  {
78  value,
82  gradient,
86  hessian
87  };
88 
89 
90 
111  template <EvaluatorVariant variant,
112  int dim,
113  int n_rows,
114  int n_columns,
115  typename Number,
116  typename Number2 = Number>
118  {};
119 
120 
121 
139  template <int dim,
140  int n_rows,
141  int n_columns,
142  typename Number,
143  typename Number2>
145  dim,
146  n_rows,
147  n_columns,
148  Number,
149  Number2>
150  {
151  static constexpr unsigned int n_rows_of_product =
152  Utilities::pow(n_rows, dim);
153  static constexpr unsigned int n_columns_of_product =
154  Utilities::pow(n_columns, dim);
155 
161  : shape_values(nullptr)
162  , shape_gradients(nullptr)
163  , shape_hessians(nullptr)
164  {}
165 
170  const AlignedVector<Number2> &shape_gradients,
171  const AlignedVector<Number2> &shape_hessians,
172  const unsigned int dummy1 = 0,
173  const unsigned int dummy2 = 0)
174  : shape_values(shape_values.begin())
175  , shape_gradients(shape_gradients.begin())
176  , shape_hessians(shape_hessians.begin())
177  {
178  // We can enter this function either for the apply() path that has
179  // n_rows * n_columns entries or for the apply_face() path that only has
180  // n_rows * 3 entries in the array. Since we cannot decide about the use
181  // we must allow for both here.
182  Assert(shape_values.size() == 0 ||
183  shape_values.size() == n_rows * n_columns ||
184  shape_values.size() == 3 * n_rows,
185  ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
186  Assert(shape_gradients.size() == 0 ||
187  shape_gradients.size() == n_rows * n_columns,
188  ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
189  Assert(shape_hessians.size() == 0 ||
190  shape_hessians.size() == n_rows * n_columns,
191  ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
192  (void)dummy1;
193  (void)dummy2;
194  }
195 
196  template <int direction, bool contract_over_rows, bool add>
197  void
198  values(const Number in[], Number out[]) const
199  {
200  apply<direction, contract_over_rows, add>(shape_values, in, out);
201  }
202 
203  template <int direction, bool contract_over_rows, bool add>
204  void
205  gradients(const Number in[], Number out[]) const
206  {
207  apply<direction, contract_over_rows, add>(shape_gradients, in, out);
208  }
209 
210  template <int direction, bool contract_over_rows, bool add>
211  void
212  hessians(const Number in[], Number out[]) const
213  {
214  apply<direction, contract_over_rows, add>(shape_hessians, in, out);
215  }
216 
217  template <int direction, bool contract_over_rows, bool add>
218  void
219  values_one_line(const Number in[], Number out[]) const
220  {
221  Assert(shape_values != nullptr, ExcNotInitialized());
222  apply<direction, contract_over_rows, add, true>(shape_values, in, out);
223  }
224 
225  template <int direction, bool contract_over_rows, bool add>
226  void
227  gradients_one_line(const Number in[], Number out[]) const
228  {
229  Assert(shape_gradients != nullptr, ExcNotInitialized());
230  apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
231  }
232 
233  template <int direction, bool contract_over_rows, bool add>
234  void
235  hessians_one_line(const Number in[], Number out[]) const
236  {
237  Assert(shape_hessians != nullptr, ExcNotInitialized());
238  apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
239  }
240 
265  template <int direction,
266  bool contract_over_rows,
267  bool add,
268  bool one_line = false>
269  static void
270  apply(const Number2 *DEAL_II_RESTRICT shape_data,
271  const Number * in,
272  Number * out);
273 
308  template <int face_direction,
309  bool contract_onto_face,
310  bool add,
311  int max_derivative,
312  bool lex_faces = false>
313  void
314  apply_face(const Number *DEAL_II_RESTRICT in,
315  Number *DEAL_II_RESTRICT out) const;
316 
317  const Number2 *shape_values;
318  const Number2 *shape_gradients;
319  const Number2 *shape_hessians;
320  };
321 
322 
323 
324  template <int dim,
325  int n_rows,
326  int n_columns,
327  typename Number,
328  typename Number2>
329  template <int direction, bool contract_over_rows, bool add, bool one_line>
330  inline void
332  dim,
333  n_rows,
334  n_columns,
335  Number,
336  Number2>::apply(const Number2 *DEAL_II_RESTRICT
337  shape_data,
338  const Number * in,
339  Number * out)
340  {
341  static_assert(one_line == false || direction == dim - 1,
342  "Single-line evaluation only works for direction=dim-1.");
343  Assert(shape_data != nullptr,
344  ExcMessage(
345  "The given array shape_data must not be the null pointer!"));
346  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
347  in != out,
348  ExcMessage("In-place operation only supported for "
349  "n_rows==n_columns or single-line interpolation"));
350  AssertIndexRange(direction, dim);
351  constexpr int mm = contract_over_rows ? n_rows : n_columns,
352  nn = contract_over_rows ? n_columns : n_rows;
353 
354  constexpr int stride = Utilities::pow(n_columns, direction);
355  constexpr int n_blocks1 = one_line ? 1 : stride;
356  constexpr int n_blocks2 =
357  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
358 
359  for (int i2 = 0; i2 < n_blocks2; ++i2)
360  {
361  for (int i1 = 0; i1 < n_blocks1; ++i1)
362  {
363  Number x[mm];
364  for (int i = 0; i < mm; ++i)
365  x[i] = in[stride * i];
366  for (int col = 0; col < nn; ++col)
367  {
368  Number2 val0;
369  if (contract_over_rows == true)
370  val0 = shape_data[col];
371  else
372  val0 = shape_data[col * n_columns];
373  Number res0 = val0 * x[0];
374  for (int i = 1; i < mm; ++i)
375  {
376  if (contract_over_rows == true)
377  val0 = shape_data[i * n_columns + col];
378  else
379  val0 = shape_data[col * n_columns + i];
380  res0 += val0 * x[i];
381  }
382  if (add == false)
383  out[stride * col] = res0;
384  else
385  out[stride * col] += res0;
386  }
387 
388  if (one_line == false)
389  {
390  ++in;
391  ++out;
392  }
393  }
394  if (one_line == false)
395  {
396  in += stride * (mm - 1);
397  out += stride * (nn - 1);
398  }
399  }
400  }
401 
402 
403 
404  template <int dim,
405  int n_rows,
406  int n_columns,
407  typename Number,
408  typename Number2>
409  template <int face_direction,
410  bool contract_onto_face,
411  bool add,
412  int max_derivative,
413  bool lex_faces>
414  inline void
416  dim,
417  n_rows,
418  n_columns,
419  Number,
420  Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
421  Number *DEAL_II_RESTRICT
422  out) const
423  {
424  Assert(dim > 0 && (lex_faces || dim < 4),
425  ExcMessage("Only dim=1,2,3 supported"));
426  static_assert(max_derivative >= 0 && max_derivative < 3,
427  "Only derivative orders 0-2 implemented");
428  Assert(shape_values != nullptr,
429  ExcMessage(
430  "The given array shape_values must not be the null pointer."));
431 
432  constexpr int n_blocks1 =
433  lex_faces ? ::Utilities::pow<unsigned int>(n_rows, face_direction) :
434  (dim > 1 ? n_rows : 1);
435  constexpr int n_blocks2 =
436  lex_faces ? ::Utilities::pow<unsigned int>(
437  n_rows, std::max(dim - face_direction - 1, 0)) :
438  (dim > 2 ? n_rows : 1);
439 
440  AssertIndexRange(face_direction, dim);
441  constexpr int stride = Utilities::pow(n_rows, face_direction);
442  constexpr int out_stride = Utilities::pow(n_rows, dim - 1);
443  const Number *DEAL_II_RESTRICT shape_values = this->shape_values;
444 
445  for (int i2 = 0; i2 < n_blocks2; ++i2)
446  {
447  for (int i1 = 0; i1 < n_blocks1; ++i1)
448  {
449  if (contract_onto_face == true)
450  {
451  Number res0 = shape_values[0] * in[0];
452  Number res1, res2;
453  if (max_derivative > 0)
454  res1 = shape_values[n_rows] * in[0];
455  if (max_derivative > 1)
456  res2 = shape_values[2 * n_rows] * in[0];
457  for (int ind = 1; ind < n_rows; ++ind)
458  {
459  res0 += shape_values[ind] * in[stride * ind];
460  if (max_derivative > 0)
461  res1 += shape_values[ind + n_rows] * in[stride * ind];
462  if (max_derivative > 1)
463  res2 += shape_values[ind + 2 * n_rows] * in[stride * ind];
464  }
465  if (add == false)
466  {
467  out[0] = res0;
468  if (max_derivative > 0)
469  out[out_stride] = res1;
470  if (max_derivative > 1)
471  out[2 * out_stride] = res2;
472  }
473  else
474  {
475  out[0] += res0;
476  if (max_derivative > 0)
477  out[out_stride] += res1;
478  if (max_derivative > 1)
479  out[2 * out_stride] += res2;
480  }
481  }
482  else
483  {
484  for (int col = 0; col < n_rows; ++col)
485  {
486  if (add == false)
487  out[col * stride] = shape_values[col] * in[0];
488  else
489  out[col * stride] += shape_values[col] * in[0];
490  if (max_derivative > 0)
491  out[col * stride] +=
492  shape_values[col + n_rows] * in[out_stride];
493  if (max_derivative > 1)
494  out[col * stride] +=
495  shape_values[col + 2 * n_rows] * in[2 * out_stride];
496  }
497  }
498 
499  if (lex_faces)
500  {
501  ++out;
502  ++in;
503  }
504  else
505  // increment: in regular case, just go to the next point in
506  // x-direction. If we are at the end of one chunk in x-dir, need
507  // to jump over to the next layer in z-direction
508  switch (face_direction)
509  {
510  case 0:
511  in += contract_onto_face ? n_rows : 1;
512  out += contract_onto_face ? 1 : n_rows;
513  break;
514  case 1:
515  ++in;
516  ++out;
517  // faces 2 and 3 in 3D use local coordinate system zx, which
518  // is the other way around compared to the tensor
519  // product. Need to take that into account.
520  if (dim == 3)
521  {
522  if (contract_onto_face)
523  out += n_rows - 1;
524  else
525  in += n_rows - 1;
526  }
527  break;
528  case 2:
529  ++in;
530  ++out;
531  break;
532  default:
533  Assert(false, ExcNotImplemented());
534  }
535  }
536  if (lex_faces)
537  {
538  if (contract_onto_face)
539  in += (::Utilities::pow(n_rows, face_direction + 1) -
540  n_blocks1);
541  else
542  out += (::Utilities::pow(n_rows, face_direction + 1) -
543  n_blocks1);
544  }
545  else if (face_direction == 1 && dim == 3)
546  {
547  // adjust for local coordinate system zx
548  if (contract_onto_face)
549  {
550  in += n_rows * (n_rows - 1);
551  out -= n_rows * n_rows - 1;
552  }
553  else
554  {
555  out += n_rows * (n_rows - 1);
556  in -= n_rows * n_rows - 1;
557  }
558  }
559  }
560  }
561 
562 
563 
577  template <int dim, typename Number, typename Number2>
578  struct EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>
579  {
580  static constexpr unsigned int n_rows_of_product =
582  static constexpr unsigned int n_columns_of_product =
584 
590  : shape_values(nullptr)
591  , shape_gradients(nullptr)
592  , shape_hessians(nullptr)
593  , n_rows(numbers::invalid_unsigned_int)
594  , n_columns(numbers::invalid_unsigned_int)
595  {}
596 
601  const AlignedVector<Number2> &shape_gradients,
602  const AlignedVector<Number2> &shape_hessians,
603  const unsigned int n_rows,
604  const unsigned int n_columns)
605  : shape_values(shape_values.begin())
606  , shape_gradients(shape_gradients.begin())
607  , shape_hessians(shape_hessians.begin())
608  , n_rows(n_rows)
609  , n_columns(n_columns)
610  {
611  // We can enter this function either for the apply() path that has
612  // n_rows * n_columns entries or for the apply_face() path that only has
613  // n_rows * 3 entries in the array. Since we cannot decide about the use
614  // we must allow for both here.
615  Assert(shape_values.size() == 0 ||
616  shape_values.size() == n_rows * n_columns ||
617  shape_values.size() == n_rows * 3,
618  ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
619  Assert(shape_gradients.size() == 0 ||
620  shape_gradients.size() == n_rows * n_columns,
621  ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
622  Assert(shape_hessians.size() == 0 ||
623  shape_hessians.size() == n_rows * n_columns,
624  ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
625  }
626 
627  template <int direction, bool contract_over_rows, bool add>
628  void
629  values(const Number *in, Number *out) const
630  {
631  apply<direction, contract_over_rows, add>(shape_values, in, out);
632  }
633 
634  template <int direction, bool contract_over_rows, bool add>
635  void
636  gradients(const Number *in, Number *out) const
637  {
638  apply<direction, contract_over_rows, add>(shape_gradients, in, out);
639  }
640 
641  template <int direction, bool contract_over_rows, bool add>
642  void
643  hessians(const Number *in, Number *out) const
644  {
645  apply<direction, contract_over_rows, add>(shape_hessians, in, out);
646  }
647 
648  template <int direction, bool contract_over_rows, bool add>
649  void
650  values_one_line(const Number in[], Number out[]) const
651  {
652  Assert(shape_values != nullptr, ExcNotInitialized());
653  apply<direction, contract_over_rows, add, true>(shape_values, in, out);
654  }
655 
656  template <int direction, bool contract_over_rows, bool add>
657  void
658  gradients_one_line(const Number in[], Number out[]) const
659  {
660  Assert(shape_gradients != nullptr, ExcNotInitialized());
661  apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
662  }
663 
664  template <int direction, bool contract_over_rows, bool add>
665  void
666  hessians_one_line(const Number in[], Number out[]) const
667  {
668  Assert(shape_hessians != nullptr, ExcNotInitialized());
669  apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
670  }
671 
672  template <int direction,
673  bool contract_over_rows,
674  bool add,
675  bool one_line = false>
676  void
677  apply(const Number2 *DEAL_II_RESTRICT shape_data,
678  const Number * in,
679  Number * out) const;
680 
681  template <int face_direction,
682  bool contract_onto_face,
683  bool add,
684  int max_derivative,
685  bool lex_faces = false>
686  void
687  apply_face(const Number *DEAL_II_RESTRICT in,
688  Number *DEAL_II_RESTRICT out) const;
689 
690  const Number2 * shape_values;
691  const Number2 * shape_gradients;
692  const Number2 * shape_hessians;
693  const unsigned int n_rows;
694  const unsigned int n_columns;
695  };
696 
697 
698 
699  template <int dim, typename Number, typename Number2>
700  template <int direction, bool contract_over_rows, bool add, bool one_line>
701  inline void
703  const Number2 *DEAL_II_RESTRICT shape_data,
704  const Number * in,
705  Number * out) const
706  {
707  static_assert(one_line == false || direction == dim - 1,
708  "Single-line evaluation only works for direction=dim-1.");
709  Assert(shape_data != nullptr,
710  ExcMessage(
711  "The given array shape_data must not be the null pointer!"));
712  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
713  in != out,
714  ExcMessage("In-place operation only supported for "
715  "n_rows==n_columns or single-line interpolation"));
716  AssertIndexRange(direction, dim);
717  const int mm = contract_over_rows ? n_rows : n_columns,
718  nn = contract_over_rows ? n_columns : n_rows;
719 
720  const int stride =
721  direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
722  const int n_blocks1 = one_line ? 1 : stride;
723  const int n_blocks2 = direction >= dim - 1 ?
724  1 :
725  Utilities::fixed_power<dim - direction - 1>(n_rows);
726  Assert(n_rows <= 128, ExcNotImplemented());
727 
728  // specialization for n_rows = 2 that manually unrolls the innermost loop
729  // to make the operation perform better (not completely as good as the
730  // templated one, but much better than the generic version down below,
731  // because the loop over col can be more effectively unrolled by the
732  // compiler)
733  if (contract_over_rows && n_rows == 2)
734  {
735  const Number2 *shape_data_1 = shape_data + n_columns;
736  for (int i2 = 0; i2 < n_blocks2; ++i2)
737  {
738  for (int i1 = 0; i1 < n_blocks1; ++i1)
739  {
740  const Number x0 = in[0], x1 = in[stride];
741  for (int col = 0; col < nn; ++col)
742  {
743  const Number result =
744  shape_data[col] * x0 + shape_data_1[col] * x1;
745  if (add == false)
746  out[stride * col] = result;
747  else
748  out[stride * col] += result;
749  }
750 
751  if (one_line == false)
752  {
753  ++in;
754  ++out;
755  }
756  }
757  if (one_line == false)
758  {
759  in += stride * (mm - 1);
760  out += stride * (nn - 1);
761  }
762  }
763  }
764  // specialization for n = 3
765  else if (contract_over_rows && n_rows == 3)
766  {
767  const Number2 *shape_data_1 = shape_data + n_columns;
768  const Number2 *shape_data_2 = shape_data + 2 * n_columns;
769  for (int i2 = 0; i2 < n_blocks2; ++i2)
770  {
771  for (int i1 = 0; i1 < n_blocks1; ++i1)
772  {
773  const Number x0 = in[0], x1 = in[stride], x2 = in[2 * stride];
774  for (int col = 0; col < nn; ++col)
775  {
776  const Number result = shape_data[col] * x0 +
777  shape_data_1[col] * x1 +
778  shape_data_2[col] * x2;
779  if (add == false)
780  out[stride * col] = result;
781  else
782  out[stride * col] += result;
783  }
784 
785  if (one_line == false)
786  {
787  ++in;
788  ++out;
789  }
790  }
791  if (one_line == false)
792  {
793  in += stride * (mm - 1);
794  out += stride * (nn - 1);
795  }
796  }
797  }
798  // general loop for all other cases
799  else
800  for (int i2 = 0; i2 < n_blocks2; ++i2)
801  {
802  for (int i1 = 0; i1 < n_blocks1; ++i1)
803  {
804  Number x[129];
805  for (int i = 0; i < mm; ++i)
806  x[i] = in[stride * i];
807  for (int col = 0; col < nn; ++col)
808  {
809  Number2 val0;
810  if (contract_over_rows == true)
811  val0 = shape_data[col];
812  else
813  val0 = shape_data[col * n_columns];
814  Number res0 = val0 * x[0];
815  for (int i = 1; i < mm; ++i)
816  {
817  if (contract_over_rows == true)
818  val0 = shape_data[i * n_columns + col];
819  else
820  val0 = shape_data[col * n_columns + i];
821  res0 += val0 * x[i];
822  }
823  if (add == false)
824  out[stride * col] = res0;
825  else
826  out[stride * col] += res0;
827  }
828 
829  if (one_line == false)
830  {
831  ++in;
832  ++out;
833  }
834  }
835  if (one_line == false)
836  {
837  in += stride * (mm - 1);
838  out += stride * (nn - 1);
839  }
840  }
841  }
842 
843 
844 
845  template <int dim, typename Number, typename Number2>
846  template <int face_direction,
847  bool contract_onto_face,
848  bool add,
849  int max_derivative,
850  bool lex_faces>
851  inline void
853  apply_face(const Number *DEAL_II_RESTRICT in,
854  Number *DEAL_II_RESTRICT out) const
855  {
856  static_assert(lex_faces == false, "Not implemented yet.");
857 
858  Assert(shape_values != nullptr,
859  ExcMessage(
860  "The given array shape_data must not be the null pointer!"));
861  static_assert(dim > 0 && dim < 4, "Only dim=1,2,3 supported");
862  const int n_blocks1 = dim > 1 ? n_rows : 1;
863  const int n_blocks2 = dim > 2 ? n_rows : 1;
864 
865  AssertIndexRange(face_direction, dim);
866  const int stride =
867  face_direction > 0 ? Utilities::fixed_power<face_direction>(n_rows) : 1;
868  const int out_stride =
869  dim > 1 ? Utilities::fixed_power<dim - 1>(n_rows) : 1;
870 
871  for (int i2 = 0; i2 < n_blocks2; ++i2)
872  {
873  for (int i1 = 0; i1 < n_blocks1; ++i1)
874  {
875  if (contract_onto_face == true)
876  {
877  Number res0 = shape_values[0] * in[0];
878  Number res1, res2;
879  if (max_derivative > 0)
880  res1 = shape_values[n_rows] * in[0];
881  if (max_derivative > 1)
882  res2 = shape_values[2 * n_rows] * in[0];
883  for (unsigned int ind = 1; ind < n_rows; ++ind)
884  {
885  res0 += shape_values[ind] * in[stride * ind];
886  if (max_derivative > 0)
887  res1 += shape_values[ind + n_rows] * in[stride * ind];
888  if (max_derivative > 1)
889  res2 += shape_values[ind + 2 * n_rows] * in[stride * ind];
890  }
891  if (add == false)
892  {
893  out[0] = res0;
894  if (max_derivative > 0)
895  out[out_stride] = res1;
896  if (max_derivative > 1)
897  out[2 * out_stride] = res2;
898  }
899  else
900  {
901  out[0] += res0;
902  if (max_derivative > 0)
903  out[out_stride] += res1;
904  if (max_derivative > 1)
905  out[2 * out_stride] += res2;
906  }
907  }
908  else
909  {
910  for (unsigned int col = 0; col < n_rows; ++col)
911  {
912  if (add == false)
913  out[col * stride] = shape_values[col] * in[0];
914  else
915  out[col * stride] += shape_values[col] * in[0];
916  if (max_derivative > 0)
917  out[col * stride] +=
918  shape_values[col + n_rows] * in[out_stride];
919  if (max_derivative > 1)
920  out[col * stride] +=
921  shape_values[col + 2 * n_rows] * in[2 * out_stride];
922  }
923  }
924 
925  // increment: in regular case, just go to the next point in
926  // x-direction. If we are at the end of one chunk in x-dir, need
927  // to jump over to the next layer in z-direction
928  switch (face_direction)
929  {
930  case 0:
931  in += contract_onto_face ? n_rows : 1;
932  out += contract_onto_face ? 1 : n_rows;
933  break;
934  case 1:
935  ++in;
936  ++out;
937  // faces 2 and 3 in 3D use local coordinate system zx, which
938  // is the other way around compared to the tensor
939  // product. Need to take that into account.
940  if (dim == 3)
941  {
942  if (contract_onto_face)
943  out += n_rows - 1;
944  else
945  in += n_rows - 1;
946  }
947  break;
948  case 2:
949  ++in;
950  ++out;
951  break;
952  default:
953  Assert(false, ExcNotImplemented());
954  }
955  }
956  if (face_direction == 1 && dim == 3)
957  {
958  // adjust for local coordinate system zx
959  if (contract_onto_face)
960  {
961  in += n_rows * (n_rows - 1);
962  out -= n_rows * n_rows - 1;
963  }
964  else
965  {
966  out += n_rows * (n_rows - 1);
967  in -= n_rows * n_rows - 1;
968  }
969  }
970  }
971  }
972 
973 
974 
995  template <int dim,
996  int n_rows,
997  int n_columns,
998  typename Number,
999  typename Number2>
1001  dim,
1002  n_rows,
1003  n_columns,
1004  Number,
1005  Number2>
1006  {
1007  static constexpr unsigned int n_rows_of_product =
1008  Utilities::pow(n_rows, dim);
1009  static constexpr unsigned int n_columns_of_product =
1010  Utilities::pow(n_columns, dim);
1011 
1016  const AlignedVector<Number2> &shape_gradients,
1017  const AlignedVector<Number2> &shape_hessians,
1018  const unsigned int dummy1 = 0,
1019  const unsigned int dummy2 = 0)
1020  : shape_values(shape_values.begin())
1021  , shape_gradients(shape_gradients.begin())
1022  , shape_hessians(shape_hessians.begin())
1023  {
1024  Assert(shape_values.size() == 0 ||
1025  shape_values.size() == n_rows * n_columns,
1026  ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1027  Assert(shape_gradients.size() == 0 ||
1028  shape_gradients.size() == n_rows * n_columns,
1029  ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
1030  Assert(shape_hessians.size() == 0 ||
1031  shape_hessians.size() == n_rows * n_columns,
1032  ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
1033  (void)dummy1;
1034  (void)dummy2;
1035  }
1036 
1037  template <int direction, bool contract_over_rows, bool add>
1038  void
1039  values(const Number in[], Number out[]) const;
1040 
1041  template <int direction, bool contract_over_rows, bool add>
1042  void
1043  gradients(const Number in[], Number out[]) const;
1044 
1045  template <int direction, bool contract_over_rows, bool add>
1046  void
1047  hessians(const Number in[], Number out[]) const;
1048 
1049  const Number2 *shape_values;
1050  const Number2 *shape_gradients;
1051  const Number2 *shape_hessians;
1052  };
1053 
1054 
1055 
1056  // In this case, the 1D shape values read (sorted lexicographically, rows
1057  // run over 1D dofs, columns over quadrature points):
1058  // Q2 --> [ 0.687 0 -0.087 ]
1059  // [ 0.4 1 0.4 ]
1060  // [-0.087 0 0.687 ]
1061  // Q3 --> [ 0.66 0.003 0.002 0.049 ]
1062  // [ 0.521 1.005 -0.01 -0.230 ]
1063  // [-0.230 -0.01 1.005 0.521 ]
1064  // [ 0.049 0.002 0.003 0.66 ]
1065  // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ]
1066  // [ 0.608 1.059 0 0.039 0.176 ]
1067  // [-0.409 -0.113 1 -0.113 -0.409 ]
1068  // [ 0.176 0.039 0 1.059 0.608 ]
1069  // [-0.032 -0.007 0 0.022 0.658 ]
1070  //
1071  // In these matrices, we want to use avoid computations involving zeros and
1072  // ones and in addition use the symmetry in entries to reduce the number of
1073  // read operations.
1074  template <int dim,
1075  int n_rows,
1076  int n_columns,
1077  typename Number,
1078  typename Number2>
1079  template <int direction, bool contract_over_rows, bool add>
1080  inline void
1082  dim,
1083  n_rows,
1084  n_columns,
1085  Number,
1086  Number2>::values(const Number in[], Number out[]) const
1087  {
1088  Assert(shape_values != nullptr, ExcNotInitialized());
1089  AssertIndexRange(direction, dim);
1090  constexpr int mm = contract_over_rows ? n_rows : n_columns,
1091  nn = contract_over_rows ? n_columns : n_rows;
1092  constexpr int n_cols = nn / 2;
1093  constexpr int mid = mm / 2;
1094 
1095  constexpr int stride = Utilities::pow(n_columns, direction);
1096  constexpr int n_blocks1 = stride;
1097  constexpr int n_blocks2 =
1098  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1099 
1100  for (int i2 = 0; i2 < n_blocks2; ++i2)
1101  {
1102  for (int i1 = 0; i1 < n_blocks1; ++i1)
1103  {
1104  for (int col = 0; col < n_cols; ++col)
1105  {
1106  Number2 val0, val1;
1107  Number in0, in1, res0, res1;
1108  if (contract_over_rows == true)
1109  {
1110  val0 = shape_values[col];
1111  val1 = shape_values[nn - 1 - col];
1112  }
1113  else
1114  {
1115  val0 = shape_values[col * n_columns];
1116  val1 = shape_values[(col + 1) * n_columns - 1];
1117  }
1118  if (mid > 0)
1119  {
1120  in0 = in[0];
1121  in1 = in[stride * (mm - 1)];
1122  res0 = val0 * in0;
1123  res1 = val1 * in0;
1124  res0 += val1 * in1;
1125  res1 += val0 * in1;
1126  for (int ind = 1; ind < mid; ++ind)
1127  {
1128  if (contract_over_rows == true)
1129  {
1130  val0 = shape_values[ind * n_columns + col];
1131  val1 = shape_values[ind * n_columns + nn - 1 - col];
1132  }
1133  else
1134  {
1135  val0 = shape_values[col * n_columns + ind];
1136  val1 =
1137  shape_values[(col + 1) * n_columns - 1 - ind];
1138  }
1139  in0 = in[stride * ind];
1140  in1 = in[stride * (mm - 1 - ind)];
1141  res0 += val0 * in0;
1142  res1 += val1 * in0;
1143  res0 += val1 * in1;
1144  res1 += val0 * in1;
1145  }
1146  }
1147  else
1148  res0 = res1 = Number();
1149  if (contract_over_rows == true)
1150  {
1151  if (mm % 2 == 1)
1152  {
1153  val0 = shape_values[mid * n_columns + col];
1154  in1 = val0 * in[stride * mid];
1155  res0 += in1;
1156  res1 += in1;
1157  }
1158  }
1159  else
1160  {
1161  if (mm % 2 == 1 && nn % 2 == 0)
1162  {
1163  val0 = shape_values[col * n_columns + mid];
1164  in1 = val0 * in[stride * mid];
1165  res0 += in1;
1166  res1 += in1;
1167  }
1168  }
1169  if (add == false)
1170  {
1171  out[stride * col] = res0;
1172  out[stride * (nn - 1 - col)] = res1;
1173  }
1174  else
1175  {
1176  out[stride * col] += res0;
1177  out[stride * (nn - 1 - col)] += res1;
1178  }
1179  }
1180  if (contract_over_rows == true && nn % 2 == 1 && mm % 2 == 1)
1181  {
1182  if (add == false)
1183  out[stride * n_cols] = in[stride * mid];
1184  else
1185  out[stride * n_cols] += in[stride * mid];
1186  }
1187  else if (contract_over_rows == true && nn % 2 == 1)
1188  {
1189  Number res0;
1190  Number2 val0 = shape_values[n_cols];
1191  if (mid > 0)
1192  {
1193  res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1194  for (int ind = 1; ind < mid; ++ind)
1195  {
1196  val0 = shape_values[ind * n_columns + n_cols];
1197  res0 += val0 * (in[stride * ind] +
1198  in[stride * (mm - 1 - ind)]);
1199  }
1200  }
1201  else
1202  res0 = Number();
1203  if (add == false)
1204  out[stride * n_cols] = res0;
1205  else
1206  out[stride * n_cols] += res0;
1207  }
1208  else if (contract_over_rows == false && nn % 2 == 1)
1209  {
1210  Number res0;
1211  if (mid > 0)
1212  {
1213  Number2 val0 = shape_values[n_cols * n_columns];
1214  res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1215  for (int ind = 1; ind < mid; ++ind)
1216  {
1217  val0 = shape_values[n_cols * n_columns + ind];
1218  Number in1 = val0 * (in[stride * ind] +
1219  in[stride * (mm - 1 - ind)]);
1220  res0 += in1;
1221  }
1222  if (mm % 2)
1223  res0 += in[stride * mid];
1224  }
1225  else
1226  res0 = in[0];
1227  if (add == false)
1228  out[stride * n_cols] = res0;
1229  else
1230  out[stride * n_cols] += res0;
1231  }
1232 
1233  ++in;
1234  ++out;
1235  }
1236  in += stride * (mm - 1);
1237  out += stride * (nn - 1);
1238  }
1239  }
1240 
1241 
1242 
1243  // For the specialized loop used for the gradient computation in
1244  // here, the 1D shape values read (sorted lexicographically, rows
1245  // run over 1D dofs, columns over quadrature points):
1246  // Q2 --> [-2.549 -1 0.549 ]
1247  // [ 3.098 0 -3.098 ]
1248  // [-0.549 1 2.549 ]
1249  // Q3 --> [-4.315 -1.03 0.5 -0.44 ]
1250  // [ 6.07 -1.44 -2.97 2.196 ]
1251  // [-2.196 2.97 1.44 -6.07 ]
1252  // [ 0.44 -0.5 1.03 4.315 ]
1253  // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ]
1254  // [10.111 -2.76 -2.667 2.066 -2.306 ]
1255  // [-5.688 5.773 0 -5.773 5.688 ]
1256  // [ 2.306 -2.066 2.667 2.76 -10.111 ]
1257  // [-0.413 0.353 -0.333 -0.353 0.413 ]
1258  //
1259  // In these matrices, we want to use avoid computations involving
1260  // zeros and ones and in addition use the symmetry in entries to
1261  // reduce the number of read operations.
1262  template <int dim,
1263  int n_rows,
1264  int n_columns,
1265  typename Number,
1266  typename Number2>
1267  template <int direction, bool contract_over_rows, bool add>
1268  inline void
1270  dim,
1271  n_rows,
1272  n_columns,
1273  Number,
1274  Number2>::gradients(const Number in[],
1275  Number out[]) const
1276  {
1277  Assert(shape_gradients != nullptr, ExcNotInitialized());
1278  AssertIndexRange(direction, dim);
1279  constexpr int mm = contract_over_rows ? n_rows : n_columns,
1280  nn = contract_over_rows ? n_columns : n_rows;
1281  constexpr int n_cols = nn / 2;
1282  constexpr int mid = mm / 2;
1283 
1284  constexpr int stride = Utilities::pow(n_columns, direction);
1285  constexpr int n_blocks1 = stride;
1286  constexpr int n_blocks2 =
1287  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1288 
1289  for (int i2 = 0; i2 < n_blocks2; ++i2)
1290  {
1291  for (int i1 = 0; i1 < n_blocks1; ++i1)
1292  {
1293  for (int col = 0; col < n_cols; ++col)
1294  {
1295  Number2 val0, val1;
1296  Number in0, in1, res0, res1;
1297  if (contract_over_rows == true)
1298  {
1299  val0 = shape_gradients[col];
1300  val1 = shape_gradients[nn - 1 - col];
1301  }
1302  else
1303  {
1304  val0 = shape_gradients[col * n_columns];
1305  val1 = shape_gradients[(nn - col - 1) * n_columns];
1306  }
1307  if (mid > 0)
1308  {
1309  in0 = in[0];
1310  in1 = in[stride * (mm - 1)];
1311  res0 = val0 * in0;
1312  res1 = val1 * in0;
1313  res0 -= val1 * in1;
1314  res1 -= val0 * in1;
1315  for (int ind = 1; ind < mid; ++ind)
1316  {
1317  if (contract_over_rows == true)
1318  {
1319  val0 = shape_gradients[ind * n_columns + col];
1320  val1 =
1321  shape_gradients[ind * n_columns + nn - 1 - col];
1322  }
1323  else
1324  {
1325  val0 = shape_gradients[col * n_columns + ind];
1326  val1 =
1327  shape_gradients[(nn - col - 1) * n_columns + ind];
1328  }
1329  in0 = in[stride * ind];
1330  in1 = in[stride * (mm - 1 - ind)];
1331  res0 += val0 * in0;
1332  res1 += val1 * in0;
1333  res0 -= val1 * in1;
1334  res1 -= val0 * in1;
1335  }
1336  }
1337  else
1338  res0 = res1 = Number();
1339  if (mm % 2 == 1)
1340  {
1341  if (contract_over_rows == true)
1342  val0 = shape_gradients[mid * n_columns + col];
1343  else
1344  val0 = shape_gradients[col * n_columns + mid];
1345  in1 = val0 * in[stride * mid];
1346  res0 += in1;
1347  res1 -= in1;
1348  }
1349  if (add == false)
1350  {
1351  out[stride * col] = res0;
1352  out[stride * (nn - 1 - col)] = res1;
1353  }
1354  else
1355  {
1356  out[stride * col] += res0;
1357  out[stride * (nn - 1 - col)] += res1;
1358  }
1359  }
1360  if (nn % 2 == 1)
1361  {
1362  Number2 val0;
1363  Number res0;
1364  if (contract_over_rows == true)
1365  val0 = shape_gradients[n_cols];
1366  else
1367  val0 = shape_gradients[n_cols * n_columns];
1368  res0 = val0 * (in[0] - in[stride * (mm - 1)]);
1369  for (int ind = 1; ind < mid; ++ind)
1370  {
1371  if (contract_over_rows == true)
1372  val0 = shape_gradients[ind * n_columns + n_cols];
1373  else
1374  val0 = shape_gradients[n_cols * n_columns + ind];
1375  Number in1 =
1376  val0 * (in[stride * ind] - in[stride * (mm - 1 - ind)]);
1377  res0 += in1;
1378  }
1379  if (add == false)
1380  out[stride * n_cols] = res0;
1381  else
1382  out[stride * n_cols] += res0;
1383  }
1384 
1385  ++in;
1386  ++out;
1387  }
1388  in += stride * (mm - 1);
1389  out += stride * (nn - 1);
1390  }
1391  }
1392 
1393 
1394 
1395  // evaluates the given shape data in 1d-3d using the tensor product
1396  // form assuming the symmetries of unit cell shape hessians for
1397  // finite elements in FEEvaluation
1398  template <int dim,
1399  int n_rows,
1400  int n_columns,
1401  typename Number,
1402  typename Number2>
1403  template <int direction, bool contract_over_rows, bool add>
1404  inline void
1406  dim,
1407  n_rows,
1408  n_columns,
1409  Number,
1410  Number2>::hessians(const Number in[],
1411  Number out[]) const
1412  {
1413  Assert(shape_hessians != nullptr, ExcNotInitialized());
1414  AssertIndexRange(direction, dim);
1415  constexpr int mm = contract_over_rows ? n_rows : n_columns;
1416  constexpr int nn = contract_over_rows ? n_columns : n_rows;
1417  constexpr int n_cols = nn / 2;
1418  constexpr int mid = mm / 2;
1419 
1420  constexpr int stride = Utilities::pow(n_columns, direction);
1421  constexpr int n_blocks1 = stride;
1422  constexpr int n_blocks2 =
1423  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1424 
1425  for (int i2 = 0; i2 < n_blocks2; ++i2)
1426  {
1427  for (int i1 = 0; i1 < n_blocks1; ++i1)
1428  {
1429  for (int col = 0; col < n_cols; ++col)
1430  {
1431  Number2 val0, val1;
1432  Number in0, in1, res0, res1;
1433  if (contract_over_rows == true)
1434  {
1435  val0 = shape_hessians[col];
1436  val1 = shape_hessians[nn - 1 - col];
1437  }
1438  else
1439  {
1440  val0 = shape_hessians[col * n_columns];
1441  val1 = shape_hessians[(col + 1) * n_columns - 1];
1442  }
1443  if (mid > 0)
1444  {
1445  in0 = in[0];
1446  in1 = in[stride * (mm - 1)];
1447  res0 = val0 * in0;
1448  res1 = val1 * in0;
1449  res0 += val1 * in1;
1450  res1 += val0 * in1;
1451  for (int ind = 1; ind < mid; ++ind)
1452  {
1453  if (contract_over_rows == true)
1454  {
1455  val0 = shape_hessians[ind * n_columns + col];
1456  val1 =
1457  shape_hessians[ind * n_columns + nn - 1 - col];
1458  }
1459  else
1460  {
1461  val0 = shape_hessians[col * n_columns + ind];
1462  val1 =
1463  shape_hessians[(col + 1) * n_columns - 1 - ind];
1464  }
1465  in0 = in[stride * ind];
1466  in1 = in[stride * (mm - 1 - ind)];
1467  res0 += val0 * in0;
1468  res1 += val1 * in0;
1469  res0 += val1 * in1;
1470  res1 += val0 * in1;
1471  }
1472  }
1473  else
1474  res0 = res1 = Number();
1475  if (mm % 2 == 1)
1476  {
1477  if (contract_over_rows == true)
1478  val0 = shape_hessians[mid * n_columns + col];
1479  else
1480  val0 = shape_hessians[col * n_columns + mid];
1481  in1 = val0 * in[stride * mid];
1482  res0 += in1;
1483  res1 += in1;
1484  }
1485  if (add == false)
1486  {
1487  out[stride * col] = res0;
1488  out[stride * (nn - 1 - col)] = res1;
1489  }
1490  else
1491  {
1492  out[stride * col] += res0;
1493  out[stride * (nn - 1 - col)] += res1;
1494  }
1495  }
1496  if (nn % 2 == 1)
1497  {
1498  Number2 val0;
1499  Number res0;
1500  if (contract_over_rows == true)
1501  val0 = shape_hessians[n_cols];
1502  else
1503  val0 = shape_hessians[n_cols * n_columns];
1504  if (mid > 0)
1505  {
1506  res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1507  for (int ind = 1; ind < mid; ++ind)
1508  {
1509  if (contract_over_rows == true)
1510  val0 = shape_hessians[ind * n_columns + n_cols];
1511  else
1512  val0 = shape_hessians[n_cols * n_columns + ind];
1513  Number in1 = val0 * (in[stride * ind] +
1514  in[stride * (mm - 1 - ind)]);
1515  res0 += in1;
1516  }
1517  }
1518  else
1519  res0 = Number();
1520  if (mm % 2 == 1)
1521  {
1522  if (contract_over_rows == true)
1523  val0 = shape_hessians[mid * n_columns + n_cols];
1524  else
1525  val0 = shape_hessians[n_cols * n_columns + mid];
1526  res0 += val0 * in[stride * mid];
1527  }
1528  if (add == false)
1529  out[stride * n_cols] = res0;
1530  else
1531  out[stride * n_cols] += res0;
1532  }
1533 
1534  ++in;
1535  ++out;
1536  }
1537  in += stride * (mm - 1);
1538  out += stride * (nn - 1);
1539  }
1540  }
1541 
1542 
1543 
1575  template <int dim,
1576  int n_rows,
1577  int n_columns,
1578  typename Number,
1579  typename Number2>
1581  dim,
1582  n_rows,
1583  n_columns,
1584  Number,
1585  Number2>
1586  {
1587  static constexpr unsigned int n_rows_of_product =
1588  Utilities::pow(n_rows, dim);
1589  static constexpr unsigned int n_columns_of_product =
1590  Utilities::pow(n_columns, dim);
1591 
1598  : shape_values(nullptr)
1599  , shape_gradients(nullptr)
1600  , shape_hessians(nullptr)
1601  {}
1602 
1608  : shape_values(shape_values.begin())
1609  , shape_gradients(nullptr)
1610  , shape_hessians(nullptr)
1611  {
1612  AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
1613  }
1614 
1620  const AlignedVector<Number2> &shape_gradients,
1621  const AlignedVector<Number2> &shape_hessians,
1622  const unsigned int dummy1 = 0,
1623  const unsigned int dummy2 = 0)
1624  : shape_values(shape_values.begin())
1625  , shape_gradients(shape_gradients.begin())
1626  , shape_hessians(shape_hessians.begin())
1627  {
1628  // In this function, we allow for dummy pointers if some of values,
1629  // gradients or hessians should not be computed
1630  if (!shape_values.empty())
1631  AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
1632  if (!shape_gradients.empty())
1633  AssertDimension(shape_gradients.size(), n_rows * ((n_columns + 1) / 2));
1634  if (!shape_hessians.empty())
1635  AssertDimension(shape_hessians.size(), n_rows * ((n_columns + 1) / 2));
1636  (void)dummy1;
1637  (void)dummy2;
1638  }
1639 
1640  template <int direction, bool contract_over_rows, bool add>
1641  void
1642  values(const Number in[], Number out[]) const
1643  {
1644  Assert(shape_values != nullptr, ExcNotInitialized());
1645  apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
1646  }
1647 
1648  template <int direction, bool contract_over_rows, bool add>
1649  void
1650  gradients(const Number in[], Number out[]) const
1651  {
1652  Assert(shape_gradients != nullptr, ExcNotInitialized());
1653  apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
1654  }
1655 
1656  template <int direction, bool contract_over_rows, bool add>
1657  void
1658  hessians(const Number in[], Number out[]) const
1659  {
1660  Assert(shape_hessians != nullptr, ExcNotInitialized());
1661  apply<direction, contract_over_rows, add, 2>(shape_hessians, in, out);
1662  }
1663 
1664  template <int direction, bool contract_over_rows, bool add>
1665  void
1666  values_one_line(const Number in[], Number out[]) const
1667  {
1668  Assert(shape_values != nullptr, ExcNotInitialized());
1669  apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
1670  }
1671 
1672  template <int direction, bool contract_over_rows, bool add>
1673  void
1674  gradients_one_line(const Number in[], Number out[]) const
1675  {
1676  Assert(shape_gradients != nullptr, ExcNotInitialized());
1677  apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
1678  in,
1679  out);
1680  }
1681 
1682  template <int direction, bool contract_over_rows, bool add>
1683  void
1684  hessians_one_line(const Number in[], Number out[]) const
1685  {
1686  Assert(shape_hessians != nullptr, ExcNotInitialized());
1687  apply<direction, contract_over_rows, add, 2, true>(shape_hessians,
1688  in,
1689  out);
1690  }
1691 
1720  template <int direction,
1721  bool contract_over_rows,
1722  bool add,
1723  int type,
1724  bool one_line = false>
1725  static void
1726  apply(const Number2 *DEAL_II_RESTRICT shape_data,
1727  const Number * in,
1728  Number * out);
1729 
1730  const Number2 *shape_values;
1731  const Number2 *shape_gradients;
1732  const Number2 *shape_hessians;
1733  };
1734 
1735 
1736 
1737  template <int dim,
1738  int n_rows,
1739  int n_columns,
1740  typename Number,
1741  typename Number2>
1742  template <int direction,
1743  bool contract_over_rows,
1744  bool add,
1745  int type,
1746  bool one_line>
1747  inline void
1749  dim,
1750  n_rows,
1751  n_columns,
1752  Number,
1753  Number2>::apply(const Number2 *DEAL_II_RESTRICT shapes,
1754  const Number * in,
1755  Number * out)
1756  {
1757  static_assert(type < 3, "Only three variants type=0,1,2 implemented");
1758  static_assert(one_line == false || direction == dim - 1,
1759  "Single-line evaluation only works for direction=dim-1.");
1760  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
1761  in != out,
1762  ExcMessage("In-place operation only supported for "
1763  "n_rows==n_columns or single-line interpolation"));
1764 
1765  // We cannot statically assert that direction is less than dim, so must do
1766  // an additional dynamic check
1767  AssertIndexRange(direction, dim);
1768 
1769  constexpr int nn = contract_over_rows ? n_columns : n_rows;
1770  constexpr int mm = contract_over_rows ? n_rows : n_columns;
1771  constexpr int n_cols = nn / 2;
1772  constexpr int mid = mm / 2;
1773 
1774  constexpr int stride = Utilities::pow(n_columns, direction);
1775  constexpr int n_blocks1 = one_line ? 1 : stride;
1776  constexpr int n_blocks2 =
1777  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1778 
1779  constexpr int offset = (n_columns + 1) / 2;
1780 
1781  // this code may look very inefficient at first sight due to the many
1782  // different cases with if's at the innermost loop part, but all of the
1783  // conditionals can be evaluated at compile time because they are
1784  // templates, so the compiler should optimize everything away
1785  for (int i2 = 0; i2 < n_blocks2; ++i2)
1786  {
1787  for (int i1 = 0; i1 < n_blocks1; ++i1)
1788  {
1789  Number xp[mid > 0 ? mid : 1], xm[mid > 0 ? mid : 1];
1790  for (int i = 0; i < mid; ++i)
1791  {
1792  if (contract_over_rows == true && type == 1)
1793  {
1794  xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
1795  xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
1796  }
1797  else
1798  {
1799  xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
1800  xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
1801  }
1802  }
1803  Number xmid = in[stride * mid];
1804  for (int col = 0; col < n_cols; ++col)
1805  {
1806  Number r0, r1;
1807  if (mid > 0)
1808  {
1809  if (contract_over_rows == true)
1810  {
1811  r0 = shapes[col] * xp[0];
1812  r1 = shapes[(n_rows - 1) * offset + col] * xm[0];
1813  }
1814  else
1815  {
1816  r0 = shapes[col * offset] * xp[0];
1817  r1 = shapes[(n_rows - 1 - col) * offset] * xm[0];
1818  }
1819  for (int ind = 1; ind < mid; ++ind)
1820  {
1821  if (contract_over_rows == true)
1822  {
1823  r0 += shapes[ind * offset + col] * xp[ind];
1824  r1 += shapes[(n_rows - 1 - ind) * offset + col] *
1825  xm[ind];
1826  }
1827  else
1828  {
1829  r0 += shapes[col * offset + ind] * xp[ind];
1830  r1 += shapes[(n_rows - 1 - col) * offset + ind] *
1831  xm[ind];
1832  }
1833  }
1834  }
1835  else
1836  r0 = r1 = Number();
1837  if (mm % 2 == 1 && contract_over_rows == true)
1838  {
1839  if (type == 1)
1840  r1 += shapes[mid * offset + col] * xmid;
1841  else
1842  r0 += shapes[mid * offset + col] * xmid;
1843  }
1844  else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0 || mm == 3))
1845  r0 += shapes[col * offset + mid] * xmid;
1846 
1847  if (add == false)
1848  {
1849  out[stride * col] = r0 + r1;
1850  if (type == 1 && contract_over_rows == false)
1851  out[stride * (nn - 1 - col)] = r1 - r0;
1852  else
1853  out[stride * (nn - 1 - col)] = r0 - r1;
1854  }
1855  else
1856  {
1857  out[stride * col] += r0 + r1;
1858  if (type == 1 && contract_over_rows == false)
1859  out[stride * (nn - 1 - col)] += r1 - r0;
1860  else
1861  out[stride * (nn - 1 - col)] += r0 - r1;
1862  }
1863  }
1864  if (type == 0 && contract_over_rows == true && nn % 2 == 1 &&
1865  mm % 2 == 1 && mm > 3)
1866  {
1867  if (add == false)
1868  out[stride * n_cols] = shapes[mid * offset + n_cols] * xmid;
1869  else
1870  out[stride * n_cols] += shapes[mid * offset + n_cols] * xmid;
1871  }
1872  else if (contract_over_rows == true && nn % 2 == 1)
1873  {
1874  Number r0;
1875  if (mid > 0)
1876  {
1877  r0 = shapes[n_cols] * xp[0];
1878  for (int ind = 1; ind < mid; ++ind)
1879  r0 += shapes[ind * offset + n_cols] * xp[ind];
1880  }
1881  else
1882  r0 = Number();
1883  if (type != 1 && mm % 2 == 1)
1884  r0 += shapes[mid * offset + n_cols] * xmid;
1885 
1886  if (add == false)
1887  out[stride * n_cols] = r0;
1888  else
1889  out[stride * n_cols] += r0;
1890  }
1891  else if (contract_over_rows == false && nn % 2 == 1)
1892  {
1893  Number r0;
1894  if (mid > 0)
1895  {
1896  if (type == 1)
1897  {
1898  r0 = shapes[n_cols * offset] * xm[0];
1899  for (int ind = 1; ind < mid; ++ind)
1900  r0 += shapes[n_cols * offset + ind] * xm[ind];
1901  }
1902  else
1903  {
1904  r0 = shapes[n_cols * offset] * xp[0];
1905  for (int ind = 1; ind < mid; ++ind)
1906  r0 += shapes[n_cols * offset + ind] * xp[ind];
1907  }
1908  }
1909  else
1910  r0 = Number();
1911 
1912  if ((type == 0 || type == 2) && mm % 2 == 1)
1913  r0 += shapes[n_cols * offset + mid] * xmid;
1914 
1915  if (add == false)
1916  out[stride * n_cols] = r0;
1917  else
1918  out[stride * n_cols] += r0;
1919  }
1920  if (one_line == false)
1921  {
1922  in += 1;
1923  out += 1;
1924  }
1925  }
1926  if (one_line == false)
1927  {
1928  in += stride * (mm - 1);
1929  out += stride * (nn - 1);
1930  }
1931  }
1932  }
1933 
1934 
1935 
1964  template <int dim,
1965  int n_rows,
1966  int n_columns,
1967  typename Number,
1968  typename Number2>
1970  dim,
1971  n_rows,
1972  n_columns,
1973  Number,
1974  Number2>
1975  {
1976  static constexpr unsigned int n_rows_of_product =
1977  Utilities::pow(n_rows, dim);
1978  static constexpr unsigned int n_columns_of_product =
1979  Utilities::pow(n_columns, dim);
1980 
1987  : shape_values(nullptr)
1988  , shape_gradients(nullptr)
1989  , shape_hessians(nullptr)
1990  {}
1991 
1997  : shape_values(shape_values.begin())
1998  , shape_gradients(nullptr)
1999  , shape_hessians(nullptr)
2000  {}
2001 
2007  const AlignedVector<Number2> &shape_gradients,
2008  const AlignedVector<Number2> &shape_hessians,
2009  const unsigned int dummy1 = 0,
2010  const unsigned int dummy2 = 0)
2011  : shape_values(shape_values.begin())
2012  , shape_gradients(shape_gradients.begin())
2013  , shape_hessians(shape_hessians.begin())
2014  {
2015  (void)dummy1;
2016  (void)dummy2;
2017  }
2018 
2019  template <int direction, bool contract_over_rows, bool add>
2020  void
2021  values(const Number in[], Number out[]) const
2022  {
2023  Assert(shape_values != nullptr, ExcNotInitialized());
2024  apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
2025  }
2026 
2027  template <int direction, bool contract_over_rows, bool add>
2028  void
2029  gradients(const Number in[], Number out[]) const
2030  {
2031  Assert(shape_gradients != nullptr, ExcNotInitialized());
2032  apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
2033  }
2034 
2035  template <int direction, bool contract_over_rows, bool add>
2036  void
2037  hessians(const Number in[], Number out[]) const
2038  {
2039  Assert(shape_hessians != nullptr, ExcNotInitialized());
2040  apply<direction, contract_over_rows, add, 0>(shape_hessians, in, out);
2041  }
2042 
2043  template <int direction, bool contract_over_rows, bool add>
2044  void
2045  values_one_line(const Number in[], Number out[]) const
2046  {
2047  Assert(shape_values != nullptr, ExcNotInitialized());
2048  apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
2049  }
2050 
2051  template <int direction, bool contract_over_rows, bool add>
2052  void
2053  gradients_one_line(const Number in[], Number out[]) const
2054  {
2055  Assert(shape_gradients != nullptr, ExcNotInitialized());
2056  apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
2057  in,
2058  out);
2059  }
2060 
2061  template <int direction, bool contract_over_rows, bool add>
2062  void
2063  hessians_one_line(const Number in[], Number out[]) const
2064  {
2065  Assert(shape_hessians != nullptr, ExcNotInitialized());
2066  apply<direction, contract_over_rows, add, 0, true>(shape_hessians,
2067  in,
2068  out);
2069  }
2070 
2098  template <int direction,
2099  bool contract_over_rows,
2100  bool add,
2101  int type,
2102  bool one_line = false>
2103  static void
2104  apply(const Number2 *DEAL_II_RESTRICT shape_data,
2105  const Number * in,
2106  Number * out);
2107 
2108  const Number2 *shape_values;
2109  const Number2 *shape_gradients;
2110  const Number2 *shape_hessians;
2111  };
2112 
2113 
2114 
2115  template <int dim,
2116  int n_rows,
2117  int n_columns,
2118  typename Number,
2119  typename Number2>
2120  template <int direction,
2121  bool contract_over_rows,
2122  bool add,
2123  int type,
2124  bool one_line>
2125  inline void
2127  dim,
2128  n_rows,
2129  n_columns,
2130  Number,
2131  Number2>::apply(const Number2 *DEAL_II_RESTRICT shapes,
2132  const Number * in,
2133  Number * out)
2134  {
2135  static_assert(one_line == false || direction == dim - 1,
2136  "Single-line evaluation only works for direction=dim-1.");
2137  static_assert(
2138  type == 0 || type == 1,
2139  "Only types 0 and 1 implemented for evaluate_symmetric_hierarchical.");
2140  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
2141  in != out,
2142  ExcMessage("In-place operation only supported for "
2143  "n_rows==n_columns or single-line interpolation"));
2144 
2145  // We cannot statically assert that direction is less than dim, so must do
2146  // an additional dynamic check
2147  AssertIndexRange(direction, dim);
2148 
2149  constexpr int nn = contract_over_rows ? n_columns : n_rows;
2150  constexpr int mm = contract_over_rows ? n_rows : n_columns;
2151  constexpr int n_cols = nn / 2;
2152  constexpr int mid = mm / 2;
2153 
2154  constexpr int stride = Utilities::pow(n_columns, direction);
2155  constexpr int n_blocks1 = one_line ? 1 : stride;
2156  constexpr int n_blocks2 =
2157  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
2158 
2159  // this code may look very inefficient at first sight due to the many
2160  // different cases with if's at the innermost loop part, but all of the
2161  // conditionals can be evaluated at compile time because they are
2162  // templates, so the compiler should optimize everything away
2163  for (int i2 = 0; i2 < n_blocks2; ++i2)
2164  {
2165  for (int i1 = 0; i1 < n_blocks1; ++i1)
2166  {
2167  if (contract_over_rows)
2168  {
2169  Number x[mm];
2170  for (unsigned int i = 0; i < mm; ++i)
2171  x[i] = in[stride * i];
2172  for (unsigned int col = 0; col < n_cols; ++col)
2173  {
2174  Number r0, r1;
2175  if (mid > 0)
2176  {
2177  r0 = shapes[col] * x[0];
2178  r1 = shapes[col + n_columns] * x[1];
2179  for (unsigned int ind = 1; ind < mid; ++ind)
2180  {
2181  r0 +=
2182  shapes[col + 2 * ind * n_columns] * x[2 * ind];
2183  r1 += shapes[col + (2 * ind + 1) * n_columns] *
2184  x[2 * ind + 1];
2185  }
2186  }
2187  else
2188  r0 = r1 = Number();
2189  if (mm % 2 == 1)
2190  r0 += shapes[col + (mm - 1) * n_columns] * x[mm - 1];
2191  if (add == false)
2192  {
2193  out[stride * col] = r0 + r1;
2194  if (type == 1)
2195  out[stride * (nn - 1 - col)] = r1 - r0;
2196  else
2197  out[stride * (nn - 1 - col)] = r0 - r1;
2198  }
2199  else
2200  {
2201  out[stride * col] += r0 + r1;
2202  if (type == 1)
2203  out[stride * (nn - 1 - col)] += r1 - r0;
2204  else
2205  out[stride * (nn - 1 - col)] += r0 - r1;
2206  }
2207  }
2208  if (nn % 2 == 1)
2209  {
2210  Number r0;
2211  const unsigned int shift = type == 1 ? 1 : 0;
2212  if (mid > 0)
2213  {
2214  r0 = shapes[n_cols + shift * n_columns] * x[shift];
2215  for (unsigned int ind = 1; ind < mid; ++ind)
2216  r0 += shapes[n_cols + (2 * ind + shift) * n_columns] *
2217  x[2 * ind + shift];
2218  }
2219  else
2220  r0 = 0;
2221  if (type != 1 && mm % 2 == 1)
2222  r0 += shapes[n_cols + (mm - 1) * n_columns] * x[mm - 1];
2223  if (add == false)
2224  out[stride * n_cols] = r0;
2225  else
2226  out[stride * n_cols] += r0;
2227  }
2228  }
2229  else
2230  {
2231  Number xp[mid + 1], xm[mid > 0 ? mid : 1];
2232  for (int i = 0; i < mid; ++i)
2233  if (type == 0)
2234  {
2235  xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
2236  xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
2237  }
2238  else
2239  {
2240  xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
2241  xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
2242  }
2243  if (mm % 2 == 1)
2244  xp[mid] = in[stride * mid];
2245  for (unsigned int col = 0; col < n_cols; ++col)
2246  {
2247  Number r0, r1;
2248  if (mid > 0)
2249  {
2250  r0 = shapes[2 * col * n_columns] * xp[0];
2251  r1 = shapes[(2 * col + 1) * n_columns] * xm[0];
2252  for (unsigned int ind = 1; ind < mid; ++ind)
2253  {
2254  r0 += shapes[2 * col * n_columns + ind] * xp[ind];
2255  r1 +=
2256  shapes[(2 * col + 1) * n_columns + ind] * xm[ind];
2257  }
2258  }
2259  else
2260  r0 = r1 = Number();
2261  if (mm % 2 == 1)
2262  {
2263  if (type == 1)
2264  r1 +=
2265  shapes[(2 * col + 1) * n_columns + mid] * xp[mid];
2266  else
2267  r0 += shapes[2 * col * n_columns + mid] * xp[mid];
2268  }
2269  if (add == false)
2270  {
2271  out[stride * (2 * col)] = r0;
2272  out[stride * (2 * col + 1)] = r1;
2273  }
2274  else
2275  {
2276  out[stride * (2 * col)] += r0;
2277  out[stride * (2 * col + 1)] += r1;
2278  }
2279  }
2280  if (nn % 2 == 1)
2281  {
2282  Number r0;
2283  if (mid > 0)
2284  {
2285  r0 = shapes[(nn - 1) * n_columns] * xp[0];
2286  for (unsigned int ind = 1; ind < mid; ++ind)
2287  r0 += shapes[(nn - 1) * n_columns + ind] * xp[ind];
2288  }
2289  else
2290  r0 = Number();
2291  if (mm % 2 == 1 && type == 0)
2292  r0 += shapes[(nn - 1) * n_columns + mid] * xp[mid];
2293  if (add == false)
2294  out[stride * (nn - 1)] = r0;
2295  else
2296  out[stride * (nn - 1)] += r0;
2297  }
2298  }
2299  if (one_line == false)
2300  {
2301  in += 1;
2302  out += 1;
2303  }
2304  }
2305  if (one_line == false)
2306  {
2307  in += stride * (mm - 1);
2308  out += stride * (nn - 1);
2309  }
2310  }
2311  }
2312 
2313 } // end of namespace internal
2314 
2315 
2317 
2318 #endif
#define DEAL_II_RESTRICT
Definition: config.h:95
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1568
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)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1636
static ::ExceptionBase & ExcNotInitialized()
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:461
T fixed_power(const T t)
Definition: utilities.h:1045
static ::ExceptionBase & ExcMessage(std::string arg1)
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)
#define Assert(cond, exc)
Definition: exceptions.h:1411
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
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)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
VectorType::value_type * begin(VectorType &V)
size_type size() const
static ::ExceptionBase & ExcNotImplemented()
bool empty() const
T max(const T &t, const MPI_Comm &mpi_communicator)
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:815