Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
22 #include <deal.II/base/aligned_vector.h>
23 #include <deal.II/base/utilities.h>
24 
25 
26 DEAL_II_NAMESPACE_OPEN
27 
28 
29 
30 namespace internal
31 {
36  enum EvaluatorVariant
37  {
42  evaluate_general,
48  evaluate_symmetric,
54  evaluate_evenodd,
65  evaluate_symmetric_hierarchical
66  };
67 
68 
69 
90  template <EvaluatorVariant variant,
91  int dim,
92  int n_rows,
93  int n_columns,
94  typename Number,
95  typename Number2 = Number>
96  struct EvaluatorTensorProduct
97  {};
98 
99 
100 
118  template <int dim,
119  int n_rows,
120  int n_columns,
121  typename Number,
122  typename Number2>
123  struct EvaluatorTensorProduct<evaluate_general,
124  dim,
125  n_rows,
126  n_columns,
127  Number,
128  Number2>
129  {
130  static constexpr unsigned int n_rows_of_product =
131  Utilities::pow(n_rows, dim);
132  static constexpr unsigned int n_columns_of_product =
133  Utilities::pow(n_columns, dim);
134 
139  EvaluatorTensorProduct()
140  : shape_values(nullptr)
141  , shape_gradients(nullptr)
142  , shape_hessians(nullptr)
143  {}
144 
148  EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
149  const AlignedVector<Number2> &shape_gradients,
150  const AlignedVector<Number2> &shape_hessians,
151  const unsigned int dummy1 = 0,
152  const unsigned int dummy2 = 0)
153  : shape_values(shape_values.begin())
154  , shape_gradients(shape_gradients.begin())
155  , shape_hessians(shape_hessians.begin())
156  {
157  // We can enter this function either for the apply() path that has
158  // n_rows * n_columns entries or for the apply_face() path that only has
159  // n_rows * 3 entries in the array. Since we cannot decide about the use
160  // we must allow for both here.
161  Assert(shape_values.size() == 0 ||
162  shape_values.size() == n_rows * n_columns ||
163  shape_values.size() == 3 * n_rows,
164  ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
165  Assert(shape_gradients.size() == 0 ||
166  shape_gradients.size() == n_rows * n_columns,
167  ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
168  Assert(shape_hessians.size() == 0 ||
169  shape_hessians.size() == n_rows * n_columns,
170  ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
171  (void)dummy1;
172  (void)dummy2;
173  }
174 
175  template <int direction, bool contract_over_rows, bool add>
176  void
177  values(const Number in[], Number out[]) const
178  {
179  apply<direction, contract_over_rows, add>(shape_values, in, out);
180  }
181 
182  template <int direction, bool contract_over_rows, bool add>
183  void
184  gradients(const Number in[], Number out[]) const
185  {
186  apply<direction, contract_over_rows, add>(shape_gradients, in, out);
187  }
188 
189  template <int direction, bool contract_over_rows, bool add>
190  void
191  hessians(const Number in[], Number out[]) const
192  {
193  apply<direction, contract_over_rows, add>(shape_hessians, in, out);
194  }
195 
196  template <int direction, bool contract_over_rows, bool add>
197  void
198  values_one_line(const Number in[], Number out[]) const
199  {
200  Assert(shape_values != nullptr, ExcNotInitialized());
201  apply<direction, contract_over_rows, add, true>(shape_values, in, out);
202  }
203 
204  template <int direction, bool contract_over_rows, bool add>
205  void
206  gradients_one_line(const Number in[], Number out[]) const
207  {
208  Assert(shape_gradients != nullptr, ExcNotInitialized());
209  apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
210  }
211 
212  template <int direction, bool contract_over_rows, bool add>
213  void
214  hessians_one_line(const Number in[], Number out[]) const
215  {
216  Assert(shape_hessians != nullptr, ExcNotInitialized());
217  apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
218  }
219 
244  template <int direction,
245  bool contract_over_rows,
246  bool add,
247  bool one_line = false>
248  static void
249  apply(const Number2 *DEAL_II_RESTRICT shape_data,
250  const Number * in,
251  Number * out);
252 
282  template <int face_direction,
283  bool contract_onto_face,
284  bool add,
285  int max_derivative>
286  void
287  apply_face(const Number *DEAL_II_RESTRICT in,
288  Number *DEAL_II_RESTRICT out) const;
289 
290  const Number2 *shape_values;
291  const Number2 *shape_gradients;
292  const Number2 *shape_hessians;
293  };
294 
295 
296 
297  template <int dim,
298  int n_rows,
299  int n_columns,
300  typename Number,
301  typename Number2>
302  template <int direction, bool contract_over_rows, bool add, bool one_line>
303  inline void
304  EvaluatorTensorProduct<evaluate_general,
305  dim,
306  n_rows,
307  n_columns,
308  Number,
309  Number2>::apply(const Number2 *DEAL_II_RESTRICT
310  shape_data,
311  const Number * in,
312  Number * out)
313  {
314  static_assert(one_line == false || direction == dim - 1,
315  "Single-line evaluation only works for direction=dim-1.");
316  Assert(shape_data != nullptr,
317  ExcMessage(
318  "The given array shape_data must not be the null pointer!"));
319  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
320  in != out,
321  ExcMessage("In-place operation only supported for "
322  "n_rows==n_columns or single-line interpolation"));
323  AssertIndexRange(direction, dim);
324  constexpr int mm = contract_over_rows ? n_rows : n_columns,
325  nn = contract_over_rows ? n_columns : n_rows;
326 
327  constexpr int stride = Utilities::pow(n_columns, direction);
328  constexpr int n_blocks1 = one_line ? 1 : stride;
329  constexpr int n_blocks2 =
330  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
331 
332  for (int i2 = 0; i2 < n_blocks2; ++i2)
333  {
334  for (int i1 = 0; i1 < n_blocks1; ++i1)
335  {
336  Number x[mm];
337  for (int i = 0; i < mm; ++i)
338  x[i] = in[stride * i];
339  for (int col = 0; col < nn; ++col)
340  {
341  Number2 val0;
342  if (contract_over_rows == true)
343  val0 = shape_data[col];
344  else
345  val0 = shape_data[col * n_columns];
346  Number res0 = val0 * x[0];
347  for (int i = 1; i < mm; ++i)
348  {
349  if (contract_over_rows == true)
350  val0 = shape_data[i * n_columns + col];
351  else
352  val0 = shape_data[col * n_columns + i];
353  res0 += val0 * x[i];
354  }
355  if (add == false)
356  out[stride * col] = res0;
357  else
358  out[stride * col] += res0;
359  }
360 
361  if (one_line == false)
362  {
363  ++in;
364  ++out;
365  }
366  }
367  if (one_line == false)
368  {
369  in += stride * (mm - 1);
370  out += stride * (nn - 1);
371  }
372  }
373  }
374 
375 
376 
377  template <int dim,
378  int n_rows,
379  int n_columns,
380  typename Number,
381  typename Number2>
382  template <int face_direction,
383  bool contract_onto_face,
384  bool add,
385  int max_derivative>
386  inline void
387  EvaluatorTensorProduct<evaluate_general,
388  dim,
389  n_rows,
390  n_columns,
391  Number,
392  Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
393  Number *DEAL_II_RESTRICT
394  out) const
395  {
396  static_assert(dim > 0 && dim < 4, "Only dim=1,2,3 supported");
397  static_assert(max_derivative >= 0 && max_derivative < 3,
398  "Only derivative orders 0-2 implemented");
399  Assert(shape_values != nullptr,
400  ExcMessage(
401  "The given array shape_values must not be the null pointer."));
402 
403  constexpr int n_blocks1 = dim > 1 ? n_rows : 1;
404  constexpr int n_blocks2 = dim > 2 ? n_rows : 1;
405 
406  AssertIndexRange(face_direction, dim);
407  constexpr int stride = Utilities::pow(n_rows, face_direction);
408  constexpr int out_stride = Utilities::pow(n_rows, dim - 1);
409  const Number *DEAL_II_RESTRICT shape_values = this->shape_values;
410 
411  for (int i2 = 0; i2 < n_blocks2; ++i2)
412  {
413  for (int i1 = 0; i1 < n_blocks1; ++i1)
414  {
415  if (contract_onto_face == true)
416  {
417  Number res0 = shape_values[0] * in[0];
418  Number res1, res2;
419  if (max_derivative > 0)
420  res1 = shape_values[n_rows] * in[0];
421  if (max_derivative > 1)
422  res2 = shape_values[2 * n_rows] * in[0];
423  for (int ind = 1; ind < n_rows; ++ind)
424  {
425  res0 += shape_values[ind] * in[stride * ind];
426  if (max_derivative > 0)
427  res1 += shape_values[ind + n_rows] * in[stride * ind];
428  if (max_derivative > 1)
429  res2 += shape_values[ind + 2 * n_rows] * in[stride * ind];
430  }
431  if (add == false)
432  {
433  out[0] = res0;
434  if (max_derivative > 0)
435  out[out_stride] = res1;
436  if (max_derivative > 1)
437  out[2 * out_stride] = res2;
438  }
439  else
440  {
441  out[0] += res0;
442  if (max_derivative > 0)
443  out[out_stride] += res1;
444  if (max_derivative > 1)
445  out[2 * out_stride] += res2;
446  }
447  }
448  else
449  {
450  for (int col = 0; col < n_rows; ++col)
451  {
452  if (add == false)
453  out[col * stride] = shape_values[col] * in[0];
454  else
455  out[col * stride] += shape_values[col] * in[0];
456  if (max_derivative > 0)
457  out[col * stride] +=
458  shape_values[col + n_rows] * in[out_stride];
459  if (max_derivative > 1)
460  out[col * stride] +=
461  shape_values[col + 2 * n_rows] * in[2 * out_stride];
462  }
463  }
464 
465  // increment: in regular case, just go to the next point in
466  // x-direction. If we are at the end of one chunk in x-dir, need
467  // to jump over to the next layer in z-direction
468  switch (face_direction)
469  {
470  case 0:
471  in += contract_onto_face ? n_rows : 1;
472  out += contract_onto_face ? 1 : n_rows;
473  break;
474  case 1:
475  ++in;
476  ++out;
477  // faces 2 and 3 in 3D use local coordinate system zx, which
478  // is the other way around compared to the tensor
479  // product. Need to take that into account.
480  if (dim == 3)
481  {
482  if (contract_onto_face)
483  out += n_rows - 1;
484  else
485  in += n_rows - 1;
486  }
487  break;
488  case 2:
489  ++in;
490  ++out;
491  break;
492  default:
493  Assert(false, ExcNotImplemented());
494  }
495  }
496  if (face_direction == 1 && dim == 3)
497  {
498  // adjust for local coordinate system zx
499  if (contract_onto_face)
500  {
501  in += n_rows * (n_rows - 1);
502  out -= n_rows * n_rows - 1;
503  }
504  else
505  {
506  out += n_rows * (n_rows - 1);
507  in -= n_rows * n_rows - 1;
508  }
509  }
510  }
511  }
512 
513 
514 
528  template <int dim, typename Number, typename Number2>
529  struct EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>
530  {
531  static constexpr unsigned int n_rows_of_product =
532  numbers::invalid_unsigned_int;
533  static constexpr unsigned int n_columns_of_product =
534  numbers::invalid_unsigned_int;
535 
540  EvaluatorTensorProduct()
541  : shape_values(nullptr)
542  , shape_gradients(nullptr)
543  , shape_hessians(nullptr)
544  , n_rows(numbers::invalid_unsigned_int)
545  , n_columns(numbers::invalid_unsigned_int)
546  {}
547 
551  EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
552  const AlignedVector<Number2> &shape_gradients,
553  const AlignedVector<Number2> &shape_hessians,
554  const unsigned int n_rows,
555  const unsigned int n_columns)
556  : shape_values(shape_values.begin())
557  , shape_gradients(shape_gradients.begin())
558  , shape_hessians(shape_hessians.begin())
559  , n_rows(n_rows)
560  , n_columns(n_columns)
561  {
562  // We can enter this function either for the apply() path that has
563  // n_rows * n_columns entries or for the apply_face() path that only has
564  // n_rows * 3 entries in the array. Since we cannot decide about the use
565  // we must allow for both here.
566  Assert(shape_values.size() == 0 ||
567  shape_values.size() == n_rows * n_columns ||
568  shape_values.size() == n_rows * 3,
569  ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
570  Assert(shape_gradients.size() == 0 ||
571  shape_gradients.size() == n_rows * n_columns,
572  ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
573  Assert(shape_hessians.size() == 0 ||
574  shape_hessians.size() == n_rows * n_columns,
575  ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
576  }
577 
578  template <int direction, bool contract_over_rows, bool add>
579  void
580  values(const Number *in, Number *out) const
581  {
582  apply<direction, contract_over_rows, add>(shape_values, in, out);
583  }
584 
585  template <int direction, bool contract_over_rows, bool add>
586  void
587  gradients(const Number *in, Number *out) const
588  {
589  apply<direction, contract_over_rows, add>(shape_gradients, in, out);
590  }
591 
592  template <int direction, bool contract_over_rows, bool add>
593  void
594  hessians(const Number *in, Number *out) const
595  {
596  apply<direction, contract_over_rows, add>(shape_hessians, in, out);
597  }
598 
599  template <int direction, bool contract_over_rows, bool add>
600  void
601  values_one_line(const Number in[], Number out[]) const
602  {
603  Assert(shape_values != nullptr, ExcNotInitialized());
604  apply<direction, contract_over_rows, add, true>(shape_values, in, out);
605  }
606 
607  template <int direction, bool contract_over_rows, bool add>
608  void
609  gradients_one_line(const Number in[], Number out[]) const
610  {
611  Assert(shape_gradients != nullptr, ExcNotInitialized());
612  apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
613  }
614 
615  template <int direction, bool contract_over_rows, bool add>
616  void
617  hessians_one_line(const Number in[], Number out[]) const
618  {
619  Assert(shape_hessians != nullptr, ExcNotInitialized());
620  apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
621  }
622 
623  template <int direction,
624  bool contract_over_rows,
625  bool add,
626  bool one_line = false>
627  void
628  apply(const Number2 *DEAL_II_RESTRICT shape_data,
629  const Number * in,
630  Number * out) const;
631 
632  template <int face_direction,
633  bool contract_onto_face,
634  bool add,
635  int max_derivative>
636  void
637  apply_face(const Number *DEAL_II_RESTRICT in,
638  Number *DEAL_II_RESTRICT out) const;
639 
640  const Number2 * shape_values;
641  const Number2 * shape_gradients;
642  const Number2 * shape_hessians;
643  const unsigned int n_rows;
644  const unsigned int n_columns;
645  };
646 
647 
648 
649  template <int dim, typename Number, typename Number2>
650  template <int direction, bool contract_over_rows, bool add, bool one_line>
651  inline void
652  EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>::apply(
653  const Number2 *DEAL_II_RESTRICT shape_data,
654  const Number * in,
655  Number * out) const
656  {
657  static_assert(one_line == false || direction == dim - 1,
658  "Single-line evaluation only works for direction=dim-1.");
659  Assert(shape_data != nullptr,
660  ExcMessage(
661  "The given array shape_data must not be the null pointer!"));
662  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
663  in != out,
664  ExcMessage("In-place operation only supported for "
665  "n_rows==n_columns or single-line interpolation"));
666  AssertIndexRange(direction, dim);
667  const int mm = contract_over_rows ? n_rows : n_columns,
668  nn = contract_over_rows ? n_columns : n_rows;
669 
670  const int stride =
671  direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
672  const int n_blocks1 = one_line ? 1 : stride;
673  const int n_blocks2 = direction >= dim - 1 ?
674  1 :
675  Utilities::fixed_power<dim - direction - 1>(n_rows);
676  Assert(n_rows <= 128, ExcNotImplemented());
677 
678  // specialization for n_rows = 2 that manually unrolls the innermost loop
679  // to make the operation perform better (not completely as good as the
680  // templated one, but much better than the generic version down below,
681  // because the loop over col can be more effectively unrolled by the
682  // compiler)
683  if (contract_over_rows && n_rows == 2)
684  {
685  const Number2 *shape_data_1 = shape_data + n_columns;
686  for (int i2 = 0; i2 < n_blocks2; ++i2)
687  {
688  for (int i1 = 0; i1 < n_blocks1; ++i1)
689  {
690  const Number x0 = in[0], x1 = in[stride];
691  for (int col = 0; col < nn; ++col)
692  {
693  const Number result =
694  shape_data[col] * x0 + shape_data_1[col] * x1;
695  if (add == false)
696  out[stride * col] = result;
697  else
698  out[stride * col] += result;
699  }
700 
701  if (one_line == false)
702  {
703  ++in;
704  ++out;
705  }
706  }
707  if (one_line == false)
708  {
709  in += stride * (mm - 1);
710  out += stride * (nn - 1);
711  }
712  }
713  }
714  // specialization for n = 3
715  else if (contract_over_rows && n_rows == 3)
716  {
717  const Number2 *shape_data_1 = shape_data + n_columns;
718  const Number2 *shape_data_2 = shape_data + 2 * n_columns;
719  for (int i2 = 0; i2 < n_blocks2; ++i2)
720  {
721  for (int i1 = 0; i1 < n_blocks1; ++i1)
722  {
723  const Number x0 = in[0], x1 = in[stride], x2 = in[2 * stride];
724  for (int col = 0; col < nn; ++col)
725  {
726  const Number result = shape_data[col] * x0 +
727  shape_data_1[col] * x1 +
728  shape_data_2[col] * x2;
729  if (add == false)
730  out[stride * col] = result;
731  else
732  out[stride * col] += result;
733  }
734 
735  if (one_line == false)
736  {
737  ++in;
738  ++out;
739  }
740  }
741  if (one_line == false)
742  {
743  in += stride * (mm - 1);
744  out += stride * (nn - 1);
745  }
746  }
747  }
748  // general loop for all other cases
749  else
750  for (int i2 = 0; i2 < n_blocks2; ++i2)
751  {
752  for (int i1 = 0; i1 < n_blocks1; ++i1)
753  {
754  Number x[129];
755  for (int i = 0; i < mm; ++i)
756  x[i] = in[stride * i];
757  for (int col = 0; col < nn; ++col)
758  {
759  Number2 val0;
760  if (contract_over_rows == true)
761  val0 = shape_data[col];
762  else
763  val0 = shape_data[col * n_columns];
764  Number res0 = val0 * x[0];
765  for (int i = 1; i < mm; ++i)
766  {
767  if (contract_over_rows == true)
768  val0 = shape_data[i * n_columns + col];
769  else
770  val0 = shape_data[col * n_columns + i];
771  res0 += val0 * x[i];
772  }
773  if (add == false)
774  out[stride * col] = res0;
775  else
776  out[stride * col] += res0;
777  }
778 
779  if (one_line == false)
780  {
781  ++in;
782  ++out;
783  }
784  }
785  if (one_line == false)
786  {
787  in += stride * (mm - 1);
788  out += stride * (nn - 1);
789  }
790  }
791  }
792 
793 
794 
795  template <int dim, typename Number, typename Number2>
796  template <int face_direction,
797  bool contract_onto_face,
798  bool add,
799  int max_derivative>
800  inline void
801  EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>::
802  apply_face(const Number *DEAL_II_RESTRICT in,
803  Number *DEAL_II_RESTRICT out) const
804  {
805  Assert(shape_values != nullptr,
806  ExcMessage(
807  "The given array shape_data must not be the null pointer!"));
808  static_assert(dim > 0 && dim < 4, "Only dim=1,2,3 supported");
809  const int n_blocks1 = dim > 1 ? n_rows : 1;
810  const int n_blocks2 = dim > 2 ? n_rows : 1;
811 
812  AssertIndexRange(face_direction, dim);
813  const int stride =
814  face_direction > 0 ? Utilities::fixed_power<face_direction>(n_rows) : 1;
815  const int out_stride =
816  dim > 1 ? Utilities::fixed_power<dim - 1>(n_rows) : 1;
817 
818  for (int i2 = 0; i2 < n_blocks2; ++i2)
819  {
820  for (int i1 = 0; i1 < n_blocks1; ++i1)
821  {
822  if (contract_onto_face == true)
823  {
824  Number res0 = shape_values[0] * in[0];
825  Number res1, res2;
826  if (max_derivative > 0)
827  res1 = shape_values[n_rows] * in[0];
828  if (max_derivative > 1)
829  res2 = shape_values[2 * n_rows] * in[0];
830  for (unsigned int ind = 1; ind < n_rows; ++ind)
831  {
832  res0 += shape_values[ind] * in[stride * ind];
833  if (max_derivative > 0)
834  res1 += shape_values[ind + n_rows] * in[stride * ind];
835  if (max_derivative > 1)
836  res2 += shape_values[ind + 2 * n_rows] * in[stride * ind];
837  }
838  if (add == false)
839  {
840  out[0] = res0;
841  if (max_derivative > 0)
842  out[out_stride] = res1;
843  if (max_derivative > 1)
844  out[2 * out_stride] = res2;
845  }
846  else
847  {
848  out[0] += res0;
849  if (max_derivative > 0)
850  out[out_stride] += res1;
851  if (max_derivative > 1)
852  out[2 * out_stride] += res2;
853  }
854  }
855  else
856  {
857  for (unsigned int col = 0; col < n_rows; ++col)
858  {
859  if (add == false)
860  out[col * stride] = shape_values[col] * in[0];
861  else
862  out[col * stride] += shape_values[col] * in[0];
863  if (max_derivative > 0)
864  out[col * stride] +=
865  shape_values[col + n_rows] * in[out_stride];
866  if (max_derivative > 1)
867  out[col * stride] +=
868  shape_values[col + 2 * n_rows] * in[2 * out_stride];
869  }
870  }
871 
872  // increment: in regular case, just go to the next point in
873  // x-direction. If we are at the end of one chunk in x-dir, need
874  // to jump over to the next layer in z-direction
875  switch (face_direction)
876  {
877  case 0:
878  in += contract_onto_face ? n_rows : 1;
879  out += contract_onto_face ? 1 : n_rows;
880  break;
881  case 1:
882  ++in;
883  ++out;
884  // faces 2 and 3 in 3D use local coordinate system zx, which
885  // is the other way around compared to the tensor
886  // product. Need to take that into account.
887  if (dim == 3)
888  {
889  if (contract_onto_face)
890  out += n_rows - 1;
891  else
892  in += n_rows - 1;
893  }
894  break;
895  case 2:
896  ++in;
897  ++out;
898  break;
899  default:
900  Assert(false, ExcNotImplemented());
901  }
902  }
903  if (face_direction == 1 && dim == 3)
904  {
905  // adjust for local coordinate system zx
906  if (contract_onto_face)
907  {
908  in += n_rows * (n_rows - 1);
909  out -= n_rows * n_rows - 1;
910  }
911  else
912  {
913  out += n_rows * (n_rows - 1);
914  in -= n_rows * n_rows - 1;
915  }
916  }
917  }
918  }
919 
920 
921 
942  template <int dim,
943  int n_rows,
944  int n_columns,
945  typename Number,
946  typename Number2>
947  struct EvaluatorTensorProduct<evaluate_symmetric,
948  dim,
949  n_rows,
950  n_columns,
951  Number,
952  Number2>
953  {
954  static constexpr unsigned int n_rows_of_product =
955  Utilities::pow(n_rows, dim);
956  static constexpr unsigned int n_columns_of_product =
957  Utilities::pow(n_columns, dim);
958 
962  EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
963  const AlignedVector<Number2> &shape_gradients,
964  const AlignedVector<Number2> &shape_hessians,
965  const unsigned int dummy1 = 0,
966  const unsigned int dummy2 = 0)
967  : shape_values(shape_values.begin())
968  , shape_gradients(shape_gradients.begin())
969  , shape_hessians(shape_hessians.begin())
970  {
971  Assert(shape_values.size() == 0 ||
972  shape_values.size() == n_rows * n_columns,
973  ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
974  Assert(shape_gradients.size() == 0 ||
975  shape_gradients.size() == n_rows * n_columns,
976  ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
977  Assert(shape_hessians.size() == 0 ||
978  shape_hessians.size() == n_rows * n_columns,
979  ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
980  (void)dummy1;
981  (void)dummy2;
982  }
983 
984  template <int direction, bool contract_over_rows, bool add>
985  void
986  values(const Number in[], Number out[]) const;
987 
988  template <int direction, bool contract_over_rows, bool add>
989  void
990  gradients(const Number in[], Number out[]) const;
991 
992  template <int direction, bool contract_over_rows, bool add>
993  void
994  hessians(const Number in[], Number out[]) const;
995 
996  const Number2 *shape_values;
997  const Number2 *shape_gradients;
998  const Number2 *shape_hessians;
999  };
1000 
1001 
1002 
1003  // In this case, the 1D shape values read (sorted lexicographically, rows
1004  // run over 1D dofs, columns over quadrature points):
1005  // Q2 --> [ 0.687 0 -0.087 ]
1006  // [ 0.4 1 0.4 ]
1007  // [-0.087 0 0.687 ]
1008  // Q3 --> [ 0.66 0.003 0.002 0.049 ]
1009  // [ 0.521 1.005 -0.01 -0.230 ]
1010  // [-0.230 -0.01 1.005 0.521 ]
1011  // [ 0.049 0.002 0.003 0.66 ]
1012  // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ]
1013  // [ 0.608 1.059 0 0.039 0.176 ]
1014  // [-0.409 -0.113 1 -0.113 -0.409 ]
1015  // [ 0.176 0.039 0 1.059 0.608 ]
1016  // [-0.032 -0.007 0 0.022 0.658 ]
1017  //
1018  // In these matrices, we want to use avoid computations involving zeros and
1019  // ones and in addition use the symmetry in entries to reduce the number of
1020  // read operations.
1021  template <int dim,
1022  int n_rows,
1023  int n_columns,
1024  typename Number,
1025  typename Number2>
1026  template <int direction, bool contract_over_rows, bool add>
1027  inline void
1028  EvaluatorTensorProduct<evaluate_symmetric,
1029  dim,
1030  n_rows,
1031  n_columns,
1032  Number,
1033  Number2>::values(const Number in[], Number out[]) const
1034  {
1035  Assert(shape_values != nullptr, ExcNotInitialized());
1036  AssertIndexRange(direction, dim);
1037  constexpr int mm = contract_over_rows ? n_rows : n_columns,
1038  nn = contract_over_rows ? n_columns : n_rows;
1039  constexpr int n_cols = nn / 2;
1040  constexpr int mid = mm / 2;
1041 
1042  constexpr int stride = Utilities::pow(n_columns, direction);
1043  constexpr int n_blocks1 = stride;
1044  constexpr int n_blocks2 =
1045  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1046 
1047  for (int i2 = 0; i2 < n_blocks2; ++i2)
1048  {
1049  for (int i1 = 0; i1 < n_blocks1; ++i1)
1050  {
1051  for (int col = 0; col < n_cols; ++col)
1052  {
1053  Number2 val0, val1;
1054  Number in0, in1, res0, res1;
1055  if (contract_over_rows == true)
1056  {
1057  val0 = shape_values[col];
1058  val1 = shape_values[nn - 1 - col];
1059  }
1060  else
1061  {
1062  val0 = shape_values[col * n_columns];
1063  val1 = shape_values[(col + 1) * n_columns - 1];
1064  }
1065  if (mid > 0)
1066  {
1067  in0 = in[0];
1068  in1 = in[stride * (mm - 1)];
1069  res0 = val0 * in0;
1070  res1 = val1 * in0;
1071  res0 += val1 * in1;
1072  res1 += val0 * in1;
1073  for (int ind = 1; ind < mid; ++ind)
1074  {
1075  if (contract_over_rows == true)
1076  {
1077  val0 = shape_values[ind * n_columns + col];
1078  val1 = shape_values[ind * n_columns + nn - 1 - col];
1079  }
1080  else
1081  {
1082  val0 = shape_values[col * n_columns + ind];
1083  val1 =
1084  shape_values[(col + 1) * n_columns - 1 - ind];
1085  }
1086  in0 = in[stride * ind];
1087  in1 = in[stride * (mm - 1 - ind)];
1088  res0 += val0 * in0;
1089  res1 += val1 * in0;
1090  res0 += val1 * in1;
1091  res1 += val0 * in1;
1092  }
1093  }
1094  else
1095  res0 = res1 = Number();
1096  if (contract_over_rows == true)
1097  {
1098  if (mm % 2 == 1)
1099  {
1100  val0 = shape_values[mid * n_columns + col];
1101  in1 = val0 * in[stride * mid];
1102  res0 += in1;
1103  res1 += in1;
1104  }
1105  }
1106  else
1107  {
1108  if (mm % 2 == 1 && nn % 2 == 0)
1109  {
1110  val0 = shape_values[col * n_columns + mid];
1111  in1 = val0 * in[stride * mid];
1112  res0 += in1;
1113  res1 += in1;
1114  }
1115  }
1116  if (add == false)
1117  {
1118  out[stride * col] = res0;
1119  out[stride * (nn - 1 - col)] = res1;
1120  }
1121  else
1122  {
1123  out[stride * col] += res0;
1124  out[stride * (nn - 1 - col)] += res1;
1125  }
1126  }
1127  if (contract_over_rows == true && nn % 2 == 1 && mm % 2 == 1)
1128  {
1129  if (add == false)
1130  out[stride * n_cols] = in[stride * mid];
1131  else
1132  out[stride * n_cols] += in[stride * mid];
1133  }
1134  else if (contract_over_rows == true && nn % 2 == 1)
1135  {
1136  Number res0;
1137  Number2 val0 = shape_values[n_cols];
1138  if (mid > 0)
1139  {
1140  res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1141  for (int ind = 1; ind < mid; ++ind)
1142  {
1143  val0 = shape_values[ind * n_columns + n_cols];
1144  res0 += val0 * (in[stride * ind] +
1145  in[stride * (mm - 1 - ind)]);
1146  }
1147  }
1148  else
1149  res0 = Number();
1150  if (add == false)
1151  out[stride * n_cols] = res0;
1152  else
1153  out[stride * n_cols] += res0;
1154  }
1155  else if (contract_over_rows == false && nn % 2 == 1)
1156  {
1157  Number res0;
1158  if (mid > 0)
1159  {
1160  Number2 val0 = shape_values[n_cols * n_columns];
1161  res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1162  for (int ind = 1; ind < mid; ++ind)
1163  {
1164  val0 = shape_values[n_cols * n_columns + ind];
1165  Number in1 = val0 * (in[stride * ind] +
1166  in[stride * (mm - 1 - ind)]);
1167  res0 += in1;
1168  }
1169  if (mm % 2)
1170  res0 += in[stride * mid];
1171  }
1172  else
1173  res0 = in[0];
1174  if (add == false)
1175  out[stride * n_cols] = res0;
1176  else
1177  out[stride * n_cols] += res0;
1178  }
1179 
1180  ++in;
1181  ++out;
1182  }
1183  in += stride * (mm - 1);
1184  out += stride * (nn - 1);
1185  }
1186  }
1187 
1188 
1189 
1190  // For the specialized loop used for the gradient computation in
1191  // here, the 1D shape values read (sorted lexicographically, rows
1192  // run over 1D dofs, columns over quadrature points):
1193  // Q2 --> [-2.549 -1 0.549 ]
1194  // [ 3.098 0 -3.098 ]
1195  // [-0.549 1 2.549 ]
1196  // Q3 --> [-4.315 -1.03 0.5 -0.44 ]
1197  // [ 6.07 -1.44 -2.97 2.196 ]
1198  // [-2.196 2.97 1.44 -6.07 ]
1199  // [ 0.44 -0.5 1.03 4.315 ]
1200  // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ]
1201  // [10.111 -2.76 -2.667 2.066 -2.306 ]
1202  // [-5.688 5.773 0 -5.773 5.688 ]
1203  // [ 2.306 -2.066 2.667 2.76 -10.111 ]
1204  // [-0.413 0.353 -0.333 -0.353 0.413 ]
1205  //
1206  // In these matrices, we want to use avoid computations involving
1207  // zeros and ones and in addition use the symmetry in entries to
1208  // reduce the number of read operations.
1209  template <int dim,
1210  int n_rows,
1211  int n_columns,
1212  typename Number,
1213  typename Number2>
1214  template <int direction, bool contract_over_rows, bool add>
1215  inline void
1216  EvaluatorTensorProduct<evaluate_symmetric,
1217  dim,
1218  n_rows,
1219  n_columns,
1220  Number,
1221  Number2>::gradients(const Number in[],
1222  Number out[]) const
1223  {
1224  Assert(shape_gradients != nullptr, ExcNotInitialized());
1225  AssertIndexRange(direction, dim);
1226  constexpr int mm = contract_over_rows ? n_rows : n_columns,
1227  nn = contract_over_rows ? n_columns : n_rows;
1228  constexpr int n_cols = nn / 2;
1229  constexpr int mid = mm / 2;
1230 
1231  constexpr int stride = Utilities::pow(n_columns, direction);
1232  constexpr int n_blocks1 = stride;
1233  constexpr int n_blocks2 =
1234  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1235 
1236  for (int i2 = 0; i2 < n_blocks2; ++i2)
1237  {
1238  for (int i1 = 0; i1 < n_blocks1; ++i1)
1239  {
1240  for (int col = 0; col < n_cols; ++col)
1241  {
1242  Number2 val0, val1;
1243  Number in0, in1, res0, res1;
1244  if (contract_over_rows == true)
1245  {
1246  val0 = shape_gradients[col];
1247  val1 = shape_gradients[nn - 1 - col];
1248  }
1249  else
1250  {
1251  val0 = shape_gradients[col * n_columns];
1252  val1 = shape_gradients[(nn - col - 1) * n_columns];
1253  }
1254  if (mid > 0)
1255  {
1256  in0 = in[0];
1257  in1 = in[stride * (mm - 1)];
1258  res0 = val0 * in0;
1259  res1 = val1 * in0;
1260  res0 -= val1 * in1;
1261  res1 -= val0 * in1;
1262  for (int ind = 1; ind < mid; ++ind)
1263  {
1264  if (contract_over_rows == true)
1265  {
1266  val0 = shape_gradients[ind * n_columns + col];
1267  val1 =
1268  shape_gradients[ind * n_columns + nn - 1 - col];
1269  }
1270  else
1271  {
1272  val0 = shape_gradients[col * n_columns + ind];
1273  val1 =
1274  shape_gradients[(nn - col - 1) * n_columns + ind];
1275  }
1276  in0 = in[stride * ind];
1277  in1 = in[stride * (mm - 1 - ind)];
1278  res0 += val0 * in0;
1279  res1 += val1 * in0;
1280  res0 -= val1 * in1;
1281  res1 -= val0 * in1;
1282  }
1283  }
1284  else
1285  res0 = res1 = Number();
1286  if (mm % 2 == 1)
1287  {
1288  if (contract_over_rows == true)
1289  val0 = shape_gradients[mid * n_columns + col];
1290  else
1291  val0 = shape_gradients[col * n_columns + mid];
1292  in1 = val0 * in[stride * mid];
1293  res0 += in1;
1294  res1 -= in1;
1295  }
1296  if (add == false)
1297  {
1298  out[stride * col] = res0;
1299  out[stride * (nn - 1 - col)] = res1;
1300  }
1301  else
1302  {
1303  out[stride * col] += res0;
1304  out[stride * (nn - 1 - col)] += res1;
1305  }
1306  }
1307  if (nn % 2 == 1)
1308  {
1309  Number2 val0;
1310  Number res0;
1311  if (contract_over_rows == true)
1312  val0 = shape_gradients[n_cols];
1313  else
1314  val0 = shape_gradients[n_cols * n_columns];
1315  res0 = val0 * (in[0] - in[stride * (mm - 1)]);
1316  for (int ind = 1; ind < mid; ++ind)
1317  {
1318  if (contract_over_rows == true)
1319  val0 = shape_gradients[ind * n_columns + n_cols];
1320  else
1321  val0 = shape_gradients[n_cols * n_columns + ind];
1322  Number in1 =
1323  val0 * (in[stride * ind] - in[stride * (mm - 1 - ind)]);
1324  res0 += in1;
1325  }
1326  if (add == false)
1327  out[stride * n_cols] = res0;
1328  else
1329  out[stride * n_cols] += res0;
1330  }
1331 
1332  ++in;
1333  ++out;
1334  }
1335  in += stride * (mm - 1);
1336  out += stride * (nn - 1);
1337  }
1338  }
1339 
1340 
1341 
1342  // evaluates the given shape data in 1d-3d using the tensor product
1343  // form assuming the symmetries of unit cell shape hessians for
1344  // finite elements in FEEvaluation
1345  template <int dim,
1346  int n_rows,
1347  int n_columns,
1348  typename Number,
1349  typename Number2>
1350  template <int direction, bool contract_over_rows, bool add>
1351  inline void
1352  EvaluatorTensorProduct<evaluate_symmetric,
1353  dim,
1354  n_rows,
1355  n_columns,
1356  Number,
1357  Number2>::hessians(const Number in[],
1358  Number out[]) const
1359  {
1360  Assert(shape_hessians != nullptr, ExcNotInitialized());
1361  AssertIndexRange(direction, dim);
1362  constexpr int mm = contract_over_rows ? n_rows : n_columns;
1363  constexpr int nn = contract_over_rows ? n_columns : n_rows;
1364  constexpr int n_cols = nn / 2;
1365  constexpr int mid = mm / 2;
1366 
1367  constexpr int stride = Utilities::pow(n_columns, direction);
1368  constexpr int n_blocks1 = stride;
1369  constexpr int n_blocks2 =
1370  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1371 
1372  for (int i2 = 0; i2 < n_blocks2; ++i2)
1373  {
1374  for (int i1 = 0; i1 < n_blocks1; ++i1)
1375  {
1376  for (int col = 0; col < n_cols; ++col)
1377  {
1378  Number2 val0, val1;
1379  Number in0, in1, res0, res1;
1380  if (contract_over_rows == true)
1381  {
1382  val0 = shape_hessians[col];
1383  val1 = shape_hessians[nn - 1 - col];
1384  }
1385  else
1386  {
1387  val0 = shape_hessians[col * n_columns];
1388  val1 = shape_hessians[(col + 1) * n_columns - 1];
1389  }
1390  if (mid > 0)
1391  {
1392  in0 = in[0];
1393  in1 = in[stride * (mm - 1)];
1394  res0 = val0 * in0;
1395  res1 = val1 * in0;
1396  res0 += val1 * in1;
1397  res1 += val0 * in1;
1398  for (int ind = 1; ind < mid; ++ind)
1399  {
1400  if (contract_over_rows == true)
1401  {
1402  val0 = shape_hessians[ind * n_columns + col];
1403  val1 =
1404  shape_hessians[ind * n_columns + nn - 1 - col];
1405  }
1406  else
1407  {
1408  val0 = shape_hessians[col * n_columns + ind];
1409  val1 =
1410  shape_hessians[(col + 1) * n_columns - 1 - ind];
1411  }
1412  in0 = in[stride * ind];
1413  in1 = in[stride * (mm - 1 - ind)];
1414  res0 += val0 * in0;
1415  res1 += val1 * in0;
1416  res0 += val1 * in1;
1417  res1 += val0 * in1;
1418  }
1419  }
1420  else
1421  res0 = res1 = Number();
1422  if (mm % 2 == 1)
1423  {
1424  if (contract_over_rows == true)
1425  val0 = shape_hessians[mid * n_columns + col];
1426  else
1427  val0 = shape_hessians[col * n_columns + mid];
1428  in1 = val0 * in[stride * mid];
1429  res0 += in1;
1430  res1 += in1;
1431  }
1432  if (add == false)
1433  {
1434  out[stride * col] = res0;
1435  out[stride * (nn - 1 - col)] = res1;
1436  }
1437  else
1438  {
1439  out[stride * col] += res0;
1440  out[stride * (nn - 1 - col)] += res1;
1441  }
1442  }
1443  if (nn % 2 == 1)
1444  {
1445  Number2 val0;
1446  Number res0;
1447  if (contract_over_rows == true)
1448  val0 = shape_hessians[n_cols];
1449  else
1450  val0 = shape_hessians[n_cols * n_columns];
1451  if (mid > 0)
1452  {
1453  res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1454  for (int ind = 1; ind < mid; ++ind)
1455  {
1456  if (contract_over_rows == true)
1457  val0 = shape_hessians[ind * n_columns + n_cols];
1458  else
1459  val0 = shape_hessians[n_cols * n_columns + ind];
1460  Number in1 = val0 * (in[stride * ind] +
1461  in[stride * (mm - 1 - ind)]);
1462  res0 += in1;
1463  }
1464  }
1465  else
1466  res0 = Number();
1467  if (mm % 2 == 1)
1468  {
1469  if (contract_over_rows == true)
1470  val0 = shape_hessians[mid * n_columns + n_cols];
1471  else
1472  val0 = shape_hessians[n_cols * n_columns + mid];
1473  res0 += val0 * in[stride * mid];
1474  }
1475  if (add == false)
1476  out[stride * n_cols] = res0;
1477  else
1478  out[stride * n_cols] += res0;
1479  }
1480 
1481  ++in;
1482  ++out;
1483  }
1484  in += stride * (mm - 1);
1485  out += stride * (nn - 1);
1486  }
1487  }
1488 
1489 
1490 
1522  template <int dim,
1523  int n_rows,
1524  int n_columns,
1525  typename Number,
1526  typename Number2>
1527  struct EvaluatorTensorProduct<evaluate_evenodd,
1528  dim,
1529  n_rows,
1530  n_columns,
1531  Number,
1532  Number2>
1533  {
1534  static constexpr unsigned int n_rows_of_product =
1535  Utilities::pow(n_rows, dim);
1536  static constexpr unsigned int n_columns_of_product =
1537  Utilities::pow(n_columns, dim);
1538 
1544  EvaluatorTensorProduct()
1545  : shape_values(nullptr)
1546  , shape_gradients(nullptr)
1547  , shape_hessians(nullptr)
1548  {}
1549 
1554  EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values)
1555  : shape_values(shape_values.begin())
1556  , shape_gradients(nullptr)
1557  , shape_hessians(nullptr)
1558  {
1559  AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
1560  }
1561 
1566  EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
1567  const AlignedVector<Number2> &shape_gradients,
1568  const AlignedVector<Number2> &shape_hessians,
1569  const unsigned int dummy1 = 0,
1570  const unsigned int dummy2 = 0)
1571  : shape_values(shape_values.begin())
1572  , shape_gradients(shape_gradients.begin())
1573  , shape_hessians(shape_hessians.begin())
1574  {
1575  // In this function, we allow for dummy pointers if some of values,
1576  // gradients or hessians should not be computed
1577  if (!shape_values.empty())
1578  AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
1579  if (!shape_gradients.empty())
1580  AssertDimension(shape_gradients.size(), n_rows * ((n_columns + 1) / 2));
1581  if (!shape_hessians.empty())
1582  AssertDimension(shape_hessians.size(), n_rows * ((n_columns + 1) / 2));
1583  (void)dummy1;
1584  (void)dummy2;
1585  }
1586 
1587  template <int direction, bool contract_over_rows, bool add>
1588  void
1589  values(const Number in[], Number out[]) const
1590  {
1591  Assert(shape_values != nullptr, ExcNotInitialized());
1592  apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
1593  }
1594 
1595  template <int direction, bool contract_over_rows, bool add>
1596  void
1597  gradients(const Number in[], Number out[]) const
1598  {
1599  Assert(shape_gradients != nullptr, ExcNotInitialized());
1600  apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
1601  }
1602 
1603  template <int direction, bool contract_over_rows, bool add>
1604  void
1605  hessians(const Number in[], Number out[]) const
1606  {
1607  Assert(shape_hessians != nullptr, ExcNotInitialized());
1608  apply<direction, contract_over_rows, add, 2>(shape_hessians, in, out);
1609  }
1610 
1611  template <int direction, bool contract_over_rows, bool add>
1612  void
1613  values_one_line(const Number in[], Number out[]) const
1614  {
1615  Assert(shape_values != nullptr, ExcNotInitialized());
1616  apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
1617  }
1618 
1619  template <int direction, bool contract_over_rows, bool add>
1620  void
1621  gradients_one_line(const Number in[], Number out[]) const
1622  {
1623  Assert(shape_gradients != nullptr, ExcNotInitialized());
1624  apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
1625  in,
1626  out);
1627  }
1628 
1629  template <int direction, bool contract_over_rows, bool add>
1630  void
1631  hessians_one_line(const Number in[], Number out[]) const
1632  {
1633  Assert(shape_hessians != nullptr, ExcNotInitialized());
1634  apply<direction, contract_over_rows, add, 2, true>(shape_hessians,
1635  in,
1636  out);
1637  }
1638 
1667  template <int direction,
1668  bool contract_over_rows,
1669  bool add,
1670  int type,
1671  bool one_line = false>
1672  static void
1673  apply(const Number2 *DEAL_II_RESTRICT shape_data,
1674  const Number * in,
1675  Number * out);
1676 
1677  const Number2 *shape_values;
1678  const Number2 *shape_gradients;
1679  const Number2 *shape_hessians;
1680  };
1681 
1682 
1683 
1684  template <int dim,
1685  int n_rows,
1686  int n_columns,
1687  typename Number,
1688  typename Number2>
1689  template <int direction,
1690  bool contract_over_rows,
1691  bool add,
1692  int type,
1693  bool one_line>
1694  inline void
1695  EvaluatorTensorProduct<evaluate_evenodd,
1696  dim,
1697  n_rows,
1698  n_columns,
1699  Number,
1700  Number2>::apply(const Number2 *DEAL_II_RESTRICT shapes,
1701  const Number * in,
1702  Number * out)
1703  {
1704  static_assert(type < 3, "Only three variants type=0,1,2 implemented");
1705  static_assert(one_line == false || direction == dim - 1,
1706  "Single-line evaluation only works for direction=dim-1.");
1707  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
1708  in != out,
1709  ExcMessage("In-place operation only supported for "
1710  "n_rows==n_columns or single-line interpolation"));
1711 
1712  // We cannot statically assert that direction is less than dim, so must do
1713  // an additional dynamic check
1714  AssertIndexRange(direction, dim);
1715 
1716  constexpr int nn = contract_over_rows ? n_columns : n_rows;
1717  constexpr int mm = contract_over_rows ? n_rows : n_columns;
1718  constexpr int n_cols = nn / 2;
1719  constexpr int mid = mm / 2;
1720 
1721  constexpr int stride = Utilities::pow(n_columns, direction);
1722  constexpr int n_blocks1 = one_line ? 1 : stride;
1723  constexpr int n_blocks2 =
1724  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1725 
1726  constexpr int offset = (n_columns + 1) / 2;
1727 
1728  // this code may look very inefficient at first sight due to the many
1729  // different cases with if's at the innermost loop part, but all of the
1730  // conditionals can be evaluated at compile time because they are
1731  // templates, so the compiler should optimize everything away
1732  for (int i2 = 0; i2 < n_blocks2; ++i2)
1733  {
1734  for (int i1 = 0; i1 < n_blocks1; ++i1)
1735  {
1736  Number xp[mid > 0 ? mid : 1], xm[mid > 0 ? mid : 1];
1737  for (int i = 0; i < mid; ++i)
1738  {
1739  if (contract_over_rows == true && type == 1)
1740  {
1741  xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
1742  xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
1743  }
1744  else
1745  {
1746  xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
1747  xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
1748  }
1749  }
1750  Number xmid = in[stride * mid];
1751  for (int col = 0; col < n_cols; ++col)
1752  {
1753  Number r0, r1;
1754  if (mid > 0)
1755  {
1756  if (contract_over_rows == true)
1757  {
1758  r0 = shapes[col] * xp[0];
1759  r1 = shapes[(n_rows - 1) * offset + col] * xm[0];
1760  }
1761  else
1762  {
1763  r0 = shapes[col * offset] * xp[0];
1764  r1 = shapes[(n_rows - 1 - col) * offset] * xm[0];
1765  }
1766  for (int ind = 1; ind < mid; ++ind)
1767  {
1768  if (contract_over_rows == true)
1769  {
1770  r0 += shapes[ind * offset + col] * xp[ind];
1771  r1 += shapes[(n_rows - 1 - ind) * offset + col] *
1772  xm[ind];
1773  }
1774  else
1775  {
1776  r0 += shapes[col * offset + ind] * xp[ind];
1777  r1 += shapes[(n_rows - 1 - col) * offset + ind] *
1778  xm[ind];
1779  }
1780  }
1781  }
1782  else
1783  r0 = r1 = Number();
1784  if (mm % 2 == 1 && contract_over_rows == true)
1785  {
1786  if (type == 1)
1787  r1 += shapes[mid * offset + col] * xmid;
1788  else
1789  r0 += shapes[mid * offset + col] * xmid;
1790  }
1791  else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0 || mm == 3))
1792  r0 += shapes[col * offset + mid] * xmid;
1793 
1794  if (add == false)
1795  {
1796  out[stride * col] = r0 + r1;
1797  if (type == 1 && contract_over_rows == false)
1798  out[stride * (nn - 1 - col)] = r1 - r0;
1799  else
1800  out[stride * (nn - 1 - col)] = r0 - r1;
1801  }
1802  else
1803  {
1804  out[stride * col] += r0 + r1;
1805  if (type == 1 && contract_over_rows == false)
1806  out[stride * (nn - 1 - col)] += r1 - r0;
1807  else
1808  out[stride * (nn - 1 - col)] += r0 - r1;
1809  }
1810  }
1811  if (type == 0 && contract_over_rows == true && nn % 2 == 1 &&
1812  mm % 2 == 1 && mm > 3)
1813  {
1814  if (add == false)
1815  out[stride * n_cols] = shapes[mid * offset + n_cols] * xmid;
1816  else
1817  out[stride * n_cols] += shapes[mid * offset + n_cols] * xmid;
1818  }
1819  else if (contract_over_rows == true && nn % 2 == 1)
1820  {
1821  Number r0;
1822  if (mid > 0)
1823  {
1824  r0 = shapes[n_cols] * xp[0];
1825  for (int ind = 1; ind < mid; ++ind)
1826  r0 += shapes[ind * offset + n_cols] * xp[ind];
1827  }
1828  else
1829  r0 = Number();
1830  if (type != 1 && mm % 2 == 1)
1831  r0 += shapes[mid * offset + n_cols] * xmid;
1832 
1833  if (add == false)
1834  out[stride * n_cols] = r0;
1835  else
1836  out[stride * n_cols] += r0;
1837  }
1838  else if (contract_over_rows == false && nn % 2 == 1)
1839  {
1840  Number r0;
1841  if (mid > 0)
1842  {
1843  if (type == 1)
1844  {
1845  r0 = shapes[n_cols * offset] * xm[0];
1846  for (int ind = 1; ind < mid; ++ind)
1847  r0 += shapes[n_cols * offset + ind] * xm[ind];
1848  }
1849  else
1850  {
1851  r0 = shapes[n_cols * offset] * xp[0];
1852  for (int ind = 1; ind < mid; ++ind)
1853  r0 += shapes[n_cols * offset + ind] * xp[ind];
1854  }
1855  }
1856  else
1857  r0 = Number();
1858 
1859  if ((type == 0 || type == 2) && mm % 2 == 1)
1860  r0 += shapes[n_cols * offset + mid] * xmid;
1861 
1862  if (add == false)
1863  out[stride * n_cols] = r0;
1864  else
1865  out[stride * n_cols] += r0;
1866  }
1867  if (one_line == false)
1868  {
1869  in += 1;
1870  out += 1;
1871  }
1872  }
1873  if (one_line == false)
1874  {
1875  in += stride * (mm - 1);
1876  out += stride * (nn - 1);
1877  }
1878  }
1879  }
1880 
1881 
1882 
1911  template <int dim,
1912  int n_rows,
1913  int n_columns,
1914  typename Number,
1915  typename Number2>
1916  struct EvaluatorTensorProduct<evaluate_symmetric_hierarchical,
1917  dim,
1918  n_rows,
1919  n_columns,
1920  Number,
1921  Number2>
1922  {
1923  static constexpr unsigned int n_rows_of_product =
1924  Utilities::pow(n_rows, dim);
1925  static constexpr unsigned int n_columns_of_product =
1926  Utilities::pow(n_columns, dim);
1927 
1933  EvaluatorTensorProduct()
1934  : shape_values(nullptr)
1935  , shape_gradients(nullptr)
1936  , shape_hessians(nullptr)
1937  {}
1938 
1943  EvaluatorTensorProduct(const AlignedVector<Number> &shape_values)
1944  : shape_values(shape_values.begin())
1945  , shape_gradients(nullptr)
1946  , shape_hessians(nullptr)
1947  {}
1948 
1953  EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
1954  const AlignedVector<Number2> &shape_gradients,
1955  const AlignedVector<Number2> &shape_hessians,
1956  const unsigned int dummy1 = 0,
1957  const unsigned int dummy2 = 0)
1958  : shape_values(shape_values.begin())
1959  , shape_gradients(shape_gradients.begin())
1960  , shape_hessians(shape_hessians.begin())
1961  {
1962  (void)dummy1;
1963  (void)dummy2;
1964  }
1965 
1966  template <int direction, bool contract_over_rows, bool add>
1967  void
1968  values(const Number in[], Number out[]) const
1969  {
1970  Assert(shape_values != nullptr, ExcNotInitialized());
1971  apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
1972  }
1973 
1974  template <int direction, bool contract_over_rows, bool add>
1975  void
1976  gradients(const Number in[], Number out[]) const
1977  {
1978  Assert(shape_gradients != nullptr, ExcNotInitialized());
1979  apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
1980  }
1981 
1982  template <int direction, bool contract_over_rows, bool add>
1983  void
1984  hessians(const Number in[], Number out[]) const
1985  {
1986  Assert(shape_hessians != nullptr, ExcNotInitialized());
1987  apply<direction, contract_over_rows, add, 0>(shape_hessians, in, out);
1988  }
1989 
1990  template <int direction, bool contract_over_rows, bool add>
1991  void
1992  values_one_line(const Number in[], Number out[]) const
1993  {
1994  Assert(shape_values != nullptr, ExcNotInitialized());
1995  apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
1996  }
1997 
1998  template <int direction, bool contract_over_rows, bool add>
1999  void
2000  gradients_one_line(const Number in[], Number out[]) const
2001  {
2002  Assert(shape_gradients != nullptr, ExcNotInitialized());
2003  apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
2004  in,
2005  out);
2006  }
2007 
2008  template <int direction, bool contract_over_rows, bool add>
2009  void
2010  hessians_one_line(const Number in[], Number out[]) const
2011  {
2012  Assert(shape_hessians != nullptr, ExcNotInitialized());
2013  apply<direction, contract_over_rows, add, 0, true>(shape_hessians,
2014  in,
2015  out);
2016  }
2017 
2045  template <int direction,
2046  bool contract_over_rows,
2047  bool add,
2048  int type,
2049  bool one_line = false>
2050  static void
2051  apply(const Number2 *DEAL_II_RESTRICT shape_data,
2052  const Number * in,
2053  Number * out);
2054 
2055  const Number2 *shape_values;
2056  const Number2 *shape_gradients;
2057  const Number2 *shape_hessians;
2058  };
2059 
2060 
2061 
2062  template <int dim,
2063  int n_rows,
2064  int n_columns,
2065  typename Number,
2066  typename Number2>
2067  template <int direction,
2068  bool contract_over_rows,
2069  bool add,
2070  int type,
2071  bool one_line>
2072  inline void
2073  EvaluatorTensorProduct<evaluate_symmetric_hierarchical,
2074  dim,
2075  n_rows,
2076  n_columns,
2077  Number,
2078  Number2>::apply(const Number2 *DEAL_II_RESTRICT shapes,
2079  const Number * in,
2080  Number * out)
2081  {
2082  static_assert(one_line == false || direction == dim - 1,
2083  "Single-line evaluation only works for direction=dim-1.");
2084  static_assert(
2085  type == 0 || type == 1,
2086  "Only types 0 and 1 implemented for evaluate_symmetric_hierarchical.");
2087  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
2088  in != out,
2089  ExcMessage("In-place operation only supported for "
2090  "n_rows==n_columns or single-line interpolation"));
2091 
2092  // We cannot statically assert that direction is less than dim, so must do
2093  // an additional dynamic check
2094  AssertIndexRange(direction, dim);
2095 
2096  constexpr int nn = contract_over_rows ? n_columns : n_rows;
2097  constexpr int mm = contract_over_rows ? n_rows : n_columns;
2098  constexpr int n_cols = nn / 2;
2099  constexpr int mid = mm / 2;
2100 
2101  constexpr int stride = Utilities::pow(n_columns, direction);
2102  constexpr int n_blocks1 = one_line ? 1 : stride;
2103  constexpr int n_blocks2 =
2104  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
2105 
2106  // this code may look very inefficient at first sight due to the many
2107  // different cases with if's at the innermost loop part, but all of the
2108  // conditionals can be evaluated at compile time because they are
2109  // templates, so the compiler should optimize everything away
2110  for (int i2 = 0; i2 < n_blocks2; ++i2)
2111  {
2112  for (int i1 = 0; i1 < n_blocks1; ++i1)
2113  {
2114  if (contract_over_rows)
2115  {
2116  Number x[mm];
2117  for (unsigned int i = 0; i < mm; ++i)
2118  x[i] = in[stride * i];
2119  for (unsigned int col = 0; col < n_cols; ++col)
2120  {
2121  Number r0, r1;
2122  if (mid > 0)
2123  {
2124  r0 = shapes[col] * x[0];
2125  r1 = shapes[col + n_columns] * x[1];
2126  for (unsigned int ind = 1; ind < mid; ++ind)
2127  {
2128  r0 +=
2129  shapes[col + 2 * ind * n_columns] * x[2 * ind];
2130  r1 += shapes[col + (2 * ind + 1) * n_columns] *
2131  x[2 * ind + 1];
2132  }
2133  }
2134  else
2135  r0 = r1 = Number();
2136  if (mm % 2 == 1)
2137  r0 += shapes[col + (mm - 1) * n_columns] * x[mm - 1];
2138  if (add == false)
2139  {
2140  out[stride * col] = r0 + r1;
2141  if (type == 1)
2142  out[stride * (nn - 1 - col)] = r1 - r0;
2143  else
2144  out[stride * (nn - 1 - col)] = r0 - r1;
2145  }
2146  else
2147  {
2148  out[stride * col] += r0 + r1;
2149  if (type == 1)
2150  out[stride * (nn - 1 - col)] += r1 - r0;
2151  else
2152  out[stride * (nn - 1 - col)] += r0 - r1;
2153  }
2154  }
2155  if (nn % 2 == 1)
2156  {
2157  Number r0;
2158  const unsigned int shift = type == 1 ? 1 : 0;
2159  if (mid > 0)
2160  {
2161  r0 = shapes[n_cols + shift * n_columns] * x[shift];
2162  for (unsigned int ind = 1; ind < mid; ++ind)
2163  r0 += shapes[n_cols + (2 * ind + shift) * n_columns] *
2164  x[2 * ind + shift];
2165  }
2166  else
2167  r0 = 0;
2168  if (type != 1 && mm % 2 == 1)
2169  r0 += shapes[n_cols + (mm - 1) * n_columns] * x[mm - 1];
2170  if (add == false)
2171  out[stride * n_cols] = r0;
2172  else
2173  out[stride * n_cols] += r0;
2174  }
2175  }
2176  else
2177  {
2178  Number xp[mid + 1], xm[mid > 0 ? mid : 1];
2179  for (int i = 0; i < mid; ++i)
2180  if (type == 0)
2181  {
2182  xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
2183  xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
2184  }
2185  else
2186  {
2187  xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
2188  xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
2189  }
2190  if (mm % 2 == 1)
2191  xp[mid] = in[stride * mid];
2192  for (unsigned int col = 0; col < n_cols; ++col)
2193  {
2194  Number r0, r1;
2195  if (mid > 0)
2196  {
2197  r0 = shapes[2 * col * n_columns] * xp[0];
2198  r1 = shapes[(2 * col + 1) * n_columns] * xm[0];
2199  for (unsigned int ind = 1; ind < mid; ++ind)
2200  {
2201  r0 += shapes[2 * col * n_columns + ind] * xp[ind];
2202  r1 +=
2203  shapes[(2 * col + 1) * n_columns + ind] * xm[ind];
2204  }
2205  }
2206  else
2207  r0 = r1 = Number();
2208  if (mm % 2 == 1)
2209  {
2210  if (type == 1)
2211  r1 +=
2212  shapes[(2 * col + 1) * n_columns + mid] * xp[mid];
2213  else
2214  r0 += shapes[2 * col * n_columns + mid] * xp[mid];
2215  }
2216  if (add == false)
2217  {
2218  out[stride * (2 * col)] = r0;
2219  out[stride * (2 * col + 1)] = r1;
2220  }
2221  else
2222  {
2223  out[stride * (2 * col)] += r0;
2224  out[stride * (2 * col + 1)] += r1;
2225  }
2226  }
2227  if (nn % 2 == 1)
2228  {
2229  Number r0;
2230  if (mid > 0)
2231  {
2232  r0 = shapes[(nn - 1) * n_columns] * xp[0];
2233  for (unsigned int ind = 1; ind < mid; ++ind)
2234  r0 += shapes[(nn - 1) * n_columns + ind] * xp[ind];
2235  }
2236  else
2237  r0 = Number();
2238  if (mm % 2 == 1 && type == 0)
2239  r0 += shapes[(nn - 1) * n_columns + mid] * xp[mid];
2240  if (add == false)
2241  out[stride * (nn - 1)] = r0;
2242  else
2243  out[stride * (nn - 1)] += r0;
2244  }
2245  }
2246  if (one_line == false)
2247  {
2248  in += 1;
2249  out += 1;
2250  }
2251  }
2252  if (one_line == false)
2253  {
2254  in += stride * (mm - 1);
2255  out += stride * (nn - 1);
2256  }
2257  }
2258  }
2259 
2260 } // end of namespace internal
2261 
2262 
2263 DEAL_II_NAMESPACE_CLOSE
2264 
2265 #endif
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::n_rows
const unsigned int n_rows
Definition: tensor_product_kernels.h:643
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::values_one_line
void values_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1613
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:476
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::gradients
void gradients(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1976
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::values
void values(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1589
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:1678
internal::EvaluatorVariant
EvaluatorVariant
Definition: tensor_product_kernels.h:36
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:1544
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:998
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::hessians
void hessians(const Number *in, Number *out) const
Definition: tensor_product_kernels.h:594
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:1679
internal::EvaluatorTensorProduct
Definition: tensor_product_kernels.h:96
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:1677
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::hessians_one_line
void hessians_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:617
utilities.h
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::values
void values(const Number *in, Number *out) const
Definition: tensor_product_kernels.h:580
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:291
AlignedVector::size
size_type size() const
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number > &shape_values)
Definition: tensor_product_kernels.h:1943
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:2056
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::hessians
void hessians(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1984
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition: tensor_product_kernels.h:148
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::values_one_line
void values_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:601
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:2057
internal::evaluate_symmetric_hierarchical
@ evaluate_symmetric_hierarchical
Definition: tensor_product_kernels.h:65
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:540
internal::evaluate_symmetric
@ evaluate_symmetric
Definition: tensor_product_kernels.h:48
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int n_rows, const unsigned int n_columns)
Definition: tensor_product_kernels.h:551
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:292
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:139
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::hessians_one_line
void hessians_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:214
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::values
void values(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1968
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::values_one_line
void values_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1992
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition: trilinos_sparse_matrix.cc:51
Utilities::fixed_power
T fixed_power(const T t)
Definition: utilities.h:1072
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values)
Definition: tensor_product_kernels.h:1554
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::gradients
void gradients(const Number *in, Number *out) const
Definition: tensor_product_kernels.h:587
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::gradients_one_line
void gradients_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:206
AlignedVector
Definition: aligned_vector.h:61
std_cxx17::apply
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std_cxx14::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:40
internal::evaluate_general
@ evaluate_general
Definition: tensor_product_kernels.h:42
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:290
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::hessians
void hessians(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1605
numbers
Definition: numbers.h:207
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::hessians_one_line
void hessians_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1631
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::gradients
void gradients(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1597
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:640
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:641
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::gradients_one_line
void gradients_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:609
DEAL_II_RESTRICT
#define DEAL_II_RESTRICT
Definition: config.h:100
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition: tensor_product_kernels.h:962
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::shape_gradients
const Number2 * shape_gradients
Definition: tensor_product_kernels.h:997
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
aligned_vector.h
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::hessians_one_line
void hessians_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:2010
GridTools::shift
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:817
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition: tensor_product_kernels.h:1953
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::values
void values(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:177
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::values_one_line
void values_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:198
config.h
internal::EvaluatorTensorProduct< evaluate_symmetric, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:996
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct()
Definition: tensor_product_kernels.h:1933
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::gradients
void gradients(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:184
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::gradients_one_line
void gradients_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:2000
internal
Definition: aligned_vector.h:369
AlignedVector::empty
bool empty() const
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
internal::evaluate_evenodd
@ evaluate_evenodd
Definition: tensor_product_kernels.h:54
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::gradients_one_line
void gradients_one_line(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:1621
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::shape_hessians
const Number2 * shape_hessians
Definition: tensor_product_kernels.h:642
internal::EvaluatorTensorProduct< evaluate_general, dim, 0, 0, Number, Number2 >::n_columns
const unsigned int n_columns
Definition: tensor_product_kernels.h:644
internal::EvaluatorTensorProduct< evaluate_symmetric_hierarchical, dim, n_rows, n_columns, Number, Number2 >::shape_values
const Number2 * shape_values
Definition: tensor_product_kernels.h:2055
internal::EvaluatorTensorProduct< evaluate_general, dim, n_rows, n_columns, Number, Number2 >::hessians
void hessians(const Number in[], Number out[]) const
Definition: tensor_product_kernels.h:191
internal::EvaluatorTensorProduct< evaluate_evenodd, dim, n_rows, n_columns, Number, Number2 >::EvaluatorTensorProduct
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
Definition: tensor_product_kernels.h:1566
Generated by   doxygen 1.8.17