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