Reference documentation for deal.II version GIT b6bf1e606d 2022-08-11 15:25:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
evaluation_kernels_hanging_nodes.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_evaluation_kernels_hanging_nodes_h
18 #define dealii_matrix_free_evaluation_kernels_hanging_nodes_h
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/ndarray.h>
24 
27 
28 
30 
31 #ifdef DEBUG
32 # define DEAL_II_ALWAYS_INLINE_RELEASE
33 #else
34 # define DEAL_II_ALWAYS_INLINE_RELEASE DEAL_II_ALWAYS_INLINE
35 #endif
36 
37 
38 
39 namespace internal
40 {
42  {
43  scalar,
45  };
46 
47 
48 
53  enum class VectorizationTypes
54  {
58  index,
62  group,
67  mask,
71  sorted
72  };
73 
74 
75 
77  int dim,
78  int fe_degree,
79  typename Number>
81 
82 
83 
84  template <int dim, int fe_degree, typename Number>
87  dim,
88  fe_degree,
89  Number>
90  {
91  private:
92  template <int structdim, unsigned int direction, bool transpose>
93  static void
94  interpolate(const unsigned int offset,
95  const unsigned int outer_stride,
96  const unsigned int given_degree,
97  const Number mask_weight,
98  const Number mask_write,
99  const Number *DEAL_II_RESTRICT weights,
100  Number *DEAL_II_RESTRICT values)
101  {
102  static constexpr unsigned int max_n_points_1D = 40;
103 
104  static_assert(structdim == 1 || structdim == 2,
105  "Only 1D and 2D interpolation implemented");
106  Number temp[fe_degree != -1 ? fe_degree + 1 : max_n_points_1D];
107 
108  const unsigned int points =
109  (fe_degree != -1 ? fe_degree : given_degree) + 1;
110 
111  AssertIndexRange(points, max_n_points_1D);
112 
113  const unsigned int stride = Utilities::pow(points, direction);
114 
115  const unsigned int end_of_outer_loop = structdim == 1 ? 2 : points - 1;
116  for (unsigned int g = 1; g < end_of_outer_loop; ++g)
117  {
118  const unsigned int my_offset =
119  offset + (structdim > 1 ? g * outer_stride : 0);
120 
121  // extract values to interpolate
122  for (unsigned int k = 0; k < points; ++k)
123  temp[k] = values[my_offset + k * stride];
124 
125  // perform interpolation point by point and write back
126  for (unsigned int k = 0; k < points / 2; ++k)
127  {
128  const unsigned int kmirror = points - 1 - k;
129  Number sum0 = Number(), sum1 = Number(), sum2 = Number(),
130  sum3 = Number();
131  for (unsigned int h = 0; h < points; ++h)
132  {
133  const unsigned int hmirror = points - 1 - h;
134  // load from both sides of the interpolation matrix to
135  // reflect symmetry between the two subfaces along that
136  // direction
137  const Number w0 = weights[(transpose ? 1 : points) * kmirror +
138  (transpose ? points : 1) * hmirror];
139  const Number w1 = weights[(transpose ? 1 : points) * k +
140  (transpose ? points : 1) * h];
141  sum0 += temp[h] * w0;
142  sum1 += temp[h] * w1;
143  sum2 += temp[hmirror] * w1;
144  sum3 += temp[hmirror] * w0;
145  }
146  values[my_offset + k * stride] =
147  temp[k] +
148  mask_write * (sum0 + mask_weight * (sum1 - sum0) - temp[k]);
149  values[my_offset + kmirror * stride] =
150  temp[kmirror] +
151  mask_write *
152  (sum2 + mask_weight * (sum3 - sum2) - temp[kmirror]);
153  }
154 
155  // cleanup case
156  if (points % 2)
157  {
158  const unsigned int k = points / 2;
159  Number sum0 = temp[k] * weights[(transpose ? 1 : points) * k +
160  (transpose ? points : 1) * k],
161  sum1 = sum0;
162  for (unsigned int h = 0; h < points / 2; ++h)
163  {
164  const unsigned int hmirror = points - 1 - h;
165  const Number w0 = weights[(transpose ? 1 : points) * k +
166  (transpose ? points : 1) * hmirror];
167  const Number w1 = weights[(transpose ? 1 : points) * k +
168  (transpose ? points : 1) * h];
169  sum0 += temp[h] * w0;
170  sum0 += temp[hmirror] * w1;
171  sum1 += temp[h] * w1;
172  sum1 += temp[hmirror] * w0;
173  }
174  values[my_offset + k * stride] =
175  temp[k] +
176  mask_write * (sum0 + mask_weight * (sum1 - sum0) - temp[k]);
177  }
178  }
179  }
180 
181  public:
182  template <bool transpose>
183  static void
185  const unsigned int n_components,
186  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
188  Number::size()> & constraint_mask,
189  Number * values)
190  {
191  const unsigned int given_degree =
192  fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
193 
194  const Number *DEAL_II_RESTRICT weights =
195  shape_info.data.front().subface_interpolation_matrices[0].data();
196 
197  const unsigned int points = given_degree + 1;
198  const unsigned int n_dofs = shape_info.dofs_per_component_on_cell;
199 
200  if (dim == 2)
201  {
202  ::ndarray<Number, 2> mask_weights = {};
203  ::ndarray<Number, 2, 2> mask_write = {};
204  ::ndarray<bool, 2, 2> do_face = {};
205 
206  for (unsigned int v = 0; v < Number::size(); ++v)
207  {
208  const auto kind = constraint_mask[v];
209  const bool subcell_x = (kind >> 0) & 1;
210  const bool subcell_y = (kind >> 1) & 1;
211  const bool face_x = (kind >> 3) & 1;
212  const bool face_y = (kind >> 4) & 1;
213 
214  if (face_y)
215  {
216  const unsigned int side = !subcell_y;
217  mask_write[0][side][v] = 1;
218  do_face[0][side] = true;
219  mask_weights[0][v] = subcell_x;
220  }
221 
222  if (face_x)
223  {
224  const unsigned int side = !subcell_x;
225  mask_write[1][side][v] = 1;
226  do_face[1][side] = true;
227  mask_weights[1][v] = subcell_y;
228  }
229  }
230 
231  // x direction
232  {
233  const std::array<unsigned int, 2> offsets = {
234  {0, (points - 1) * points}};
235  for (unsigned int c = 0; c < n_components; ++c)
236  for (unsigned int face = 0; face < 2; ++face)
237  if (do_face[0][face])
238  interpolate<1, 0, transpose>(offsets[face],
239  0,
240  given_degree,
241  mask_weights[0],
242  mask_write[0][face],
243  weights,
244  values + c * n_dofs);
245  }
246 
247  // y direction
248  {
249  const std::array<unsigned int, 2> offsets = {{0, points - 1}};
250  for (unsigned int c = 0; c < n_components; ++c)
251  for (unsigned int face = 0; face < 2; ++face)
252  if (do_face[1][face])
253  interpolate<1, 1, transpose>(offsets[face],
254  0,
255  given_degree,
256  mask_weights[1],
257  mask_write[1][face],
258  weights,
259  values + c * n_dofs);
260  }
261  }
262  else if (dim == 3)
263  {
264  const unsigned int p0 = 0;
265  const unsigned int p1 = points - 1;
266  const unsigned int p2 = points * points - points;
267  const unsigned int p3 = points * points - 1;
268  const unsigned int p4 = points * points * points - points * points;
269  const unsigned int p5 =
270  points * points * points - points * points + points - 1;
271  const unsigned int p6 = points * points * points - points;
272 
273  ::ndarray<bool, 3, 4> process_edge = {};
274  ::ndarray<bool, 3, 4> process_face = {};
275  ::ndarray<Number, 3, 4> mask_edge = {};
276  ::ndarray<Number, 3, 4> mask_face = {};
277  ::ndarray<Number, 3> mask_weights = {};
278 
279  for (unsigned int v = 0; v < Number::size(); ++v)
280  {
281  const auto kind = constraint_mask[v];
282 
283  const bool subcell_x = (kind >> 0) & 1;
284  const bool subcell_y = (kind >> 1) & 1;
285  const bool subcell_z = (kind >> 2) & 1;
286  const bool face_x = ((kind >> 3) & 1) ? (kind >> 5) & 1 : 0;
287  const bool face_y = ((kind >> 3) & 1) ? (kind >> 6) & 1 : 0;
288  const bool face_z = ((kind >> 3) & 1) ? (kind >> 7) & 1 : 0;
289  const bool edge_x = ((kind >> 4) & 1) ? (kind >> 5) & 1 : 0;
290  const bool edge_y = ((kind >> 4) & 1) ? (kind >> 6) & 1 : 0;
291  const bool edge_z = ((kind >> 4) & 1) ? (kind >> 7) & 1 : 0;
292 
293  if (subcell_x)
294  mask_weights[0][v] = 1;
295  if (subcell_y)
296  mask_weights[1][v] = 1;
297  if (subcell_z)
298  mask_weights[2][v] = 1;
299 
300  if (face_x)
301  {
302  const unsigned int side = !subcell_x;
303 
304  mask_face[1][side][v] = process_face[1][side] = true;
305  mask_edge[1][side][v] = process_edge[1][side] = true;
306  mask_edge[1][2 + side][v] = process_edge[1][2 + side] = true;
307  mask_face[2][side][v] = process_face[2][side] = true;
308  mask_edge[2][side][v] = process_edge[2][side] = true;
309  mask_edge[2][2 + side][v] = process_edge[2][2 + side] = true;
310  }
311  if (face_y)
312  {
313  const unsigned int side = !subcell_y;
314 
315  mask_face[0][side][v] = process_face[0][side] = true;
316  mask_edge[0][side][v] = process_edge[0][side] = true;
317  mask_edge[0][2 + side][v] = process_edge[0][2 + side] = true;
318  mask_face[2][2 + side][v] = process_face[2][2 + side] = true;
319  mask_edge[2][2 * side][v] = process_edge[2][2 * side] = true;
320  mask_edge[2][2 * side + 1][v] =
321  process_edge[2][2 * side + 1] = true;
322  }
323  if (face_z)
324  {
325  const unsigned int side = !subcell_z;
326 
327  mask_face[0][2 + side][v] = process_face[0][2 + side] = true;
328  mask_edge[0][2 * side][v] = process_edge[0][2 * side] = true;
329  mask_edge[0][2 * side + 1][v] =
330  process_edge[0][2 * side + 1] = true;
331  mask_face[1][2 + side][v] = process_face[1][2 + side] = true;
332  mask_edge[1][2 * side][v] = process_edge[1][2 * side] = true;
333  mask_edge[1][2 * side + 1][v] =
334  process_edge[1][2 * side + 1] = true;
335  }
336  if (edge_x)
337  {
338  const unsigned int index = (!subcell_z) * 2 + (!subcell_y);
339  mask_edge[0][index][v] = process_edge[0][index] = true;
340  }
341  if (edge_y)
342  {
343  const unsigned int index = (!subcell_z) * 2 + (!subcell_x);
344  mask_edge[1][index][v] = process_edge[1][index] = true;
345  }
346  if (edge_z)
347  {
348  const unsigned int index = (!subcell_y) * 2 + (!subcell_x);
349  mask_edge[2][index][v] = process_edge[2][index] = true;
350  }
351  }
352 
353  // direction 0:
354  if (given_degree > 1)
355  {
356  const std::array<unsigned int, 4> face_offsets = {
357  {p0, p2, p0, p4}};
358  const std::array<unsigned int, 2> outer_strides = {
359  {points * points, points}};
360  for (unsigned int c = 0; c < n_components; ++c)
361  for (unsigned int face = 0; face < 4; ++face)
362  if (process_face[0][face])
363  interpolate<2, 0, transpose>(face_offsets[face],
364  outer_strides[face / 2],
365  given_degree,
366  mask_weights[0],
367  mask_face[0][face],
368  weights,
369  values + c * n_dofs);
370  }
371  {
372  const std::array<unsigned int, 4> edge_offsets = {{p0, p2, p4, p6}};
373  for (unsigned int c = 0; c < n_components; ++c)
374  for (unsigned int edge = 0; edge < 4; ++edge)
375  if (process_edge[0][edge])
376  interpolate<1, 0, transpose>(edge_offsets[edge],
377  0,
378  given_degree,
379  mask_weights[0],
380  mask_edge[0][edge],
381  weights,
382  values + c * n_dofs);
383  }
384 
385  // direction 1:
386  if (given_degree > 1)
387  {
388  const std::array<unsigned int, 4> face_offsets = {
389  {p0, p1, p0, p4}};
390  const std::array<unsigned int, 2> outer_strides = {
391  {points * points, 1}};
392  for (unsigned int c = 0; c < n_components; ++c)
393  for (unsigned int face = 0; face < 4; ++face)
394  if (process_face[1][face])
395  interpolate<2, 1, transpose>(face_offsets[face],
396  outer_strides[face / 2],
397  given_degree,
398  mask_weights[1],
399  mask_face[1][face],
400  weights,
401  values + c * n_dofs);
402  }
403 
404  {
405  const std::array<unsigned int, 4> edge_offsets = {{p0, p1, p4, p5}};
406  for (unsigned int c = 0; c < n_components; ++c)
407  for (unsigned int edge = 0; edge < 4; ++edge)
408  if (process_edge[1][edge])
409  interpolate<1, 1, transpose>(edge_offsets[edge],
410  0,
411  given_degree,
412  mask_weights[1],
413  mask_edge[1][edge],
414  weights,
415  values + c * n_dofs);
416  }
417 
418  // direction 2:
419  if (given_degree > 1)
420  {
421  const std::array<unsigned int, 4> face_offsets = {
422  {p0, p1, p0, p2}};
423  const std::array<unsigned int, 2> outer_strides = {{points, 1}};
424  for (unsigned int c = 0; c < n_components; ++c)
425  for (unsigned int face = 0; face < 4; ++face)
426  if (process_face[2][face])
427  interpolate<2, 2, transpose>(face_offsets[face],
428  outer_strides[face / 2],
429  given_degree,
430  mask_weights[2],
431  mask_face[2][face],
432  weights,
433  values + c * n_dofs);
434  }
435 
436  {
437  const std::array<unsigned int, 4> edge_offsets = {{p0, p1, p2, p3}};
438  for (unsigned int c = 0; c < n_components; ++c)
439  for (unsigned int edge = 0; edge < 4; ++edge)
440  if (process_edge[2][edge])
441  interpolate<1, 2, transpose>(edge_offsets[edge],
442  0,
443  given_degree,
444  mask_weights[2],
445  mask_edge[2][edge],
446  weights,
447  values + c * n_dofs);
448  }
449  }
450  else
451  {
452  Assert(false, ExcNotImplemented());
453  }
454  }
455  };
456 
457  template <typename T1, VectorizationTypes VT>
458  struct Trait;
459 
460  template <typename T1>
462  {
463  using value_type = typename T1::value_type;
464  using index_type = unsigned int;
466 
467  template <typename T>
468  static inline const std::array<AlignedVector<interpolation_type>, 2> &
469  get_interpolation_matrix(const T &shape_info)
470  {
471  return shape_info.data.front().subface_interpolation_matrices_scalar;
472  }
473 
474  static inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
476  T1::size()> mask,
478  T1::size()> mask_new,
479  const unsigned int v)
480  {
481  (void)mask;
482  (void)mask_new;
483  return v;
484  }
485 
486  static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
487  do_break(unsigned int v,
489  {
490  (void)v;
491  (void)kind;
492  return false;
493  }
494 
495  static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
496  do_continue(unsigned int v,
498  {
499  (void)v;
500  return kind ==
502  }
503 
504  static inline DEAL_II_ALWAYS_INLINE_RELEASE
508  T1::size()> mask)
509  {
510  return mask;
511  }
512 
513  static inline DEAL_II_ALWAYS_INLINE_RELEASE typename T1::value_type
514  get_value(const typename T1::value_type &value, const index_type &i)
515  {
516  (void)i;
517  return value;
518  }
519 
520  static inline DEAL_II_ALWAYS_INLINE_RELEASE typename T1::value_type
521  get_value(const T1 &value, const index_type &i)
522  {
523  return value[i];
524  }
525 
526  static inline DEAL_II_ALWAYS_INLINE_RELEASE void
527  set_value(T1 & result,
528  const typename T1::value_type &value,
529  const index_type & i)
530  {
531  result[i] = value;
532  }
533  };
534 
535  template <typename T1>
537  {
538  using value_type = T1;
539  using index_type = std::pair<T1, T1>;
540  using interpolation_type = T1;
541 
542  template <typename T>
543  static inline const std::array<AlignedVector<T1>, 2> &
544  get_interpolation_matrix(const T &shape_info)
545  {
546  return shape_info.data.front().subface_interpolation_matrices;
547  }
548 
549  static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
550  do_break(unsigned int v,
552  {
553  (void)v;
554  (void)kind;
555  return false;
556  }
557 
558  static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
559  do_continue(unsigned int v,
561  {
562  (void)v;
563  return kind ==
565  }
566 
567  static inline DEAL_II_ALWAYS_INLINE_RELEASE index_type
569  T1::size()> mask,
571  T1::size()> mask_new,
572  const unsigned int v)
573  {
574  (void)mask;
575  (void)mask_new;
576  T1 result = 0.0;
577  result[v] = 1.0;
578  return {result, T1(1.0) - result};
579  }
580 
581  static inline DEAL_II_ALWAYS_INLINE_RELEASE
585  T1::size()> mask)
586  {
587  return mask;
588  }
589 
590  static inline DEAL_II_ALWAYS_INLINE_RELEASE T1
591  get_value(const T1 &value, const index_type &)
592  {
593  return value;
594  }
595 
596  static inline DEAL_II_ALWAYS_INLINE_RELEASE void
597  set_value(T1 &result, const T1 &value, const index_type &i)
598  {
599  result = result * i.second + value * i.first;
600  }
601  };
602 
603  template <typename T1>
605  {
606  using value_type = T1;
607  using index_type = std::pair<T1, T1>;
608  using interpolation_type = T1;
609 
610  template <typename T>
611  static inline const std::array<AlignedVector<T1>, 2> &
612  get_interpolation_matrix(const T &shape_info)
613  {
614  return shape_info.data.front().subface_interpolation_matrices;
615  }
616 
617  static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
618  do_break(unsigned int v,
620  {
621  (void)v;
622  return kind ==
624  }
625 
626  static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
627  do_continue(unsigned int v,
629  {
630  (void)v;
631  return kind ==
633  }
634 
635  static inline DEAL_II_ALWAYS_INLINE_RELEASE index_type
637  T1::size()> mask,
639  T1::size()> mask_new,
640  const unsigned int v)
641  {
642  T1 result;
643 
644  for (unsigned int i = 0; i < T1::size(); ++i)
645  result[i] = mask_new[v] == mask[i];
646 
647  return {result, T1(1.0) - result};
648  }
649 
650  static inline DEAL_II_ALWAYS_INLINE_RELEASE
654  T1::size()> mask)
655  {
656  auto new_mask = mask;
657 
658  std::sort(new_mask.begin(), new_mask.end());
659  std::fill(std::unique(new_mask.begin(), new_mask.end()),
660  new_mask.end(),
662 
663  return new_mask;
664  }
665 
666  static inline DEAL_II_ALWAYS_INLINE_RELEASE T1
667  get_value(const T1 &value, const index_type &)
668  {
669  return value;
670  }
671 
672  static inline DEAL_II_ALWAYS_INLINE_RELEASE void
673  set_value(T1 &result, const T1 &value, const index_type &i)
674  {
675  result = result * i.second + value * i.first;
676  }
677  };
678 
679  template <typename T1>
681  {
682  using value_type = T1;
683  using index_type = T1;
684  using interpolation_type = T1;
685 
686  template <typename T>
687  static inline const std::array<AlignedVector<T1>, 2> &
688  get_interpolation_matrix(const T &shape_info)
689  {
690  return shape_info.data.front().subface_interpolation_matrices;
691  }
692 
693  static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
694  do_break(unsigned int v,
696  {
697  (void)kind;
698  return v > 0;
699  }
700 
701  static inline DEAL_II_ALWAYS_INLINE_RELEASE bool
702  do_continue(unsigned int v,
704  {
705  (void)kind;
706 
707  Assert(false, ExcInternalError());
708 
709  return v > 0; // should not be called
710  }
711 
712  static inline DEAL_II_ALWAYS_INLINE_RELEASE T1
714  T1::size()> mask,
716  T1::size()> mask_new,
717  const unsigned int v)
718  {
719  (void)mask;
720  (void)mask_new;
721  (void)v;
722  return 1.0; // return something since not used
723  }
724 
725  static inline DEAL_II_ALWAYS_INLINE_RELEASE
729  T1::size()> mask)
730  {
731  return mask;
732  }
733 
734  static inline DEAL_II_ALWAYS_INLINE_RELEASE T1
735  get_value(const T1 &value, const index_type &)
736  {
737  return value;
738  }
739 
740  static inline DEAL_II_ALWAYS_INLINE_RELEASE void
741  set_value(T1 &result, const T1 &value, const index_type &i)
742  {
743  (void)i;
744  result = value;
745  }
746  };
747 
748 
749 
750  template <typename T,
751  typename Number,
752  VectorizationTypes VectorizationType,
753  int fe_degree,
754  bool transpose>
756  {
757  public:
760  const T & t,
761  const unsigned int & given_degree,
762  const bool & type_x,
763  const bool & type_y,
764  const bool & type_z,
766  const std::array<
770  Number *values)
771  : t(t)
773  , type_x(type_x)
774  , type_y(type_y)
775  , type_z(type_z)
776  , v(v)
778  , values(values)
779  {}
780 
781  template <unsigned int direction, unsigned int d, bool skip_borders>
782  static inline DEAL_II_ALWAYS_INLINE_RELEASE void
784  const unsigned int dof_offset,
785  const unsigned int given_degree,
788  *DEAL_II_RESTRICT weight,
789  Number *DEAL_II_RESTRICT values)
790  {
791  static constexpr unsigned int max_n_points_1D = 40;
792 
794  temp[fe_degree != -1 ? (fe_degree + 1) : max_n_points_1D];
795 
796  const unsigned int points =
797  (fe_degree != -1 ? fe_degree : given_degree) + 1;
798 
799  AssertIndexRange(given_degree, max_n_points_1D);
800 
801  const unsigned int stride = fe_degree != -1 ?
802  Utilities::pow(fe_degree + 1, direction) :
803  Utilities::pow(given_degree + 1, direction);
804 
805  // direction side0 side1 side2
806  // 0 - p^2 p
807  // 1 p^2 - 1
808  // 2 p - 1
809  const unsigned int stride2 =
810  ((direction == 0 && d == 1) || (direction == 1 && d == 0)) ?
811  (points * points) :
812  (((direction == 0 && d == 2) || (direction == 2 && d == 0)) ? points :
813  1);
814 
815  for (unsigned int g = (skip_borders ? 1 : 0);
816  g < points - (skip_borders ? 1 : 0);
817  ++g)
818  {
819  // copy result back
820  for (unsigned int k = 0; k < points; ++k)
822  values[dof_offset + k * stride + stride2 * g], v);
823 
824  // perform interpolation point by point
825  for (unsigned int k = 0; k < points; ++k)
826  {
828  weight[(transpose ? 1 : points) * k], v) *
829  temp[0];
830  for (unsigned int h = 1; h < points; ++h)
832  weight[(transpose ? 1 : points) * k +
833  (transpose ? points : 1) * h],
834  v) *
835  temp[h];
837  values[dof_offset + k * stride + stride2 * g], sum, v);
838  }
839  }
840  }
841 
842  template <unsigned int direction>
843  static inline DEAL_II_ALWAYS_INLINE_RELEASE void
845  const unsigned int p,
846  const unsigned int given_degree,
849  *DEAL_II_RESTRICT weight,
850  Number *DEAL_II_RESTRICT values)
851  {
852  static constexpr unsigned int max_n_points_1D = 40;
853 
855  temp[fe_degree != -1 ? (fe_degree + 1) : max_n_points_1D];
856 
857  const unsigned int points =
858  (fe_degree != -1 ? fe_degree : given_degree) + 1;
859 
860  AssertIndexRange(given_degree, max_n_points_1D);
861 
862  const unsigned int stride = fe_degree != -1 ?
863  Utilities::pow(fe_degree + 1, direction) :
864  Utilities::pow(given_degree + 1, direction);
865 
866  // copy result back
867  for (unsigned int k = 0; k < points; ++k)
868  temp[k] =
870  v);
871 
872  // perform interpolation point by point
873  for (unsigned int k = 0; k < points; ++k)
874  {
876  weight[(transpose ? 1 : points) * k], v) *
877  temp[0];
878  for (unsigned int h = 1; h < points; ++h)
880  weight[(transpose ? 1 : points) * k +
881  (transpose ? points : 1) * h],
882  v) *
883  temp[h];
885  sum,
886  v);
887  }
888  }
889 
890  template <bool do_x, bool do_y, bool do_z>
892  process_edge() const
893  {
894  if (do_x)
895  interpolate_3D_edge<0>(t.line(0, type_y, type_z),
896  given_degree,
897  v,
899  values);
900 
901  if (do_y)
902  interpolate_3D_edge<1>(t.line(1, type_x, type_z),
903  given_degree,
904  v,
906  values);
907 
908  if (do_z)
909  interpolate_3D_edge<2>(t.line(2, type_x, type_y),
910  given_degree,
911  v,
913  values);
914  }
915 
916  template <bool do_x, bool do_y, bool do_z>
919  {
920  static_assert((do_x && !do_y && !do_z) || (!do_x && do_y && !do_z) ||
921  (!do_x && !do_y && do_z),
922  "Only one face can be chosen.");
923 
924  static const unsigned int direction = do_x ? 0 : (do_y ? 1 : 2);
925  const bool type = do_x ? type_x : (do_y ? type_y : type_z);
926 
927  if (!do_x)
928  interpolate_3D_face<0, direction, false>(
929  t.face(direction, type),
930  given_degree,
931  v,
933  values);
934 
935  if (!do_y)
936  interpolate_3D_face<1, direction, false>(
937  t.face(direction, type),
938  given_degree,
939  v,
941  values);
942 
943  if (!do_z)
944  interpolate_3D_face<2, direction, false>(
945  t.face(direction, type),
946  given_degree,
947  v,
949  values);
950  }
951 
952  template <bool do_x, bool do_y, bool do_z>
955  {
956  static_assert(((do_x && !do_y && !do_z) || (!do_x && do_y && !do_z) ||
957  (!do_x && !do_y && do_z)) == false,
958  "Only one face can be chosen.");
959 
960  // direction 0
961  {
962  const auto inpterolation_matrix =
964 
965  // faces
966  if (do_y && given_degree > 1)
967  interpolate_3D_face<0, 1, true>(
968  t.face(1, type_y), given_degree, v, inpterolation_matrix, values);
969 
970  if (do_z && given_degree > 1)
971  interpolate_3D_face<0, 2, true>(
972  t.face(2, type_z), given_degree, v, inpterolation_matrix, values);
973 
974  // direction 0 -> edges
975  interpolate_3D_edge<0>((do_x && do_y && !do_z) ?
976  (t.lines_plane(0, type_x, type_y, 0)) :
977  ((do_x && !do_y && do_z) ?
978  (t.lines_plane(1, type_x, type_z, 0)) :
979  (t.lines(0, type_y, type_z, 0))),
980  given_degree,
981  v,
982  inpterolation_matrix,
983  values);
984 
985 
986  interpolate_3D_edge<0>((do_x && do_y && !do_z) ?
987  (t.lines_plane(0, type_x, type_y, 1)) :
988  ((do_x && !do_y && do_z) ?
989  (t.lines_plane(1, type_x, type_z, 1)) :
990  (t.lines(0, type_y, type_z, 1))),
991  given_degree,
992  v,
993  inpterolation_matrix,
994  values);
995 
996  if (do_y && do_z)
997  interpolate_3D_edge<0>(t.lines(0, type_y, type_z, 2),
998  given_degree,
999  v,
1000  inpterolation_matrix,
1001  values);
1002  }
1003 
1004  // direction 1
1005  {
1006  const auto inpterolation_matrix =
1007  interpolation_matrices[!type_y].data();
1008 
1009  // faces
1010  if (do_x && given_degree > 1)
1011  interpolate_3D_face<1, 0, true>(
1012  t.face(0, type_x), given_degree, v, inpterolation_matrix, values);
1013 
1014  if (do_z && given_degree > 1)
1015  interpolate_3D_face<1, 2, true>(
1016  t.face(2, type_z), given_degree, v, inpterolation_matrix, values);
1017 
1018  // lines
1019  interpolate_3D_edge<1>((do_x && do_y && !do_z) ?
1020  (t.lines_plane(0, type_x, type_y, 2)) :
1021  ((!do_x && do_y && do_z) ?
1022  (t.lines_plane(2, type_y, type_z, 0)) :
1023  (t.lines(1, type_x, type_z, 0))),
1024  given_degree,
1025  v,
1026  inpterolation_matrix,
1027  values);
1028 
1029  interpolate_3D_edge<1>((do_x && do_y && !do_z) ?
1030  (t.lines_plane(0, type_x, type_y, 3)) :
1031  ((!do_x && do_y && do_z) ?
1032  (t.lines_plane(2, type_y, type_z, 1)) :
1033  (t.lines(1, type_x, type_z, 1))),
1034  given_degree,
1035  v,
1036  inpterolation_matrix,
1037  values);
1038 
1039  if (do_x && do_z)
1040  interpolate_3D_edge<1>(t.lines(1, type_x, type_z, 2),
1041  given_degree,
1042  v,
1043  inpterolation_matrix,
1044  values);
1045  }
1046 
1047  // direction 2 -> faces
1048  {
1049  const auto inpterolation_matrix =
1050  interpolation_matrices[!type_z].data();
1051 
1052  if (do_x && given_degree > 1)
1053  interpolate_3D_face<2, 0, true>(
1054  t.face(0, type_x), given_degree, v, inpterolation_matrix, values);
1055 
1056  if (do_y && given_degree > 1)
1057  interpolate_3D_face<2, 1, true>(
1058  t.face(1, type_y), given_degree, v, inpterolation_matrix, values);
1059 
1060  // direction 2 -> edges
1061  interpolate_3D_edge<2>((do_x && !do_y && do_z) ?
1062  (t.lines_plane(1, type_x, type_z, 2)) :
1063  ((!do_x && do_y && do_z) ?
1064  (t.lines_plane(2, type_y, type_z, 2)) :
1065  (t.lines(2, type_x, type_y, 0))),
1066  given_degree,
1067  v,
1068  inpterolation_matrix,
1069  values);
1070 
1071  interpolate_3D_edge<2>((do_x && !do_y && do_z) ?
1072  (t.lines_plane(1, type_x, type_z, 3)) :
1073  ((!do_x && do_y && do_z) ?
1074  (t.lines_plane(2, type_y, type_z, 3)) :
1075  (t.lines(2, type_x, type_y, 1))),
1076  given_degree,
1077  v,
1078  inpterolation_matrix,
1079  values);
1080 
1081  if (do_x && do_y)
1082  interpolate_3D_edge<2>(t.lines(2, type_x, type_y, 2),
1083  given_degree,
1084  v,
1085  inpterolation_matrix,
1086  values);
1087  }
1088  }
1089 
1090  private:
1091  const T & t;
1092  const unsigned int & given_degree;
1093  const bool & type_x;
1094  const bool & type_y;
1095  const bool & type_z;
1097  const std::array<
1098  AlignedVector<
1101  Number *values;
1102  };
1103 
1107  enum class HelperType
1108  {
1113  constant,
1118  dynamic
1119  };
1120 
1121  template <HelperType helper_type,
1122  typename Number,
1123  VectorizationTypes VectorizationType,
1124  int fe_degree,
1125  bool transpose>
1126  class Helper;
1127 
1128  template <typename Number,
1129  VectorizationTypes VectorizationType,
1130  int fe_degree,
1131  bool transpose>
1133  Number,
1134  VectorizationType,
1135  fe_degree,
1136  transpose> : public HelperBase<Helper<HelperType::dynamic,
1137  Number,
1138  VectorizationType,
1139  fe_degree,
1140  transpose>,
1141  Number,
1142  VectorizationType,
1143  fe_degree,
1144  transpose>
1145  {
1146  public:
1148  Helper(const unsigned int &given_degree,
1149  const bool & type_x,
1150  const bool & type_y,
1151  const bool & type_z,
1153  const std::array<
1154  AlignedVector<
1156  2> & interpolation_matrices,
1157  Number *values)
1159  Number,
1160  VectorizationType,
1161  fe_degree,
1162  transpose>,
1163  Number,
1164  VectorizationType,
1165  fe_degree,
1166  transpose>(*this,
1167  given_degree,
1168  type_x,
1169  type_y,
1170  type_z,
1171  v,
1172  interpolation_matrices,
1173  values)
1174  , points(given_degree + 1)
1175  {
1176  static_assert(fe_degree == -1, "Only working for fe_degree = -1.");
1177  }
1178 
1179  const unsigned int points;
1180 
1181  inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1182  line(unsigned int i, unsigned int j, unsigned int k) const
1183  {
1184  return line_array[i][j][k];
1185  }
1186 
1187  inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1188  face(unsigned int i, unsigned int j) const
1189  {
1190  return face_array[i][j];
1191  }
1192 
1193  inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1194  lines_plane(unsigned int i,
1195  unsigned int j,
1196  unsigned int k,
1197  unsigned int l) const
1198  {
1199  return lines_plane_array[i][j][k][l];
1200  }
1201 
1202  inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1203  lines(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
1204  {
1205  return lines_array[i][j][k][l];
1206  }
1207 
1208  private:
1209  const ::ndarray<unsigned int, 3, 2, 2> line_array = {
1210  {{{{{points * points * points - points, points *points - points}},
1211  {{points * points * points - points * points, 0}}}},
1212  {{{{points * points * points - points * points + points - 1,
1213  points - 1}},
1214  {{points * points * points - points * points, 0}}}},
1215  {{{{points * points - 1, points - 1}},
1216  {{points * points - points, 0}}}}}};
1217 
1218  const ::ndarray<unsigned int, 3, 2> face_array = {
1219  {{{points - 1, 0}},
1220  {{points * points - points, 0}},
1221  {{points * points * points - points * points, 0}}}};
1222 
1223  const ::ndarray<unsigned int, 3, 2, 2, 4> lines_plane_array = {
1224  {{{{{{{points * points - points,
1225  points *points *points - points,
1226  points - 1,
1227  points *points *points - points *points + points - 1}},
1228  {{0,
1229  points *points *points - points *points,
1230  points - 1,
1231  points *points *points - points *points + points - 1}}}},
1232  {{{{points * points - points,
1233  points *points *points - points,
1234  0,
1235  points *points *points - points *points}},
1236  {{0,
1237  points *points *points - points *points,
1238  0,
1239  points *points *points - points *points}}}}}},
1240  {{{{{{points * points * points - points * points,
1241  points *points *points - points,
1242  points - 1,
1243  points *points - 1}},
1244  {{0, points *points - points, points - 1, points *points - 1}}}},
1245  {{{{points * points * points - points * points,
1246  points *points *points - points,
1247  0,
1248  points *points - points}},
1249  {{0, points *points - points, 0, points *points - points}}}}}},
1250  {{{{{{points * points * points - points * points,
1251  points *points *points - points *points + points - 1,
1252  points * points - points,
1253  points * points - 1}},
1254  {{0, points - 1, points *points - points, points *points - 1}}}},
1255  {{{{points * points * points - points * points,
1256  points *points *points - points *points + points - 1,
1257  0,
1258  points - 1}},
1259  {{0, points - 1, 0, points - 1}}}}}}}};
1260 
1261  const ::ndarray<unsigned int, 3, 2, 2, 3> lines_array = {
1262  {{{{{{{points * points - points,
1263  points *points *points - points *points,
1264  points *points *points - points}},
1265  {{0, points *points - points, points *points *points - points}}}},
1266  {{{{0,
1267  points *points *points - points *points,
1268  points *points *points - points}},
1269  {{0,
1270  points *points - points,
1271  points *points *points - points *points}}}}}},
1272  {{{{{{points - 1,
1273  points *points *points - points *points,
1274  points *points *points - points *points + points - 1}},
1275  {{0,
1276  points - 1,
1277  points *points *points - points *points + points - 1}}}},
1278  {{{{0,
1279  points *points *points - points *points,
1280  points *points *points - points *points + points - 1}},
1281  {{0, points - 1, points *points *points - points *points}}}}}},
1282  {{{{{{points - 1, points *points - points, points *points - 1}},
1283  {{0, points - 1, points *points - 1}}}},
1284  {{{{0, points *points - points, points *points - 1}},
1285  {{0, points - 1, points *points - points}}}}}}}};
1286  };
1287 
1288  template <typename Number,
1289  VectorizationTypes VectorizationType,
1290  int fe_degree,
1291  bool transpose>
1293  Number,
1294  VectorizationType,
1295  fe_degree,
1296  transpose> : public HelperBase<Helper<HelperType::constant,
1297  Number,
1298  VectorizationType,
1299  fe_degree,
1300  transpose>,
1301  Number,
1302  VectorizationType,
1303  fe_degree,
1304  transpose>
1305  {
1306  public:
1308  Helper(const unsigned int &given_degree,
1309  const bool & type_x,
1310  const bool & type_y,
1311  const bool & type_z,
1313  const std::array<
1314  AlignedVector<
1316  2> & interpolation_matrices,
1317  Number *values)
1319  Number,
1320  VectorizationType,
1321  fe_degree,
1322  transpose>,
1323  Number,
1324  VectorizationType,
1325  fe_degree,
1326  transpose>(*this,
1327  given_degree,
1328  type_x,
1329  type_y,
1330  type_z,
1331  v,
1332  interpolation_matrices,
1333  values)
1334  {
1335  static_assert(fe_degree != -1, "Only working for fe_degree != -1.");
1336  }
1337 
1338 
1339  inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1340  line(unsigned int i, unsigned int j, unsigned int k) const
1341  {
1342  static constexpr unsigned int points = fe_degree + 1;
1343 
1344  static constexpr ::ndarray<unsigned int, 3, 2, 2> line_array = {
1345  {{{{{points * points * points - points, points * points - points}},
1346  {{points * points * points - points * points, 0}}}},
1347  {{{{points * points * points - points * points + points - 1,
1348  points - 1}},
1349  {{points * points * points - points * points, 0}}}},
1350  {{{{points * points - 1, points - 1}},
1351  {{points * points - points, 0}}}}}};
1352 
1353  return line_array[i][j][k];
1354  }
1355 
1356  inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1357  face(unsigned int i, unsigned int j) const
1358  {
1359  static constexpr unsigned int points = fe_degree + 1;
1360 
1361  static constexpr ::ndarray<unsigned int, 3, 2> face_array = {
1362  {{{points - 1, 0}},
1363  {{points * points - points, 0}},
1364  {{points * points * points - points * points, 0}}}};
1365 
1366  return face_array[i][j];
1367  }
1368 
1369  inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1370  lines_plane(unsigned int i,
1371  unsigned int j,
1372  unsigned int k,
1373  unsigned int l) const
1374  {
1375  static constexpr unsigned int points = fe_degree + 1;
1376 
1377  static constexpr ::ndarray<unsigned int, 3, 2, 2, 4>
1378  lines_plane_array = {
1379  {{{{{{{points * points - points,
1380  points * points * points - points,
1381  points - 1,
1382  points * points * points - points * points + points - 1}},
1383  {{0,
1384  points * points * points - points * points,
1385  points - 1,
1386  points * points * points - points * points + points - 1}}}},
1387  {{{{points * points - points,
1388  points * points * points - points,
1389  0,
1390  points * points * points - points * points}},
1391  {{0,
1392  points * points * points - points * points,
1393  0,
1394  points * points * points - points * points}}}}}},
1395  {{{{{{points * points * points - points * points,
1396  points * points * points - points,
1397  points - 1,
1398  points * points - 1}},
1399  {{0,
1400  points * points - points,
1401  points - 1,
1402  points * points - 1}}}},
1403  {{{{points * points * points - points * points,
1404  points * points * points - points,
1405  0,
1406  points * points - points}},
1407  {{0, points * points - points, 0, points * points - points}}}}}},
1408  {{{{{{points * points * points - points * points,
1409  points * points * points - points * points + points - 1,
1410  points * points - points,
1411  points * points - 1}},
1412  {{0,
1413  points - 1,
1414  points * points - points,
1415  points * points - 1}}}},
1416  {{{{points * points * points - points * points,
1417  points * points * points - points * points + points - 1,
1418  0,
1419  points - 1}},
1420  {{0, points - 1, 0, points - 1}}}}}}}};
1421 
1422  return lines_plane_array[i][j][k][l];
1423  }
1424 
1425  inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int
1426  lines(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
1427  {
1428  static constexpr unsigned int points = fe_degree + 1;
1429 
1430  static constexpr ::ndarray<unsigned int, 3, 2, 2, 3> lines_array = {
1431  {{{{{{{points * points - points,
1432  points * points * points - points * points,
1433  points * points * points - points}},
1434  {{0,
1435  points * points - points,
1436  points * points * points - points}}}},
1437  {{{{0,
1438  points * points * points - points * points,
1439  points * points * points - points}},
1440  {{0,
1441  points * points - points,
1442  points * points * points - points * points}}}}}},
1443  {{{{{{points - 1,
1444  points * points * points - points * points,
1445  points * points * points - points * points + points - 1}},
1446  {{0,
1447  points - 1,
1448  points * points * points - points * points + points - 1}}}},
1449  {{{{0,
1450  points * points * points - points * points,
1451  points * points * points - points * points + points - 1}},
1452  {{0, points - 1, points * points * points - points * points}}}}}},
1453  {{{{{{points - 1, points * points - points, points * points - 1}},
1454  {{0, points - 1, points * points - 1}}}},
1455  {{{{0, points * points - points, points * points - 1}},
1456  {{0, points - 1, points * points - points}}}}}}}};
1457 
1458  return lines_array[i][j][k][l];
1459  }
1460  };
1461 
1462 
1463  template <int dim, int fe_degree, typename Number>
1466  dim,
1467  fe_degree,
1468  Number>
1469  {
1470  public:
1471  static const VectorizationTypes VectorizationType =
1473 
1474  private:
1475  template <unsigned int side, bool transpose>
1476  static inline DEAL_II_ALWAYS_INLINE_RELEASE void
1478  const unsigned int given_degree,
1481  *DEAL_II_RESTRICT weight,
1482  Number *DEAL_II_RESTRICT values)
1483  {
1484  static constexpr unsigned int max_n_points_1D = 40;
1485 
1487  temp[fe_degree != -1 ? (fe_degree + 1) : max_n_points_1D];
1488 
1489  const unsigned int points =
1490  (fe_degree != -1 ? fe_degree : given_degree) + 1;
1491 
1492  AssertIndexRange(given_degree, max_n_points_1D);
1493 
1494  const unsigned int d = side / 2; // direction
1495  const unsigned int s = side % 2; // left or right surface
1496 
1497  const unsigned int offset = ::Utilities::pow(points, d + 1);
1498  const unsigned int stride =
1499  (s == 0 ? 0 : (points - 1)) * ::Utilities::pow(points, d);
1500 
1501  const unsigned int r1 = ::Utilities::pow(points, dim - d - 1);
1502  const unsigned int r2 = ::Utilities::pow(points, d);
1503 
1504  // copy result back
1505  for (unsigned int i = 0, k = 0; i < r1; ++i)
1506  for (unsigned int j = 0; j < r2; ++j, ++k)
1508  values[i * offset + stride + j], v);
1509 
1510  // perform interpolation point by point (note: r1 * r2 ==
1511  // points^(dim-1))
1512  for (unsigned int i = 0, k = 0; i < r1; ++i)
1513  for (unsigned int j = 0; j < r2; ++j, ++k)
1514  {
1516  for (unsigned int h = 0; h < points; ++h)
1518  weight[(transpose ? 1 : points) * k +
1519  (transpose ? points : 1) * h],
1520  v) *
1521  temp[h];
1523  values[i * offset + stride + j], sum, v);
1524  }
1525  }
1526 
1527  public:
1528  template <bool transpose>
1529  static void
1531  const unsigned int n_desired_components,
1532  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1534  Number::size()> & constraint_mask,
1535  Number * values)
1536  {
1537  const unsigned int given_degree =
1538  fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
1539 
1540  const auto &interpolation_matrices =
1542 
1543  const auto constraint_mask_sorted =
1545 
1546  for (unsigned int c = 0; c < n_desired_components; ++c)
1547  {
1548  for (unsigned int v = 0; v < Number::size(); ++v)
1549  {
1550  const auto mask = constraint_mask_sorted[v];
1551 
1553  break;
1554 
1556  continue;
1557 
1558  const auto vv =
1560  constraint_mask_sorted,
1561  v);
1562 
1563  if (dim == 2) // 2D: only faces
1564  {
1565  const bool subcell_x = (mask >> 0) & 1;
1566  const bool subcell_y = (mask >> 1) & 1;
1567  const bool face_x = (mask >> 3) & 1;
1568  const bool face_y = (mask >> 4) & 1;
1569 
1570  // direction 0:
1571  if (face_y)
1572  {
1573  const auto *weights =
1574  interpolation_matrices[!subcell_x].data();
1575 
1576  if (subcell_y)
1577  interpolate_2D<2, transpose>(given_degree,
1578  vv,
1579  weights,
1580  values); // face 2
1581  else
1582  interpolate_2D<3, transpose>(given_degree,
1583  vv,
1584  weights,
1585  values); // face 3
1586  }
1587 
1588  // direction 1:
1589  if (face_x)
1590  {
1591  const auto *weights =
1592  interpolation_matrices[!subcell_y].data();
1593 
1594  if (subcell_x)
1595  interpolate_2D<0, transpose>(given_degree,
1596  vv,
1597  weights,
1598  values); // face 0
1599  else
1600  interpolate_2D<1, transpose>(given_degree,
1601  vv,
1602  weights,
1603  values); // face 1
1604  }
1605  }
1606  else if (dim == 3) // 3D faces and edges
1607  {
1608  const bool type_x = (mask >> 0) & 1;
1609  const bool type_y = (mask >> 1) & 1;
1610  const bool type_z = (mask >> 2) & 1;
1611 
1612  const auto flag_0 = (mask >> 3) & 3;
1613  const auto flag_1 = (mask >> 5) & 7;
1614  const auto faces = (flag_0 & 0b01) ? flag_1 : 0;
1615  const auto edges = (flag_0 & 0b10) ? flag_1 : 0;
1616 
1617  Helper<fe_degree == -1 ? HelperType::dynamic :
1619  Number,
1620  VectorizationType,
1621  fe_degree,
1622  transpose>
1623  helper(given_degree,
1624  type_x,
1625  type_y,
1626  type_z,
1627  vv,
1628  interpolation_matrices,
1629  values);
1630 
1631  if (faces > 0)
1632  switch (faces)
1633  {
1634  case 0:
1635  break;
1636  case 1:
1637  helper
1638  .template process_faces_fast<true, false, false>();
1639  break;
1640  case 2:
1641  helper
1642  .template process_faces_fast<false, true, false>();
1643  break;
1644  case 3:
1645  helper.template process_faces<true, true, false>();
1646  break;
1647  case 4:
1648  helper
1649  .template process_faces_fast<false, false, true>();
1650  break;
1651  case 5:
1652  helper.template process_faces<true, false, true>();
1653  break;
1654  case 6:
1655  helper.template process_faces<false, true, true>();
1656  break;
1657  case 7:
1658  helper.template process_faces<true, true, true>();
1659  break;
1660  }
1661 
1662  if (edges > 0)
1663  switch (edges)
1664  {
1665  case 0:
1666  break;
1667  case 1:
1668  helper.template process_edge<true, false, false>();
1669  break;
1670  case 2:
1671  helper.template process_edge<false, true, false>();
1672  break;
1673  case 3:
1674  helper.template process_edge<true, true, false>();
1675  break;
1676  case 4:
1677  helper.template process_edge<false, false, true>();
1678  break;
1679  case 5:
1680  helper.template process_edge<true, false, true>();
1681  break;
1682  case 6:
1683  helper.template process_edge<false, true, true>();
1684  break;
1685  case 7:
1686  helper.template process_edge<true, true, true>();
1687  break;
1688  }
1689  }
1690  else
1691  {
1692  Assert(false, ExcNotImplemented());
1693  }
1694  }
1695 
1696  values += shape_info.dofs_per_component_on_cell;
1697  }
1698  }
1699  };
1700 
1701 
1702 
1703  template <int dim, typename Number>
1705  {
1706  public:
1707  template <int fe_degree, int n_q_points_1d>
1708  static bool
1709  run(const unsigned int n_desired_components,
1710  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1711  const bool transpose,
1713  Number::size()> & c_mask,
1714  Number * values)
1715  {
1716  using RunnerType =
1718  dim,
1719  fe_degree,
1720  Number>;
1721 
1722  if (transpose)
1723  RunnerType::template run_internal<true>(n_desired_components,
1724  shape_info,
1725  c_mask,
1726  values);
1727  else
1728  RunnerType::template run_internal<false>(n_desired_components,
1729  shape_info,
1730  c_mask,
1731  values);
1732 
1733  return false;
1734  }
1735 
1736  template <int fe_degree>
1739  {
1740  return ((Number::size() > 2) && (fe_degree == -1 || fe_degree > 2)) ?
1743  }
1744  };
1745 
1746 
1747 } // end of namespace internal
1748 
1749 #undef DEAL_II_ALWAYS_INLINE_RELEASE
1750 
1751 
1753 
1754 #endif
static DEAL_II_ALWAYS_INLINE_RELEASE void interpolate_2D(const unsigned int given_degree, const typename Trait< Number, VectorizationType >::index_type v, const typename Trait< Number, VectorizationType >::interpolation_type *DEAL_II_RESTRICT weight, Number *DEAL_II_RESTRICT values)
static void run_internal(const unsigned int n_desired_components, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const std::array< MatrixFreeFunctions::compressed_constraint_kind, Number::size()> &constraint_mask, Number *values)
static void interpolate(const unsigned int offset, const unsigned int outer_stride, const unsigned int given_degree, const Number mask_weight, const Number mask_write, const Number *DEAL_II_RESTRICT weights, Number *DEAL_II_RESTRICT values)
static void run_internal(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const std::array< MatrixFreeFunctions::compressed_constraint_kind, Number::size()> &constraint_mask, Number *values)
DEAL_II_ALWAYS_INLINE_RELEASE void process_faces() const
const std::array< AlignedVector< typename Trait< Number, VectorizationType >::interpolation_type >, 2 > & interpolation_matrices
const Trait< Number, VectorizationType >::index_type & v
DEAL_II_ALWAYS_INLINE_RELEASE HelperBase(const T &t, const unsigned int &given_degree, const bool &type_x, const bool &type_y, const bool &type_z, const typename Trait< Number, VectorizationType >::index_type &v, const std::array< AlignedVector< typename Trait< Number, VectorizationType >::interpolation_type >, 2 > &interpolation_matrices, Number *values)
static DEAL_II_ALWAYS_INLINE_RELEASE void interpolate_3D_edge(const unsigned int p, const unsigned int given_degree, const typename Trait< Number, VectorizationType >::index_type v, const typename Trait< Number, VectorizationType >::interpolation_type *DEAL_II_RESTRICT weight, Number *DEAL_II_RESTRICT values)
DEAL_II_ALWAYS_INLINE_RELEASE void process_edge() const
static DEAL_II_ALWAYS_INLINE_RELEASE void interpolate_3D_face(const unsigned int dof_offset, const unsigned int given_degree, const typename Trait< Number, VectorizationType >::index_type v, const typename Trait< Number, VectorizationType >::interpolation_type *DEAL_II_RESTRICT weight, Number *DEAL_II_RESTRICT values)
DEAL_II_ALWAYS_INLINE_RELEASE void process_faces_fast() const
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int lines_plane(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
DEAL_II_ALWAYS_INLINE_RELEASE Helper(const unsigned int &given_degree, const bool &type_x, const bool &type_y, const bool &type_z, const typename Trait< Number, VectorizationType >::index_type &v, const std::array< AlignedVector< typename Trait< Number, VectorizationType >::interpolation_type >, 2 > &interpolation_matrices, Number *values)
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int lines(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int face(unsigned int i, unsigned int j) const
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int line(unsigned int i, unsigned int j, unsigned int k) const
DEAL_II_ALWAYS_INLINE_RELEASE Helper(const unsigned int &given_degree, const bool &type_x, const bool &type_y, const bool &type_z, const typename Trait< Number, VectorizationType >::index_type &v, const std::array< AlignedVector< typename Trait< Number, VectorizationType >::interpolation_type >, 2 > &interpolation_matrices, Number *values)
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int lines_plane(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int face(unsigned int i, unsigned int j) const
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int line(unsigned int i, unsigned int j, unsigned int k) const
DEAL_II_ALWAYS_INLINE_RELEASE unsigned int lines(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_RESTRICT
Definition: config.h:103
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
#define DEAL_II_ALWAYS_INLINE_RELEASE
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
void get_interpolation_matrix(const FiniteElement< dim, spacedim > &fe1, const FiniteElement< dim, spacedim > &fe2, FullMatrix< number > &interpolation_matrix)
static const char T
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm &mpi_communicator)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:462
std::uint8_t compressed_constraint_kind
Definition: dof_info.h:83
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108
static constexpr FEEvaluationImplHangingNodesRunnerTypes used_runner_type()
static bool run(const unsigned int n_desired_components, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const bool transpose, const std::array< MatrixFreeFunctions::compressed_constraint_kind, Number::size()> &c_mask, Number *values)
std::vector< UnivariateShapeData< Number > > data
Definition: shape_info.h:423
static const std::array< AlignedVector< T1 >, 2 > & get_interpolation_matrix(const T &shape_info)
static DEAL_II_ALWAYS_INLINE_RELEASE index_type create(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask, const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask_new, const unsigned int v)
static DEAL_II_ALWAYS_INLINE_RELEASE T1 get_value(const T1 &value, const index_type &)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_continue(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)
static DEAL_II_ALWAYS_INLINE_RELEASE void set_value(T1 &result, const T1 &value, const index_type &i)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_break(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)
static DEAL_II_ALWAYS_INLINE_RELEASE std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> create_mask(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_break(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)
static DEAL_II_ALWAYS_INLINE_RELEASE T1::value_type get_value(const T1 &value, const index_type &i)
static DEAL_II_ALWAYS_INLINE_RELEASE void set_value(T1 &result, const typename T1::value_type &value, const index_type &i)
static DEAL_II_ALWAYS_INLINE_RELEASE unsigned int create(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask, const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask_new, const unsigned int v)
static const std::array< AlignedVector< interpolation_type >, 2 > & get_interpolation_matrix(const T &shape_info)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_continue(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)
static DEAL_II_ALWAYS_INLINE_RELEASE T1::value_type get_value(const typename T1::value_type &value, const index_type &i)
static DEAL_II_ALWAYS_INLINE_RELEASE std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> create_mask(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask)
static const std::array< AlignedVector< T1 >, 2 > & get_interpolation_matrix(const T &shape_info)
static DEAL_II_ALWAYS_INLINE_RELEASE index_type create(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask, const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask_new, const unsigned int v)
static DEAL_II_ALWAYS_INLINE_RELEASE std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> create_mask(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask)
static DEAL_II_ALWAYS_INLINE_RELEASE T1 get_value(const T1 &value, const index_type &)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_continue(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)
static DEAL_II_ALWAYS_INLINE_RELEASE void set_value(T1 &result, const T1 &value, const index_type &i)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_break(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)
static DEAL_II_ALWAYS_INLINE_RELEASE T1 create(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask, const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask_new, const unsigned int v)
static DEAL_II_ALWAYS_INLINE_RELEASE std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> create_mask(const std::array< MatrixFreeFunctions::compressed_constraint_kind, T1::size()> mask)
static DEAL_II_ALWAYS_INLINE_RELEASE T1 get_value(const T1 &value, const index_type &)
static DEAL_II_ALWAYS_INLINE_RELEASE void set_value(T1 &result, const T1 &value, const index_type &i)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_break(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)
static const std::array< AlignedVector< T1 >, 2 > & get_interpolation_matrix(const T &shape_info)
static DEAL_II_ALWAYS_INLINE_RELEASE bool do_continue(unsigned int v, const MatrixFreeFunctions::compressed_constraint_kind &kind)