Reference documentation for deal.II version GIT 3c37cfefb2 2023-03-24 04:00:03+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tensor_product_kernels.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_tensor_product_kernels_h
18 #define dealii_matrix_free_tensor_product_kernels_h
19 
20 #include <deal.II/base/config.h>
21 
23 #include <deal.II/base/ndarray.h>
25 #include <deal.II/base/utilities.h>
26 
27 
29 
30 
31 
32 namespace internal
33 {
39  {
72  };
73 
74 
75 
79  enum class EvaluatorQuantity
80  {
84  value,
88  gradient,
92  hessian
93  };
94 
95 
96 
117  template <EvaluatorVariant variant,
118  int dim,
119  int n_rows,
120  int n_columns,
121  typename Number,
122  typename Number2 = Number>
124  {};
125 
148  template <EvaluatorVariant variant,
149  int dim,
150  int n_rows,
151  int n_columns,
152  typename Number,
153  int normal_dir,
154  typename Number2 = Number>
156  {};
157 
158 
159 
177  template <int dim,
178  int n_rows,
179  int n_columns,
180  typename Number,
181  typename Number2>
183  dim,
184  n_rows,
185  n_columns,
186  Number,
187  Number2>
188  {
189  static constexpr unsigned int n_rows_of_product =
190  Utilities::pow(n_rows, dim);
191  static constexpr unsigned int n_columns_of_product =
192  Utilities::pow(n_columns, dim);
193 
199  : shape_values(nullptr)
200  , shape_gradients(nullptr)
201  , shape_hessians(nullptr)
202  {}
203 
208  const AlignedVector<Number2> &shape_gradients,
209  const AlignedVector<Number2> &shape_hessians,
210  const unsigned int dummy1 = 0,
211  const unsigned int dummy2 = 0)
212  : shape_values(shape_values.begin())
213  , shape_gradients(shape_gradients.begin())
214  , shape_hessians(shape_hessians.begin())
215  {
216  // We can enter this function either for the apply() path that has
217  // n_rows * n_columns entries or for the apply_face() path that only has
218  // n_rows * 3 entries in the array. Since we cannot decide about the use
219  // we must allow for both here.
220  Assert(shape_values.size() == 0 ||
221  shape_values.size() == n_rows * n_columns ||
222  shape_values.size() == 3 * n_rows,
223  ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
224  Assert(shape_gradients.size() == 0 ||
225  shape_gradients.size() == n_rows * n_columns,
226  ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
227  Assert(shape_hessians.size() == 0 ||
228  shape_hessians.size() == n_rows * n_columns,
229  ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
230  (void)dummy1;
231  (void)dummy2;
232  }
233 
234  template <int direction, bool contract_over_rows, bool add>
235  void
236  values(const Number in[], Number out[]) const
237  {
238  apply<direction, contract_over_rows, add>(shape_values, in, out);
239  }
240 
241  template <int direction, bool contract_over_rows, bool add>
242  void
243  gradients(const Number in[], Number out[]) const
244  {
245  apply<direction, contract_over_rows, add>(shape_gradients, in, out);
246  }
247 
248  template <int direction, bool contract_over_rows, bool add>
249  void
250  hessians(const Number in[], Number out[]) const
251  {
252  apply<direction, contract_over_rows, add>(shape_hessians, in, out);
253  }
254 
255  template <int direction, bool contract_over_rows, bool add>
256  void
257  values_one_line(const Number in[], Number out[]) const
258  {
259  Assert(shape_values != nullptr, ExcNotInitialized());
260  apply<direction, contract_over_rows, add, true>(shape_values, in, out);
261  }
262 
263  template <int direction, bool contract_over_rows, bool add>
264  void
265  gradients_one_line(const Number in[], Number out[]) const
266  {
267  Assert(shape_gradients != nullptr, ExcNotInitialized());
268  apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
269  }
270 
271  template <int direction, bool contract_over_rows, bool add>
272  void
273  hessians_one_line(const Number in[], Number out[]) const
274  {
275  Assert(shape_hessians != nullptr, ExcNotInitialized());
276  apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
277  }
278 
303  template <int direction,
304  bool contract_over_rows,
305  bool add,
306  bool one_line = false>
307  static void
308  apply(const Number2 *DEAL_II_RESTRICT shape_data,
309  const Number * in,
310  Number * out);
311 
341  template <int face_direction,
342  bool contract_onto_face,
343  bool add,
344  int max_derivative>
345  void
346  apply_face(const Number *DEAL_II_RESTRICT in,
347  Number *DEAL_II_RESTRICT out) const;
348 
349  private:
350  const Number2 *shape_values;
351  const Number2 *shape_gradients;
352  const Number2 *shape_hessians;
353  };
354 
355 
356 
357  template <int dim,
358  int n_rows,
359  int n_columns,
360  typename Number,
361  typename Number2>
362  template <int direction, bool contract_over_rows, bool add, bool one_line>
363  inline void
365  dim,
366  n_rows,
367  n_columns,
368  Number,
369  Number2>::apply(const Number2 *DEAL_II_RESTRICT
370  shape_data,
371  const Number *in,
372  Number * out)
373  {
374  static_assert(one_line == false || direction == dim - 1,
375  "Single-line evaluation only works for direction=dim-1.");
376  Assert(shape_data != nullptr,
377  ExcMessage(
378  "The given array shape_data must not be the null pointer!"));
379  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
380  in != out,
381  ExcMessage("In-place operation only supported for "
382  "n_rows==n_columns or single-line interpolation"));
383  AssertIndexRange(direction, dim);
384  constexpr int mm = contract_over_rows ? n_rows : n_columns,
385  nn = contract_over_rows ? n_columns : n_rows;
386 
387  constexpr int stride = Utilities::pow(n_columns, direction);
388  constexpr int n_blocks1 = one_line ? 1 : stride;
389  constexpr int n_blocks2 =
390  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
391 
392  for (int i2 = 0; i2 < n_blocks2; ++i2)
393  {
394  for (int i1 = 0; i1 < n_blocks1; ++i1)
395  {
396  Number x[mm];
397  for (int i = 0; i < mm; ++i)
398  x[i] = in[stride * i];
399  for (int col = 0; col < nn; ++col)
400  {
401  Number2 val0;
402  if (contract_over_rows == true)
403  val0 = shape_data[col];
404  else
405  val0 = shape_data[col * n_columns];
406  Number res0 = val0 * x[0];
407  for (int i = 1; i < mm; ++i)
408  {
409  if (contract_over_rows == true)
410  val0 = shape_data[i * n_columns + col];
411  else
412  val0 = shape_data[col * n_columns + i];
413  res0 += val0 * x[i];
414  }
415  if (add)
416  out[stride * col] += res0;
417  else
418  out[stride * col] = res0;
419  }
420 
421  if (one_line == false)
422  {
423  ++in;
424  ++out;
425  }
426  }
427  if (one_line == false)
428  {
429  in += stride * (mm - 1);
430  out += stride * (nn - 1);
431  }
432  }
433  }
434 
435 
436 
437  template <int dim,
438  int n_rows,
439  int n_columns,
440  typename Number,
441  typename Number2>
442  template <int face_direction,
443  bool contract_onto_face,
444  bool add,
445  int max_derivative>
446  inline void
448  dim,
449  n_rows,
450  n_columns,
451  Number,
452  Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
453  Number *DEAL_II_RESTRICT
454  out) const
455  {
456  Assert(dim > 0, ExcMessage("Only dim=1,2,3 supported"));
457  static_assert(max_derivative >= 0 && max_derivative < 3,
458  "Only derivative orders 0-2 implemented");
459  Assert(shape_values != nullptr,
460  ExcMessage(
461  "The given array shape_values must not be the null pointer."));
462 
463  constexpr int n_blocks1 = (dim > 1 ? n_rows : 1);
464  constexpr int n_blocks2 = (dim > 2 ? n_rows : 1);
465 
466  AssertIndexRange(face_direction, dim);
467  constexpr int in_stride = Utilities::pow(n_rows, face_direction);
468  constexpr int out_stride = Utilities::pow(n_rows, dim - 1);
469  const Number *DEAL_II_RESTRICT shape_values = this->shape_values;
470 
471  for (int i2 = 0; i2 < n_blocks2; ++i2)
472  {
473  for (int i1 = 0; i1 < n_blocks1; ++i1)
474  {
475  if (contract_onto_face == true)
476  {
477  Number res0 = shape_values[0] * in[0];
478  Number res1, res2;
479  if (max_derivative > 0)
480  res1 = shape_values[n_rows] * in[0];
481  if (max_derivative > 1)
482  res2 = shape_values[2 * n_rows] * in[0];
483  for (int ind = 1; ind < n_rows; ++ind)
484  {
485  res0 += shape_values[ind] * in[in_stride * ind];
486  if (max_derivative > 0)
487  res1 += shape_values[ind + n_rows] * in[in_stride * ind];
488  if (max_derivative > 1)
489  res2 +=
490  shape_values[ind + 2 * n_rows] * in[in_stride * ind];
491  }
492  if (add)
493  {
494  out[0] += res0;
495  if (max_derivative > 0)
496  out[out_stride] += res1;
497  if (max_derivative > 1)
498  out[2 * out_stride] += res2;
499  }
500  else
501  {
502  out[0] = res0;
503  if (max_derivative > 0)
504  out[out_stride] = res1;
505  if (max_derivative > 1)
506  out[2 * out_stride] = res2;
507  }
508  }
509  else
510  {
511  for (int col = 0; col < n_rows; ++col)
512  {
513  if (add)
514  out[col * in_stride] += shape_values[col] * in[0];
515  else
516  out[col * in_stride] = shape_values[col] * in[0];
517  if (max_derivative > 0)
518  out[col * in_stride] +=
519  shape_values[col + n_rows] * in[out_stride];
520  if (max_derivative > 1)
521  out[col * in_stride] +=
522  shape_values[col + 2 * n_rows] * in[2 * out_stride];
523  }
524  }
525 
526  // increment: in regular case, just go to the next point in
527  // x-direction. If we are at the end of one chunk in x-dir, need
528  // to jump over to the next layer in z-direction
529  switch (face_direction)
530  {
531  case 0:
532  in += contract_onto_face ? n_rows : 1;
533  out += contract_onto_face ? 1 : n_rows;
534  break;
535  case 1:
536  ++in;
537  ++out;
538  // faces 2 and 3 in 3d use local coordinate system zx, which
539  // is the other way around compared to the tensor
540  // product. Need to take that into account.
541  if (dim == 3)
542  {
543  if (contract_onto_face)
544  out += n_rows - 1;
545  else
546  in += n_rows - 1;
547  }
548  break;
549  case 2:
550  ++in;
551  ++out;
552  break;
553  default:
554  Assert(false, ExcNotImplemented());
555  }
556  }
557 
558  // adjust for local coordinate system zx
559  if (face_direction == 1 && dim == 3)
560  {
561  if (contract_onto_face)
562  {
563  in += n_rows * (n_rows - 1);
564  out -= n_rows * n_rows - 1;
565  }
566  else
567  {
568  out += n_rows * (n_rows - 1);
569  in -= n_rows * n_rows - 1;
570  }
571  }
572  }
573  }
574 
575 
576 
590  template <int dim, typename Number, typename Number2>
591  struct EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>
592  {
593  static constexpr unsigned int n_rows_of_product =
595  static constexpr unsigned int n_columns_of_product =
597 
603  : shape_values(nullptr)
604  , shape_gradients(nullptr)
605  , shape_hessians(nullptr)
606  , n_rows(numbers::invalid_unsigned_int)
607  , n_columns(numbers::invalid_unsigned_int)
608  {}
609 
614  const AlignedVector<Number2> &shape_gradients,
615  const AlignedVector<Number2> &shape_hessians,
616  const unsigned int n_rows,
617  const unsigned int n_columns)
618  : shape_values(shape_values.begin())
619  , shape_gradients(shape_gradients.begin())
620  , shape_hessians(shape_hessians.begin())
621  , n_rows(n_rows)
622  , n_columns(n_columns)
623  {
624  // We can enter this function either for the apply() path that has
625  // n_rows * n_columns entries or for the apply_face() path that only has
626  // n_rows * 3 entries in the array. Since we cannot decide about the use
627  // we must allow for both here.
628  Assert(shape_values.size() == 0 ||
629  shape_values.size() == n_rows * n_columns ||
630  shape_values.size() == n_rows * 3,
631  ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
632  Assert(shape_gradients.size() == 0 ||
633  shape_gradients.size() == n_rows * n_columns,
634  ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
635  Assert(shape_hessians.size() == 0 ||
636  shape_hessians.size() == n_rows * n_columns,
637  ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
638  }
639 
643  EvaluatorTensorProduct(const Number2 * shape_values,
644  const Number2 * shape_gradients,
645  const Number2 * shape_hessians,
646  const unsigned int n_rows,
647  const unsigned int n_columns)
648  : shape_values(shape_values)
649  , shape_gradients(shape_gradients)
650  , shape_hessians(shape_hessians)
651  , n_rows(n_rows)
652  , n_columns(n_columns)
653  {}
654 
655  template <int direction, bool contract_over_rows, bool add>
656  void
657  values(const Number *in, Number *out) const
658  {
659  apply<direction, contract_over_rows, add>(shape_values, in, out);
660  }
661 
662  template <int direction, bool contract_over_rows, bool add>
663  void
664  gradients(const Number *in, Number *out) const
665  {
666  apply<direction, contract_over_rows, add>(shape_gradients, in, out);
667  }
668 
669  template <int direction, bool contract_over_rows, bool add>
670  void
671  hessians(const Number *in, Number *out) const
672  {
673  apply<direction, contract_over_rows, add>(shape_hessians, in, out);
674  }
675 
676  template <int direction, bool contract_over_rows, bool add>
677  void
678  values_one_line(const Number in[], Number out[]) const
679  {
680  Assert(shape_values != nullptr, ExcNotInitialized());
681  apply<direction, contract_over_rows, add, true>(shape_values, in, out);
682  }
683 
684  template <int direction, bool contract_over_rows, bool add>
685  void
686  gradients_one_line(const Number in[], Number out[]) const
687  {
688  Assert(shape_gradients != nullptr, ExcNotInitialized());
689  apply<direction, contract_over_rows, add, true>(shape_gradients, in, out);
690  }
691 
692  template <int direction, bool contract_over_rows, bool add>
693  void
694  hessians_one_line(const Number in[], Number out[]) const
695  {
696  Assert(shape_hessians != nullptr, ExcNotInitialized());
697  apply<direction, contract_over_rows, add, true>(shape_hessians, in, out);
698  }
699 
700  template <int direction,
701  bool contract_over_rows,
702  bool add,
703  bool one_line = false>
704  void
705  apply(const Number2 *DEAL_II_RESTRICT shape_data,
706  const Number * in,
707  Number * out) const;
708 
709  template <int face_direction,
710  bool contract_onto_face,
711  bool add,
712  int max_derivative>
713  void
714  apply_face(const Number *DEAL_II_RESTRICT in,
715  Number *DEAL_II_RESTRICT out) const;
716 
717  const Number2 * shape_values;
718  const Number2 * shape_gradients;
719  const Number2 * shape_hessians;
720  const unsigned int n_rows;
721  const unsigned int n_columns;
722  };
723 
724 
725 
726  template <int dim, typename Number, typename Number2>
727  template <int direction, bool contract_over_rows, bool add, bool one_line>
728  inline void
730  const Number2 *DEAL_II_RESTRICT shape_data,
731  const Number * in,
732  Number * out) const
733  {
734  static_assert(one_line == false || direction == dim - 1,
735  "Single-line evaluation only works for direction=dim-1.");
736  Assert(shape_data != nullptr,
737  ExcMessage(
738  "The given array shape_data must not be the null pointer!"));
739  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
740  in != out,
741  ExcMessage("In-place operation only supported for "
742  "n_rows==n_columns or single-line interpolation"));
743  AssertIndexRange(direction, dim);
744  const int mm = contract_over_rows ? n_rows : n_columns,
745  nn = contract_over_rows ? n_columns : n_rows;
746 
747  const int stride =
748  direction == 0 ? 1 : Utilities::fixed_power<direction>(n_columns);
749  const int n_blocks1 = one_line ? 1 : stride;
750  const int n_blocks2 = direction >= dim - 1 ?
751  1 :
752  Utilities::fixed_power<dim - direction - 1>(n_rows);
753  Assert(n_rows <= 128, ExcNotImplemented());
754 
755  // specialization for n_rows = 2 that manually unrolls the innermost loop
756  // to make the operation perform better (not completely as good as the
757  // templated one, but much better than the generic version down below,
758  // because the loop over col can be more effectively unrolled by the
759  // compiler)
760  if (contract_over_rows && n_rows == 2)
761  {
762  const Number2 *shape_data_1 = shape_data + n_columns;
763  for (int i2 = 0; i2 < n_blocks2; ++i2)
764  {
765  for (int i1 = 0; i1 < n_blocks1; ++i1)
766  {
767  const Number x0 = in[0], x1 = in[stride];
768  for (int col = 0; col < nn; ++col)
769  {
770  const Number result =
771  shape_data[col] * x0 + shape_data_1[col] * x1;
772  if (add)
773  out[stride * col] += result;
774  else
775  out[stride * col] = result;
776  }
777 
778  if (one_line == false)
779  {
780  ++in;
781  ++out;
782  }
783  }
784  if (one_line == false)
785  {
786  in += stride * (mm - 1);
787  out += stride * (nn - 1);
788  }
789  }
790  }
791  // specialization for n = 3
792  else if (contract_over_rows && n_rows == 3)
793  {
794  const Number2 *shape_data_1 = shape_data + n_columns;
795  const Number2 *shape_data_2 = shape_data + 2 * n_columns;
796  for (int i2 = 0; i2 < n_blocks2; ++i2)
797  {
798  for (int i1 = 0; i1 < n_blocks1; ++i1)
799  {
800  const Number x0 = in[0], x1 = in[stride], x2 = in[2 * stride];
801  for (int col = 0; col < nn; ++col)
802  {
803  const Number result = shape_data[col] * x0 +
804  shape_data_1[col] * x1 +
805  shape_data_2[col] * x2;
806  if (add)
807  out[stride * col] += result;
808  else
809  out[stride * col] = result;
810  }
811 
812  if (one_line == false)
813  {
814  ++in;
815  ++out;
816  }
817  }
818  if (one_line == false)
819  {
820  in += stride * (mm - 1);
821  out += stride * (nn - 1);
822  }
823  }
824  }
825  // general loop for all other cases
826  else
827  for (int i2 = 0; i2 < n_blocks2; ++i2)
828  {
829  for (int i1 = 0; i1 < n_blocks1; ++i1)
830  {
831  Number x[129];
832  for (int i = 0; i < mm; ++i)
833  x[i] = in[stride * i];
834  for (int col = 0; col < nn; ++col)
835  {
836  Number2 val0;
837  if (contract_over_rows == true)
838  val0 = shape_data[col];
839  else
840  val0 = shape_data[col * n_columns];
841  Number res0 = val0 * x[0];
842  for (int i = 1; i < mm; ++i)
843  {
844  if (contract_over_rows == true)
845  val0 = shape_data[i * n_columns + col];
846  else
847  val0 = shape_data[col * n_columns + i];
848  res0 += val0 * x[i];
849  }
850  if (add)
851  out[stride * col] += res0;
852  else
853  out[stride * col] = res0;
854  }
855 
856  if (one_line == false)
857  {
858  ++in;
859  ++out;
860  }
861  }
862  if (one_line == false)
863  {
864  in += stride * (mm - 1);
865  out += stride * (nn - 1);
866  }
867  }
868  }
869 
870 
871 
872  template <int dim, typename Number, typename Number2>
873  template <int face_direction,
874  bool contract_onto_face,
875  bool add,
876  int max_derivative>
877  inline void
879  apply_face(const Number *DEAL_II_RESTRICT in,
880  Number *DEAL_II_RESTRICT out) const
881  {
882  Assert(shape_values != nullptr,
883  ExcMessage(
884  "The given array shape_data must not be the null pointer!"));
885  static_assert(dim > 0 && dim < 4, "Only dim=1,2,3 supported");
886  const int n_blocks1 = dim > 1 ? n_rows : 1;
887  const int n_blocks2 = dim > 2 ? n_rows : 1;
888 
889  AssertIndexRange(face_direction, dim);
890  const int in_stride =
891  face_direction > 0 ? Utilities::fixed_power<face_direction>(n_rows) : 1;
892  const int out_stride =
893  dim > 1 ? Utilities::fixed_power<dim - 1>(n_rows) : 1;
894 
895  for (int i2 = 0; i2 < n_blocks2; ++i2)
896  {
897  for (int i1 = 0; i1 < n_blocks1; ++i1)
898  {
899  if (contract_onto_face == true)
900  {
901  Number res0 = shape_values[0] * in[0];
902  Number res1, res2;
903  if (max_derivative > 0)
904  res1 = shape_values[n_rows] * in[0];
905  if (max_derivative > 1)
906  res2 = shape_values[2 * n_rows] * in[0];
907  for (unsigned int ind = 1; ind < n_rows; ++ind)
908  {
909  res0 += shape_values[ind] * in[in_stride * ind];
910  if (max_derivative > 0)
911  res1 += shape_values[ind + n_rows] * in[in_stride * ind];
912  if (max_derivative > 1)
913  res2 +=
914  shape_values[ind + 2 * n_rows] * in[in_stride * ind];
915  }
916  if (add)
917  {
918  out[0] += res0;
919  if (max_derivative > 0)
920  out[out_stride] += res1;
921  if (max_derivative > 1)
922  out[2 * out_stride] += res2;
923  }
924  else
925  {
926  out[0] = res0;
927  if (max_derivative > 0)
928  out[out_stride] = res1;
929  if (max_derivative > 1)
930  out[2 * out_stride] = res2;
931  }
932  }
933  else
934  {
935  for (unsigned int col = 0; col < n_rows; ++col)
936  {
937  if (add)
938  out[col * in_stride] += shape_values[col] * in[0];
939  else
940  out[col * in_stride] = shape_values[col] * in[0];
941  if (max_derivative > 0)
942  out[col * in_stride] +=
943  shape_values[col + n_rows] * in[out_stride];
944  if (max_derivative > 1)
945  out[col * in_stride] +=
946  shape_values[col + 2 * n_rows] * in[2 * out_stride];
947  }
948  }
949 
950  // increment: in regular case, just go to the next point in
951  // x-direction. If we are at the end of one chunk in x-dir, need
952  // to jump over to the next layer in z-direction
953  switch (face_direction)
954  {
955  case 0:
956  in += contract_onto_face ? n_rows : 1;
957  out += contract_onto_face ? 1 : n_rows;
958  break;
959  case 1:
960  ++in;
961  ++out;
962  // faces 2 and 3 in 3d use local coordinate system zx, which
963  // is the other way around compared to the tensor
964  // product. Need to take that into account.
965  if (dim == 3)
966  {
967  if (contract_onto_face)
968  out += n_rows - 1;
969  else
970  in += n_rows - 1;
971  }
972  break;
973  case 2:
974  ++in;
975  ++out;
976  break;
977  default:
978  Assert(false, ExcNotImplemented());
979  }
980  }
981  if (face_direction == 1 && dim == 3)
982  {
983  // adjust for local coordinate system zx
984  if (contract_onto_face)
985  {
986  in += n_rows * (n_rows - 1);
987  out -= n_rows * n_rows - 1;
988  }
989  else
990  {
991  out += n_rows * (n_rows - 1);
992  in -= n_rows * n_rows - 1;
993  }
994  }
995  }
996  }
997 
998 
999 
1020  template <int dim,
1021  int n_rows,
1022  int n_columns,
1023  typename Number,
1024  typename Number2>
1026  dim,
1027  n_rows,
1028  n_columns,
1029  Number,
1030  Number2>
1031  {
1032  static constexpr unsigned int n_rows_of_product =
1033  Utilities::pow(n_rows, dim);
1034  static constexpr unsigned int n_columns_of_product =
1035  Utilities::pow(n_columns, dim);
1036 
1041  const AlignedVector<Number2> &shape_gradients,
1042  const AlignedVector<Number2> &shape_hessians,
1043  const unsigned int dummy1 = 0,
1044  const unsigned int dummy2 = 0)
1045  : shape_values(shape_values.begin())
1046  , shape_gradients(shape_gradients.begin())
1047  , shape_hessians(shape_hessians.begin())
1048  {
1049  Assert(shape_values.size() == 0 ||
1050  shape_values.size() == n_rows * n_columns,
1051  ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
1052  Assert(shape_gradients.size() == 0 ||
1053  shape_gradients.size() == n_rows * n_columns,
1054  ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
1055  Assert(shape_hessians.size() == 0 ||
1056  shape_hessians.size() == n_rows * n_columns,
1057  ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
1058  (void)dummy1;
1059  (void)dummy2;
1060  }
1061 
1062  template <int direction, bool contract_over_rows, bool add>
1063  void
1064  values(const Number in[], Number out[]) const;
1065 
1066  template <int direction, bool contract_over_rows, bool add>
1067  void
1068  gradients(const Number in[], Number out[]) const;
1069 
1070  template <int direction, bool contract_over_rows, bool add>
1071  void
1072  hessians(const Number in[], Number out[]) const;
1073 
1074  private:
1075  const Number2 *shape_values;
1076  const Number2 *shape_gradients;
1077  const Number2 *shape_hessians;
1078  };
1079 
1080 
1081 
1082  // In this case, the 1d shape values read (sorted lexicographically, rows
1083  // run over 1d dofs, columns over quadrature points):
1084  // Q2 --> [ 0.687 0 -0.087 ]
1085  // [ 0.4 1 0.4 ]
1086  // [-0.087 0 0.687 ]
1087  // Q3 --> [ 0.66 0.003 0.002 0.049 ]
1088  // [ 0.521 1.005 -0.01 -0.230 ]
1089  // [-0.230 -0.01 1.005 0.521 ]
1090  // [ 0.049 0.002 0.003 0.66 ]
1091  // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ]
1092  // [ 0.608 1.059 0 0.039 0.176 ]
1093  // [-0.409 -0.113 1 -0.113 -0.409 ]
1094  // [ 0.176 0.039 0 1.059 0.608 ]
1095  // [-0.032 -0.007 0 0.022 0.658 ]
1096  //
1097  // In these matrices, we want to use avoid computations involving zeros and
1098  // ones and in addition use the symmetry in entries to reduce the number of
1099  // read operations.
1100  template <int dim,
1101  int n_rows,
1102  int n_columns,
1103  typename Number,
1104  typename Number2>
1105  template <int direction, bool contract_over_rows, bool add>
1106  inline void
1108  dim,
1109  n_rows,
1110  n_columns,
1111  Number,
1112  Number2>::values(const Number in[], Number out[]) const
1113  {
1114  Assert(shape_values != nullptr, ExcNotInitialized());
1115  AssertIndexRange(direction, dim);
1116  constexpr int mm = contract_over_rows ? n_rows : n_columns,
1117  nn = contract_over_rows ? n_columns : n_rows;
1118  constexpr int n_cols = nn / 2;
1119  constexpr int mid = mm / 2;
1120 
1121  constexpr int stride = Utilities::pow(n_columns, direction);
1122  constexpr int n_blocks1 = stride;
1123  constexpr int n_blocks2 =
1124  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1125 
1126  for (int i2 = 0; i2 < n_blocks2; ++i2)
1127  {
1128  for (int i1 = 0; i1 < n_blocks1; ++i1)
1129  {
1130  for (int col = 0; col < n_cols; ++col)
1131  {
1132  Number2 val0, val1;
1133  Number in0, in1, res0, res1;
1134  if (contract_over_rows == true)
1135  {
1136  val0 = shape_values[col];
1137  val1 = shape_values[nn - 1 - col];
1138  }
1139  else
1140  {
1141  val0 = shape_values[col * n_columns];
1142  val1 = shape_values[(col + 1) * n_columns - 1];
1143  }
1144  if (mid > 0)
1145  {
1146  in0 = in[0];
1147  in1 = in[stride * (mm - 1)];
1148  res0 = val0 * in0;
1149  res1 = val1 * in0;
1150  res0 += val1 * in1;
1151  res1 += val0 * in1;
1152  for (int ind = 1; ind < mid; ++ind)
1153  {
1154  if (contract_over_rows == true)
1155  {
1156  val0 = shape_values[ind * n_columns + col];
1157  val1 = shape_values[ind * n_columns + nn - 1 - col];
1158  }
1159  else
1160  {
1161  val0 = shape_values[col * n_columns + ind];
1162  val1 =
1163  shape_values[(col + 1) * n_columns - 1 - ind];
1164  }
1165  in0 = in[stride * ind];
1166  in1 = in[stride * (mm - 1 - ind)];
1167  res0 += val0 * in0;
1168  res1 += val1 * in0;
1169  res0 += val1 * in1;
1170  res1 += val0 * in1;
1171  }
1172  }
1173  else
1174  res0 = res1 = Number();
1175  if (contract_over_rows == true)
1176  {
1177  if (mm % 2 == 1)
1178  {
1179  val0 = shape_values[mid * n_columns + col];
1180  in1 = val0 * in[stride * mid];
1181  res0 += in1;
1182  res1 += in1;
1183  }
1184  }
1185  else
1186  {
1187  if (mm % 2 == 1 && nn % 2 == 0)
1188  {
1189  val0 = shape_values[col * n_columns + mid];
1190  in1 = val0 * in[stride * mid];
1191  res0 += in1;
1192  res1 += in1;
1193  }
1194  }
1195  if (add)
1196  {
1197  out[stride * col] += res0;
1198  out[stride * (nn - 1 - col)] += res1;
1199  }
1200  else
1201  {
1202  out[stride * col] = res0;
1203  out[stride * (nn - 1 - col)] = res1;
1204  }
1205  }
1206  if (contract_over_rows == true && nn % 2 == 1 && mm % 2 == 1)
1207  {
1208  if (add)
1209  out[stride * n_cols] += in[stride * mid];
1210  else
1211  out[stride * n_cols] = in[stride * mid];
1212  }
1213  else if (contract_over_rows == true && nn % 2 == 1)
1214  {
1215  Number res0;
1216  Number2 val0 = shape_values[n_cols];
1217  if (mid > 0)
1218  {
1219  res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1220  for (int ind = 1; ind < mid; ++ind)
1221  {
1222  val0 = shape_values[ind * n_columns + n_cols];
1223  res0 += val0 * (in[stride * ind] +
1224  in[stride * (mm - 1 - ind)]);
1225  }
1226  }
1227  else
1228  res0 = Number();
1229  if (add)
1230  out[stride * n_cols] += res0;
1231  else
1232  out[stride * n_cols] = res0;
1233  }
1234  else if (contract_over_rows == false && nn % 2 == 1)
1235  {
1236  Number res0;
1237  if (mid > 0)
1238  {
1239  Number2 val0 = shape_values[n_cols * n_columns];
1240  res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1241  for (int ind = 1; ind < mid; ++ind)
1242  {
1243  val0 = shape_values[n_cols * n_columns + ind];
1244  Number in1 = val0 * (in[stride * ind] +
1245  in[stride * (mm - 1 - ind)]);
1246  res0 += in1;
1247  }
1248  if (mm % 2)
1249  res0 += in[stride * mid];
1250  }
1251  else
1252  res0 = in[0];
1253  if (add)
1254  out[stride * n_cols] += res0;
1255  else
1256  out[stride * n_cols] = res0;
1257  }
1258 
1259  ++in;
1260  ++out;
1261  }
1262  in += stride * (mm - 1);
1263  out += stride * (nn - 1);
1264  }
1265  }
1266 
1267 
1268 
1269  // For the specialized loop used for the gradient computation in
1270  // here, the 1d shape values read (sorted lexicographically, rows
1271  // run over 1d dofs, columns over quadrature points):
1272  // Q2 --> [-2.549 -1 0.549 ]
1273  // [ 3.098 0 -3.098 ]
1274  // [-0.549 1 2.549 ]
1275  // Q3 --> [-4.315 -1.03 0.5 -0.44 ]
1276  // [ 6.07 -1.44 -2.97 2.196 ]
1277  // [-2.196 2.97 1.44 -6.07 ]
1278  // [ 0.44 -0.5 1.03 4.315 ]
1279  // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ]
1280  // [10.111 -2.76 -2.667 2.066 -2.306 ]
1281  // [-5.688 5.773 0 -5.773 5.688 ]
1282  // [ 2.306 -2.066 2.667 2.76 -10.111 ]
1283  // [-0.413 0.353 -0.333 -0.353 0.413 ]
1284  //
1285  // In these matrices, we want to use avoid computations involving
1286  // zeros and ones and in addition use the symmetry in entries to
1287  // reduce the number of read operations.
1288  template <int dim,
1289  int n_rows,
1290  int n_columns,
1291  typename Number,
1292  typename Number2>
1293  template <int direction, bool contract_over_rows, bool add>
1294  inline void
1296  dim,
1297  n_rows,
1298  n_columns,
1299  Number,
1300  Number2>::gradients(const Number in[],
1301  Number out[]) const
1302  {
1303  Assert(shape_gradients != nullptr, ExcNotInitialized());
1304  AssertIndexRange(direction, dim);
1305  constexpr int mm = contract_over_rows ? n_rows : n_columns,
1306  nn = contract_over_rows ? n_columns : n_rows;
1307  constexpr int n_cols = nn / 2;
1308  constexpr int mid = mm / 2;
1309 
1310  constexpr int stride = Utilities::pow(n_columns, direction);
1311  constexpr int n_blocks1 = stride;
1312  constexpr int n_blocks2 =
1313  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1314 
1315  for (int i2 = 0; i2 < n_blocks2; ++i2)
1316  {
1317  for (int i1 = 0; i1 < n_blocks1; ++i1)
1318  {
1319  for (int col = 0; col < n_cols; ++col)
1320  {
1321  Number2 val0, val1;
1322  Number in0, in1, res0, res1;
1323  if (contract_over_rows == true)
1324  {
1325  val0 = shape_gradients[col];
1326  val1 = shape_gradients[nn - 1 - col];
1327  }
1328  else
1329  {
1330  val0 = shape_gradients[col * n_columns];
1331  val1 = shape_gradients[(nn - col - 1) * n_columns];
1332  }
1333  if (mid > 0)
1334  {
1335  in0 = in[0];
1336  in1 = in[stride * (mm - 1)];
1337  res0 = val0 * in0;
1338  res1 = val1 * in0;
1339  res0 -= val1 * in1;
1340  res1 -= val0 * in1;
1341  for (int ind = 1; ind < mid; ++ind)
1342  {
1343  if (contract_over_rows == true)
1344  {
1345  val0 = shape_gradients[ind * n_columns + col];
1346  val1 =
1347  shape_gradients[ind * n_columns + nn - 1 - col];
1348  }
1349  else
1350  {
1351  val0 = shape_gradients[col * n_columns + ind];
1352  val1 =
1353  shape_gradients[(nn - col - 1) * n_columns + ind];
1354  }
1355  in0 = in[stride * ind];
1356  in1 = in[stride * (mm - 1 - ind)];
1357  res0 += val0 * in0;
1358  res1 += val1 * in0;
1359  res0 -= val1 * in1;
1360  res1 -= val0 * in1;
1361  }
1362  }
1363  else
1364  res0 = res1 = Number();
1365  if (mm % 2 == 1)
1366  {
1367  if (contract_over_rows == true)
1368  val0 = shape_gradients[mid * n_columns + col];
1369  else
1370  val0 = shape_gradients[col * n_columns + mid];
1371  in1 = val0 * in[stride * mid];
1372  res0 += in1;
1373  res1 -= in1;
1374  }
1375  if (add)
1376  {
1377  out[stride * col] += res0;
1378  out[stride * (nn - 1 - col)] += res1;
1379  }
1380  else
1381  {
1382  out[stride * col] = res0;
1383  out[stride * (nn - 1 - col)] = res1;
1384  }
1385  }
1386  if (nn % 2 == 1)
1387  {
1388  Number2 val0;
1389  Number res0;
1390  if (contract_over_rows == true)
1391  val0 = shape_gradients[n_cols];
1392  else
1393  val0 = shape_gradients[n_cols * n_columns];
1394  res0 = val0 * (in[0] - in[stride * (mm - 1)]);
1395  for (int ind = 1; ind < mid; ++ind)
1396  {
1397  if (contract_over_rows == true)
1398  val0 = shape_gradients[ind * n_columns + n_cols];
1399  else
1400  val0 = shape_gradients[n_cols * n_columns + ind];
1401  Number in1 =
1402  val0 * (in[stride * ind] - in[stride * (mm - 1 - ind)]);
1403  res0 += in1;
1404  }
1405  if (add)
1406  out[stride * n_cols] += res0;
1407  else
1408  out[stride * n_cols] = res0;
1409  }
1410 
1411  ++in;
1412  ++out;
1413  }
1414  in += stride * (mm - 1);
1415  out += stride * (nn - 1);
1416  }
1417  }
1418 
1419 
1420 
1421  // evaluates the given shape data in 1d-3d using the tensor product
1422  // form assuming the symmetries of unit cell shape hessians for
1423  // finite elements in FEEvaluation
1424  template <int dim,
1425  int n_rows,
1426  int n_columns,
1427  typename Number,
1428  typename Number2>
1429  template <int direction, bool contract_over_rows, bool add>
1430  inline void
1432  dim,
1433  n_rows,
1434  n_columns,
1435  Number,
1436  Number2>::hessians(const Number in[],
1437  Number out[]) const
1438  {
1439  Assert(shape_hessians != nullptr, ExcNotInitialized());
1440  AssertIndexRange(direction, dim);
1441  constexpr int mm = contract_over_rows ? n_rows : n_columns;
1442  constexpr int nn = contract_over_rows ? n_columns : n_rows;
1443  constexpr int n_cols = nn / 2;
1444  constexpr int mid = mm / 2;
1445 
1446  constexpr int stride = Utilities::pow(n_columns, direction);
1447  constexpr int n_blocks1 = stride;
1448  constexpr int n_blocks2 =
1449  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1450 
1451  for (int i2 = 0; i2 < n_blocks2; ++i2)
1452  {
1453  for (int i1 = 0; i1 < n_blocks1; ++i1)
1454  {
1455  for (int col = 0; col < n_cols; ++col)
1456  {
1457  Number2 val0, val1;
1458  Number in0, in1, res0, res1;
1459  if (contract_over_rows == true)
1460  {
1461  val0 = shape_hessians[col];
1462  val1 = shape_hessians[nn - 1 - col];
1463  }
1464  else
1465  {
1466  val0 = shape_hessians[col * n_columns];
1467  val1 = shape_hessians[(col + 1) * n_columns - 1];
1468  }
1469  if (mid > 0)
1470  {
1471  in0 = in[0];
1472  in1 = in[stride * (mm - 1)];
1473  res0 = val0 * in0;
1474  res1 = val1 * in0;
1475  res0 += val1 * in1;
1476  res1 += val0 * in1;
1477  for (int ind = 1; ind < mid; ++ind)
1478  {
1479  if (contract_over_rows == true)
1480  {
1481  val0 = shape_hessians[ind * n_columns + col];
1482  val1 =
1483  shape_hessians[ind * n_columns + nn - 1 - col];
1484  }
1485  else
1486  {
1487  val0 = shape_hessians[col * n_columns + ind];
1488  val1 =
1489  shape_hessians[(col + 1) * n_columns - 1 - ind];
1490  }
1491  in0 = in[stride * ind];
1492  in1 = in[stride * (mm - 1 - ind)];
1493  res0 += val0 * in0;
1494  res1 += val1 * in0;
1495  res0 += val1 * in1;
1496  res1 += val0 * in1;
1497  }
1498  }
1499  else
1500  res0 = res1 = Number();
1501  if (mm % 2 == 1)
1502  {
1503  if (contract_over_rows == true)
1504  val0 = shape_hessians[mid * n_columns + col];
1505  else
1506  val0 = shape_hessians[col * n_columns + mid];
1507  in1 = val0 * in[stride * mid];
1508  res0 += in1;
1509  res1 += in1;
1510  }
1511  if (add)
1512  {
1513  out[stride * col] += res0;
1514  out[stride * (nn - 1 - col)] += res1;
1515  }
1516  else
1517  {
1518  out[stride * col] = res0;
1519  out[stride * (nn - 1 - col)] = res1;
1520  }
1521  }
1522  if (nn % 2 == 1)
1523  {
1524  Number2 val0;
1525  Number res0;
1526  if (contract_over_rows == true)
1527  val0 = shape_hessians[n_cols];
1528  else
1529  val0 = shape_hessians[n_cols * n_columns];
1530  if (mid > 0)
1531  {
1532  res0 = val0 * (in[0] + in[stride * (mm - 1)]);
1533  for (int ind = 1; ind < mid; ++ind)
1534  {
1535  if (contract_over_rows == true)
1536  val0 = shape_hessians[ind * n_columns + n_cols];
1537  else
1538  val0 = shape_hessians[n_cols * n_columns + ind];
1539  Number in1 = val0 * (in[stride * ind] +
1540  in[stride * (mm - 1 - ind)]);
1541  res0 += in1;
1542  }
1543  }
1544  else
1545  res0 = Number();
1546  if (mm % 2 == 1)
1547  {
1548  if (contract_over_rows == true)
1549  val0 = shape_hessians[mid * n_columns + n_cols];
1550  else
1551  val0 = shape_hessians[n_cols * n_columns + mid];
1552  res0 += val0 * in[stride * mid];
1553  }
1554  if (add)
1555  out[stride * n_cols] += res0;
1556  else
1557  out[stride * n_cols] = res0;
1558  }
1559 
1560  ++in;
1561  ++out;
1562  }
1563  in += stride * (mm - 1);
1564  out += stride * (nn - 1);
1565  }
1566  }
1567 
1568 
1569 
1570  template <int dim,
1571  int n_rows_static,
1572  int n_columns_static,
1573  typename Number,
1574  typename Number2,
1575  int direction,
1576  bool contract_over_rows,
1577  bool add,
1578  int type,
1579  bool one_line>
1580  inline void
1581  even_odd_apply(const int n_rows_in,
1582  const int n_columns_in,
1583  const Number2 *DEAL_II_RESTRICT shapes,
1584  const Number * in,
1585  Number * out)
1586  {
1587  static_assert(type < 3, "Only three variants type=0,1,2 implemented");
1588  static_assert(one_line == false || direction == dim - 1,
1589  "Single-line evaluation only works for direction=dim-1.");
1590 
1591  const int n_rows = n_rows_static == -1 ? n_rows_in : n_rows_static;
1592  const int n_columns =
1593  n_columns_static == -1 ? n_columns_in : n_columns_static;
1594 
1595  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
1596  in != out,
1597  ExcMessage("In-place operation only supported for "
1598  "n_rows==n_columns or single-line interpolation"));
1599 
1600  // We cannot statically assert that direction is less than dim, so must do
1601  // an additional dynamic check
1602  AssertIndexRange(direction, dim);
1603 
1604  const int nn = contract_over_rows ? n_columns : n_rows;
1605  const int mm = contract_over_rows ? n_rows : n_columns;
1606  constexpr int mm_static =
1607  contract_over_rows ? n_rows_static : n_columns_static;
1608  const int n_cols = nn / 2;
1609  const int mid = mm / 2;
1610  constexpr int mid_static = mm_static / 2;
1611  constexpr int max_mid = 15; // for non-templated execution
1612 
1613  Assert((n_rows_static != -1 && n_columns_static != -1) || mid <= max_mid,
1614  ExcNotImplemented());
1615 
1616  const int stride = Utilities::pow(n_columns, direction);
1617  const int n_blocks1 = one_line ? 1 : stride;
1618  const int n_blocks2 =
1619  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
1620 
1621  const int offset = (n_columns + 1) / 2;
1622 
1623  // this code may look very inefficient at first sight due to the many
1624  // different cases with if's at the innermost loop part, but all of the
1625  // conditionals can be evaluated at compile time because they are
1626  // templates, so the compiler should optimize everything away
1627  for (int i2 = 0; i2 < n_blocks2; ++i2)
1628  {
1629  for (int i1 = 0; i1 < n_blocks1; ++i1)
1630  {
1631  constexpr unsigned int mid_size =
1632  (n_rows_static == -1 || n_columns_static == -1) ?
1633  max_mid :
1634  (mid_static > 0 ? mid_static : 1);
1635  Number xp[mid_size], xm[mid_size];
1636  for (int i = 0; i < mid; ++i)
1637  {
1638  if (contract_over_rows == true && type == 1)
1639  {
1640  xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
1641  xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
1642  }
1643  else
1644  {
1645  xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
1646  xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
1647  }
1648  }
1649  Number xmid = in[stride * mid];
1650  for (int col = 0; col < n_cols; ++col)
1651  {
1652  Number r0, r1;
1653  if (mid > 0)
1654  {
1655  if (contract_over_rows == true)
1656  {
1657  r0 = shapes[col] * xp[0];
1658  r1 = shapes[(n_rows - 1) * offset + col] * xm[0];
1659  }
1660  else
1661  {
1662  r0 = shapes[col * offset] * xp[0];
1663  r1 = shapes[(n_rows - 1 - col) * offset] * xm[0];
1664  }
1665  for (int ind = 1; ind < mid; ++ind)
1666  {
1667  if (contract_over_rows == true)
1668  {
1669  r0 += shapes[ind * offset + col] * xp[ind];
1670  r1 += shapes[(n_rows - 1 - ind) * offset + col] *
1671  xm[ind];
1672  }
1673  else
1674  {
1675  r0 += shapes[col * offset + ind] * xp[ind];
1676  r1 += shapes[(n_rows - 1 - col) * offset + ind] *
1677  xm[ind];
1678  }
1679  }
1680  }
1681  else
1682  r0 = r1 = Number();
1683  if (mm % 2 == 1 && contract_over_rows == true)
1684  {
1685  if (type == 1)
1686  r1 += shapes[mid * offset + col] * xmid;
1687  else
1688  r0 += shapes[mid * offset + col] * xmid;
1689  }
1690  else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0 || mm == 3))
1691  r0 += shapes[col * offset + mid] * xmid;
1692 
1693  if (add)
1694  {
1695  out[stride * col] += r0 + r1;
1696  if (type == 1 && contract_over_rows == false)
1697  out[stride * (nn - 1 - col)] += r1 - r0;
1698  else
1699  out[stride * (nn - 1 - col)] += r0 - r1;
1700  }
1701  else
1702  {
1703  out[stride * col] = r0 + r1;
1704  if (type == 1 && contract_over_rows == false)
1705  out[stride * (nn - 1 - col)] = r1 - r0;
1706  else
1707  out[stride * (nn - 1 - col)] = r0 - r1;
1708  }
1709  }
1710  if (type == 0 && contract_over_rows == true && nn % 2 == 1 &&
1711  mm % 2 == 1 && mm > 3)
1712  {
1713  if (add)
1714  out[stride * n_cols] += shapes[mid * offset + n_cols] * xmid;
1715  else
1716  out[stride * n_cols] = shapes[mid * offset + n_cols] * xmid;
1717  }
1718  else if (contract_over_rows == true && nn % 2 == 1)
1719  {
1720  Number r0;
1721  if (mid > 0)
1722  {
1723  r0 = shapes[n_cols] * xp[0];
1724  for (int ind = 1; ind < mid; ++ind)
1725  r0 += shapes[ind * offset + n_cols] * xp[ind];
1726  }
1727  else
1728  r0 = Number();
1729  if (type != 1 && mm % 2 == 1)
1730  r0 += shapes[mid * offset + n_cols] * xmid;
1731 
1732  if (add)
1733  out[stride * n_cols] += r0;
1734  else
1735  out[stride * n_cols] = r0;
1736  }
1737  else if (contract_over_rows == false && nn % 2 == 1)
1738  {
1739  Number r0;
1740  if (mid > 0)
1741  {
1742  if (type == 1)
1743  {
1744  r0 = shapes[n_cols * offset] * xm[0];
1745  for (int ind = 1; ind < mid; ++ind)
1746  r0 += shapes[n_cols * offset + ind] * xm[ind];
1747  }
1748  else
1749  {
1750  r0 = shapes[n_cols * offset] * xp[0];
1751  for (int ind = 1; ind < mid; ++ind)
1752  r0 += shapes[n_cols * offset + ind] * xp[ind];
1753  }
1754  }
1755  else
1756  r0 = Number();
1757 
1758  if ((type == 0 || type == 2) && mm % 2 == 1)
1759  r0 += shapes[n_cols * offset + mid] * xmid;
1760 
1761  if (add)
1762  out[stride * n_cols] += r0;
1763  else
1764  out[stride * n_cols] = r0;
1765  }
1766  if (one_line == false)
1767  {
1768  in += 1;
1769  out += 1;
1770  }
1771  }
1772  if (one_line == false)
1773  {
1774  in += stride * (mm - 1);
1775  out += stride * (nn - 1);
1776  }
1777  }
1778  }
1779 
1780 
1781 
1813  template <int dim,
1814  int n_rows,
1815  int n_columns,
1816  typename Number,
1817  typename Number2>
1819  dim,
1820  n_rows,
1821  n_columns,
1822  Number,
1823  Number2>
1824  {
1825  static constexpr unsigned int n_rows_of_product =
1826  Utilities::pow(n_rows, dim);
1827  static constexpr unsigned int n_columns_of_product =
1828  Utilities::pow(n_columns, dim);
1829 
1836  : shape_values(nullptr)
1837  , shape_gradients(nullptr)
1838  , shape_hessians(nullptr)
1839  {}
1840 
1846  : shape_values(shape_values.begin())
1847  , shape_gradients(nullptr)
1848  , shape_hessians(nullptr)
1849  {
1850  AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
1851  }
1852 
1858  const AlignedVector<Number2> &shape_gradients,
1859  const AlignedVector<Number2> &shape_hessians,
1860  const unsigned int dummy1 = 0,
1861  const unsigned int dummy2 = 0)
1862  : shape_values(shape_values.begin())
1863  , shape_gradients(shape_gradients.begin())
1864  , shape_hessians(shape_hessians.begin())
1865  {
1866  // In this function, we allow for dummy pointers if some of values,
1867  // gradients or hessians should not be computed
1868  if (!shape_values.empty())
1869  AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
1870  if (!shape_gradients.empty())
1871  AssertDimension(shape_gradients.size(), n_rows * ((n_columns + 1) / 2));
1872  if (!shape_hessians.empty())
1873  AssertDimension(shape_hessians.size(), n_rows * ((n_columns + 1) / 2));
1874  (void)dummy1;
1875  (void)dummy2;
1876  }
1877 
1878  template <int direction, bool contract_over_rows, bool add>
1879  void
1880  values(const Number in[], Number out[]) const
1881  {
1882  Assert(shape_values != nullptr, ExcNotInitialized());
1883  apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
1884  }
1885 
1886  template <int direction, bool contract_over_rows, bool add>
1887  void
1888  gradients(const Number in[], Number out[]) const
1889  {
1890  Assert(shape_gradients != nullptr, ExcNotInitialized());
1891  apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
1892  }
1893 
1894  template <int direction, bool contract_over_rows, bool add>
1895  void
1896  hessians(const Number in[], Number out[]) const
1897  {
1898  Assert(shape_hessians != nullptr, ExcNotInitialized());
1899  apply<direction, contract_over_rows, add, 2>(shape_hessians, in, out);
1900  }
1901 
1902  template <int direction, bool contract_over_rows, bool add>
1903  void
1904  values_one_line(const Number in[], Number out[]) const
1905  {
1906  Assert(shape_values != nullptr, ExcNotInitialized());
1907  apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
1908  }
1909 
1910  template <int direction, bool contract_over_rows, bool add>
1911  void
1912  gradients_one_line(const Number in[], Number out[]) const
1913  {
1914  Assert(shape_gradients != nullptr, ExcNotInitialized());
1915  apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
1916  in,
1917  out);
1918  }
1919 
1920  template <int direction, bool contract_over_rows, bool add>
1921  void
1922  hessians_one_line(const Number in[], Number out[]) const
1923  {
1924  Assert(shape_hessians != nullptr, ExcNotInitialized());
1925  apply<direction, contract_over_rows, add, 2, true>(shape_hessians,
1926  in,
1927  out);
1928  }
1929 
1958  template <int direction,
1959  bool contract_over_rows,
1960  bool add,
1961  int type,
1962  bool one_line = false>
1963  static void
1964  apply(const Number2 *DEAL_II_RESTRICT shape_data,
1965  const Number * in,
1966  Number * out)
1967  {
1968  even_odd_apply<dim,
1969  n_rows,
1970  n_columns,
1971  Number,
1972  Number2,
1973  direction,
1974  contract_over_rows,
1975  add,
1976  type,
1977  one_line>(n_rows, n_columns, shape_data, in, out);
1978  }
1979 
1980  private:
1981  const Number2 *shape_values;
1982  const Number2 *shape_gradients;
1983  const Number2 *shape_hessians;
1984  };
1985 
1986 
1993  template <int dim, typename Number, typename Number2>
1994  struct EvaluatorTensorProduct<evaluate_evenodd, dim, 0, 0, Number, Number2>
1995  {
1997  : shape_values(nullptr)
1998  , shape_gradients(nullptr)
1999  , shape_hessians(nullptr)
2000  , n_rows(numbers::invalid_unsigned_int)
2001  , n_columns(numbers::invalid_unsigned_int)
2002  {}
2003 
2005  const unsigned int n_rows = 0,
2006  const unsigned int n_columns = 0)
2007  : shape_values(shape_values.begin())
2008  , shape_gradients(nullptr)
2009  , shape_hessians(nullptr)
2010  , n_rows(n_rows)
2011  , n_columns(n_columns)
2012  {
2013  AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
2014  }
2015 
2017  const AlignedVector<Number2> &shape_gradients,
2018  const AlignedVector<Number2> &shape_hessians,
2019  const unsigned int n_rows = 0,
2020  const unsigned int n_columns = 0)
2021  : shape_values(shape_values.begin())
2022  , shape_gradients(shape_gradients.begin())
2023  , shape_hessians(shape_hessians.begin())
2024  , n_rows(n_rows)
2025  , n_columns(n_columns)
2026  {
2027  if (!shape_values.empty())
2028  AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
2029  if (!shape_gradients.empty())
2030  AssertDimension(shape_gradients.size(), n_rows * ((n_columns + 1) / 2));
2031  if (!shape_hessians.empty())
2032  AssertDimension(shape_hessians.size(), n_rows * ((n_columns + 1) / 2));
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, 2>(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, 2, true>(shape_hessians,
2083  in,
2084  out);
2085  }
2086 
2087  template <int direction,
2088  bool contract_over_rows,
2089  bool add,
2090  int type,
2091  bool one_line = false>
2092  void
2093  apply(const Number2 *DEAL_II_RESTRICT shape_data,
2094  const Number * in,
2095  Number * out) const
2096  {
2097  even_odd_apply<dim,
2098  -1,
2099  -1,
2100  Number,
2101  Number2,
2102  direction,
2103  contract_over_rows,
2104  add,
2105  type,
2106  one_line>(n_rows, n_columns, shape_data, in, out);
2107  }
2108 
2109  private:
2110  const Number2 * shape_values;
2111  const Number2 * shape_gradients;
2112  const Number2 * shape_hessians;
2113  const unsigned int n_rows;
2114  const unsigned int n_columns;
2115  };
2116 
2117 
2118 
2147  template <int dim,
2148  int n_rows,
2149  int n_columns,
2150  typename Number,
2151  typename Number2>
2153  dim,
2154  n_rows,
2155  n_columns,
2156  Number,
2157  Number2>
2158  {
2159  static constexpr unsigned int n_rows_of_product =
2160  Utilities::pow(n_rows, dim);
2161  static constexpr unsigned int n_columns_of_product =
2162  Utilities::pow(n_columns, dim);
2163 
2170  : shape_values(nullptr)
2171  , shape_gradients(nullptr)
2172  , shape_hessians(nullptr)
2173  {}
2174 
2180  : shape_values(shape_values.begin())
2181  , shape_gradients(nullptr)
2182  , shape_hessians(nullptr)
2183  {}
2184 
2190  const AlignedVector<Number2> &shape_gradients,
2191  const AlignedVector<Number2> &shape_hessians,
2192  const unsigned int dummy1 = 0,
2193  const unsigned int dummy2 = 0)
2194  : shape_values(shape_values.begin())
2195  , shape_gradients(shape_gradients.begin())
2196  , shape_hessians(shape_hessians.begin())
2197  {
2198  (void)dummy1;
2199  (void)dummy2;
2200  }
2201 
2202  template <int direction, bool contract_over_rows, bool add>
2203  void
2204  values(const Number in[], Number out[]) const
2205  {
2206  Assert(shape_values != nullptr, ExcNotInitialized());
2207  apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
2208  }
2209 
2210  template <int direction, bool contract_over_rows, bool add>
2211  void
2212  gradients(const Number in[], Number out[]) const
2213  {
2214  Assert(shape_gradients != nullptr, ExcNotInitialized());
2215  apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
2216  }
2217 
2218  template <int direction, bool contract_over_rows, bool add>
2219  void
2220  hessians(const Number in[], Number out[]) const
2221  {
2222  Assert(shape_hessians != nullptr, ExcNotInitialized());
2223  apply<direction, contract_over_rows, add, 0>(shape_hessians, in, out);
2224  }
2225 
2226  template <int direction, bool contract_over_rows, bool add>
2227  void
2228  values_one_line(const Number in[], Number out[]) const
2229  {
2230  Assert(shape_values != nullptr, ExcNotInitialized());
2231  apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
2232  }
2233 
2234  template <int direction, bool contract_over_rows, bool add>
2235  void
2236  gradients_one_line(const Number in[], Number out[]) const
2237  {
2238  Assert(shape_gradients != nullptr, ExcNotInitialized());
2239  apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
2240  in,
2241  out);
2242  }
2243 
2244  template <int direction, bool contract_over_rows, bool add>
2245  void
2246  hessians_one_line(const Number in[], Number out[]) const
2247  {
2248  Assert(shape_hessians != nullptr, ExcNotInitialized());
2249  apply<direction, contract_over_rows, add, 0, true>(shape_hessians,
2250  in,
2251  out);
2252  }
2253 
2281  template <int direction,
2282  bool contract_over_rows,
2283  bool add,
2284  int type,
2285  bool one_line = false>
2286  static void
2287  apply(const Number2 *DEAL_II_RESTRICT shape_data,
2288  const Number * in,
2289  Number * out);
2290 
2291  private:
2292  const Number2 *shape_values;
2293  const Number2 *shape_gradients;
2294  const Number2 *shape_hessians;
2295  };
2296 
2297 
2298 
2299  template <int dim,
2300  int n_rows,
2301  int n_columns,
2302  typename Number,
2303  typename Number2>
2304  template <int direction,
2305  bool contract_over_rows,
2306  bool add,
2307  int type,
2308  bool one_line>
2309  inline void
2311  dim,
2312  n_rows,
2313  n_columns,
2314  Number,
2315  Number2>::apply(const Number2 *DEAL_II_RESTRICT shapes,
2316  const Number * in,
2317  Number * out)
2318  {
2319  static_assert(one_line == false || direction == dim - 1,
2320  "Single-line evaluation only works for direction=dim-1.");
2321  static_assert(
2322  type == 0 || type == 1,
2323  "Only types 0 and 1 implemented for evaluate_symmetric_hierarchical.");
2324  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
2325  in != out,
2326  ExcMessage("In-place operation only supported for "
2327  "n_rows==n_columns or single-line interpolation"));
2328 
2329  // We cannot statically assert that direction is less than dim, so must do
2330  // an additional dynamic check
2331  AssertIndexRange(direction, dim);
2332 
2333  constexpr int nn = contract_over_rows ? n_columns : n_rows;
2334  constexpr int mm = contract_over_rows ? n_rows : n_columns;
2335  constexpr int n_cols = nn / 2;
2336  constexpr int mid = mm / 2;
2337 
2338  constexpr int stride = Utilities::pow(n_columns, direction);
2339  constexpr int n_blocks1 = one_line ? 1 : stride;
2340  constexpr int n_blocks2 =
2341  Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
2342 
2343  // this code may look very inefficient at first sight due to the many
2344  // different cases with if's at the innermost loop part, but all of the
2345  // conditionals can be evaluated at compile time because they are
2346  // templates, so the compiler should optimize everything away
2347  for (int i2 = 0; i2 < n_blocks2; ++i2)
2348  {
2349  for (int i1 = 0; i1 < n_blocks1; ++i1)
2350  {
2351  if (contract_over_rows)
2352  {
2353  Number x[mm];
2354  for (unsigned int i = 0; i < mm; ++i)
2355  x[i] = in[stride * i];
2356  for (unsigned int col = 0; col < n_cols; ++col)
2357  {
2358  Number r0, r1;
2359  if (mid > 0)
2360  {
2361  r0 = shapes[col] * x[0];
2362  r1 = shapes[col + n_columns] * x[1];
2363  for (unsigned int ind = 1; ind < mid; ++ind)
2364  {
2365  r0 +=
2366  shapes[col + 2 * ind * n_columns] * x[2 * ind];
2367  r1 += shapes[col + (2 * ind + 1) * n_columns] *
2368  x[2 * ind + 1];
2369  }
2370  }
2371  else
2372  r0 = r1 = Number();
2373  if (mm % 2 == 1)
2374  r0 += shapes[col + (mm - 1) * n_columns] * x[mm - 1];
2375  if (add)
2376  {
2377  out[stride * col] += r0 + r1;
2378  if (type == 1)
2379  out[stride * (nn - 1 - col)] += r1 - r0;
2380  else
2381  out[stride * (nn - 1 - col)] += r0 - r1;
2382  }
2383  else
2384  {
2385  out[stride * col] = r0 + r1;
2386  if (type == 1)
2387  out[stride * (nn - 1 - col)] = r1 - r0;
2388  else
2389  out[stride * (nn - 1 - col)] = r0 - r1;
2390  }
2391  }
2392  if (nn % 2 == 1)
2393  {
2394  Number r0;
2395  const unsigned int shift = type == 1 ? 1 : 0;
2396  if (mid > 0)
2397  {
2398  r0 = shapes[n_cols + shift * n_columns] * x[shift];
2399  for (unsigned int ind = 1; ind < mid; ++ind)
2400  r0 += shapes[n_cols + (2 * ind + shift) * n_columns] *
2401  x[2 * ind + shift];
2402  }
2403  else
2404  r0 = 0;
2405  if (type != 1 && mm % 2 == 1)
2406  r0 += shapes[n_cols + (mm - 1) * n_columns] * x[mm - 1];
2407  if (add)
2408  out[stride * n_cols] += r0;
2409  else
2410  out[stride * n_cols] = r0;
2411  }
2412  }
2413  else
2414  {
2415  Number xp[mid + 1], xm[mid > 0 ? mid : 1];
2416  for (int i = 0; i < mid; ++i)
2417  if (type == 0)
2418  {
2419  xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
2420  xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
2421  }
2422  else
2423  {
2424  xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
2425  xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
2426  }
2427  if (mm % 2 == 1)
2428  xp[mid] = in[stride * mid];
2429  for (unsigned int col = 0; col < n_cols; ++col)
2430  {
2431  Number r0, r1;
2432  if (mid > 0)
2433  {
2434  r0 = shapes[2 * col * n_columns] * xp[0];
2435  r1 = shapes[(2 * col + 1) * n_columns] * xm[0];
2436  for (unsigned int ind = 1; ind < mid; ++ind)
2437  {
2438  r0 += shapes[2 * col * n_columns + ind] * xp[ind];
2439  r1 +=
2440  shapes[(2 * col + 1) * n_columns + ind] * xm[ind];
2441  }
2442  }
2443  else
2444  r0 = r1 = Number();
2445  if (mm % 2 == 1)
2446  {
2447  if (type == 1)
2448  r1 +=
2449  shapes[(2 * col + 1) * n_columns + mid] * xp[mid];
2450  else
2451  r0 += shapes[2 * col * n_columns + mid] * xp[mid];
2452  }
2453  if (add)
2454  {
2455  out[stride * (2 * col)] += r0;
2456  out[stride * (2 * col + 1)] += r1;
2457  }
2458  else
2459  {
2460  out[stride * (2 * col)] = r0;
2461  out[stride * (2 * col + 1)] = r1;
2462  }
2463  }
2464  if (nn % 2 == 1)
2465  {
2466  Number r0;
2467  if (mid > 0)
2468  {
2469  r0 = shapes[(nn - 1) * n_columns] * xp[0];
2470  for (unsigned int ind = 1; ind < mid; ++ind)
2471  r0 += shapes[(nn - 1) * n_columns + ind] * xp[ind];
2472  }
2473  else
2474  r0 = Number();
2475  if (mm % 2 == 1 && type == 0)
2476  r0 += shapes[(nn - 1) * n_columns + mid] * xp[mid];
2477  if (add)
2478  out[stride * (nn - 1)] += r0;
2479  else
2480  out[stride * (nn - 1)] = r0;
2481  }
2482  }
2483  if (one_line == false)
2484  {
2485  in += 1;
2486  out += 1;
2487  }
2488  }
2489  if (one_line == false)
2490  {
2491  in += stride * (mm - 1);
2492  out += stride * (nn - 1);
2493  }
2494  }
2495  }
2496 
2497 
2498 
2518  template <int dim,
2519  int n_rows,
2520  int n_columns,
2521  typename Number,
2522  int normal_dir,
2523  typename Number2>
2525  dim,
2526  n_rows,
2527  n_columns,
2528  Number,
2529  normal_dir,
2530  Number2>
2531  {
2532  static constexpr unsigned int n_rows_of_product =
2534  static constexpr unsigned int n_columns_of_product =
2536 
2542  : shape_values(nullptr)
2543  , shape_gradients(nullptr)
2544  , shape_hessians(nullptr)
2545  {}
2546 
2551  const AlignedVector<Number2> &shape_values,
2552  const AlignedVector<Number2> &shape_gradients,
2553  const AlignedVector<Number2> &shape_hessians,
2554  const unsigned int dummy1 = 0,
2555  const unsigned int dummy2 = 0)
2556  : shape_values(shape_values.begin())
2557  , shape_gradients(shape_gradients.begin())
2558  , shape_hessians(shape_hessians.begin())
2559  {
2560  // We can enter this function either for the apply() path that has
2561  // n_rows * n_columns entries or for the apply_face() path that only has
2562  // n_rows * 3 entries in the array. Since we cannot decide about the use
2563  // we must allow for both here.
2564  Assert(shape_values.size() == 0 ||
2565  shape_values.size() == n_rows * n_columns ||
2566  shape_values.size() == 3 * n_rows,
2567  ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
2568  Assert(shape_gradients.size() == 0 ||
2569  shape_gradients.size() == n_rows * n_columns,
2570  ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
2571  Assert(shape_hessians.size() == 0 ||
2572  shape_hessians.size() == n_rows * n_columns,
2573  ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
2574  (void)dummy1;
2575  (void)dummy2;
2576  }
2577 
2578  template <int direction, bool contract_over_rows, bool add>
2579  void
2580  values(const Number in[], Number out[]) const
2581  {
2582  apply<direction, contract_over_rows, add>(shape_values, in, out);
2583  }
2584 
2585  template <int direction, bool contract_over_rows, bool add>
2586  void
2587  gradients(const Number in[], Number out[]) const
2588  {
2589  apply<direction, contract_over_rows, add>(shape_gradients, in, out);
2590  }
2591 
2592  template <int direction, bool contract_over_rows, bool add>
2593  void
2594  hessians(const Number in[], Number out[]) const
2595  {
2596  apply<direction, contract_over_rows, add>(shape_hessians, in, out);
2597  }
2598 
2627  template <int direction,
2628  bool contract_over_rows,
2629  bool add,
2630  bool one_line = false>
2631  static void
2632  apply(const Number2 *DEAL_II_RESTRICT shape_data,
2633  const Number * in,
2634  Number * out);
2635 
2636  template <int face_direction,
2637  bool contract_onto_face,
2638  bool add,
2639  int max_derivative>
2640  void
2641  apply_face(const Number *DEAL_II_RESTRICT in,
2642  Number *DEAL_II_RESTRICT out) const;
2643 
2644  private:
2645  const Number2 *shape_values;
2646  const Number2 *shape_gradients;
2647  const Number2 *shape_hessians;
2648  };
2649 
2650  template <int dim,
2651  int n_rows,
2652  int n_columns,
2653  typename Number,
2654  int normal_dir,
2655  typename Number2>
2656  template <int direction, bool contract_over_rows, bool add, bool one_line>
2657  inline void
2660  dim,
2661  n_rows,
2662  n_columns,
2663  Number,
2664  normal_dir,
2665  Number2>::apply(const Number2 *DEAL_II_RESTRICT shape_data,
2666  const Number * in,
2667  Number * out)
2668  {
2669  static_assert(one_line == false || direction == dim - 1,
2670  "Single-line evaluation only works for direction=dim-1.");
2671  Assert(shape_data != nullptr,
2672  ExcMessage(
2673  "The given array shape_data must not be the null pointer!"));
2674  Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
2675  in != out,
2676  ExcMessage("In-place operation only supported for "
2677  "n_rows==n_columns or single-line interpolation"));
2678  AssertIndexRange(direction, dim);
2679  constexpr int mm = contract_over_rows ? n_rows : n_columns,
2680  nn = contract_over_rows ? n_columns : n_rows;
2681 
2682  constexpr int stride = Utilities::pow(n_columns, direction);
2683  constexpr int n_blocks1 = one_line ? 1 : stride;
2684 
2685  // The number of blocks depend on both direction and dimension.
2686  constexpr int n_blocks2 =
2687  (dim - direction - 1 == 0) ?
2688  1 :
2689  ((direction == normal_dir) ?
2690  Utilities::pow((n_rows - 1),
2691  (direction >= dim) ? 0 : dim - direction - 1) :
2692  (((direction < normal_dir) ? (n_rows + 1) : n_rows) *
2693  ((dim - direction == 3) ? n_rows : 1)));
2694 
2695  for (int i2 = 0; i2 < n_blocks2; ++i2)
2696  {
2697  for (int i1 = 0; i1 < n_blocks1; ++i1)
2698  {
2699  Number x[mm];
2700  for (int i = 0; i < mm; ++i)
2701  x[i] = in[stride * i];
2702 
2703  for (int col = 0; col < nn; ++col)
2704  {
2705  Number2 val0;
2706 
2707  if (contract_over_rows)
2708  val0 = shape_data[col];
2709  else
2710  val0 = shape_data[col * n_columns];
2711 
2712  Number res0 = val0 * x[0];
2713  for (int i = 1; i < mm; ++i)
2714  {
2715  if (contract_over_rows)
2716  val0 = shape_data[i * n_columns + col];
2717  else
2718  val0 = shape_data[col * n_columns + i];
2719 
2720  res0 += val0 * x[i];
2721  }
2722  if (add)
2723  out[stride * col] += res0;
2724 
2725  else
2726  out[stride * col] = res0;
2727  }
2728 
2729  if (one_line == false)
2730  {
2731  ++in;
2732  ++out;
2733  }
2734  }
2735  if (one_line == false)
2736  {
2737  in += stride * (mm - 1);
2738  out += stride * (nn - 1);
2739  }
2740  }
2741  }
2742 
2743  template <int dim,
2744  int n_rows,
2745  int n_columns,
2746  typename Number,
2747  int normal_dir,
2748  typename Number2>
2749  template <int face_direction,
2750  bool contract_onto_face,
2751  bool add,
2752  int max_derivative>
2753  inline void
2756  dim,
2757  n_rows,
2758  n_columns,
2759  Number,
2760  normal_dir,
2761  Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
2762  Number *DEAL_II_RESTRICT out) const
2763  {
2764  Assert(dim > 1 && dim < 4, ExcMessage("Only dim=2,3 supported"));
2765  static_assert(max_derivative >= 0 && max_derivative < 3,
2766  "Only derivative orders 0-2 implemented");
2767  Assert(shape_values != nullptr,
2768  ExcMessage(
2769  "The given array shape_values must not be the null pointer."));
2770 
2771  // Determine the number of blocks depending on the face and normaldirection,
2772  // as well as dimension.
2773  constexpr int n_blocks1 = (face_direction == normal_dir) ? (n_rows - 1) :
2774  ((face_direction == 0 && normal_dir == 2) ||
2775  (face_direction == 1 && normal_dir == 2) ||
2776  (face_direction == 2 && normal_dir == 1)) ?
2777  n_rows :
2778  (n_rows + 1);
2779  constexpr int n_blocks2 = (dim == 2) ?
2780  1 :
2781  ((face_direction == normal_dir) ?
2782  (n_rows - 1) :
2783  (((face_direction == 0 && normal_dir == 1) ||
2784  (face_direction == 1 && normal_dir == 0) ||
2785  (face_direction == 2 && normal_dir == 0)) ?
2786  n_rows :
2787  (n_rows + 1)));
2788 
2789  AssertIndexRange(face_direction, dim);
2790 
2791  constexpr int in_stride =
2792  (face_direction == normal_dir) ?
2793  Utilities::pow(n_rows - 1, face_direction) :
2794  ((face_direction == 0) ?
2795  1 :
2796  ((face_direction == 2) ?
2797  n_rows * (n_rows + 1) :
2798  ((face_direction == 1 && normal_dir == 0) ? (n_rows + 1) :
2799  n_rows)));
2800  constexpr int out_stride = n_blocks1 * n_blocks2;
2801 
2802  const Number *DEAL_II_RESTRICT shape_values = this->shape_values;
2803 
2804  for (int i2 = 0; i2 < n_blocks2; ++i2)
2805  {
2806  for (int i1 = 0; i1 < n_blocks1; ++i1)
2807  {
2808  if (contract_onto_face == true)
2809  {
2810  Number res0 = shape_values[0] * in[0];
2811  Number res1, res2;
2812 
2813  if (max_derivative > 0)
2814  res1 = shape_values[n_rows] * in[0];
2815 
2816  if (max_derivative > 1)
2817  res2 = shape_values[2 * n_rows] * in[0];
2818 
2819  for (int ind = 1; ind < n_rows; ++ind)
2820  {
2821  res0 += shape_values[ind] * in[in_stride * ind];
2822  if (max_derivative > 0)
2823  res1 += shape_values[ind + n_rows] * in[in_stride * ind];
2824 
2825  if (max_derivative > 1)
2826  res2 +=
2827  shape_values[ind + 2 * n_rows] * in[in_stride * ind];
2828  }
2829  if (add)
2830  {
2831  out[0] += res0;
2832 
2833  if (max_derivative > 0)
2834  out[out_stride] += res1;
2835 
2836  if (max_derivative > 1)
2837  out[2 * out_stride] += res2;
2838  }
2839  else
2840  {
2841  out[0] = res0;
2842 
2843  if (max_derivative > 0)
2844  out[out_stride] = res1;
2845 
2846  if (max_derivative > 1)
2847  out[2 * out_stride] = res2;
2848  }
2849  }
2850  else
2851  {
2852  for (int col = 0; col < n_rows; ++col)
2853  {
2854  if (add)
2855  out[col * in_stride] += shape_values[col] * in[0];
2856  else
2857  out[col * in_stride] = shape_values[col] * in[0];
2858 
2859  if (max_derivative > 0)
2860  out[col * in_stride] +=
2861  shape_values[col + n_rows] * in[out_stride];
2862 
2863  if (max_derivative > 1)
2864  out[col * in_stride] +=
2865  shape_values[col + 2 * n_rows] * in[2 * out_stride];
2866  }
2867  }
2868 
2869  // increment: in regular case, just go to the next point in
2870  // x-direction. If we are at the end of one chunk in x-dir, need
2871  // to jump over to the next layer in z-direction
2872  switch (face_direction)
2873  {
2874  case 0:
2875  in += contract_onto_face ? n_rows : 1;
2876  out += contract_onto_face ? 1 : n_rows;
2877  break;
2878 
2879  case 1:
2880  ++in;
2881  ++out;
2882  // faces 2 and 3 in 3d use local coordinate system zx, which
2883  // is the other way around compared to the tensor
2884  // product. Need to take that into account.
2885  if (dim == 3)
2886  {
2887  if (normal_dir == 0)
2888  {
2889  if (contract_onto_face)
2890  out += n_rows - 1;
2891  else
2892  in += n_rows - 1;
2893  }
2894  if (normal_dir == 1)
2895  {
2896  if (contract_onto_face)
2897  out += n_rows - 2;
2898  else
2899  in += n_rows - 2;
2900  }
2901  if (normal_dir == 2)
2902  {
2903  if (contract_onto_face)
2904  out += n_rows;
2905  else
2906  in += n_rows;
2907  }
2908  }
2909  break;
2910 
2911  case 2:
2912  ++in;
2913  ++out;
2914  break;
2915 
2916  default:
2917  Assert(false, ExcNotImplemented());
2918  }
2919  }
2920  if (face_direction == 1 && dim == 3)
2921  {
2922  // adjust for local coordinate system zx
2923  if (contract_onto_face)
2924  {
2925  if (normal_dir == 0)
2926  {
2927  in += (n_rows + 1) * (n_rows - 1);
2928  out -= n_rows * (n_rows + 1) - 1;
2929  }
2930  if (normal_dir == 1)
2931  {
2932  in += (n_rows - 1) * (n_rows - 1);
2933  out -= (n_rows - 1) * (n_rows - 1) - 1;
2934  }
2935  if (normal_dir == 2)
2936  {
2937  in += (n_rows - 1) * (n_rows);
2938  out -= (n_rows) * (n_rows + 1) - 1;
2939  }
2940  }
2941  else
2942  {
2943  if (normal_dir == 0)
2944  {
2945  out += (n_rows + 1) * (n_rows - 1);
2946  in -= n_rows * (n_rows + 1) - 1;
2947  }
2948  if (normal_dir == 1)
2949  {
2950  out += (n_rows - 1) * (n_rows - 1);
2951  in -= (n_rows - 1) * (n_rows - 1) - 1;
2952  }
2953  if (normal_dir == 2)
2954  {
2955  out += (n_rows - 1) * (n_rows);
2956  in -= (n_rows) * (n_rows + 1) - 1;
2957  }
2958  }
2959  }
2960  }
2961  }
2962 
2963 
2964 
2971  template <typename Number, typename Number2>
2973  {
2975  };
2976 
2977  template <int dim, typename Number, typename Number2>
2978  struct ProductTypeNoPoint<Point<dim, Number>, Number2>
2979  {
2981  };
2982 
2983 
2984 
3019  template <int dim, typename Number, typename Number2>
3020  inline std::pair<
3024  const std::vector<Polynomials::Polynomial<double>> &poly,
3025  const std::vector<Number> & values,
3026  const Point<dim, Number2> & p,
3027  const bool d_linear = false,
3028  const std::vector<unsigned int> & renumber = {})
3029  {
3030  static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
3031 
3032  using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
3033 
3034  // use `int` type for this variable and the loops below to inform the
3035  // compiler that the loops below will never overflow, which allows it to
3036  // generate more optimized code for the variable loop bounds in the
3037  // present context
3038  const int n_shapes = poly.size();
3039  AssertDimension(Utilities::pow(n_shapes, dim), values.size());
3040  Assert(renumber.empty() || renumber.size() == values.size(),
3041  ExcDimensionMismatch(renumber.size(), values.size()));
3042 
3043  // shortcut for linear interpolation to speed up evaluation
3044  if (d_linear)
3045  {
3046  AssertDimension(poly.size(), 2);
3047  for (unsigned int i = 0; i < renumber.size(); ++i)
3048  AssertDimension(renumber[i], i);
3049 
3050  if (dim == 1)
3051  {
3052  Tensor<1, dim, Number3> derivative;
3053  derivative[0] = values[1] - values[0];
3054  return std::make_pair((1. - p[0]) * values[0] + p[0] * values[1],
3055  derivative);
3056  }
3057  else if (dim == 2)
3058  {
3059  const Number2 x0 = 1. - p[0], x1 = p[0];
3060  const Number3 tmp0 = x0 * values[0] + x1 * values[1];
3061  const Number3 tmp1 = x0 * values[2] + x1 * values[3];
3062  const Number3 mapped = (1. - p[1]) * tmp0 + p[1] * tmp1;
3063  Tensor<1, dim, Number3> derivative;
3064  derivative[0] = (1. - p[1]) * (values[1] - values[0]) +
3065  p[1] * (values[3] - values[2]);
3066  derivative[1] = tmp1 - tmp0;
3067  return std::make_pair(mapped, derivative);
3068  }
3069  else if (dim == 3)
3070  {
3071  const Number2 x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1],
3072  z0 = 1. - p[2], z1 = p[2];
3073  const Number3 tmp0 = x0 * values[0] + x1 * values[1];
3074  const Number3 tmp1 = x0 * values[2] + x1 * values[3];
3075  const Number3 tmpy0 = y0 * tmp0 + y1 * tmp1;
3076  const Number3 tmp2 = x0 * values[4] + x1 * values[5];
3077  const Number3 tmp3 = x0 * values[6] + x1 * values[7];
3078  const Number3 tmpy1 = y0 * tmp2 + y1 * tmp3;
3079  const Number3 mapped = z0 * tmpy0 + z1 * tmpy1;
3080  Tensor<1, dim, Number3> derivative;
3081  derivative[2] = tmpy1 - tmpy0;
3082  derivative[1] = z0 * (tmp1 - tmp0) + z1 * (tmp3 - tmp2);
3083  derivative[0] =
3084  z0 *
3085  (y0 * (values[1] - values[0]) + y1 * (values[3] - values[2])) +
3086  z1 *
3087  (y0 * (values[5] - values[4]) + y1 * (values[7] - values[6]));
3088  return std::make_pair(mapped, derivative);
3089  }
3090  }
3091 
3092  AssertIndexRange(n_shapes, 200);
3094 
3095  // Evaluate 1d polynomials and their derivatives
3096  std::array<Number2, dim> point;
3097  for (unsigned int d = 0; d < dim; ++d)
3098  point[d] = p[d];
3099  for (int i = 0; i < n_shapes; ++i)
3100  poly[i].values_of_array(point, 1, &shapes[i][0]);
3101 
3102  // Go through the tensor product of shape functions and interpolate
3103  // with optimal algorithm
3104  std::pair<Number3, Tensor<1, dim, Number3>> result = {};
3105  for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
3106  {
3107  Number3 value_y = {}, deriv_x = {}, deriv_y = {};
3108  for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
3109  {
3110  // Interpolation + derivative x direction
3111  Number3 value = {}, deriv = {};
3112 
3113  // Distinguish the inner loop based on whether we have a
3114  // renumbering or not
3115  if (renumber.empty())
3116  for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3117  {
3118  value += shapes[i0][0][0] * values[i];
3119  deriv += shapes[i0][1][0] * values[i];
3120  }
3121  else
3122  for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3123  {
3124  value += shapes[i0][0][0] * values[renumber[i]];
3125  deriv += shapes[i0][1][0] * values[renumber[i]];
3126  }
3127 
3128  // Interpolation + derivative in y direction
3129  if (dim > 1)
3130  {
3131  value_y += value * shapes[i1][0][1];
3132  deriv_x += deriv * shapes[i1][0][1];
3133  deriv_y += value * shapes[i1][1][1];
3134  }
3135  else
3136  {
3137  result.first = value;
3138  result.second[0] = deriv;
3139  }
3140  }
3141  if (dim == 3)
3142  {
3143  // Interpolation + derivative in z direction
3144  result.first += value_y * shapes[i2][0][2];
3145  result.second[0] += deriv_x * shapes[i2][0][2];
3146  result.second[1] += deriv_y * shapes[i2][0][2];
3147  result.second[2] += value_y * shapes[i2][1][2];
3148  }
3149  else if (dim == 2)
3150  {
3151  result.first = value_y;
3152  result.second[0] = deriv_x;
3153  result.second[1] = deriv_y;
3154  }
3155  }
3156 
3157  return result;
3158  }
3159 
3160 
3161 
3162  template <int dim, typename Number, typename Number2>
3165  const std::vector<Polynomials::Polynomial<double>> &poly,
3166  const std::vector<Number> & values,
3167  const Point<dim, Number2> & p,
3168  const std::vector<unsigned int> & renumber = {})
3169  {
3170  static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
3171 
3172  using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
3173 
3174  // use `int` type for this variable and the loops below to inform the
3175  // compiler that the loops below will never overflow, which allows it to
3176  // generate more optimized code for the variable loop bounds in the
3177  // present context
3178  const int n_shapes = poly.size();
3179  AssertDimension(Utilities::pow(n_shapes, dim), values.size());
3180  Assert(renumber.empty() || renumber.size() == values.size(),
3181  ExcDimensionMismatch(renumber.size(), values.size()));
3182 
3183  AssertIndexRange(n_shapes, 200);
3185 
3186  // Evaluate 1d polynomials and their derivatives
3187  std::array<Number2, dim> point;
3188  for (unsigned int d = 0; d < dim; ++d)
3189  point[d] = p[d];
3190  for (int i = 0; i < n_shapes; ++i)
3191  poly[i].values_of_array(point, 2, &shapes[i][0]);
3192 
3193  // Go through the tensor product of shape functions and interpolate
3194  // with optimal algorithm
3196  for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
3197  {
3198  Number3 value_y = {}, deriv_x = {}, deriv_y = {}, deriv_xx = {},
3199  deriv_xy = {}, deriv_yy = {};
3200  for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
3201  {
3202  // Interpolation + derivative x direction
3203  Number3 value = {}, deriv_1 = {}, deriv_2 = {};
3204 
3205  // Distinguish the inner loop based on whether we have a
3206  // renumbering or not
3207  if (renumber.empty())
3208  for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3209  {
3210  value += shapes[i0][0][0] * values[i];
3211  deriv_1 += shapes[i0][1][0] * values[i];
3212  deriv_2 += shapes[i0][2][0] * values[i];
3213  }
3214  else
3215  for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3216  {
3217  value += shapes[i0][0][0] * values[renumber[i]];
3218  deriv_1 += shapes[i0][1][0] * values[renumber[i]];
3219  deriv_2 += shapes[i0][2][0] * values[renumber[i]];
3220  }
3221 
3222  // Interpolation + derivative in y direction
3223  if (dim > 1)
3224  {
3225  if (dim > 2)
3226  {
3227  value_y += value * shapes[i1][0][1];
3228  deriv_x += deriv_1 * shapes[i1][0][1];
3229  deriv_y += value * shapes[i1][1][1];
3230  }
3231  deriv_xx += deriv_2 * shapes[i1][0][1];
3232  deriv_xy += deriv_1 * shapes[i1][1][1];
3233  deriv_yy += value * shapes[i1][2][1];
3234  }
3235  else
3236  {
3237  result[0][0] = deriv_2;
3238  }
3239  }
3240  if (dim == 3)
3241  {
3242  // Interpolation + derivative in z direction
3243  result[0][0] += deriv_xx * shapes[i2][0][2];
3244  result[0][1] += deriv_xy * shapes[i2][0][2];
3245  result[0][2] += deriv_x * shapes[i2][1][2];
3246  result[1][1] += deriv_yy * shapes[i2][0][2];
3247  result[1][2] += deriv_y * shapes[i2][1][2];
3248  result[2][2] += value_y * shapes[i2][2][2];
3249  }
3250  else if (dim == 2)
3251  {
3252  result[0][0] = deriv_xx;
3253  result[1][0] = deriv_xy;
3254  result[1][1] = deriv_yy;
3255  }
3256  }
3257 
3258  return result;
3259  }
3260 
3261 
3262 
3266  template <int dim, typename Number, typename Number2>
3267  inline void
3269  const std::vector<Polynomials::Polynomial<double>> &poly,
3270  const Number2 & value,
3271  const Tensor<1, dim, Number2> & gradient,
3272  const Point<dim, Number> & p,
3274  const std::vector<unsigned int> & renumber = {})
3275  {
3276  static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
3277 
3278  // as in evaluate, use `int` type to produce better code in this context
3279  const int n_shapes = poly.size();
3280  AssertDimension(Utilities::pow(n_shapes, dim), values.size());
3281  Assert(renumber.empty() || renumber.size() == values.size(),
3282  ExcDimensionMismatch(renumber.size(), values.size()));
3283 
3284  AssertIndexRange(n_shapes, 200);
3286 
3287  // Evaluate 1d polynomials and their derivatives
3288  std::array<Number, dim> point;
3289  for (unsigned int d = 0; d < dim; ++d)
3290  point[d] = p[d];
3291  for (int i = 0; i < n_shapes; ++i)
3292  poly[i].values_of_array(point, 1, &shapes[i][0]);
3293 
3294  // Implement the transpose of the function above
3295  for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
3296  {
3297  const Number2 test_value_z =
3298  dim > 2 ?
3299  (value * shapes[i2][0][2] + gradient[2] * shapes[i2][1][2]) :
3300  value;
3301  const Number2 test_grad_x =
3302  dim > 2 ? gradient[0] * shapes[i2][0][2] : gradient[0];
3303  const Number2 test_grad_y = dim > 2 ?
3304  gradient[1] * shapes[i2][0][2] :
3305  (dim > 1 ? gradient[1] : Number2());
3306  for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
3307  {
3308  const Number2 test_value_y = dim > 1 ?
3309  (test_value_z * shapes[i1][0][1] +
3310  test_grad_y * shapes[i1][1][1]) :
3311  test_value_z;
3312  const Number2 test_grad_xy =
3313  dim > 1 ? test_grad_x * shapes[i1][0][1] : test_grad_x;
3314  if (renumber.empty())
3315  for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3316  values[i] += shapes[i0][0][0] * test_value_y +
3317  shapes[i0][1][0] * test_grad_xy;
3318  else
3319  for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
3320  values[renumber[i]] += shapes[i0][0][0] * test_value_y +
3321  shapes[i0][1][0] * test_grad_xy;
3322  }
3323  }
3324  }
3325 
3326 
3327  template <int dim, int n_points_1d_template, typename Number>
3328  inline void
3329  weight_fe_q_dofs_by_entity(const Number * weights,
3330  const unsigned int n_components,
3331  const int n_points_1d_non_template,
3332  Number * data)
3333  {
3334  const int n_points_1d = n_points_1d_template != -1 ?
3335  n_points_1d_template :
3336  n_points_1d_non_template;
3337 
3338  Assert(n_points_1d > 0, ExcNotImplemented());
3339  Assert(n_points_1d < 100, ExcNotImplemented());
3340 
3341  unsigned int compressed_index[100];
3342  compressed_index[0] = 0;
3343  for (int i = 1; i < n_points_1d - 1; ++i)
3344  compressed_index[i] = 1;
3345  compressed_index[n_points_1d - 1] = 2;
3346 
3347  for (unsigned int c = 0; c < n_components; ++c)
3348  for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3349  for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
3350  {
3351  const unsigned int shift =
3352  9 * compressed_index[k] + 3 * compressed_index[j];
3353  data[0] *= weights[shift];
3354  // loop bound as int avoids compiler warnings in case n_points_1d
3355  // == 1 (polynomial degree 0)
3356  for (int i = 1; i < n_points_1d - 1; ++i)
3357  data[i] *= weights[shift + 1];
3358  data[n_points_1d - 1] *= weights[shift + 2];
3359  data += n_points_1d;
3360  }
3361  }
3362 
3363 
3364  template <int dim, int n_points_1d_template, typename Number>
3365  inline void
3366  weight_fe_q_dofs_by_entity_shifted(const Number * weights,
3367  const unsigned int n_components,
3368  const int n_points_1d_non_template,
3369  Number * data)
3370  {
3371  const int n_points_1d = n_points_1d_template != -1 ?
3372  n_points_1d_template :
3373  n_points_1d_non_template;
3374 
3375  Assert((n_points_1d % 2) == 1,
3376  ExcMessage("The function can only with add number of points"));
3377  Assert(n_points_1d > 0, ExcNotImplemented());
3378  Assert(n_points_1d < 100, ExcNotImplemented());
3379 
3380  const unsigned int n_inside_1d = n_points_1d / 2;
3381 
3382  unsigned int compressed_index[100];
3383 
3384  unsigned int c = 0;
3385  for (int i = 0; i < n_inside_1d; ++i)
3386  compressed_index[c++] = 0;
3387  compressed_index[c++] = 1;
3388  for (int i = 0; i < n_inside_1d; ++i)
3389  compressed_index[c++] = 2;
3390 
3391  for (unsigned int c = 0; c < n_components; ++c)
3392  for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3393  for (int j = 0; j < (dim > 1 ? n_points_1d : 1); ++j)
3394  {
3395  const unsigned int shift =
3396  9 * compressed_index[k] + 3 * compressed_index[j];
3397 
3398  unsigned int c = 0;
3399  for (int i = 0; i < n_inside_1d; ++i)
3400  data[c++] *= weights[shift];
3401  data[c++] *= weights[shift + 1];
3402  for (int i = 0; i < n_inside_1d; ++i)
3403  data[c++] *= weights[shift + 2];
3404  data += n_points_1d;
3405  }
3406  }
3407 
3408 
3409  template <int dim, int n_points_1d_template, typename Number>
3410  inline bool
3412  const unsigned int n_components,
3413  const int n_points_1d_non_template,
3414  Number * weights)
3415  {
3416  const int n_points_1d = n_points_1d_template != -1 ?
3417  n_points_1d_template :
3418  n_points_1d_non_template;
3419 
3420  Assert(n_points_1d > 0, ExcNotImplemented());
3421  Assert(n_points_1d < 100, ExcNotImplemented());
3422 
3423  unsigned int compressed_index[100];
3424  compressed_index[0] = 0;
3425  for (int i = 1; i < n_points_1d - 1; ++i)
3426  compressed_index[i] = 1;
3427  compressed_index[n_points_1d - 1] = 2;
3428 
3429  // Insert the number data into a storage position for weight,
3430  // ensuring that the array has either not been touched before
3431  // or the previous content is the same. In case the previous
3432  // content has a different value, we exit this function and
3433  // signal to outer functions that the compression was not possible.
3434  const auto check_and_set = [](Number &weight, const Number &data) {
3435  if (weight == Number(-1.0) || weight == data)
3436  {
3437  weight = data;
3438  return true; // success for the entry
3439  }
3440 
3441  return false; // failure for the entry
3442  };
3443 
3444  for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
3445  weights[c] = Number(-1.0);
3446 
3447  for (unsigned int c = 0; c < n_components; ++c)
3448  for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3449  for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
3450  ++j, data += n_points_1d)
3451  {
3452  const unsigned int shift =
3453  9 * compressed_index[k] + 3 * compressed_index[j];
3454 
3455  if (!check_and_set(weights[shift], data[0]))
3456  return false; // failure
3457 
3458  for (int i = 1; i < n_points_1d - 1; ++i)
3459  if (!check_and_set(weights[shift + 1], data[i]))
3460  return false; // failure
3461 
3462  if (!check_and_set(weights[shift + 2], data[n_points_1d - 1]))
3463  return false; // failure
3464  }
3465 
3466  return true; // success
3467  }
3468 
3469 
3470  template <int dim, int n_points_1d_template, typename Number>
3471  inline bool
3473  const Number * data,
3474  const unsigned int n_components,
3475  const int n_points_1d_non_template,
3476  Number * weights)
3477  {
3478  const int n_points_1d = n_points_1d_template != -1 ?
3479  n_points_1d_template :
3480  n_points_1d_non_template;
3481 
3482  Assert((n_points_1d % 2) == 1,
3483  ExcMessage("The function can only with add number of points"));
3484  Assert(n_points_1d > 0, ExcNotImplemented());
3485  Assert(n_points_1d < 100, ExcNotImplemented());
3486 
3487  const unsigned int n_inside_1d = n_points_1d / 2;
3488 
3489  unsigned int compressed_index[100];
3490 
3491  unsigned int c = 0;
3492  for (int i = 0; i < n_inside_1d; ++i)
3493  compressed_index[c++] = 0;
3494  compressed_index[c++] = 1;
3495  for (int i = 0; i < n_inside_1d; ++i)
3496  compressed_index[c++] = 2;
3497 
3498  // Insert the number data into a storage position for weight,
3499  // ensuring that the array has either not been touched before
3500  // or the previous content is the same. In case the previous
3501  // content has a different value, we exit this function and
3502  // signal to outer functions that the compression was not possible.
3503  const auto check_and_set = [](Number &weight, const Number &data) {
3504  if (weight == Number(-1.0) || weight == data)
3505  {
3506  weight = data;
3507  return true; // success for the entry
3508  }
3509 
3510  return false; // failure for the entry
3511  };
3512 
3513  for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
3514  weights[c] = Number(-1.0);
3515 
3516  for (unsigned int comp = 0; comp < n_components; ++comp)
3517  for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
3518  for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
3519  ++j, data += n_points_1d)
3520  {
3521  const unsigned int shift =
3522  9 * compressed_index[k] + 3 * compressed_index[j];
3523 
3524  unsigned int c = 0;
3525 
3526  for (int i = 0; i < n_inside_1d; ++i)
3527  if (!check_and_set(weights[shift], data[c++]))
3528  return false; // failure
3529 
3530  if (!check_and_set(weights[shift + 1], data[c++]))
3531  return false; // failure
3532 
3533  for (int i = 0; i < n_inside_1d; ++i)
3534  if (!check_and_set(weights[shift + 2], data[c++]))
3535  return false; // failure
3536  }
3537 
3538  return true; // success
3539  }
3540 
3541 
3542 } // end of namespace internal
3543 
3544 
3546 
3547 #endif
bool empty() const
size_type size() const
Definition: point.h:110
Definition: tensor.h:516
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_RESTRICT
Definition: config.h:109
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
static ::ExceptionBase & ExcMessage(std::string arg1)
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2131
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:189
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:447
T fixed_power(const T t)
Definition: utilities.h:966
void integrate_add_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double >> &poly, const Number2 &value, const Tensor< 1, dim, Number2 > &gradient, const Point< dim, Number > &p, AlignedVector< Number2 > &values, const std::vector< unsigned int > &renumber={})
void weight_fe_q_dofs_by_entity(const Number *weights, const unsigned int n_components, const int n_points_1d_non_template, Number *data)
SymmetricTensor< 2, dim, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_hessian(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const std::vector< unsigned int > &renumber={})
std::pair< typename ProductTypeNoPoint< Number, Number2 >::type, Tensor< 1, dim, typename ProductTypeNoPoint< Number, Number2 >::type > > evaluate_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
void even_odd_apply(const int n_rows_in, const int n_columns_in, const Number2 *DEAL_II_RESTRICT shapes, const Number *in, Number *out)
void weight_fe_q_dofs_by_entity_shifted(const Number *weights, const unsigned int n_components, const int n_points_1d_non_template, Number *data)
bool compute_weights_fe_q_dofs_by_entity(const Number *data, const unsigned int n_components, const int n_points_1d_non_template, Number *weights)
bool compute_weights_fe_q_dofs_by_entity_shifted(const Number *data, const unsigned int n_components, const int n_points_1d_non_template, Number *weights)
static const unsigned int invalid_unsigned_int
Definition: types.h:213
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:36
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
EvaluatorTensorProductAnisotropic(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int n_rows=0, const unsigned int n_columns=0)
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const unsigned int n_rows=0, const unsigned int n_columns=0)
void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out) const
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)
static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, Number *out)
EvaluatorTensorProduct(const Number2 *shape_values, const Number2 *shape_gradients, const Number2 *shape_hessians, const unsigned int n_rows, const unsigned int n_columns)
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int n_rows, const unsigned int n_columns)
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
EvaluatorTensorProduct(const AlignedVector< Number2 > &shape_values, const AlignedVector< Number2 > &shape_gradients, const AlignedVector< Number2 > &shape_hessians, const unsigned int dummy1=0, const unsigned int dummy2=0)
typename ProductType< Tensor< 1, dim, Number >, Number2 >::type type
typename ProductType< Number, Number2 >::type type