Reference documentation for deal.II version 9.3.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
evaluation_kernels.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2021 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_h
18 #define dealii_matrix_free_evaluation_kernels_h
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/utilities.h>
24 
30 
31 
33 
34 
35 // forward declaration
36 template <int, typename, bool, typename>
38 
39 
40 
41 namespace internal
42 {
43  // Select evaluator type from element shape function type
44  template <MatrixFreeFunctions::ElementType element, bool is_long>
46  {};
47 
48  template <bool is_long>
49  struct EvaluatorSelector<MatrixFreeFunctions::tensor_general, is_long>
50  {
51  static const EvaluatorVariant variant = evaluate_general;
52  };
53 
54  template <>
55  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, false>
56  {
57  static const EvaluatorVariant variant = evaluate_symmetric;
58  };
59 
60  template <>
61  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, true>
62  {
63  static const EvaluatorVariant variant = evaluate_evenodd;
64  };
65 
66  template <bool is_long>
67  struct EvaluatorSelector<MatrixFreeFunctions::truncated_tensor, is_long>
68  {
69  static const EvaluatorVariant variant = evaluate_general;
70  };
71 
72  template <>
73  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,
74  false>
75  {
76  static const EvaluatorVariant variant = evaluate_general;
77  };
78 
79  template <>
80  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0, true>
81  {
82  static const EvaluatorVariant variant = evaluate_evenodd;
83  };
84 
85  template <bool is_long>
86  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_collocation,
87  is_long>
88  {
89  static const EvaluatorVariant variant = evaluate_evenodd;
90  };
91 
92 
93 
110  template <MatrixFreeFunctions::ElementType type,
111  int dim,
112  int fe_degree,
113  int n_q_points_1d,
114  typename Number>
116  {
117  static void
118  evaluate(const unsigned int n_components,
119  const EvaluationFlags::EvaluationFlags evaluation_flag,
120  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
121  const Number * values_dofs_actual,
122  Number * values_quad,
123  Number * gradients_quad,
124  Number * hessians_quad,
125  Number * scratch_data);
126 
127  static void
128  integrate(const unsigned int n_components,
129  const EvaluationFlags::EvaluationFlags integration_flag,
130  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
131  Number * values_dofs_actual,
132  Number * values_quad,
133  Number * gradients_quad,
134  Number * scratch_data,
135  const bool add_into_values_array);
136  };
137 
138 
139 
144  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
145  struct FEEvaluationImpl<MatrixFreeFunctions::tensor_none,
146  dim,
147  fe_degree,
148  n_q_points_1d,
149  Number>
150  {
151  static void
152  evaluate(const unsigned int n_components,
153  const EvaluationFlags::EvaluationFlags evaluation_flag,
154  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
155  const Number * values_dofs_actual,
156  Number * values_quad,
157  Number * gradients_quad,
158  Number * hessians_quad,
159  Number * scratch_data);
160 
161  static void
162  integrate(const unsigned int n_components,
163  const EvaluationFlags::EvaluationFlags integration_flag,
164  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
165  Number * values_dofs_actual,
166  Number * values_quad,
167  Number * gradients_quad,
168  Number * scratch_data,
169  const bool add_into_values_array);
170  };
171 
172 
173 
174  template <MatrixFreeFunctions::ElementType type,
175  int dim,
176  int fe_degree,
177  int n_q_points_1d,
178  typename Number>
179  inline void
181  const unsigned int n_components,
182  const EvaluationFlags::EvaluationFlags evaluation_flag,
183  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
184  const Number * values_dofs_actual,
185  Number * values_quad,
186  Number * gradients_quad,
187  Number * hessians_quad,
188  Number * scratch_data)
189  {
190  if (evaluation_flag == EvaluationFlags::nothing)
191  return;
192 
193  const EvaluatorVariant variant =
195  using Eval = EvaluatorTensorProduct<variant,
196  dim,
197  fe_degree + 1,
198  n_q_points_1d,
199  Number>;
200  Eval eval(variant == evaluate_evenodd ?
201  shape_info.data.front().shape_values_eo :
202  shape_info.data.front().shape_values,
203  variant == evaluate_evenodd ?
204  shape_info.data.front().shape_gradients_eo :
205  shape_info.data.front().shape_gradients,
206  variant == evaluate_evenodd ?
207  shape_info.data.front().shape_hessians_eo :
208  shape_info.data.front().shape_hessians,
209  shape_info.data.front().fe_degree + 1,
210  shape_info.data.front().n_q_points_1d);
211 
212  const unsigned int temp_size =
213  Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
214  0 :
215  (Eval::n_rows_of_product > Eval::n_columns_of_product ?
216  Eval::n_rows_of_product :
217  Eval::n_columns_of_product);
218  Number *temp1 = scratch_data;
219  Number *temp2;
220  if (temp_size == 0)
221  {
222  temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
223  shape_info.data.front().fe_degree + 1),
224  Utilities::fixed_power<dim>(
225  shape_info.data.front().n_q_points_1d));
226  }
227  else
228  {
229  temp2 = temp1 + temp_size;
230  }
231 
232  const unsigned int n_q_points =
233  temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
234  const unsigned int dofs_per_comp =
236  Utilities::fixed_power<dim>(shape_info.data.front().fe_degree + 1) :
237  shape_info.dofs_per_component_on_cell;
238  const Number *values_dofs = values_dofs_actual;
240  {
241  Number *values_dofs_tmp =
242  scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
243  shape_info.n_q_points));
244  const int degree =
245  fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
246  for (unsigned int c = 0; c < n_components; ++c)
247  for (int i = 0, count_p = 0, count_q = 0;
248  i < (dim > 2 ? degree + 1 : 1);
249  ++i)
250  {
251  for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
252  {
253  for (int k = 0; k < degree + 1 - j - i;
254  ++k, ++count_p, ++count_q)
255  values_dofs_tmp[c * dofs_per_comp + count_q] =
256  values_dofs_actual
257  [c * shape_info.dofs_per_component_on_cell + count_p];
258  for (int k = degree + 1 - j - i; k < degree + 1;
259  ++k, ++count_q)
260  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
261  }
262  for (int j = degree + 1 - i; j < degree + 1; ++j)
263  for (int k = 0; k < degree + 1; ++k, ++count_q)
264  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
265  }
266  values_dofs = values_dofs_tmp;
267  }
268 
269  switch (dim)
270  {
271  case 1:
272  for (unsigned int c = 0; c < n_components; c++)
273  {
274  if (evaluation_flag & EvaluationFlags::values)
275  eval.template values<0, true, false>(values_dofs, values_quad);
276  if (evaluation_flag & EvaluationFlags::gradients)
277  eval.template gradients<0, true, false>(values_dofs,
278  gradients_quad);
279  if (evaluation_flag & EvaluationFlags::hessians)
280  eval.template hessians<0, true, false>(values_dofs,
281  hessians_quad);
282 
283  // advance the next component in 1D array
284  values_dofs += dofs_per_comp;
285  values_quad += n_q_points;
286  gradients_quad += n_q_points;
287  hessians_quad += n_q_points;
288  }
289  break;
290 
291  case 2:
292  for (unsigned int c = 0; c < n_components; c++)
293  {
294  // grad x
295  if (evaluation_flag & EvaluationFlags::gradients)
296  {
297  eval.template gradients<0, true, false>(values_dofs, temp1);
298  eval.template values<1, true, false>(temp1, gradients_quad);
299  }
300  if (evaluation_flag & EvaluationFlags::hessians)
301  {
302  // grad xy
303  if (!(evaluation_flag & EvaluationFlags::gradients))
304  eval.template gradients<0, true, false>(values_dofs, temp1);
305  eval.template gradients<1, true, false>(temp1,
306  hessians_quad +
307  2 * n_q_points);
308 
309  // grad xx
310  eval.template hessians<0, true, false>(values_dofs, temp1);
311  eval.template values<1, true, false>(temp1, hessians_quad);
312  }
313 
314  // grad y
315  eval.template values<0, true, false>(values_dofs, temp1);
316  if (evaluation_flag & EvaluationFlags::gradients)
317  eval.template gradients<1, true, false>(temp1,
318  gradients_quad +
319  n_q_points);
320 
321  // grad yy
322  if (evaluation_flag & EvaluationFlags::hessians)
323  eval.template hessians<1, true, false>(temp1,
324  hessians_quad +
325  n_q_points);
326 
327  // val: can use values applied in x
328  if (evaluation_flag & EvaluationFlags::values)
329  eval.template values<1, true, false>(temp1, values_quad);
330 
331  // advance to the next component in 1D array
332  values_dofs += dofs_per_comp;
333  values_quad += n_q_points;
334  gradients_quad += 2 * n_q_points;
335  hessians_quad += 3 * n_q_points;
336  }
337  break;
338 
339  case 3:
340  for (unsigned int c = 0; c < n_components; c++)
341  {
342  if (evaluation_flag & EvaluationFlags::gradients)
343  {
344  // grad x
345  eval.template gradients<0, true, false>(values_dofs, temp1);
346  eval.template values<1, true, false>(temp1, temp2);
347  eval.template values<2, true, false>(temp2, gradients_quad);
348  }
349 
350  if (evaluation_flag & EvaluationFlags::hessians)
351  {
352  // grad xz
353  if (!(evaluation_flag & EvaluationFlags::gradients))
354  {
355  eval.template gradients<0, true, false>(values_dofs,
356  temp1);
357  eval.template values<1, true, false>(temp1, temp2);
358  }
359  eval.template gradients<2, true, false>(temp2,
360  hessians_quad +
361  4 * n_q_points);
362 
363  // grad xy
364  eval.template gradients<1, true, false>(temp1, temp2);
365  eval.template values<2, true, false>(temp2,
366  hessians_quad +
367  3 * n_q_points);
368 
369  // grad xx
370  eval.template hessians<0, true, false>(values_dofs, temp1);
371  eval.template values<1, true, false>(temp1, temp2);
372  eval.template values<2, true, false>(temp2, hessians_quad);
373  }
374 
375  // grad y
376  eval.template values<0, true, false>(values_dofs, temp1);
377  if (evaluation_flag & EvaluationFlags::gradients)
378  {
379  eval.template gradients<1, true, false>(temp1, temp2);
380  eval.template values<2, true, false>(temp2,
381  gradients_quad +
382  n_q_points);
383  }
384 
385  if (evaluation_flag & EvaluationFlags::hessians)
386  {
387  // grad yz
388  if (!(evaluation_flag & EvaluationFlags::gradients))
389  eval.template gradients<1, true, false>(temp1, temp2);
390  eval.template gradients<2, true, false>(temp2,
391  hessians_quad +
392  5 * n_q_points);
393 
394  // grad yy
395  eval.template hessians<1, true, false>(temp1, temp2);
396  eval.template values<2, true, false>(temp2,
397  hessians_quad +
398  n_q_points);
399  }
400 
401  // grad z: can use the values applied in x direction stored in
402  // temp1
403  eval.template values<1, true, false>(temp1, temp2);
404  if (evaluation_flag & EvaluationFlags::gradients)
405  eval.template gradients<2, true, false>(temp2,
406  gradients_quad +
407  2 * n_q_points);
408 
409  // grad zz: can use the values applied in x and y direction stored
410  // in temp2
411  if (evaluation_flag & EvaluationFlags::hessians)
412  eval.template hessians<2, true, false>(temp2,
413  hessians_quad +
414  2 * n_q_points);
415 
416  // val: can use the values applied in x & y direction stored in
417  // temp2
418  if (evaluation_flag & EvaluationFlags::values)
419  eval.template values<2, true, false>(temp2, values_quad);
420 
421  // advance to the next component in 1D array
422  values_dofs += dofs_per_comp;
423  values_quad += n_q_points;
424  gradients_quad += 3 * n_q_points;
425  hessians_quad += 6 * n_q_points;
426  }
427  break;
428 
429  default:
430  AssertThrow(false, ExcNotImplemented());
431  }
432 
433  // case additional dof for FE_Q_DG0: add values; gradients and second
434  // derivatives evaluate to zero
436  (evaluation_flag & EvaluationFlags::values))
437  {
438  values_quad -= n_components * n_q_points;
439  values_dofs -= n_components * dofs_per_comp;
440  for (unsigned int c = 0; c < n_components; ++c)
441  for (unsigned int q = 0; q < shape_info.n_q_points; ++q)
442  values_quad[c * shape_info.n_q_points + q] +=
443  values_dofs[(c + 1) * shape_info.dofs_per_component_on_cell - 1];
444  }
445  }
446 
447 
448 
449  template <MatrixFreeFunctions::ElementType type,
450  int dim,
451  int fe_degree,
452  int n_q_points_1d,
453  typename Number>
454  inline void
456  const unsigned int n_components,
457  const EvaluationFlags::EvaluationFlags integration_flag,
458  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
459  Number * values_dofs_actual,
460  Number * values_quad,
461  Number * gradients_quad,
462  Number * scratch_data,
463  const bool add_into_values_array)
464  {
465  const EvaluatorVariant variant =
467  using Eval = EvaluatorTensorProduct<variant,
468  dim,
469  fe_degree + 1,
470  n_q_points_1d,
471  Number>;
472  Eval eval(variant == evaluate_evenodd ?
473  shape_info.data.front().shape_values_eo :
474  shape_info.data.front().shape_values,
475  variant == evaluate_evenodd ?
476  shape_info.data.front().shape_gradients_eo :
477  shape_info.data.front().shape_gradients,
478  variant == evaluate_evenodd ?
479  shape_info.data.front().shape_hessians_eo :
480  shape_info.data.front().shape_hessians,
481  shape_info.data.front().fe_degree + 1,
482  shape_info.data.front().n_q_points_1d);
483 
484  const unsigned int temp_size =
485  Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
486  0 :
487  (Eval::n_rows_of_product > Eval::n_columns_of_product ?
488  Eval::n_rows_of_product :
489  Eval::n_columns_of_product);
490  Number *temp1 = scratch_data;
491  Number *temp2;
492  if (temp_size == 0)
493  {
494  temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
495  shape_info.data.front().fe_degree + 1),
496  Utilities::fixed_power<dim>(
497  shape_info.data.front().n_q_points_1d));
498  }
499  else
500  {
501  temp2 = temp1 + temp_size;
502  }
503 
504  const unsigned int n_q_points =
505  temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
506  const unsigned int dofs_per_comp =
508  Utilities::fixed_power<dim>(shape_info.data.front().fe_degree + 1) :
509  shape_info.dofs_per_component_on_cell;
510  // expand dof_values to tensor product for truncated tensor products
511  Number *values_dofs =
513  scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
514  shape_info.n_q_points)) :
515  values_dofs_actual;
516 
517  switch (dim)
518  {
519  case 1:
520  for (unsigned int c = 0; c < n_components; c++)
521  {
522  if (integration_flag & EvaluationFlags::values)
523  {
524  if (add_into_values_array == false)
525  eval.template values<0, false, false>(values_quad,
526  values_dofs);
527  else
528  eval.template values<0, false, true>(values_quad,
529  values_dofs);
530  }
531  if (integration_flag & EvaluationFlags::gradients)
532  {
533  if (integration_flag & EvaluationFlags::values ||
534  add_into_values_array == true)
535  eval.template gradients<0, false, true>(gradients_quad,
536  values_dofs);
537  else
538  eval.template gradients<0, false, false>(gradients_quad,
539  values_dofs);
540  }
541 
542  // advance to the next component in 1D array
543  values_dofs += dofs_per_comp;
544  values_quad += n_q_points;
545  gradients_quad += n_q_points;
546  }
547  break;
548 
549  case 2:
550  for (unsigned int c = 0; c < n_components; c++)
551  {
552  if ((integration_flag & EvaluationFlags::values) &&
553  !(integration_flag & EvaluationFlags::gradients))
554  {
555  eval.template values<1, false, false>(values_quad, temp1);
556  if (add_into_values_array == false)
557  eval.template values<0, false, false>(temp1, values_dofs);
558  else
559  eval.template values<0, false, true>(temp1, values_dofs);
560  }
561  if (integration_flag & EvaluationFlags::gradients)
562  {
563  eval.template gradients<1, false, false>(gradients_quad +
564  n_q_points,
565  temp1);
566  if (integration_flag & EvaluationFlags::values)
567  eval.template values<1, false, true>(values_quad, temp1);
568  if (add_into_values_array == false)
569  eval.template values<0, false, false>(temp1, values_dofs);
570  else
571  eval.template values<0, false, true>(temp1, values_dofs);
572  eval.template values<1, false, false>(gradients_quad, temp1);
573  eval.template gradients<0, false, true>(temp1, values_dofs);
574  }
575 
576  // advance to the next component in 1D array
577  values_dofs += dofs_per_comp;
578  values_quad += n_q_points;
579  gradients_quad += 2 * n_q_points;
580  }
581  break;
582 
583  case 3:
584  for (unsigned int c = 0; c < n_components; c++)
585  {
586  if ((integration_flag & EvaluationFlags::values) &&
587  !(integration_flag & EvaluationFlags::gradients))
588  {
589  eval.template values<2, false, false>(values_quad, temp1);
590  eval.template values<1, false, false>(temp1, temp2);
591  if (add_into_values_array == false)
592  eval.template values<0, false, false>(temp2, values_dofs);
593  else
594  eval.template values<0, false, true>(temp2, values_dofs);
595  }
596  if (integration_flag & EvaluationFlags::gradients)
597  {
598  eval.template gradients<2, false, false>(gradients_quad +
599  2 * n_q_points,
600  temp1);
601  if (integration_flag & EvaluationFlags::values)
602  eval.template values<2, false, true>(values_quad, temp1);
603  eval.template values<1, false, false>(temp1, temp2);
604  eval.template values<2, false, false>(gradients_quad +
605  n_q_points,
606  temp1);
607  eval.template gradients<1, false, true>(temp1, temp2);
608  if (add_into_values_array == false)
609  eval.template values<0, false, false>(temp2, values_dofs);
610  else
611  eval.template values<0, false, true>(temp2, values_dofs);
612  eval.template values<2, false, false>(gradients_quad, temp1);
613  eval.template values<1, false, false>(temp1, temp2);
614  eval.template gradients<0, false, true>(temp2, values_dofs);
615  }
616 
617  // advance to the next component in 1D array
618  values_dofs += dofs_per_comp;
619  values_quad += n_q_points;
620  gradients_quad += 3 * n_q_points;
621  }
622  break;
623 
624  default:
625  AssertThrow(false, ExcNotImplemented());
626  }
627 
628  // case FE_Q_DG0: add values, gradients and second derivatives are zero
630  {
631  values_dofs -= n_components * dofs_per_comp -
632  shape_info.dofs_per_component_on_cell + 1;
633  values_quad -= n_components * n_q_points;
634  if (integration_flag & EvaluationFlags::values)
635  for (unsigned int c = 0; c < n_components; ++c)
636  {
637  values_dofs[0] = values_quad[0];
638  for (unsigned int q = 1; q < shape_info.n_q_points; ++q)
639  values_dofs[0] += values_quad[q];
640  values_dofs += dofs_per_comp;
641  values_quad += n_q_points;
642  }
643  else
644  {
645  for (unsigned int c = 0; c < n_components; ++c)
646  values_dofs[c * shape_info.dofs_per_component_on_cell] = Number();
647  values_dofs += n_components * shape_info.dofs_per_component_on_cell;
648  }
649  }
650 
652  {
653  values_dofs -= dofs_per_comp * n_components;
654  const int degree =
655  fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
656  for (unsigned int c = 0; c < n_components; ++c)
657  for (int i = 0, count_p = 0, count_q = 0;
658  i < (dim > 2 ? degree + 1 : 1);
659  ++i)
660  {
661  for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
662  {
663  for (int k = 0; k < degree + 1 - j - i;
664  ++k, ++count_p, ++count_q)
665  values_dofs_actual[c *
666  shape_info.dofs_per_component_on_cell +
667  count_p] =
668  values_dofs[c * dofs_per_comp + count_q];
669  count_q += j + i;
670  }
671  count_q += i * (degree + 1);
672  }
673  }
674  }
675 
676 
677 
678  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
679  inline void
682  dim,
683  fe_degree,
684  n_q_points_1d,
685  Number>::evaluate(const unsigned int n_components,
686  const EvaluationFlags::EvaluationFlags evaluation_flag,
687  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
688  const Number *values_dofs_actual,
689  Number * values_quad,
690  Number * gradients_quad,
691  Number * hessians_quad,
692  Number * scratch_data)
693  {
694  (void)scratch_data;
695 
696  const unsigned int n_dofs = shape_info.dofs_per_component_on_cell;
697  const unsigned int n_q_points = shape_info.n_q_points;
698 
699  using Eval =
701 
702  if (evaluation_flag & EvaluationFlags::values)
703  {
704  const auto shape_values = shape_info.data.front().shape_values.data();
705  auto values_quad_ptr = values_quad;
706  auto values_dofs_actual_ptr = values_dofs_actual;
707 
708  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
709  for (unsigned int c = 0; c < n_components; ++c)
710  {
711  eval.template values<0, true, false>(values_dofs_actual_ptr,
712  values_quad_ptr);
713 
714  values_quad_ptr += n_q_points;
715  values_dofs_actual_ptr += n_dofs;
716  }
717  }
718 
719  if (evaluation_flag & EvaluationFlags::gradients)
720  {
721  const auto shape_gradients =
722  shape_info.data.front().shape_gradients.data();
723  auto gradients_quad_ptr = gradients_quad;
724  auto values_dofs_actual_ptr = values_dofs_actual;
725 
726  for (unsigned int c = 0; c < n_components; ++c)
727  {
728  for (unsigned int d = 0; d < dim; ++d)
729  {
730  Eval eval(nullptr,
731  shape_gradients + n_q_points * n_dofs * d,
732  nullptr,
733  n_dofs,
734  n_q_points);
735 
736  eval.template gradients<0, true, false>(values_dofs_actual_ptr,
737  gradients_quad_ptr);
738 
739  gradients_quad_ptr += n_q_points;
740  }
741  values_dofs_actual_ptr += n_dofs;
742  }
743  }
744 
745  if (evaluation_flag & EvaluationFlags::hessians)
746  {
747  Assert(false, ExcNotImplemented());
748  (void)hessians_quad;
749  }
750  }
751 
752 
753 
754  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
755  inline void
758  dim,
759  fe_degree,
760  n_q_points_1d,
761  Number>::integrate(const unsigned int n_components,
762  const EvaluationFlags::EvaluationFlags integration_flag,
763  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
764  Number * values_dofs_actual,
765  Number * values_quad,
766  Number * gradients_quad,
767  Number * scratch_data,
768  const bool add_into_values_array)
769  {
770  (void)scratch_data;
771 
772  const unsigned int n_dofs = shape_info.dofs_per_component_on_cell;
773  const unsigned int n_q_points = shape_info.n_q_points;
774 
775  using Eval =
777 
778  if (integration_flag & EvaluationFlags::values)
779  {
780  const auto shape_values = shape_info.data.front().shape_values.data();
781  auto values_quad_ptr = values_quad;
782  auto values_dofs_actual_ptr = values_dofs_actual;
783 
784  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
785  for (unsigned int c = 0; c < n_components; ++c)
786  {
787  if (add_into_values_array == false)
788  eval.template values<0, false, false>(values_quad_ptr,
789  values_dofs_actual_ptr);
790  else
791  eval.template values<0, false, true>(values_quad_ptr,
792  values_dofs_actual_ptr);
793 
794  values_quad_ptr += n_q_points;
795  values_dofs_actual_ptr += n_dofs;
796  }
797  }
798 
799  if (integration_flag & EvaluationFlags::gradients)
800  {
801  const auto shape_gradients =
802  shape_info.data.front().shape_gradients.data();
803  auto gradients_quad_ptr = gradients_quad;
804  auto values_dofs_actual_ptr = values_dofs_actual;
805 
806  for (unsigned int c = 0; c < n_components; ++c)
807  {
808  for (unsigned int d = 0; d < dim; ++d)
809  {
810  Eval eval(nullptr,
811  shape_gradients + n_q_points * n_dofs * d,
812  nullptr,
813  n_dofs,
814  n_q_points);
815 
816  if ((add_into_values_array == false &&
817  (integration_flag & EvaluationFlags::values) == false) &&
818  d == 0)
819  eval.template gradients<0, false, false>(
820  gradients_quad_ptr, values_dofs_actual_ptr);
821  else
822  eval.template gradients<0, false, true>(
823  gradients_quad_ptr, values_dofs_actual_ptr);
824 
825  gradients_quad_ptr += n_q_points;
826  }
827  values_dofs_actual_ptr += n_dofs;
828  }
829  }
830  }
831 
832 
833 
843  template <EvaluatorVariant variant,
844  EvaluatorQuantity quantity,
845  int dim,
846  int basis_size_1,
847  int basis_size_2,
848  typename Number,
849  typename Number2>
851  {
852  static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
853  "The second dimension must not be smaller than the first");
854 
877 #ifndef DEBUG
879 #endif
880  static void
882  const unsigned int n_components,
883  const AlignedVector<Number2> &transformation_matrix,
884  const Number * values_in,
885  Number * values_out,
886  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
887  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
888  {
889  Assert(
890  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
891  ExcMessage("The second dimension must not be smaller than the first"));
892 
894 
895  // we do recursion until dim==1 or dim==2 and we have
896  // basis_size_1==basis_size_2. The latter optimization increases
897  // optimization possibilities for the compiler but does only work for
898  // aliased pointers if the sizes are equal.
899  constexpr int next_dim =
900  (dim > 2 ||
901  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
902  dim - 1 :
903  dim;
904 
905  EvaluatorTensorProduct<variant,
906  dim,
907  basis_size_1,
908  (basis_size_1 == 0 ? 0 : basis_size_2),
909  Number,
910  Number2>
911  eval_val(transformation_matrix,
914  basis_size_1_variable,
915  basis_size_2_variable);
916  const unsigned int np_1 =
917  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
918  const unsigned int np_2 =
919  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
920  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
921  ExcMessage("Cannot transform with 0-point basis"));
922  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
923  ExcMessage("Cannot transform with 0-point basis"));
924 
925  // run loop backwards to ensure correctness if values_in aliases with
926  // values_out in case with basis_size_1 < basis_size_2
927  values_in = values_in + n_components * Utilities::fixed_power<dim>(np_1);
928  values_out =
929  values_out + n_components * Utilities::fixed_power<dim>(np_2);
930  for (unsigned int c = n_components; c != 0; --c)
931  {
932  values_in -= Utilities::fixed_power<dim>(np_1);
933  values_out -= Utilities::fixed_power<dim>(np_2);
934  if (next_dim < dim)
935  for (unsigned int q = np_1; q != 0; --q)
937  variant,
938  quantity,
939  next_dim,
940  basis_size_1,
941  basis_size_2,
942  Number,
943  Number2>::do_forward(1,
944  transformation_matrix,
945  values_in +
946  (q - 1) *
947  Utilities::fixed_power<next_dim>(np_1),
948  values_out +
949  (q - 1) *
950  Utilities::fixed_power<next_dim>(np_2),
951  basis_size_1_variable,
952  basis_size_2_variable);
953 
954  // the recursion stops if dim==1 or if dim==2 and
955  // basis_size_1==basis_size_2 (the latter is used because the
956  // compiler generates nicer code)
957  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
958  {
959  eval_val.template values<0, true, false>(values_in, values_out);
960  eval_val.template values<1, true, false>(values_out, values_out);
961  }
962  else if (dim == 1)
963  eval_val.template values<dim - 1, true, false>(values_in,
964  values_out);
965  else
966  eval_val.template values<dim - 1, true, false>(values_out,
967  values_out);
968  }
969  }
970 
1001 #ifndef DEBUG
1003 #endif
1004  static void
1006  const unsigned int n_components,
1007  const AlignedVector<Number2> &transformation_matrix,
1008  const bool add_into_result,
1009  Number * values_in,
1010  Number * values_out,
1011  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
1012  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
1013  {
1014  Assert(
1015  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
1016  ExcMessage("The second dimension must not be smaller than the first"));
1017  Assert(add_into_result == false || values_in != values_out,
1018  ExcMessage(
1019  "Input and output cannot alias with each other when "
1020  "adding the result of the basis change to existing data"));
1021 
1022  Assert(quantity == EvaluatorQuantity::value ||
1023  quantity == EvaluatorQuantity::hessian,
1024  ExcInternalError());
1025 
1026  constexpr int next_dim =
1027  (dim > 2 ||
1028  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
1029  dim - 1 :
1030  dim;
1031  EvaluatorTensorProduct<variant,
1032  dim,
1033  basis_size_1,
1034  (basis_size_1 == 0 ? 0 : basis_size_2),
1035  Number,
1036  Number2>
1037  eval_val(transformation_matrix,
1038  transformation_matrix,
1039  transformation_matrix,
1040  basis_size_1_variable,
1041  basis_size_2_variable);
1042  const unsigned int np_1 =
1043  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
1044  const unsigned int np_2 =
1045  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
1046  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
1047  ExcMessage("Cannot transform with 0-point basis"));
1048  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
1049  ExcMessage("Cannot transform with 0-point basis"));
1050 
1051  for (unsigned int c = 0; c < n_components; ++c)
1052  {
1053  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
1054  {
1055  if (quantity == EvaluatorQuantity::value)
1056  eval_val.template values<1, false, false>(values_in, values_in);
1057  else
1058  eval_val.template hessians<1, false, false>(values_in,
1059  values_in);
1060 
1061  if (add_into_result)
1062  {
1063  if (quantity == EvaluatorQuantity::value)
1064  eval_val.template values<0, false, true>(values_in,
1065  values_out);
1066  else
1067  eval_val.template hessians<0, false, true>(values_in,
1068  values_out);
1069  }
1070  else
1071  {
1072  if (quantity == EvaluatorQuantity::value)
1073  eval_val.template values<0, false, false>(values_in,
1074  values_out);
1075  else
1076  eval_val.template hessians<0, false, false>(values_in,
1077  values_out);
1078  }
1079  }
1080  else
1081  {
1082  if (dim == 1 && add_into_result)
1083  {
1084  if (quantity == EvaluatorQuantity::value)
1085  eval_val.template values<0, false, true>(values_in,
1086  values_out);
1087  else
1088  eval_val.template hessians<0, false, true>(values_in,
1089  values_out);
1090  }
1091  else if (dim == 1)
1092  {
1093  if (quantity == EvaluatorQuantity::value)
1094  eval_val.template values<0, false, false>(values_in,
1095  values_out);
1096  else
1097  eval_val.template hessians<0, false, false>(values_in,
1098  values_out);
1099  }
1100  else
1101  {
1102  if (quantity == EvaluatorQuantity::value)
1103  eval_val.template values<dim - 1, false, false>(values_in,
1104  values_in);
1105  else
1106  eval_val.template hessians<dim - 1, false, false>(
1107  values_in, values_in);
1108  }
1109  }
1110  if (next_dim < dim)
1111  for (unsigned int q = 0; q < np_1; ++q)
1113  quantity,
1114  next_dim,
1115  basis_size_1,
1116  basis_size_2,
1117  Number,
1118  Number2>::
1119  do_backward(1,
1120  transformation_matrix,
1121  add_into_result,
1122  values_in +
1123  q * Utilities::fixed_power<next_dim>(np_2),
1124  values_out +
1125  q * Utilities::fixed_power<next_dim>(np_1),
1126  basis_size_1_variable,
1127  basis_size_2_variable);
1128 
1129  values_in += Utilities::fixed_power<dim>(np_2);
1130  values_out += Utilities::fixed_power<dim>(np_1);
1131  }
1132  }
1133 
1154  static void
1155  do_mass(const unsigned int n_components,
1156  const AlignedVector<Number2> &transformation_matrix,
1157  const AlignedVector<Number> & coefficients,
1158  const Number * values_in,
1159  Number * scratch_data,
1160  Number * values_out)
1161  {
1162  constexpr int next_dim = dim > 1 ? dim - 1 : dim;
1163  Number * my_scratch =
1164  basis_size_1 != basis_size_2 ? scratch_data : values_out;
1165 
1166  const unsigned int size_per_component = Utilities::pow(basis_size_2, dim);
1167  Assert(coefficients.size() == size_per_component ||
1168  coefficients.size() == n_components * size_per_component,
1169  ExcDimensionMismatch(coefficients.size(), size_per_component));
1170  const unsigned int stride =
1171  coefficients.size() == size_per_component ? 0 : 1;
1172 
1173  for (unsigned int q = basis_size_1; q != 0; --q)
1175  variant,
1177  next_dim,
1178  basis_size_1,
1179  basis_size_2,
1180  Number,
1181  Number2>::do_forward(n_components,
1182  transformation_matrix,
1183  values_in +
1184  (q - 1) *
1185  Utilities::pow(basis_size_1, dim - 1),
1186  my_scratch +
1187  (q - 1) *
1188  Utilities::pow(basis_size_2, dim - 1));
1189  EvaluatorTensorProduct<variant,
1190  dim,
1191  basis_size_1,
1192  basis_size_2,
1193  Number,
1194  Number2>
1195  eval_val(transformation_matrix);
1196  const unsigned int n_inner_blocks =
1197  (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
1198  const unsigned int n_blocks = Utilities::pow(basis_size_2, dim - 1);
1199  for (unsigned int ii = 0; ii < n_blocks; ii += n_inner_blocks)
1200  for (unsigned int c = 0; c < n_components; ++c)
1201  {
1202  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1203  eval_val.template values_one_line<dim - 1, true, false>(
1204  my_scratch + i, my_scratch + i);
1205  for (unsigned int q = 0; q < basis_size_2; ++q)
1206  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1207  my_scratch[i + q * n_blocks + c * size_per_component] *=
1208  coefficients[i + q * n_blocks +
1209  c * stride * size_per_component];
1210  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1211  eval_val.template values_one_line<dim - 1, false, false>(
1212  my_scratch + i, my_scratch + i);
1213  }
1214  for (unsigned int q = 0; q < basis_size_1; ++q)
1216  variant,
1218  next_dim,
1219  basis_size_1,
1220  basis_size_2,
1221  Number,
1222  Number2>::do_backward(n_components,
1223  transformation_matrix,
1224  false,
1225  my_scratch +
1226  q * Utilities::pow(basis_size_2, dim - 1),
1227  values_out +
1228  q * Utilities::pow(basis_size_1, dim - 1));
1229  }
1230  };
1231 
1232 
1233 
1246  template <int dim, int fe_degree, typename Number>
1248  {
1249  static void
1250  evaluate(const unsigned int n_components,
1251  const EvaluationFlags::EvaluationFlags evaluation_flag,
1252  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1253  const Number * values_dofs,
1254  Number * values_quad,
1255  Number * gradients_quad,
1256  Number * hessians_quad,
1257  Number * scratch_data);
1258 
1259  static void
1260  integrate(const unsigned int n_components,
1261  const EvaluationFlags::EvaluationFlags integration_flag,
1262  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1263  Number * values_dofs,
1264  Number * values_quad,
1265  Number * gradients_quad,
1266  Number * scratch_data,
1267  const bool add_into_values_array);
1268  };
1269 
1270 
1271 
1272  template <int dim, int fe_degree, typename Number>
1273  inline void
1275  const unsigned int n_components,
1276  const EvaluationFlags::EvaluationFlags evaluation_flag,
1277  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1278  const Number * values_dofs,
1279  Number * values_quad,
1280  Number * gradients_quad,
1281  Number * hessians_quad,
1282  Number *)
1283  {
1285  shape_info.data.front().shape_gradients_collocation_eo.size(),
1286  (fe_degree + 2) / 2 * (fe_degree + 1));
1287 
1289  dim,
1290  fe_degree + 1,
1291  fe_degree + 1,
1292  Number>
1293  eval(AlignedVector<Number>(),
1294  shape_info.data.front().shape_gradients_collocation_eo,
1295  shape_info.data.front().shape_hessians_collocation_eo);
1296  constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
1297 
1298  for (unsigned int c = 0; c < n_components; c++)
1299  {
1300  if (evaluation_flag & EvaluationFlags::values)
1301  for (unsigned int i = 0; i < n_q_points; ++i)
1302  values_quad[i] = values_dofs[i];
1303  if (evaluation_flag &
1305  {
1306  eval.template gradients<0, true, false>(values_dofs,
1307  gradients_quad);
1308  if (dim > 1)
1309  eval.template gradients<1, true, false>(values_dofs,
1310  gradients_quad +
1311  n_q_points);
1312  if (dim > 2)
1313  eval.template gradients<2, true, false>(values_dofs,
1314  gradients_quad +
1315  2 * n_q_points);
1316  }
1317  if (evaluation_flag & EvaluationFlags::hessians)
1318  {
1319  eval.template hessians<0, true, false>(values_dofs, hessians_quad);
1320  if (dim > 1)
1321  {
1322  eval.template gradients<1, true, false>(gradients_quad,
1323  hessians_quad +
1324  dim * n_q_points);
1325  eval.template hessians<1, true, false>(values_dofs,
1326  hessians_quad +
1327  n_q_points);
1328  }
1329  if (dim > 2)
1330  {
1331  eval.template gradients<2, true, false>(gradients_quad,
1332  hessians_quad +
1333  4 * n_q_points);
1334  eval.template gradients<2, true, false>(
1335  gradients_quad + n_q_points, hessians_quad + 5 * n_q_points);
1336  eval.template hessians<2, true, false>(values_dofs,
1337  hessians_quad +
1338  2 * n_q_points);
1339  }
1340  hessians_quad += (dim * (dim + 1)) / 2 * n_q_points;
1341  }
1342  gradients_quad += dim * n_q_points;
1343  values_quad += n_q_points;
1344  values_dofs += n_q_points;
1345  }
1346  }
1347 
1348 
1349 
1350  template <int dim, int fe_degree, typename Number>
1351  inline void
1353  const unsigned int n_components,
1354  const EvaluationFlags::EvaluationFlags integration_flag,
1355  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1356  Number * values_dofs,
1357  Number * values_quad,
1358  Number * gradients_quad,
1359  Number *,
1360  const bool add_into_values_array)
1361  {
1363  shape_info.data.front().shape_gradients_collocation_eo.size(),
1364  (fe_degree + 2) / 2 * (fe_degree + 1));
1365 
1367  dim,
1368  fe_degree + 1,
1369  fe_degree + 1,
1370  Number>
1371  eval(AlignedVector<Number>(),
1372  shape_info.data.front().shape_gradients_collocation_eo,
1373  shape_info.data.front().shape_hessians_collocation_eo);
1374  constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
1375 
1376  for (unsigned int c = 0; c < n_components; c++)
1377  {
1378  if (integration_flag & EvaluationFlags::values)
1379  {
1380  if (add_into_values_array == false)
1381  for (unsigned int i = 0; i < n_q_points; ++i)
1382  values_dofs[i] = values_quad[i];
1383  else
1384  for (unsigned int i = 0; i < n_q_points; ++i)
1385  values_dofs[i] += values_quad[i];
1386  }
1387  if (integration_flag & EvaluationFlags::gradients)
1388  {
1389  if (integration_flag & EvaluationFlags::values ||
1390  add_into_values_array == true)
1391  eval.template gradients<0, false, true>(gradients_quad,
1392  values_dofs);
1393  else
1394  eval.template gradients<0, false, false>(gradients_quad,
1395  values_dofs);
1396  if (dim > 1)
1397  eval.template gradients<1, false, true>(gradients_quad +
1398  n_q_points,
1399  values_dofs);
1400  if (dim > 2)
1401  eval.template gradients<2, false, true>(gradients_quad +
1402  2 * n_q_points,
1403  values_dofs);
1404  }
1405  gradients_quad += dim * n_q_points;
1406  values_quad += n_q_points;
1407  values_dofs += n_q_points;
1408  }
1409  }
1410 
1411 
1412 
1423  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1425  {
1426  static void
1427  evaluate(const unsigned int n_components,
1428  const EvaluationFlags::EvaluationFlags evaluation_flag,
1429  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1430  const Number * values_dofs,
1431  Number * values_quad,
1432  Number * gradients_quad,
1433  Number * hessians_quad,
1434  Number * scratch_data);
1435 
1436  static void
1437  integrate(const unsigned int n_components,
1438  const EvaluationFlags::EvaluationFlags evaluation_flag,
1439  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1440  Number * values_dofs,
1441  Number * values_quad,
1442  Number * gradients_quad,
1443  Number * scratch_data,
1444  const bool add_into_values_array);
1445  };
1446 
1447 
1448 
1449  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1450  inline void
1452  dim,
1453  fe_degree,
1454  n_q_points_1d,
1455  Number>::evaluate(const unsigned int n_components,
1456  const EvaluationFlags::EvaluationFlags evaluation_flag,
1457  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1458  const Number * values_dofs,
1459  Number * values_quad,
1460  Number *gradients_quad,
1461  Number *hessians_quad,
1462  Number *)
1463  {
1464  Assert(n_q_points_1d > fe_degree,
1465  ExcMessage("You lose information when going to a collocation space "
1466  "of lower degree, so the evaluation results would be "
1467  "wrong. Thus, this class does not permit the desired "
1468  "operation."));
1469  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
1470 
1471  for (unsigned int c = 0; c < n_components; c++)
1472  {
1476  dim,
1477  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1478  n_q_points_1d,
1479  Number,
1480  Number>::do_forward(1,
1481  shape_info.data.front().shape_values_eo,
1482  values_dofs,
1483  values_quad);
1484 
1485  // apply derivatives in the collocation space
1486  if (evaluation_flag &
1489  1,
1490  evaluation_flag &
1492  shape_info,
1493  values_quad,
1494  nullptr,
1495  gradients_quad,
1496  hessians_quad,
1497  nullptr);
1498 
1499  values_dofs += shape_info.dofs_per_component_on_cell;
1500  values_quad += n_q_points;
1501  gradients_quad += dim * n_q_points;
1502  hessians_quad += (dim * (dim + 1)) / 2 * n_q_points;
1503  }
1504  }
1505 
1506 
1507 
1508  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1509  inline void
1511  dim,
1512  fe_degree,
1513  n_q_points_1d,
1514  Number>::integrate(const unsigned int n_components,
1515  const EvaluationFlags::EvaluationFlags integration_flag,
1516  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1517  Number *values_dofs,
1518  Number *values_quad,
1519  Number *gradients_quad,
1520  Number *,
1521  const bool add_into_values_array)
1522  {
1523  Assert(n_q_points_1d > fe_degree,
1524  ExcMessage("You lose information when going to a collocation space "
1525  "of lower degree, so the evaluation results would be "
1526  "wrong. Thus, this class does not permit the desired "
1527  "operation."));
1529  shape_info.data.front().shape_gradients_collocation_eo.size(),
1530  (n_q_points_1d + 1) / 2 * n_q_points_1d);
1531  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
1532 
1533  for (unsigned int c = 0; c < n_components; c++)
1534  {
1535  // apply derivatives in collocation space
1536  if (integration_flag & EvaluationFlags::gradients)
1538  integrate(1,
1539  integration_flag & EvaluationFlags::gradients,
1540  shape_info,
1541  values_quad,
1542  nullptr,
1543  gradients_quad,
1544  nullptr,
1545  /*add_into_values_array=*/integration_flag &
1547 
1548  // transform back to the original space
1552  dim,
1553  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1554  n_q_points_1d,
1555  Number,
1556  Number>::do_backward(1,
1557  shape_info.data.front().shape_values_eo,
1558  add_into_values_array,
1559  values_quad,
1560  values_dofs);
1561  gradients_quad += dim * n_q_points;
1562  values_quad += n_q_points;
1563  values_dofs += shape_info.dofs_per_component_on_cell;
1564  }
1565  }
1566 
1567 
1568 
1584  template <int dim, typename Number>
1586  {
1587  template <int fe_degree, int n_q_points_1d>
1588  static bool
1589  run(const unsigned int n_components,
1590  const EvaluationFlags::EvaluationFlags evaluation_flag,
1592  Number *values_dofs_actual,
1593  Number *values_quad,
1594  Number *gradients_quad,
1595  Number *hessians_quad,
1596  Number *scratch_data)
1597  {
1598  // We enable a transformation to collocation for derivatives if it gives
1599  // correct results (first condition), if it is the most efficient choice
1600  // in terms of operation counts (second condition) and if we were able to
1601  // initialize the fields in shape_info.templates.h from the polynomials
1602  // (third condition).
1603  static constexpr bool use_collocation =
1604  n_q_points_1d > fe_degree && n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
1605  n_q_points_1d < 200;
1606 
1607  if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
1608  shape_info.element_type ==
1610  {
1612  evaluate(n_components,
1613  evaluation_flag,
1614  shape_info,
1615  values_dofs_actual,
1616  values_quad,
1617  gradients_quad,
1618  hessians_quad,
1619  scratch_data);
1620  }
1621  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
1622  // shape_info.h for more details
1623  else if (fe_degree >= 0 && use_collocation &&
1624  shape_info.element_type <=
1626  {
1628  dim,
1629  fe_degree,
1630  n_q_points_1d,
1631  Number>::evaluate(n_components,
1632  evaluation_flag,
1633  shape_info,
1634  values_dofs_actual,
1635  values_quad,
1636  gradients_quad,
1637  hessians_quad,
1638  scratch_data);
1639  }
1640  else if (fe_degree >= 0 &&
1641  shape_info.element_type <=
1643  {
1646  dim,
1647  fe_degree,
1648  n_q_points_1d,
1649  Number>::evaluate(n_components,
1650  evaluation_flag,
1651  shape_info,
1652  values_dofs_actual,
1653  values_quad,
1654  gradients_quad,
1655  hessians_quad,
1656  scratch_data);
1657  }
1658  else if (shape_info.element_type ==
1660  {
1663  dim,
1664  fe_degree,
1665  n_q_points_1d,
1666  Number>::evaluate(n_components,
1667  evaluation_flag,
1668  shape_info,
1669  values_dofs_actual,
1670  values_quad,
1671  gradients_quad,
1672  hessians_quad,
1673  scratch_data);
1674  }
1675  else if (shape_info.element_type ==
1677  {
1680  dim,
1681  fe_degree,
1682  n_q_points_1d,
1683  Number>::evaluate(n_components,
1684  evaluation_flag,
1685  shape_info,
1686  values_dofs_actual,
1687  values_quad,
1688  gradients_quad,
1689  hessians_quad,
1690  scratch_data);
1691  }
1692  else if (shape_info.element_type ==
1694  {
1696  dim,
1697  fe_degree,
1698  n_q_points_1d,
1699  Number>::evaluate(n_components,
1700  evaluation_flag,
1701  shape_info,
1702  values_dofs_actual,
1703  values_quad,
1704  gradients_quad,
1705  hessians_quad,
1706  scratch_data);
1707  }
1708  else
1709  {
1712  dim,
1713  fe_degree,
1714  n_q_points_1d,
1715  Number>::evaluate(n_components,
1716  evaluation_flag,
1717  shape_info,
1718  values_dofs_actual,
1719  values_quad,
1720  gradients_quad,
1721  hessians_quad,
1722  scratch_data);
1723  }
1724 
1725  return false;
1726  }
1727  };
1728 
1729 
1730 
1746  template <int dim, typename Number>
1748  {
1749  template <int fe_degree, int n_q_points_1d>
1750  static bool
1751  run(const unsigned int n_components,
1752  const EvaluationFlags::EvaluationFlags integration_flag,
1754  Number * values_dofs_actual,
1755  Number * values_quad,
1756  Number * gradients_quad,
1757  Number * scratch_data,
1758  const bool sum_into_values_array)
1759  {
1760  // We enable a transformation to collocation for derivatives if it gives
1761  // correct results (first condition), if it is the most efficient choice
1762  // in terms of operation counts (second condition) and if we were able to
1763  // initialize the fields in shape_info.templates.h from the polynomials
1764  // (third condition).
1765  constexpr bool use_collocation = n_q_points_1d > fe_degree &&
1766  n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
1767  n_q_points_1d < 200;
1768 
1769  if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
1770  shape_info.element_type ==
1772  {
1774  integrate(n_components,
1775  integration_flag,
1776  shape_info,
1777  values_dofs_actual,
1778  values_quad,
1779  gradients_quad,
1780  scratch_data,
1781  sum_into_values_array);
1782  }
1783  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
1784  // shape_info.h for more details
1785  else if (fe_degree >= 0 && use_collocation &&
1786  shape_info.element_type <=
1788  {
1790  dim,
1791  fe_degree,
1792  n_q_points_1d,
1793  Number>::integrate(n_components,
1794  integration_flag,
1795  shape_info,
1796  values_dofs_actual,
1797  values_quad,
1798  gradients_quad,
1799  scratch_data,
1800  sum_into_values_array);
1801  }
1802  else if (fe_degree >= 0 &&
1803  shape_info.element_type <=
1805  {
1808  dim,
1809  fe_degree,
1810  n_q_points_1d,
1811  Number>::integrate(n_components,
1812  integration_flag,
1813  shape_info,
1814  values_dofs_actual,
1815  values_quad,
1816  gradients_quad,
1817  scratch_data,
1818  sum_into_values_array);
1819  }
1820  else if (shape_info.element_type ==
1822  {
1825  dim,
1826  fe_degree,
1827  n_q_points_1d,
1828  Number>::integrate(n_components,
1829  integration_flag,
1830  shape_info,
1831  values_dofs_actual,
1832  values_quad,
1833  gradients_quad,
1834  scratch_data,
1835  sum_into_values_array);
1836  }
1837  else if (shape_info.element_type ==
1839  {
1842  dim,
1843  fe_degree,
1844  n_q_points_1d,
1845  Number>::integrate(n_components,
1846  integration_flag,
1847  shape_info,
1848  values_dofs_actual,
1849  values_quad,
1850  gradients_quad,
1851  scratch_data,
1852  sum_into_values_array);
1853  }
1854  else if (shape_info.element_type ==
1856  {
1858  dim,
1859  fe_degree,
1860  n_q_points_1d,
1861  Number>::integrate(n_components,
1862  integration_flag,
1863  shape_info,
1864  values_dofs_actual,
1865  values_quad,
1866  gradients_quad,
1867  scratch_data,
1868  sum_into_values_array);
1869  }
1870  else
1871  {
1874  dim,
1875  fe_degree,
1876  n_q_points_1d,
1877  Number>::integrate(n_components,
1878  integration_flag,
1879  shape_info,
1880  values_dofs_actual,
1881  values_quad,
1882  gradients_quad,
1883  scratch_data,
1884  sum_into_values_array);
1885  }
1886 
1887  return false;
1888  }
1889  };
1890 
1891 
1892 
1893  template <bool symmetric_evaluate,
1894  int dim,
1895  int fe_degree,
1896  int n_q_points_1d,
1897  typename Number>
1899  {
1900  // We enable a transformation to collocation for derivatives if it gives
1901  // correct results (first two conditions), if it is the most efficient
1902  // choice in terms of operation counts (third condition) and if we were
1903  // able to initialize the fields in shape_info.templates.h from the
1904  // polynomials (fourth condition).
1905  static constexpr bool use_collocation =
1906  symmetric_evaluate &&
1907  n_q_points_1d > fe_degree &&n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
1908  n_q_points_1d < 200;
1909 
1910  static void
1911  evaluate_in_face(const unsigned int n_components,
1913  Number * values_dofs,
1914  Number * values_quad,
1915  Number * gradients_quad,
1916  Number * scratch_data,
1917  const bool evaluate_val,
1918  const bool evaluate_grad,
1919  const unsigned int subface_index)
1920  {
1921  const AlignedVector<Number> &val1 =
1922  symmetric_evaluate ?
1923  data.data.front().shape_values_eo :
1924  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1925  data.data.front().shape_values :
1926  data.data.front().values_within_subface[subface_index % 2]);
1927  const AlignedVector<Number> &val2 =
1928  symmetric_evaluate ?
1929  data.data.front().shape_values_eo :
1930  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1931  data.data.front().shape_values :
1932  data.data.front().values_within_subface[subface_index / 2]);
1933 
1934  const AlignedVector<Number> &grad1 =
1935  symmetric_evaluate ?
1936  data.data.front().shape_gradients_eo :
1937  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1938  data.data.front().shape_gradients :
1939  data.data.front().gradients_within_subface[subface_index % 2]);
1940  const AlignedVector<Number> &grad2 =
1941  symmetric_evaluate ?
1942  data.data.front().shape_gradients_eo :
1943  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1944  data.data.front().shape_gradients :
1945  data.data.front().gradients_within_subface[subface_index / 2]);
1946 
1947  using Eval =
1948  internal::EvaluatorTensorProduct<symmetric_evaluate ?
1951  dim - 1,
1952  fe_degree + 1,
1953  n_q_points_1d,
1954  Number>;
1955  Eval eval1(val1,
1956  grad1,
1958  data.data.front().fe_degree + 1,
1959  data.data.front().n_q_points_1d);
1960  Eval eval2(val2,
1961  grad2,
1963  data.data.front().fe_degree + 1,
1964  data.data.front().n_q_points_1d);
1965 
1966  const unsigned int size_deg =
1967  fe_degree > -1 ?
1968  Utilities::pow(fe_degree + 1, dim - 1) :
1969  (dim > 1 ?
1970  Utilities::fixed_power<dim - 1>(data.data.front().fe_degree + 1) :
1971  1);
1972 
1973  const unsigned int n_q_points = fe_degree > -1 ?
1974  Utilities::pow(n_q_points_1d, dim - 1) :
1975  data.n_q_points_face;
1976 
1977  if (evaluate_grad == false)
1978  for (unsigned int c = 0; c < n_components; ++c)
1979  {
1980  switch (dim)
1981  {
1982  case 3:
1983  eval1.template values<0, true, false>(values_dofs,
1984  values_quad);
1985  eval2.template values<1, true, false>(values_quad,
1986  values_quad);
1987  break;
1988  case 2:
1989  eval1.template values<0, true, false>(values_dofs,
1990  values_quad);
1991  break;
1992  case 1:
1993  values_quad[0] = values_dofs[0];
1994  break;
1995  default:
1996  Assert(false, ExcNotImplemented());
1997  }
1998  values_dofs += 2 * size_deg;
1999  values_quad += n_q_points;
2000  }
2001  else
2002  for (unsigned int c = 0; c < n_components; ++c)
2003  {
2004  switch (dim)
2005  {
2006  case 3:
2007  if (use_collocation)
2008  {
2009  eval1.template values<0, true, false>(values_dofs,
2010  values_quad);
2011  eval1.template values<1, true, false>(values_quad,
2012  values_quad);
2015  dim - 1,
2016  n_q_points_1d,
2017  n_q_points_1d,
2018  Number>
2019  eval_grad(
2021  data.data.front().shape_gradients_collocation_eo,
2023  eval_grad.template gradients<0, true, false>(
2024  values_quad, gradients_quad);
2025  eval_grad.template gradients<1, true, false>(
2026  values_quad, gradients_quad + n_q_points);
2027  }
2028  else
2029  {
2030  eval1.template gradients<0, true, false>(values_dofs,
2031  scratch_data);
2032  eval2.template values<1, true, false>(scratch_data,
2033  gradients_quad);
2034 
2035  eval1.template values<0, true, false>(values_dofs,
2036  scratch_data);
2037  eval2.template gradients<1, true, false>(scratch_data,
2038  gradients_quad +
2039  n_q_points);
2040 
2041  if (evaluate_val == true)
2042  eval2.template values<1, true, false>(scratch_data,
2043  values_quad);
2044  }
2045  eval1.template values<0, true, false>(values_dofs + size_deg,
2046  scratch_data);
2047  eval2.template values<1, true, false>(
2048  scratch_data, gradients_quad + (dim - 1) * n_q_points);
2049 
2050  break;
2051  case 2:
2052  eval1.template values<0, true, false>(values_dofs + size_deg,
2053  gradients_quad +
2054  (dim - 1) *
2055  n_q_points);
2056  eval1.template gradients<0, true, false>(values_dofs,
2057  gradients_quad);
2058  if (evaluate_val == true)
2059  eval1.template values<0, true, false>(values_dofs,
2060  values_quad);
2061  break;
2062  case 1:
2063  values_quad[0] = values_dofs[0];
2064  gradients_quad[0] = values_dofs[1];
2065  break;
2066  default:
2067  AssertThrow(false, ExcNotImplemented());
2068  }
2069  values_dofs += 2 * size_deg;
2070  values_quad += n_q_points;
2071  gradients_quad += dim * n_q_points;
2072  }
2073  }
2074 
2075  static void
2076  integrate_in_face(const unsigned int n_components,
2078  Number * values_dofs,
2079  Number * values_quad,
2080  Number * gradients_quad,
2081  Number * scratch_data,
2082  const bool integrate_val,
2083  const bool integrate_grad,
2084  const unsigned int subface_index)
2085  {
2086  const AlignedVector<Number> &val1 =
2087  symmetric_evaluate ?
2088  data.data.front().shape_values_eo :
2089  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
2090  data.data.front().shape_values :
2091  data.data.front().values_within_subface[subface_index % 2]);
2092  const AlignedVector<Number> &val2 =
2093  symmetric_evaluate ?
2094  data.data.front().shape_values_eo :
2095  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
2096  data.data.front().shape_values :
2097  data.data.front().values_within_subface[subface_index / 2]);
2098 
2099  const AlignedVector<Number> &grad1 =
2100  symmetric_evaluate ?
2101  data.data.front().shape_gradients_eo :
2102  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
2103  data.data.front().shape_gradients :
2104  data.data.front().gradients_within_subface[subface_index % 2]);
2105  const AlignedVector<Number> &grad2 =
2106  symmetric_evaluate ?
2107  data.data.front().shape_gradients_eo :
2108  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
2109  data.data.front().shape_gradients :
2110  data.data.front().gradients_within_subface[subface_index / 2]);
2111 
2112  using Eval =
2113  internal::EvaluatorTensorProduct<symmetric_evaluate ?
2116  dim - 1,
2117  fe_degree + 1,
2118  n_q_points_1d,
2119  Number>;
2120  Eval eval1(val1,
2121  grad1,
2122  val1,
2123  data.data.front().fe_degree + 1,
2124  data.data.front().n_q_points_1d);
2125  Eval eval2(val2,
2126  grad2,
2127  val1,
2128  data.data.front().fe_degree + 1,
2129  data.data.front().n_q_points_1d);
2130 
2131  const unsigned int size_deg =
2132  fe_degree > -1 ?
2133  Utilities::pow(fe_degree + 1, dim - 1) :
2134  (dim > 1 ?
2135  Utilities::fixed_power<dim - 1>(data.data.front().fe_degree + 1) :
2136  1);
2137 
2138  const unsigned int n_q_points = fe_degree > -1 ?
2139  Utilities::pow(n_q_points_1d, dim - 1) :
2140  data.n_q_points_face;
2141 
2142  if (integrate_grad == false)
2143  for (unsigned int c = 0; c < n_components; ++c)
2144  {
2145  switch (dim)
2146  {
2147  case 3:
2148  eval2.template values<1, false, false>(values_quad,
2149  values_quad);
2150  eval1.template values<0, false, false>(values_quad,
2151  values_dofs);
2152  break;
2153  case 2:
2154  eval1.template values<0, false, false>(values_quad,
2155  values_dofs);
2156  break;
2157  case 1:
2158  values_dofs[0] = values_quad[0];
2159  break;
2160  default:
2161  Assert(false, ExcNotImplemented());
2162  }
2163  values_dofs += 2 * size_deg;
2164  values_quad += n_q_points;
2165  }
2166  else
2167  for (unsigned int c = 0; c < n_components; ++c)
2168  {
2169  switch (dim)
2170  {
2171  case 3:
2172  eval2.template values<1, false, false>(gradients_quad +
2173  2 * n_q_points,
2174  gradients_quad +
2175  2 * n_q_points);
2176  eval1.template values<0, false, false>(
2177  gradients_quad + 2 * n_q_points, values_dofs + size_deg);
2178  if (use_collocation)
2179  {
2182  dim - 1,
2183  n_q_points_1d,
2184  n_q_points_1d,
2185  Number>
2186  eval_grad(
2188  data.data.front().shape_gradients_collocation_eo,
2190  if (integrate_val)
2191  eval_grad.template gradients<1, false, true>(
2192  gradients_quad + n_q_points, values_quad);
2193  else
2194  eval_grad.template gradients<1, false, false>(
2195  gradients_quad + n_q_points, values_quad);
2196  eval_grad.template gradients<0, false, true>(
2197  gradients_quad, values_quad);
2198  eval1.template values<1, false, false>(values_quad,
2199  values_quad);
2200  eval1.template values<0, false, false>(values_quad,
2201  values_dofs);
2202  }
2203  else
2204  {
2205  if (integrate_val)
2206  {
2207  eval2.template values<1, false, false>(values_quad,
2208  scratch_data);
2209  eval2.template gradients<1, false, true>(
2210  gradients_quad + n_q_points, scratch_data);
2211  }
2212  else
2213  eval2.template gradients<1, false, false>(
2214  gradients_quad + n_q_points, scratch_data);
2215 
2216  eval1.template values<0, false, false>(scratch_data,
2217  values_dofs);
2218  eval2.template values<1, false, false>(gradients_quad,
2219  scratch_data);
2220  eval1.template gradients<0, false, true>(scratch_data,
2221  values_dofs);
2222  }
2223  break;
2224  case 2:
2225  eval1.template values<0, false, false>(
2226  gradients_quad + n_q_points, values_dofs + size_deg);
2227  eval1.template gradients<0, false, false>(gradients_quad,
2228  values_dofs);
2229  if (integrate_val == true)
2230  eval1.template values<0, false, true>(values_quad,
2231  values_dofs);
2232  break;
2233  case 1:
2234  values_dofs[0] = values_quad[0];
2235  values_dofs[1] = gradients_quad[0];
2236  break;
2237  default:
2238  AssertThrow(false, ExcNotImplemented());
2239  }
2240  values_dofs += 2 * size_deg;
2241  values_quad += n_q_points;
2242  gradients_quad += dim * n_q_points;
2243  }
2244  }
2245  };
2246 
2247 
2248 
2249  template <int dim, int fe_degree, typename Number, bool lex_faces = false>
2251  {
2252  template <bool do_evaluate, bool add_into_output>
2253  static void
2254  interpolate(const unsigned int n_components,
2256  const Number * input,
2257  Number * output,
2258  const bool do_gradients,
2259  const unsigned int face_no)
2260  {
2261  Assert(static_cast<unsigned int>(fe_degree) ==
2262  data.data.front().fe_degree ||
2263  fe_degree == -1,
2264  ExcInternalError());
2265 
2266  interpolate_generic<do_evaluate, add_into_output>(
2267  n_components,
2268  input,
2269  output,
2270  do_gradients,
2271  face_no,
2272  data.data.front().fe_degree + 1,
2273  data.data.front().shape_data_on_face,
2275  2 * data.dofs_per_component_on_face);
2276  }
2277 
2281  template <bool do_evaluate, bool add_into_output>
2282  static void
2283  interpolate_quadrature(const unsigned int n_components,
2285  const Number * input,
2286  Number * output,
2287  const bool do_gradients,
2288  const unsigned int face_no)
2289  {
2290  Assert(static_cast<unsigned int>(fe_degree + 1) ==
2291  data.data.front().quadrature.size() ||
2292  fe_degree == -1,
2293  ExcInternalError());
2294 
2295  interpolate_generic<do_evaluate, add_into_output>(
2296  n_components,
2297  input,
2298  output,
2299  do_gradients,
2300  face_no,
2301  data.data.front().quadrature.size(),
2302  data.data.front().quadrature_data_on_face,
2303  data.n_q_points,
2304  data.n_q_points_face);
2305  }
2306 
2307  private:
2308  template <bool do_evaluate, bool add_into_output, int face_direction = 0>
2309  static void
2310  interpolate_generic(const unsigned int n_components,
2311  const Number * input,
2312  Number * output,
2313  const bool do_gradients,
2314  const unsigned int face_no,
2315  const unsigned int n_points_1d,
2316  const std::array<AlignedVector<Number>, 2> &shape_data,
2317  const unsigned int dofs_per_component_on_cell,
2318  const unsigned int dofs_per_component_on_face)
2319  {
2320  if (face_direction == face_no / 2)
2321  {
2323  dim,
2324  fe_degree + 1,
2325  0,
2326  Number>
2327  evalf(shape_data[face_no % 2],
2330  n_points_1d,
2331  0);
2332 
2333  const unsigned int in_stride = do_evaluate ?
2334  dofs_per_component_on_cell :
2335  dofs_per_component_on_face;
2336  const unsigned int out_stride = do_evaluate ?
2337  dofs_per_component_on_face :
2338  dofs_per_component_on_cell;
2339 
2340  for (unsigned int c = 0; c < n_components; c++)
2341  {
2342  if (do_gradients)
2343  evalf.template apply_face<face_direction,
2344  do_evaluate,
2345  add_into_output,
2346  1,
2347  lex_faces>(input, output);
2348  else
2349  evalf.template apply_face<face_direction,
2350  do_evaluate,
2351  add_into_output,
2352  0,
2353  lex_faces>(input, output);
2354  input += in_stride;
2355  output += out_stride;
2356  }
2357  }
2358  else if (face_direction < dim)
2359  {
2360  interpolate_generic<do_evaluate,
2361  add_into_output,
2362  std::min(face_direction + 1, dim - 1)>(
2363  n_components,
2364  input,
2365  output,
2366  do_gradients,
2367  face_no,
2368  n_points_1d,
2369  shape_data,
2370  dofs_per_component_on_cell,
2371  dofs_per_component_on_face);
2372  }
2373  }
2374  };
2375 
2376 
2377 
2378  // internal helper function for reading data; base version of different types
2379  template <typename VectorizedArrayType, typename Number2>
2380  void
2381  do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
2382  {
2383  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2384  dst[v] = src_ptr[v];
2385  }
2386 
2387 
2388 
2389  // internal helper function for reading data; specialized version where we
2390  // can use a dedicated load function
2391  template <typename Number, unsigned int width>
2392  void
2394  {
2395  dst.load(src_ptr);
2396  }
2397 
2398 
2399 
2400  // internal helper function for reading data; base version of different types
2401  template <typename VectorizedArrayType, typename Number2>
2402  void
2403  do_vectorized_gather(const Number2 * src_ptr,
2404  const unsigned int * indices,
2405  VectorizedArrayType &dst)
2406  {
2407  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2408  dst[v] = src_ptr[indices[v]];
2409  }
2410 
2411 
2412 
2413  // internal helper function for reading data; specialized version where we
2414  // can use a dedicated gather function
2415  template <typename Number, unsigned int width>
2416  void
2417  do_vectorized_gather(const Number * src_ptr,
2418  const unsigned int * indices,
2420  {
2421  dst.gather(src_ptr, indices);
2422  }
2423 
2424 
2425 
2426  // internal helper function for reading data; base version of different types
2427  template <typename VectorizedArrayType, typename Number2>
2428  void
2429  do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
2430  {
2431  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2432  dst_ptr[v] += src[v];
2433  }
2434 
2435 
2436 
2437  // internal helper function for reading data; specialized version where we
2438  // can use a dedicated load function
2439  template <typename Number, unsigned int width>
2440  void
2442  {
2444  tmp.load(dst_ptr);
2445  (tmp + src).store(dst_ptr);
2446  }
2447 
2448 
2449 
2450  // internal helper function for reading data; base version of different types
2451  template <typename VectorizedArrayType, typename Number2>
2452  void
2453  do_vectorized_scatter_add(const VectorizedArrayType src,
2454  const unsigned int * indices,
2455  Number2 * dst_ptr)
2456  {
2457  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2458  dst_ptr[indices[v]] += src[v];
2459  }
2460 
2461 
2462 
2463  // internal helper function for reading data; specialized version where we
2464  // can use a dedicated gather function
2465  template <typename Number, unsigned int width>
2466  void
2468  const unsigned int * indices,
2469  Number * dst_ptr)
2470  {
2471 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
2472  for (unsigned int v = 0; v < width; ++v)
2473  dst_ptr[indices[v]] += src[v];
2474 #else
2476  tmp.gather(dst_ptr, indices);
2477  (tmp + src).scatter(indices, dst_ptr);
2478 #endif
2479  }
2480 
2481 
2482 
2483  template <typename Number>
2484  void
2485  adjust_for_face_orientation(const unsigned int dim,
2486  const unsigned int n_components,
2487  const unsigned int face_orientation,
2488  const Table<2, unsigned int> &orientation_map,
2489  const bool integrate,
2490  const bool values,
2491  const bool gradients,
2492  const unsigned int n_q_points,
2493  Number * tmp_values,
2494  Number * values_quad,
2495  Number * gradients_quad)
2496  {
2497  Assert(face_orientation, ExcInternalError());
2498  const unsigned int *orientation = &orientation_map[face_orientation][0];
2499  for (unsigned int c = 0; c < n_components; ++c)
2500  {
2501  if (values == true)
2502  {
2503  if (integrate)
2504  for (unsigned int q = 0; q < n_q_points; ++q)
2505  tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
2506  else
2507  for (unsigned int q = 0; q < n_q_points; ++q)
2508  tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
2509  for (unsigned int q = 0; q < n_q_points; ++q)
2510  values_quad[c * n_q_points + q] = tmp_values[q];
2511  }
2512  if (gradients == true)
2513  for (unsigned int d = 0; d < dim; ++d)
2514  {
2515  if (integrate)
2516  for (unsigned int q = 0; q < n_q_points; ++q)
2517  tmp_values[q] =
2518  gradients_quad[(c * dim + d) * n_q_points + orientation[q]];
2519  else
2520  for (unsigned int q = 0; q < n_q_points; ++q)
2521  tmp_values[orientation[q]] =
2522  gradients_quad[(c * dim + d) * n_q_points + q];
2523  for (unsigned int q = 0; q < n_q_points; ++q)
2524  gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q];
2525  }
2526  }
2527  }
2528 
2529 
2530 
2531  template <int dim, typename VectorizedArrayType>
2533  {
2534  template <int fe_degree, int n_q_points_1d>
2535  static bool
2536  run(const unsigned int n_components,
2538  const VectorizedArrayType * values_array,
2539  VectorizedArrayType * values_quad,
2540  VectorizedArrayType * gradients_quad,
2541  VectorizedArrayType * scratch_data,
2542  const bool evaluate_values,
2543  const bool evaluate_gradients,
2544  const unsigned int face_no,
2545  const unsigned int subface_index,
2546  const unsigned int face_orientation,
2547  const Table<2, unsigned int> &orientation_map)
2548  {
2549  if (data.element_type == MatrixFreeFunctions::tensor_none)
2550  {
2551  const unsigned int n_dofs = data.dofs_per_component_on_cell;
2552  const unsigned int n_q_points = data.n_q_points_faces[face_no];
2553  const auto shape_info = data.data.front();
2554 
2556  1,
2557  0,
2558  0,
2559  VectorizedArrayType,
2560  VectorizedArrayType>;
2561 
2562  if (evaluate_values)
2563  {
2564  const auto shape_values =
2565  &shape_info.shape_values_face(face_no, face_orientation, 0);
2566 
2567  auto values_quad_ptr = values_quad;
2568  auto values_dofs_actual_ptr = values_array;
2569 
2570  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
2571  for (unsigned int c = 0; c < n_components; ++c)
2572  {
2573  eval.template values<0, true, false>(values_dofs_actual_ptr,
2574  values_quad_ptr);
2575 
2576  values_quad_ptr += n_q_points;
2577  values_dofs_actual_ptr += n_dofs;
2578  }
2579  }
2580 
2581  if (evaluate_gradients)
2582  {
2583  auto gradients_quad_ptr = gradients_quad;
2584  auto values_dofs_actual_ptr = values_array;
2585 
2586  std::array<const VectorizedArrayType *, dim> shape_gradients;
2587  for (unsigned int d = 0; d < dim; ++d)
2588  shape_gradients[d] = &shape_info.shape_gradients_face(
2589  face_no, face_orientation, d, 0);
2590 
2591  for (unsigned int c = 0; c < n_components; ++c)
2592  {
2593  for (unsigned int d = 0; d < dim; ++d)
2594  {
2595  Eval eval(nullptr,
2596  shape_gradients[d],
2597  nullptr,
2598  n_dofs,
2599  n_q_points);
2600 
2601  eval.template gradients<0, true, false>(
2602  values_dofs_actual_ptr, gradients_quad_ptr);
2603 
2604  gradients_quad_ptr += n_q_points;
2605  }
2606  values_dofs_actual_ptr += n_dofs;
2607  }
2608  }
2609 
2610 
2611  return true;
2612  }
2613 
2614  constexpr unsigned int static_dofs_per_face =
2615  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
2617  const unsigned int dofs_per_face =
2618  fe_degree > -1 ?
2619  static_dofs_per_face :
2620  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
2621 
2622  VectorizedArrayType *temp1 = scratch_data;
2623 
2625  template interpolate<true, false>(
2626  n_components, data, values_array, temp1, evaluate_gradients, face_no);
2627 
2628  const unsigned int n_q_points_1d_actual =
2629  fe_degree > -1 ? n_q_points_1d : 0;
2630  if (fe_degree > -1 &&
2631  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
2632  data.element_type <= MatrixFreeFunctions::tensor_symmetric)
2634  true,
2635  dim,
2636  fe_degree,
2637  n_q_points_1d_actual,
2638  VectorizedArrayType>::evaluate_in_face(n_components,
2639  data,
2640  temp1,
2641  values_quad,
2642  gradients_quad,
2643  scratch_data + 2 *
2644  n_components *
2645  dofs_per_face,
2646  evaluate_values,
2647  evaluate_gradients,
2648  subface_index);
2649  else
2651  false,
2652  dim,
2653  fe_degree,
2654  n_q_points_1d_actual,
2655  VectorizedArrayType>::evaluate_in_face(n_components,
2656  data,
2657  temp1,
2658  values_quad,
2659  gradients_quad,
2660  scratch_data + 2 *
2661  n_components *
2662  dofs_per_face,
2663  evaluate_values,
2664  evaluate_gradients,
2665  subface_index);
2666 
2667  if (face_orientation)
2669  n_components,
2670  face_orientation,
2671  orientation_map,
2672  false,
2673  evaluate_values,
2674  evaluate_gradients,
2675  data.n_q_points_face,
2676  scratch_data,
2677  values_quad,
2678  gradients_quad);
2679 
2680  return false;
2681  }
2682  };
2683 
2684 
2685 
2686  template <int dim, typename VectorizedArrayType>
2688  {
2689  template <int fe_degree, int n_q_points_1d>
2690  static bool
2691  run(const unsigned int n_components,
2693  VectorizedArrayType * values_array,
2694  VectorizedArrayType * values_quad,
2695  VectorizedArrayType * gradients_quad,
2696  VectorizedArrayType * scratch_data,
2697  const bool integrate_values,
2698  const bool integrate_gradients,
2699  const unsigned int face_no,
2700  const unsigned int subface_index,
2701  const unsigned int face_orientation,
2702  const Table<2, unsigned int> &orientation_map)
2703  {
2704  if (data.element_type == MatrixFreeFunctions::tensor_none)
2705  {
2706  const unsigned int n_dofs = data.dofs_per_component_on_cell;
2707  const unsigned int n_q_points = data.n_q_points_faces[face_no];
2708  const auto shape_info = data.data.front();
2709 
2711  1,
2712  0,
2713  0,
2714  VectorizedArrayType,
2715  VectorizedArrayType>;
2716 
2717  if (integrate_values)
2718  {
2719  const auto shape_values =
2720  &shape_info.shape_values_face(face_no, face_orientation, 0);
2721 
2722  auto values_quad_ptr = values_quad;
2723  auto values_dofs_actual_ptr = values_array;
2724 
2725  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
2726  for (unsigned int c = 0; c < n_components; ++c)
2727  {
2728  eval.template values<0, false, false>(values_quad_ptr,
2729  values_dofs_actual_ptr);
2730 
2731  values_quad_ptr += n_q_points;
2732  values_dofs_actual_ptr += n_dofs;
2733  }
2734  }
2735 
2736  if (integrate_gradients)
2737  {
2738  auto gradients_quad_ptr = gradients_quad;
2739  auto values_dofs_actual_ptr = values_array;
2740 
2741  std::array<const VectorizedArrayType *, dim> shape_gradients;
2742  for (unsigned int d = 0; d < dim; ++d)
2743  shape_gradients[d] = &shape_info.shape_gradients_face(
2744  face_no, face_orientation, d, 0);
2745 
2746  for (unsigned int c = 0; c < n_components; ++c)
2747  {
2748  for (unsigned int d = 0; d < dim; ++d)
2749  {
2750  Eval eval(nullptr,
2751  shape_gradients[d],
2752  nullptr,
2753  n_dofs,
2754  n_q_points);
2755 
2756  if ((integrate_values == false) && d == 0)
2757  eval.template gradients<0, false, false>(
2758  gradients_quad_ptr, values_dofs_actual_ptr);
2759  else
2760  eval.template gradients<0, false, true>(
2761  gradients_quad_ptr, values_dofs_actual_ptr);
2762 
2763  gradients_quad_ptr += n_q_points;
2764  }
2765  values_dofs_actual_ptr += n_dofs;
2766  }
2767  }
2768 
2769 
2770  return true;
2771  }
2772 
2773  if (face_orientation)
2775  n_components,
2776  face_orientation,
2777  orientation_map,
2778  true,
2779  integrate_values,
2780  integrate_gradients,
2781  data.n_q_points_face,
2782  scratch_data,
2783  values_quad,
2784  gradients_quad);
2785 
2786  constexpr unsigned int static_dofs_per_face =
2787  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
2789  const unsigned int dofs_per_face =
2790  fe_degree > -1 ?
2791  static_dofs_per_face :
2792  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
2793 
2794  VectorizedArrayType *temp1 = scratch_data;
2795 
2796  const unsigned int n_q_points_1d_actual =
2797  fe_degree > -1 ? n_q_points_1d : 0;
2798  if (fe_degree > -1 &&
2802  true,
2803  dim,
2804  fe_degree,
2805  n_q_points_1d_actual,
2806  VectorizedArrayType>::integrate_in_face(n_components,
2807  data,
2808  temp1,
2809  values_quad,
2810  gradients_quad,
2811  scratch_data +
2812  2 * n_components *
2813  dofs_per_face,
2814  integrate_values,
2815  integrate_gradients,
2816  subface_index);
2817  else
2819  false,
2820  dim,
2821  fe_degree,
2822  n_q_points_1d_actual,
2823  VectorizedArrayType>::integrate_in_face(n_components,
2824  data,
2825  temp1,
2826  values_quad,
2827  gradients_quad,
2828  scratch_data +
2829  2 * n_components *
2830  dofs_per_face,
2831  integrate_values,
2832  integrate_gradients,
2833  subface_index);
2834 
2836  template interpolate<false, false>(n_components,
2837  data,
2838  temp1,
2839  values_array,
2840  integrate_gradients,
2841  face_no);
2842  return false;
2843  }
2844  };
2845 
2846 
2847 
2848  template <int n_face_orientations, typename Processor>
2849  static bool
2851  {
2852  auto n_components = proc.n_components;
2853  auto integrate = proc.integrate;
2854  auto global_vector_ptr = proc.global_vector_ptr;
2855  auto &sm_ptr = proc.sm_ptr;
2856  auto &data = proc.data;
2857  auto &dof_info = proc.dof_info;
2858  auto values_quad = proc.values_quad;
2859  auto gradients_quad = proc.gradients_quad;
2860  auto scratch_data = proc.scratch_data;
2861  auto do_values = proc.do_values;
2862  auto do_gradients = proc.do_gradients;
2863  auto active_fe_index = proc.active_fe_index;
2864  auto first_selected_component = proc.first_selected_component;
2865  auto cells = proc.cells;
2866  auto face_nos = proc.face_nos;
2867  auto subface_index = proc.subface_index;
2868  auto dof_access_index = proc.dof_access_index;
2869  auto face_orientations = proc.face_orientations;
2870  auto &orientation_map = proc.orientation_map;
2871 
2872  static const int dim = Processor::dim_;
2873  static const int fe_degree = Processor::fe_degree_;
2874  using VectorizedArrayType = typename Processor::VectorizedArrayType_;
2875 
2876  using Number = typename Processor::Number_;
2877  using Number2_ = typename Processor::Number2_;
2878 
2879  const unsigned int cell = cells[0];
2880 
2881  // In the case of integration, we do not need to reshuffle the
2882  // data at the quadrature points to adjust for the face
2883  // orientation if the shape functions are nodal at the cell
2884  // boundaries (and we only requested the integration of the
2885  // values) or Hermite shape functions are used. These cases are
2886  // handled later when the values are written back into the
2887  // glrobal vector.
2888  if (integrate &&
2889  (face_orientations[0] > 0 &&
2890  (subface_index < GeometryInfo<dim>::max_children_per_cell ||
2891  !(((do_gradients == false &&
2892  data.data.front().nodal_at_cell_boundaries == true &&
2893  fe_degree > 0) ||
2894  (data.element_type ==
2896  fe_degree > 1)) &&
2897  (dof_info.index_storage_variants[dof_access_index][cell] ==
2899  interleaved_contiguous ||
2900  dof_info.index_storage_variants[dof_access_index][cell] ==
2902  interleaved_contiguous_strided ||
2903  dof_info.index_storage_variants[dof_access_index][cell] ==
2905  interleaved_contiguous_mixed_strides ||
2906  dof_info.index_storage_variants[dof_access_index][cell] ==
2908  contiguous)))))
2909  {
2910  AssertDimension(n_face_orientations, 1);
2912  n_components,
2913  face_orientations[0],
2914  orientation_map,
2915  true,
2916  do_values,
2917  do_gradients,
2918  data.n_q_points_face,
2919  scratch_data,
2920  values_quad,
2921  gradients_quad);
2922  }
2923 
2924  // we know that the gradient weights for the Hermite case on the
2925  // right (side==1) are the negative from the value at the left
2926  // (side==0), so we only read out one of them.
2927  VectorizedArrayType grad_weight =
2928  (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
2929  data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ?
2930  data.data.front()
2931  .shape_data_on_face[0][fe_degree + (integrate ?
2932  (2 - (face_nos[0] % 2)) :
2933  (1 + (face_nos[0] % 2)))] :
2934  VectorizedArrayType(0.0 /*dummy*/);
2935 
2936  constexpr unsigned int static_dofs_per_component =
2937  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) :
2939  constexpr unsigned int static_dofs_per_face =
2940  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
2942  const unsigned int dofs_per_face =
2943  fe_degree > -1 ? static_dofs_per_face :
2944  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
2945 
2946  VectorizedArrayType *temp1 = scratch_data;
2947 
2948  const unsigned int dummy = 0;
2949 
2950  // re-orientation
2951  std::array<const unsigned int *, n_face_orientations> orientation = {};
2952 
2953  if (n_face_orientations == 1)
2954  orientation[0] = (data.data.front().nodal_at_cell_boundaries == true) ?
2955  &data.face_orientations[face_orientations[0]][0] :
2956  &dummy;
2957  else
2958  {
2959  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2960  {
2961  // the loop breaks once an invalid_unsigned_int is hit for
2962  // all cases except the exterior faces in the ECL loop (where
2963  // some faces might be at the boundaries but others not)
2964  if (cells[v] == numbers::invalid_unsigned_int)
2965  continue;
2966 
2967  orientation[v] =
2968  (data.data.front().nodal_at_cell_boundaries == true) ?
2969  &data.face_orientations[face_orientations[v]][0] :
2970  &dummy;
2971  }
2972  }
2973 
2974  // face_to_cell_index_hermite
2975  std::array<const unsigned int *, n_face_orientations> index_array_hermite =
2976  {};
2977 
2978  if (n_face_orientations == 1)
2979  index_array_hermite[0] =
2980  (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
2981  data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ?
2982  &data.face_to_cell_index_hermite(face_nos[0], 0) :
2983  &dummy;
2984 
2985  if (n_face_orientations > 1 &&
2986  data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
2987  data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite)
2988  {
2989  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2990  {
2991  if (cells[v] == numbers::invalid_unsigned_int)
2992  continue;
2993 
2994  grad_weight[v] =
2995  data.data.front().shape_data_on_face
2996  [0][fe_degree + (integrate ? (2 - (face_nos[v] % 2)) :
2997  (1 + (face_nos[v] % 2)))][v];
2998 
2999  index_array_hermite[v] =
3000  &data.face_to_cell_index_hermite(face_nos[v], 0);
3001  }
3002  }
3003 
3004  // face_to_cell_index_nodal
3005  std::array<const unsigned int *, n_face_orientations> index_array_nodal =
3006  {};
3007 
3008  if (n_face_orientations == 1)
3009  index_array_nodal[0] =
3010  (data.data.front().nodal_at_cell_boundaries == true) ?
3011  &data.face_to_cell_index_nodal(face_nos[0], 0) :
3012  &dummy;
3013 
3014  if (n_face_orientations > 1 &&
3015  (data.data.front().nodal_at_cell_boundaries == true))
3016  {
3017  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3018  {
3019  if (cells[v] == numbers::invalid_unsigned_int)
3020  continue;
3021 
3022  index_array_nodal[v] =
3023  &data.face_to_cell_index_nodal(face_nos[v], 0);
3024  }
3025  }
3026 
3027  const auto reorientate = [&](const unsigned int v, const unsigned int i) {
3028  return (dim < 3 ||
3029  face_orientations[n_face_orientations == 1 ? 0 : v] == 0 ||
3030  subface_index < GeometryInfo<dim>::max_children_per_cell) ?
3031  i :
3032  orientation[v][i];
3033  };
3034 
3035  // this variable keeps track of whether we are able to directly write
3036  // the results into the result (function returns true) or not, requiring
3037  // an additional call to another function
3038  bool accesses_global_vector = true;
3039 
3040  for (unsigned int comp = 0; comp < n_components; ++comp)
3041  {
3042  if (integrate)
3043  proc.in_face_operation(temp1, comp);
3044 
3045  // we can only use the fast functions if we know the polynomial degree
3046  // as a template parameter (fe_degree != -1), and it only makes sense
3047  // to use the functions for at least linear functions for values on
3048  // the faces and quadratic functions for gradients on the faces, so
3049  // include the switch here
3050  if ((do_gradients == false &&
3051  data.data.front().nodal_at_cell_boundaries == true &&
3052  fe_degree > 0) ||
3053  (data.element_type ==
3055  fe_degree > 1))
3056  {
3057  // case 1: contiguous and interleaved indices
3058  if (n_face_orientations == 1 &&
3059  dof_info.index_storage_variants[dof_access_index][cell] ==
3061  interleaved_contiguous)
3062  {
3063  AssertDimension(n_face_orientations, 1);
3064 
3066  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
3067  VectorizedArrayType::size());
3068  Number2_ *vector_ptr =
3069  global_vector_ptr +
3070  dof_info.dof_indices_contiguous[dof_access_index]
3071  [cell *
3072  VectorizedArrayType::size()] +
3073  (dof_info
3074  .component_dof_indices_offset[active_fe_index]
3075  [first_selected_component] +
3076  comp * static_dofs_per_component) *
3077  VectorizedArrayType::size();
3078 
3079  if (fe_degree > 1 && do_gradients == true)
3080  {
3081  for (unsigned int i = 0; i < dofs_per_face; ++i)
3082  {
3083  if (n_face_orientations == 1)
3084  {
3085  const unsigned int ind1 =
3086  index_array_hermite[0][2 * i];
3087  const unsigned int ind2 =
3088  index_array_hermite[0][2 * i + 1];
3089  AssertIndexRange(ind1,
3090  data.dofs_per_component_on_cell);
3091  AssertIndexRange(ind2,
3092  data.dofs_per_component_on_cell);
3093  const unsigned int i_ = reorientate(0, i);
3094  proc.hermite_grad_vectorized(
3095  temp1[i_],
3096  temp1[i_ + dofs_per_face],
3097  vector_ptr + ind1 * VectorizedArrayType::size(),
3098  vector_ptr + ind2 * VectorizedArrayType::size(),
3099  grad_weight);
3100  }
3101  else
3102  {
3103  Assert(false, ExcNotImplemented());
3104  }
3105  }
3106  }
3107  else
3108  {
3109  for (unsigned int i = 0; i < dofs_per_face; ++i)
3110  {
3111  if (n_face_orientations == 1)
3112  {
3113  const unsigned int i_ = reorientate(0, i);
3114  const unsigned int ind = index_array_nodal[0][i];
3115  proc.value_vectorized(
3116  temp1[i_],
3117  vector_ptr + ind * VectorizedArrayType::size());
3118  }
3119  else
3120  {
3121  Assert(false, ExcNotImplemented());
3122  }
3123  }
3124  }
3125  }
3126 
3127  // case 2: contiguous and interleaved indices with fixed stride
3128  else if (n_face_orientations == 1 &&
3129  dof_info.index_storage_variants[dof_access_index][cell] ==
3131  interleaved_contiguous_strided)
3132  {
3133  AssertDimension(n_face_orientations, 1);
3134 
3136  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
3137  VectorizedArrayType::size());
3138  const unsigned int *indices =
3139  &dof_info.dof_indices_contiguous[dof_access_index]
3140  [cell *
3141  VectorizedArrayType::size()];
3142  Number2_ *vector_ptr =
3143  global_vector_ptr +
3144  (comp * static_dofs_per_component +
3145  dof_info
3146  .component_dof_indices_offset[active_fe_index]
3147  [first_selected_component]) *
3148  VectorizedArrayType::size();
3149  if (fe_degree > 1 && do_gradients == true)
3150  {
3151  for (unsigned int i = 0; i < dofs_per_face; ++i)
3152  {
3153  if (n_face_orientations == 1)
3154  {
3155  const unsigned int i_ = reorientate(0, i);
3156  const unsigned int ind1 =
3157  index_array_hermite[0][2 * i] *
3158  VectorizedArrayType::size();
3159  const unsigned int ind2 =
3160  index_array_hermite[0][2 * i + 1] *
3161  VectorizedArrayType::size();
3162  proc.hermite_grad_vectorized_indexed(
3163  temp1[i_],
3164  temp1[i_ + dofs_per_face],
3165  vector_ptr + ind1,
3166  vector_ptr + ind2,
3167  grad_weight,
3168  indices,
3169  indices);
3170  }
3171  else
3172  {
3173  Assert(false, ExcNotImplemented());
3174  }
3175  }
3176  }
3177  else
3178  {
3179  for (unsigned int i = 0; i < dofs_per_face; ++i)
3180  {
3181  if (n_face_orientations == 1)
3182  {
3183  const unsigned int i_ = reorientate(0, i);
3184  const unsigned int ind =
3185  index_array_nodal[0][i] *
3186  VectorizedArrayType::size();
3187  proc.value_vectorized_indexed(temp1[i_],
3188  vector_ptr + ind,
3189  indices);
3190  }
3191  else
3192  {
3193  Assert(false, ExcNotImplemented());
3194  }
3195  }
3196  }
3197  }
3198 
3199  // case 3: contiguous and interleaved indices with mixed stride
3200  else if (n_face_orientations == 1 &&
3201  dof_info.index_storage_variants[dof_access_index][cell] ==
3203  interleaved_contiguous_mixed_strides)
3204  {
3205  AssertDimension(n_face_orientations, 1);
3206 
3207  const unsigned int *strides =
3208  &dof_info.dof_indices_interleave_strides
3209  [dof_access_index][cell * VectorizedArrayType::size()];
3210  unsigned int indices[VectorizedArrayType::size()];
3211  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3212  indices[v] =
3213  dof_info.dof_indices_contiguous
3214  [dof_access_index]
3215  [cell * VectorizedArrayType::size() + v] +
3216  (dof_info
3217  .component_dof_indices_offset[active_fe_index]
3218  [first_selected_component] +
3219  comp * static_dofs_per_component) *
3220  strides[v];
3221  const unsigned int n_filled_lanes =
3222  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
3223 
3224  if (fe_degree > 1 && do_gradients == true)
3225  {
3226  if (n_filled_lanes == VectorizedArrayType::size())
3227  for (unsigned int i = 0; i < dofs_per_face; ++i)
3228  {
3229  if (n_face_orientations == 1)
3230  {
3231  const unsigned int i_ = reorientate(0, i);
3232  unsigned int ind1[VectorizedArrayType::size()];
3234  for (unsigned int v = 0;
3235  v < VectorizedArrayType::size();
3236  ++v)
3237  ind1[v] =
3238  indices[v] +
3239  index_array_hermite[0 /*TODO*/][2 * i] *
3240  strides[v];
3241  unsigned int ind2[VectorizedArrayType::size()];
3243  for (unsigned int v = 0;
3244  v < VectorizedArrayType::size();
3245  ++v)
3246  ind2[v] =
3247  indices[v] +
3248  index_array_hermite[0 /*TODO*/][2 * i + 1] *
3249  strides[v];
3250  proc.hermite_grad_vectorized_indexed(
3251  temp1[i_],
3252  temp1[i_ + dofs_per_face],
3253  global_vector_ptr,
3254  global_vector_ptr,
3255  grad_weight,
3256  ind1,
3257  ind2);
3258  }
3259  else
3260  {
3261  Assert(false, ExcNotImplemented());
3262  }
3263  }
3264  else
3265  {
3266  if (integrate == false)
3267  for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
3268  temp1[i] = VectorizedArrayType();
3269 
3270  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3271  for (unsigned int i = 0; i < dofs_per_face; ++i)
3272  {
3273  const unsigned int i_ =
3274  reorientate(n_face_orientations == 1 ? 0 : v,
3275  i);
3276  proc.hermite_grad(
3277  temp1[i_][v],
3278  temp1[i_ + dofs_per_face][v],
3279  global_vector_ptr
3280  [indices[v] +
3281  index_array_hermite
3282  [n_face_orientations == 1 ? 0 : v]
3283  [2 * i] *
3284  strides[v]],
3285  global_vector_ptr
3286  [indices[v] +
3287  index_array_hermite
3288  [n_face_orientations == 1 ? 0 : v]
3289  [2 * i + 1] *
3290  strides[v]],
3291  grad_weight[n_face_orientations == 1 ? 0 : v]);
3292  }
3293  }
3294  }
3295  else
3296  {
3297  if (n_filled_lanes == VectorizedArrayType::size())
3298  for (unsigned int i = 0; i < dofs_per_face; ++i)
3299  {
3300  if (n_face_orientations == 1)
3301  {
3302  unsigned int ind[VectorizedArrayType::size()];
3304  for (unsigned int v = 0;
3305  v < VectorizedArrayType::size();
3306  ++v)
3307  ind[v] = indices[v] +
3308  index_array_nodal[0][i] * strides[v];
3309  const unsigned int i_ = reorientate(0, i);
3310  proc.value_vectorized_indexed(temp1[i_],
3311  global_vector_ptr,
3312  ind);
3313  }
3314  else
3315  {
3316  Assert(false, ExcNotImplemented());
3317  }
3318  }
3319  else
3320  {
3321  if (integrate == false)
3322  for (unsigned int i = 0; i < dofs_per_face; ++i)
3323  temp1[i] = VectorizedArrayType();
3324 
3325  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3326  for (unsigned int i = 0; i < dofs_per_face; ++i)
3327  proc.value(
3328  temp1[reorientate(
3329  n_face_orientations == 1 ? 0 : v, i)][v],
3330  global_vector_ptr
3331  [indices[v] +
3332  index_array_nodal
3333  [n_face_orientations == 1 ? 0 : v][i] *
3334  strides[v]]);
3335  }
3336  }
3337  }
3338 
3339  // case 4: contiguous indices without interleaving
3340  else if (n_face_orientations > 1 ||
3341  dof_info.index_storage_variants[dof_access_index][cell] ==
3343  contiguous)
3344  {
3345  const unsigned int *indices =
3346  &dof_info.dof_indices_contiguous[dof_access_index]
3347  [cell *
3348  VectorizedArrayType::size()];
3349  Number2_ *vector_ptr =
3350  global_vector_ptr + comp * static_dofs_per_component +
3351  dof_info
3352  .component_dof_indices_offset[active_fe_index]
3353  [first_selected_component];
3354 
3355  const unsigned int n_filled_lanes =
3356  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
3357 
3358  const bool vectorization_possible =
3359  (n_face_orientations == 1) &&
3360  (n_filled_lanes == VectorizedArrayType::size()) &&
3361  (sm_ptr != nullptr);
3362 
3363  std::array<Number2_ *, VectorizedArrayType::size()>
3364  vector_ptrs = {};
3365 
3366  if (vectorization_possible == false)
3367  {
3368  if (n_face_orientations == 1)
3369  {
3370  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3371  if (sm_ptr == nullptr)
3372  {
3373  vector_ptrs[v] = vector_ptr + indices[v];
3374  }
3375  else
3376  {
3377  const auto &temp =
3378  dof_info.dof_indices_contiguous_sm
3379  [dof_access_index]
3380  [cell * VectorizedArrayType::size() + v];
3381  vector_ptrs[v] = const_cast<Number *>(
3382  sm_ptr->operator[](temp.first).data() +
3383  temp.second + comp * static_dofs_per_component +
3384  dof_info.component_dof_indices_offset
3385  [active_fe_index][first_selected_component]);
3386  }
3387  }
3388  else if (n_face_orientations == VectorizedArrayType::size())
3389  {
3390  for (unsigned int v = 0;
3391  v < VectorizedArrayType::size();
3392  ++v)
3393  if (cells[v] != numbers::invalid_unsigned_int)
3394  {
3395  if (sm_ptr == nullptr)
3396  {
3397  vector_ptrs[v] =
3398  vector_ptr +
3399  dof_info
3400  .dof_indices_contiguous[dof_access_index]
3401  [cells[v]];
3402  }
3403  else
3404  {
3405  const auto &temp =
3406  dof_info.dof_indices_contiguous_sm
3407  [dof_access_index][cells[v]];
3408  vector_ptrs[v] = const_cast<Number *>(
3409  sm_ptr->operator[](temp.first).data() +
3410  temp.second +
3411  comp * static_dofs_per_component +
3412  dof_info.component_dof_indices_offset
3413  [active_fe_index]
3414  [first_selected_component]);
3415  }
3416  }
3417  }
3418  else
3419  {
3420  Assert(false, ExcNotImplemented());
3421  }
3422  }
3423 
3424  if (do_gradients == true &&
3425  data.element_type ==
3427  {
3428  if (vectorization_possible)
3429  for (unsigned int i = 0; i < dofs_per_face; ++i)
3430  {
3431  const unsigned int ind1 =
3432  index_array_hermite[0][2 * i];
3433  const unsigned int ind2 =
3434  index_array_hermite[0][2 * i + 1];
3435  const unsigned int i_ = reorientate(0, i);
3436 
3437  proc.hermite_grad_vectorized_indexed(
3438  temp1[i_],
3439  temp1[i_ + dofs_per_face],
3440  vector_ptr + ind1,
3441  vector_ptr + ind2,
3442  grad_weight,
3443  indices,
3444  indices);
3445  }
3446  else if (n_face_orientations == 1)
3447  for (unsigned int i = 0; i < dofs_per_face; ++i)
3448  {
3449  const unsigned int ind1 =
3450  index_array_hermite[0][2 * i];
3451  const unsigned int ind2 =
3452  index_array_hermite[0][2 * i + 1];
3453  const unsigned int i_ = reorientate(0, i);
3454 
3455  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3456  proc.hermite_grad(temp1[i_][v],
3457  temp1[i_ + dofs_per_face][v],
3458  vector_ptrs[v][ind1],
3459  vector_ptrs[v][ind2],
3460  grad_weight[v]);
3461 
3462  if (integrate == false)
3463  for (unsigned int v = n_filled_lanes;
3464  v < VectorizedArrayType::size();
3465  ++v)
3466  {
3467  temp1[i_][v] = 0.0;
3468  temp1[i_ + dofs_per_face][v] = 0.0;
3469  }
3470  }
3471  else
3472  {
3473  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3474  for (unsigned int i = 0; i < dofs_per_face; ++i)
3475  proc.hermite_grad(
3476  temp1[reorientate(v, i)][v],
3477  temp1[reorientate(v, i) + dofs_per_face][v],
3478  vector_ptrs[v][index_array_hermite[v][2 * i]],
3479  vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
3480  grad_weight[v]);
3481  }
3482  }
3483  else
3484  {
3485  if (vectorization_possible)
3486  for (unsigned int i = 0; i < dofs_per_face; ++i)
3487  {
3488  const unsigned int ind = index_array_nodal[0][i];
3489  const unsigned int i_ = reorientate(0, i);
3490 
3491  proc.value_vectorized_indexed(temp1[i_],
3492  vector_ptr + ind,
3493  indices);
3494  }
3495  else if (n_face_orientations == 1)
3496  for (unsigned int i = 0; i < dofs_per_face; ++i)
3497  {
3498  const unsigned int ind = index_array_nodal[0][i];
3499  const unsigned int i_ = reorientate(0, i);
3500 
3501  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3502  proc.value(temp1[i_][v], vector_ptrs[v][ind]);
3503 
3504  if (integrate == false)
3505  for (unsigned int v = n_filled_lanes;
3506  v < VectorizedArrayType::size();
3507  ++v)
3508  temp1[i_][v] = 0.0;
3509  }
3510  else
3511  for (unsigned int i = 0; i < dofs_per_face; ++i)
3512  {
3513  for (unsigned int v = 0;
3514  v < VectorizedArrayType::size();
3515  ++v)
3516  if (cells[v] != numbers::invalid_unsigned_int)
3517  proc.value(
3518  temp1[reorientate(v, i)][v],
3519  vector_ptrs[v][index_array_nodal[v][i]]);
3520  }
3521  }
3522  }
3523  else
3524  {
3525  // case 5: default vector access
3526  // for the integrate_scatter path (integrate == true), we
3527  // need to only prepare the data in this function for all
3528  // components to later call distribute_local_to_global();
3529  // for the gather_evaluate path (integrate == false), we
3530  // instead want to leave early because we need to get the
3531  // vector data from somewhere else
3532  proc.default_operation(temp1, comp);
3533  if (integrate)
3534  accesses_global_vector = false;
3535  else
3536  return false;
3537  }
3538  }
3539  else
3540  {
3541  // case 5: default vector access
3542  proc.default_operation(temp1, comp);
3543  if (integrate)
3544  accesses_global_vector = false;
3545  else
3546  return false;
3547  }
3548 
3549  if (!integrate)
3550  proc.in_face_operation(temp1, comp);
3551  }
3552 
3553  if (!integrate &&
3554  (face_orientations[0] > 0 &&
3556  {
3557  AssertDimension(n_face_orientations, 1);
3559  n_components,
3560  face_orientations[0],
3561  orientation_map,
3562  false,
3563  do_values,
3564  do_gradients,
3565  data.n_q_points_face,
3566  scratch_data,
3567  values_quad,
3568  gradients_quad);
3569  }
3570 
3571  return accesses_global_vector;
3572  }
3573 
3574 
3575 
3576  template <int dim,
3577  typename Number,
3578  typename VectorizedArrayType,
3579  typename Number2 = Number>
3581  {
3582  template <int fe_degree, int n_q_points_1d>
3583  static bool
3584  run(const unsigned int n_components,
3585  const unsigned int n_face_orientations,
3586  const Number2 * src_ptr,
3587  const std::vector<ArrayView<const Number>> *sm_ptr,
3589  const MatrixFreeFunctions::DoFInfo & dof_info,
3590  VectorizedArrayType * values_quad,
3591  VectorizedArrayType *gradients_quad,
3592  VectorizedArrayType *scratch_data,
3593  const bool evaluate_values,
3594  const bool evaluate_gradients,
3595  const unsigned int active_fe_index,
3596  const unsigned int first_selected_component,
3597  const std::array<unsigned int, VectorizedArrayType::size()> cells,
3598  const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
3599  const unsigned int subface_index,
3600  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
3601  const std::array<unsigned int, VectorizedArrayType::size()>
3602  face_orientations,
3603  const Table<2, unsigned int> &orientation_map)
3604  {
3605  if (src_ptr == nullptr)
3606  return false;
3607 
3609  return false;
3610 
3611  (void)sm_ptr;
3612 
3613  Processor<fe_degree, n_q_points_1d> p(n_components,
3614  false,
3615  src_ptr,
3616  sm_ptr,
3617  data,
3618  dof_info,
3619  values_quad,
3620  gradients_quad,
3621  scratch_data,
3622  evaluate_values,
3623  evaluate_gradients,
3624  active_fe_index,
3625  first_selected_component,
3626  cells,
3627  face_nos,
3628  subface_index,
3629  dof_access_index,
3630  face_orientations,
3631  orientation_map);
3632 
3633  if (n_face_orientations == VectorizedArrayType::size())
3634  return fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
3635  p);
3636  else
3637  return fe_face_evaluation_process_and_io<1>(p);
3638  }
3639 
3640  private:
3641  template <int fe_degree, int n_q_points_1d>
3642  struct Processor
3643  {
3644  static const int dim_ = dim;
3645  static const int fe_degree_ = fe_degree;
3646  static const int n_q_points_1d_ = n_q_points_1d;
3647  using VectorizedArrayType_ = VectorizedArrayType;
3648  using Number_ = Number;
3649  using Number2_ = const Number2;
3650 
3652  const unsigned int n_components,
3653  const bool integrate,
3654  const Number2 * global_vector_ptr,
3655  const std::vector<ArrayView<const Number>> *sm_ptr,
3657  const MatrixFreeFunctions::DoFInfo & dof_info,
3658  VectorizedArrayType * values_quad,
3659  VectorizedArrayType *gradients_quad,
3660  VectorizedArrayType *scratch_data,
3661  const bool do_values,
3662  const bool do_gradients,
3663  const unsigned int active_fe_index,
3664  const unsigned int first_selected_component,
3665  const std::array<unsigned int, VectorizedArrayType::size()> cells,
3666  const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
3667  const unsigned int subface_index,
3668  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
3669  const std::array<unsigned int, VectorizedArrayType::size()>
3670  face_orientations,
3671  const Table<2, unsigned int> &orientation_map)
3672  : n_components(n_components)
3673  , integrate(integrate)
3674  , global_vector_ptr(global_vector_ptr)
3675  , sm_ptr(sm_ptr)
3676  , data(data)
3677  , dof_info(dof_info)
3678  , values_quad(values_quad)
3679  , gradients_quad(gradients_quad)
3680  , scratch_data(scratch_data)
3681  , do_values(do_values)
3682  , do_gradients(do_gradients)
3683  , active_fe_index(active_fe_index)
3684  , first_selected_component(first_selected_component)
3685  , cells(cells)
3686  , face_nos(face_nos)
3687  , subface_index(subface_index)
3688  , dof_access_index(dof_access_index)
3689  , face_orientations(face_orientations)
3690  , orientation_map(orientation_map)
3691  {}
3692 
3693  template <typename T0, typename T1, typename T2>
3694  void
3696  T0 & temp_2,
3697  const T1 src_ptr_1,
3698  const T1 src_ptr_2,
3699  const T2 &grad_weight)
3700  {
3701  do_vectorized_read(src_ptr_1, temp_1);
3702  do_vectorized_read(src_ptr_2, temp_2);
3703  temp_2 = grad_weight * (temp_1 - temp_2);
3704  }
3705 
3706  template <typename T1, typename T2>
3707  void
3708  value_vectorized(T1 &temp, const T2 src_ptr)
3709  {
3710  do_vectorized_read(src_ptr, temp);
3711  }
3712 
3713  template <typename T0, typename T1, typename T2, typename T3>
3714  void
3716  T0 & temp_2,
3717  const T1 src_ptr_1,
3718  const T1 src_ptr_2,
3719  const T2 &grad_weight,
3720  const T3 &indices_1,
3721  const T3 &indices_2)
3722  {
3723  do_vectorized_gather(src_ptr_1, indices_1, temp_1);
3724  do_vectorized_gather(src_ptr_2, indices_2, temp_2);
3725  temp_2 = grad_weight * (temp_1 - temp_2);
3726  }
3727 
3728  template <typename T0, typename T1, typename T2>
3729  void
3730  value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
3731  {
3732  do_vectorized_gather(src_ptr, indices, temp);
3733  }
3734 
3735  template <typename T0, typename T1, typename T2>
3736  void
3737  hermite_grad(T0 & temp_1,
3738  T0 & temp_2,
3739  const T1 &src_ptr_1,
3740  const T2 &src_ptr_2,
3741  const T2 &grad_weight)
3742  {
3743  // case 3a)
3744  temp_1 = src_ptr_1;
3745  temp_2 = grad_weight * (temp_1 - src_ptr_2);
3746  }
3747 
3748  template <typename T1, typename T2>
3749  void
3750  value(T1 &temp, const T2 &src_ptr)
3751  {
3752  // case 3b)
3753  temp = src_ptr;
3754  }
3755 
3756  template <typename T1>
3757  void
3758  default_operation(const T1 &, const unsigned int)
3759  {
3760  // case 5)
3761  }
3762 
3763  template <typename T1>
3764  void
3765  in_face_operation(T1 &temp1, const unsigned int comp)
3766  {
3767  const unsigned int dofs_per_face =
3768  fe_degree > -1 ?
3769  Utilities::pow(fe_degree + 1, dim - 1) :
3770  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
3771  const unsigned int n_q_points =
3772  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
3773  data.n_q_points_face;
3774  if (fe_degree > -1 &&
3775  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
3777  FEFaceEvaluationImpl<true,
3778  dim,
3779  fe_degree,
3780  n_q_points_1d,
3781  VectorizedArrayType>::
3782  evaluate_in_face(/* n_components */ 1,
3783  data,
3784  temp1,
3785  values_quad + comp * n_q_points,
3786  gradients_quad + comp * dim * n_q_points,
3787  scratch_data + 2 * dofs_per_face,
3788  do_values,
3789  do_gradients,
3790  subface_index);
3791  else
3792  FEFaceEvaluationImpl<false,
3793  dim,
3794  fe_degree,
3795  n_q_points_1d,
3796  VectorizedArrayType>::
3797  evaluate_in_face(/* n_components */ 1,
3798  data,
3799  temp1,
3800  values_quad + comp * n_q_points,
3801  gradients_quad + comp * dim * n_q_points,
3802  scratch_data + 2 * dofs_per_face,
3803  do_values,
3804  do_gradients,
3805  subface_index);
3806  }
3807 
3808  const unsigned int n_components;
3809  const bool integrate;
3810  const Number2 * global_vector_ptr;
3811  const std::vector<ArrayView<const Number>> *sm_ptr;
3814  VectorizedArrayType * values_quad;
3815  VectorizedArrayType * gradients_quad;
3816  VectorizedArrayType * scratch_data;
3817  const bool do_values;
3818  const bool do_gradients;
3819  const unsigned int active_fe_index;
3820  const unsigned int first_selected_component;
3821  const std::array<unsigned int, VectorizedArrayType::size()> cells;
3822  const std::array<unsigned int, VectorizedArrayType::size()> face_nos;
3823  const unsigned int subface_index;
3825  const std::array<unsigned int, VectorizedArrayType::size()>
3828  };
3829  };
3830 
3831  template <int dim,
3832  typename Number,
3833  typename VectorizedArrayType,
3834  typename Number2 = Number>
3836  {
3837  template <int fe_degree, int n_q_points_1d>
3838  static bool
3839  run(const unsigned int n_components,
3840  const unsigned int n_face_orientations,
3841  Number2 * dst_ptr,
3842  const std::vector<ArrayView<const Number2>> *sm_ptr,
3844  const MatrixFreeFunctions::DoFInfo & dof_info,
3845  VectorizedArrayType * values_array,
3846  VectorizedArrayType * values_quad,
3847  VectorizedArrayType *gradients_quad,
3848  VectorizedArrayType *scratch_data,
3849  const bool integrate_values,
3850  const bool integrate_gradients,
3851  const unsigned int active_fe_index,
3852  const unsigned int first_selected_component,
3853  const std::array<unsigned int, VectorizedArrayType::size()> cells,
3854  const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
3855  const unsigned int subface_index,
3856  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
3857  const std::array<unsigned int, VectorizedArrayType::size()>
3858  face_orientations,
3859  const Table<2, unsigned int> &orientation_map)
3860  {
3861  (void)sm_ptr;
3862 
3863  if (dst_ptr == nullptr ||
3865  {
3866  AssertDimension(n_face_orientations, 1);
3867 
3868  // for block vectors simply integrate
3870  template run<fe_degree, n_q_points_1d>(n_components,
3871  data,
3872  values_array,
3873  values_quad,
3874  gradients_quad,
3875  scratch_data,
3876  integrate_values,
3877  integrate_gradients,
3878  face_nos[0],
3879  subface_index,
3880  face_orientations[0],
3881  orientation_map);
3882 
3883  // default vector access
3884  return false;
3885  }
3886 
3887 
3888  Processor<fe_degree, n_q_points_1d> p(values_array,
3889  n_components,
3890  true,
3891  dst_ptr,
3892  sm_ptr,
3893  data,
3894  dof_info,
3895  values_quad,
3896  gradients_quad,
3897  scratch_data,
3898  integrate_values,
3899  integrate_gradients,
3900  active_fe_index,
3901  first_selected_component,
3902  cells,
3903  face_nos,
3904  subface_index,
3905  dof_access_index,
3906  face_orientations,
3907  orientation_map);
3908 
3909  if (n_face_orientations == VectorizedArrayType::size())
3910  return fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
3911  p);
3912  else
3913  return fe_face_evaluation_process_and_io<1>(p);
3914  }
3915 
3916  private:
3917  template <int fe_degree, int n_q_points_1d>
3918  struct Processor
3919  {
3920  static const int dim_ = dim;
3921  static const int fe_degree_ = fe_degree;
3922  static const int n_q_points_1d_ = n_q_points_1d;
3923  using VectorizedArrayType_ = VectorizedArrayType;
3924  using Number_ = Number;
3925  using Number2_ = Number2;
3926 
3927 
3929  VectorizedArrayType * values_array,
3930  const unsigned int n_components,
3931  const bool integrate,
3932  Number2 * global_vector_ptr,
3933  const std::vector<ArrayView<const Number>> *sm_ptr,
3935  const MatrixFreeFunctions::DoFInfo & dof_info,
3936  VectorizedArrayType * values_quad,
3937  VectorizedArrayType *gradients_quad,
3938  VectorizedArrayType *scratch_data,
3939  const bool do_values,
3940  const bool do_gradients,
3941  const unsigned int active_fe_index,
3942  const unsigned int first_selected_component,
3943  const std::array<unsigned int, VectorizedArrayType::size()> cells,
3944  const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
3945  const unsigned int subface_index,
3946  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
3947  const std::array<unsigned int, VectorizedArrayType::size()>
3948  face_orientations,
3949  const Table<2, unsigned int> &orientation_map)
3950  : values_array(values_array)
3951  , n_components(n_components)
3952  , integrate(integrate)
3953  , global_vector_ptr(global_vector_ptr)
3954  , sm_ptr(sm_ptr)
3955  , data(data)
3956  , dof_info(dof_info)
3957  , values_quad(values_quad)
3958  , gradients_quad(gradients_quad)
3959  , scratch_data(scratch_data)
3960  , do_values(do_values)
3961  , do_gradients(do_gradients)
3962  , active_fe_index(active_fe_index)
3963  , first_selected_component(first_selected_component)
3964  , cells(cells)
3965  , face_nos(face_nos)
3966  , subface_index(subface_index)
3967  , dof_access_index(dof_access_index)
3968  , face_orientations(face_orientations)
3969  , orientation_map(orientation_map)
3970  {}
3971 
3972  template <typename T0, typename T1, typename T2, typename T3, typename T4>
3973  void
3974  hermite_grad_vectorized(const T0 &temp_1,
3975  const T1 &temp_2,
3976  T2 dst_ptr_1,
3977  T3 dst_ptr_2,
3978  const T4 &grad_weight)
3979  {
3980  // case 1a)
3981  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
3982  const VectorizedArrayType grad = grad_weight * temp_2;
3983  do_vectorized_add(val, dst_ptr_1);
3984  do_vectorized_add(grad, dst_ptr_2);
3985  }
3986 
3987  template <typename T0, typename T1>
3988  void
3989  value_vectorized(const T0 &temp, T1 dst_ptr)
3990  {
3991  // case 1b)
3992  do_vectorized_add(temp, dst_ptr);
3993  }
3994 
3995  template <typename T0, typename T1, typename T2, typename T3>
3996  void
3998  const T0 &temp_2,
3999  T1 dst_ptr_1,
4000  T1 dst_ptr_2,
4001  const T2 &grad_weight,
4002  const T3 &indices_1,
4003  const T3 &indices_2)
4004  {
4005  // case 2a)
4006  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
4007  const VectorizedArrayType grad = grad_weight * temp_2;
4008  do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
4009  do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
4010  }
4011 
4012  template <typename T0, typename T1, typename T2>
4013  void
4014  value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
4015  {
4016  // case 2b)
4017  do_vectorized_scatter_add(temp, indices, dst_ptr);
4018  }
4019 
4020  template <typename T0, typename T1, typename T2>
4021  void
4022  hermite_grad(const T0 &temp_1,
4023  const T0 &temp_2,
4024  T1 & dst_ptr_1,
4025  T1 & dst_ptr_2,
4026  const T2 &grad_weight)
4027  {
4028  // case 3a)
4029  const Number val = temp_1 - grad_weight * temp_2;
4030  const Number grad = grad_weight * temp_2;
4031  dst_ptr_1 += val;
4032  dst_ptr_2 += grad;
4033  }
4034 
4035  template <typename T0, typename T1>
4036  void
4037  value(const T0 &temp, T1 &dst_ptr)
4038  {
4039  // case 3b)
4040  dst_ptr += temp;
4041  }
4042 
4043  template <typename T0>
4044  void
4045  default_operation(const T0 &temp1, const unsigned int comp)
4046  {
4047  // case 5: default vector access, must be handled separately, just do
4048  // the face-normal interpolation
4049 
4051  template interpolate<false, false>(
4052  /* n_components */ 1,
4053  data,
4054  temp1,
4055  values_array + comp * data.dofs_per_component_on_cell,
4056  do_gradients,
4057  face_nos[0]);
4058  }
4059 
4060  template <typename T0>
4061  void
4062  in_face_operation(T0 &temp1, const unsigned int comp)
4063  {
4064  const unsigned int dofs_per_face =
4065  fe_degree > -1 ?
4066  Utilities::pow(fe_degree + 1, dim - 1) :
4067  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
4068  const unsigned int n_q_points =
4069  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
4070  data.n_q_points_face;
4071  if (fe_degree > -1 &&
4072  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
4073  data.element_type <=
4076  dim,
4077  fe_degree,
4078  n_q_points_1d,
4079  VectorizedArrayType>::
4080  integrate_in_face(/* n_components */ 1,
4081  data,
4082  temp1,
4083  values_quad + comp * n_q_points,
4084  gradients_quad + dim * comp * n_q_points,
4085  scratch_data + 2 * dofs_per_face,
4086  do_values,
4087  do_gradients,
4088  subface_index);
4089  else
4091  dim,
4092  fe_degree,
4093  n_q_points_1d,
4094  VectorizedArrayType>::
4095  integrate_in_face(/* n_components */ 1,
4096  data,
4097  temp1,
4098  values_quad + comp * n_q_points,
4099  gradients_quad + dim * comp * n_q_points,
4100  scratch_data + 2 * dofs_per_face,
4101  do_values,
4102  do_gradients,
4103  subface_index);
4104  }
4105 
4106  VectorizedArrayType *values_array;
4107 
4108 
4109  const unsigned int n_components;
4110  const bool integrate;
4112  const std::vector<ArrayView<const Number>> *sm_ptr;
4115  VectorizedArrayType * values_quad;
4116  VectorizedArrayType * gradients_quad;
4117  VectorizedArrayType * scratch_data;
4118  const bool do_values;
4119  const bool do_gradients;
4120  const unsigned int active_fe_index;
4121  const unsigned int first_selected_component;
4122  const std::array<unsigned int, VectorizedArrayType::size()> cells;
4123  const std::array<unsigned int, VectorizedArrayType::size()> face_nos;
4124  const unsigned int subface_index;
4126  const std::array<unsigned int, VectorizedArrayType::size()>
4129  };
4130  };
4131 
4132 
4133 
4138  template <int dim, typename Number>
4140  {
4141  template <int fe_degree, int = 0>
4142  static bool
4143  run(const unsigned int n_components,
4144  const FEEvaluationBaseData<dim,
4145  typename Number::value_type,
4146  false,
4147  Number> &fe_eval,
4148  const Number * in_array,
4149  Number * out_array,
4150  typename std::enable_if<fe_degree != -1>::type * = nullptr)
4151  {
4152  constexpr unsigned int dofs_per_component =
4153  Utilities::pow(fe_degree + 1, dim);
4154 
4155  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
4156  Assert(fe_eval.get_shape_info().element_type <=
4158  ExcNotImplemented());
4159 
4161  dim,
4162  fe_degree + 1,
4163  fe_degree + 1,
4164  Number>
4165  evaluator(
4168  fe_eval.get_shape_info().data.front().inverse_shape_values_eo);
4169 
4170  for (unsigned int d = 0; d < n_components; ++d)
4171  {
4172  const Number *in = in_array + d * dofs_per_component;
4173  Number * out = out_array + d * dofs_per_component;
4174  // Need to select 'apply' method with hessian slot because values
4175  // assume symmetries that do not exist in the inverse shapes
4176  evaluator.template hessians<0, true, false>(in, out);
4177  if (dim > 1)
4178  evaluator.template hessians<1, true, false>(out, out);
4179  if (dim > 2)
4180  evaluator.template hessians<2, true, false>(out, out);
4181  }
4182  for (unsigned int q = 0; q < dofs_per_component; ++q)
4183  {
4184  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
4185  for (unsigned int d = 0; d < n_components; ++d)
4186  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
4187  }
4188  for (unsigned int d = 0; d < n_components; ++d)
4189  {
4190  Number *out = out_array + d * dofs_per_component;
4191  if (dim > 2)
4192  evaluator.template hessians<2, false, false>(out, out);
4193  if (dim > 1)
4194  evaluator.template hessians<1, false, false>(out, out);
4195  evaluator.template hessians<0, false, false>(out, out);
4196  }
4197  return false;
4198  }
4199 
4200  template <int fe_degree, int = 0>
4201  static bool
4202  run(const unsigned int n_components,
4203  const FEEvaluationBaseData<dim,
4204  typename Number::value_type,
4205  false,
4206  Number> &fe_eval,
4207  const Number * in_array,
4208  Number * out_array,
4209  typename std::enable_if<fe_degree == -1>::type * = nullptr)
4210  {
4211  static_assert(fe_degree == -1, "Only usable for degree -1");
4212  const unsigned int dofs_per_component =
4213  fe_eval.get_shape_info().dofs_per_component_on_cell;
4214 
4215  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
4216 
4217  internal::
4218  EvaluatorTensorProduct<internal::evaluate_general, dim, 0, 0, Number>
4219  evaluator(fe_eval.get_shape_info().data.front().inverse_shape_values,
4222  fe_eval.get_shape_info().data.front().fe_degree + 1,
4223  fe_eval.get_shape_info().data.front().fe_degree + 1);
4224 
4225  for (unsigned int d = 0; d < n_components; ++d)
4226  {
4227  const Number *in = in_array + d * dofs_per_component;
4228  Number * out = out_array + d * dofs_per_component;
4229  // Need to select 'apply' method with hessian slot because values
4230  // assume symmetries that do not exist in the inverse shapes
4231  evaluator.template values<0, true, false>(in, out);
4232  if (dim > 1)
4233  evaluator.template values<1, true, false>(out, out);
4234  if (dim > 2)
4235  evaluator.template values<2, true, false>(out, out);
4236  }
4237  for (unsigned int q = 0; q < dofs_per_component; ++q)
4238  {
4239  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
4240  for (unsigned int d = 0; d < n_components; ++d)
4241  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
4242  }
4243  for (unsigned int d = 0; d < n_components; ++d)
4244  {
4245  Number *out = out_array + d * dofs_per_component;
4246  if (dim > 2)
4247  evaluator.template values<2, false, false>(out, out);
4248  if (dim > 1)
4249  evaluator.template values<1, false, false>(out, out);
4250  evaluator.template values<0, false, false>(out, out);
4251  }
4252  return false;
4253  }
4254  };
4255 
4256 
4257 
4262  template <int dim, typename Number>
4264  {
4265  template <int fe_degree, int = 0>
4266  static bool
4267  run(const unsigned int n_desired_components,
4268  const AlignedVector<Number> &inverse_shape,
4269  const AlignedVector<Number> &inverse_coefficients,
4270  const Number * in_array,
4271  Number * out_array,
4272  typename std::enable_if<fe_degree != -1>::type * = nullptr)
4273  {
4274  constexpr unsigned int dofs_per_component =
4275  Utilities::pow(fe_degree + 1, dim);
4276  Assert(inverse_coefficients.size() > 0 &&
4277  inverse_coefficients.size() % dofs_per_component == 0,
4278  ExcMessage(
4279  "Expected diagonal to be a multiple of scalar dof per cells"));
4280  if (inverse_coefficients.size() != dofs_per_component)
4281  AssertDimension(n_desired_components * dofs_per_component,
4282  inverse_coefficients.size());
4283 
4284  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
4285 
4287  dim,
4288  fe_degree + 1,
4289  fe_degree + 1,
4290  Number>
4291  evaluator(AlignedVector<Number>(),
4293  inverse_shape);
4294 
4295  const unsigned int shift_coefficient =
4296  inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
4297  0;
4298  const Number *inv_coefficient = inverse_coefficients.data();
4299  for (unsigned int d = 0; d < n_desired_components; ++d)
4300  {
4301  const Number *in = in_array + d * dofs_per_component;
4302  Number * out = out_array + d * dofs_per_component;
4303  // Need to select 'apply' method with hessian slot because values
4304  // assume symmetries that do not exist in the inverse shapes
4305  evaluator.template hessians<0, true, false>(in, out);
4306  if (dim > 1)
4307  evaluator.template hessians<1, true, false>(out, out);
4308  if (dim > 2)
4309  evaluator.template hessians<2, true, false>(out, out);
4310 
4311  for (unsigned int q = 0; q < dofs_per_component; ++q)
4312  out[q] *= inv_coefficient[q];
4313 
4314  if (dim > 2)
4315  evaluator.template hessians<2, false, false>(out, out);
4316  if (dim > 1)
4317  evaluator.template hessians<1, false, false>(out, out);
4318  evaluator.template hessians<0, false, false>(out, out);
4319 
4320  inv_coefficient += shift_coefficient;
4321  }
4322  return false;
4323  }
4324 
4328  template <int fe_degree, int = 0>
4329  static bool
4330  run(const unsigned int,
4331  const AlignedVector<Number> &,
4332  const AlignedVector<Number> &,
4333  const Number *,
4334  Number *,
4335  typename std::enable_if<fe_degree == -1>::type * = nullptr)
4336  {
4337  static_assert(fe_degree == -1, "Only usable for degree -1");
4338  Assert(false, ExcNotImplemented());
4339  return false;
4340  }
4341  };
4342 
4343 
4344 
4349  template <int dim, typename Number>
4351  {
4352  template <int fe_degree, int = 0>
4353  static bool
4354  run(const unsigned int n_desired_components,
4355  const AlignedVector<Number> &inverse_shape,
4356  const Number * in_array,
4357  Number * out_array,
4358  typename std::enable_if<fe_degree != -1>::type * = nullptr)
4359  {
4360  constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim);
4362  dim,
4363  fe_degree + 1,
4364  fe_degree + 1,
4365  Number>
4366  evaluator(AlignedVector<Number>(),
4368  inverse_shape);
4369 
4370  for (unsigned int d = 0; d < n_desired_components; ++d)
4371  {
4372  const Number *in = in_array + d * dofs_per_cell;
4373  Number * out = out_array + d * dofs_per_cell;
4374 
4375  if (dim == 3)
4376  {
4377  evaluator.template hessians<2, false, false>(in, out);
4378  evaluator.template hessians<1, false, false>(out, out);
4379  evaluator.template hessians<0, false, false>(out, out);
4380  }
4381  if (dim == 2)
4382  {
4383  evaluator.template hessians<1, false, false>(in, out);
4384  evaluator.template hessians<0, false, false>(out, out);
4385  }
4386  if (dim == 1)
4387  evaluator.template hessians<0, false, false>(in, out);
4388  }
4389  return false;
4390  }
4391 
4392  template <int fe_degree, int = 0>
4393  static bool
4394  run(const unsigned int,
4395  const AlignedVector<Number> &,
4396  const Number *,
4397  Number *,
4398  typename std::enable_if<fe_degree == -1>::type * = nullptr)
4399  {
4400  static_assert(fe_degree == -1, "Only usable for degree -1");
4401  Assert(false, ExcNotImplemented());
4402  return false;
4403  }
4404  };
4405 
4406 } // end of namespace internal
4407 
4408 
4410 
4411 #endif
static void do_backward(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const bool add_into_result, Number *values_in, Number *values_out, const unsigned int basis_size_1_variable=numbers::invalid_unsigned_int, const unsigned int basis_size_2_variable=numbers::invalid_unsigned_int)
void value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs_actual, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool add_into_values_array)
static bool run(const unsigned int n_components, const unsigned int n_face_orientations, Number2 *dst_ptr, const std::vector< ArrayView< const Number2 >> *sm_ptr, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, const MatrixFreeFunctions::DoFInfo &dof_info, VectorizedArrayType *values_array, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool integrate_values, const bool integrate_gradients, const unsigned int active_fe_index, const unsigned int first_selected_component, const std::array< unsigned int, VectorizedArrayType::size()> cells, const std::array< unsigned int, VectorizedArrayType::size()> face_nos, const unsigned int subface_index, const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, const std::array< unsigned int, VectorizedArrayType::size()> face_orientations, const Table< 2, unsigned int > &orientation_map)
void hermite_grad_vectorized(T0 &temp_1, T0 &temp_2, const T1 src_ptr_1, const T1 src_ptr_2, const T2 &grad_weight)
const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const internal::MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs_actual, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool sum_into_values_array)
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
void adjust_for_face_orientation(const unsigned int dim, const unsigned int n_components, const unsigned int face_orientation, const Table< 2, unsigned int > &orientation_map, const bool integrate, const bool values, const bool gradients, const unsigned int n_q_points, Number *tmp_values, Number *values_quad, Number *gradients_quad)
void hermite_grad(const T0 &temp_1, const T0 &temp_2, T1 &dst_ptr_1, T1 &dst_ptr_2, const T2 &grad_weight)
void do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
static void interpolate_quadrature(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< Number > &data, const Number *input, Number *output, const bool do_gradients, const unsigned int face_no)
static void integrate_in_face(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< Number > &data, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool integrate_val, const bool integrate_grad, const unsigned int subface_index)
void do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
static void interpolate_generic(const unsigned int n_components, const Number *input, Number *output, const bool do_gradients, const unsigned int face_no, const unsigned int n_points_1d, const std::array< AlignedVector< Number >, 2 > &shape_data, const unsigned int dofs_per_component_on_cell, const unsigned int dofs_per_component_on_face)
void hermite_grad(T0 &temp_1, T0 &temp_2, const T1 &src_ptr_1, const T2 &src_ptr_2, const T2 &grad_weight)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data)
pointer data()
void hermite_grad_vectorized(const T0 &temp_1, const T1 &temp_2, T2 dst_ptr_1, T3 dst_ptr_2, const T4 &grad_weight)
static void evaluate_in_face(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< Number > &data, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool evaluate_val, const bool evaluate_grad, const unsigned int subface_index)
static bool run(const unsigned int n_components, const FEEvaluationBaseData< dim, typename Number::value_type, false, Number > &fe_eval, const Number *in_array, Number *out_array, typename std::enable_if< fe_degree !=-1 >::type *=nullptr)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
static bool run(const unsigned int, const AlignedVector< Number > &, const Number *, Number *, typename std::enable_if< fe_degree==-1 >::type *=nullptr)
std::vector< unsigned int > n_q_points_faces
Definition: shape_info.h:427
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:461
void hermite_grad_vectorized_indexed(const T0 &temp_1, const T0 &temp_2, T1 dst_ptr_1, T1 dst_ptr_2, const T2 &grad_weight, const T3 &indices_1, const T3 &indices_2)
const std::array< unsigned int, VectorizedArrayType::size()> face_orientations
T fixed_power(const T t)
Definition: utilities.h:1081
const std::vector< ArrayView< const Number > > * sm_ptr
const std::array< unsigned int, VectorizedArrayType::size()> cells
static bool run(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, const VectorizedArrayType *values_array, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool evaluate_values, const bool evaluate_gradients, const unsigned int face_no, const unsigned int subface_index, const unsigned int face_orientation, const Table< 2, unsigned int > &orientation_map)
static bool run(const unsigned int n_desired_components, const AlignedVector< Number > &inverse_shape, const AlignedVector< Number > &inverse_coefficients, const Number *in_array, Number *out_array, typename std::enable_if< fe_degree !=-1 >::type *=nullptr)
static ::ExceptionBase & ExcMessage(std::string arg1)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const Number *values_dofs_actual, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data)
void gather(const Number *base_ptr, const unsigned int *offsets)
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static bool run(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, VectorizedArrayType *values_array, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool integrate_values, const bool integrate_gradients, const unsigned int face_no, const unsigned int subface_index, const unsigned int face_orientation, const Table< 2, unsigned int > &orientation_map)
static bool run(const unsigned int n_components, const FEEvaluationBaseData< dim, typename Number::value_type, false, Number > &fe_eval, const Number *in_array, Number *out_array, typename std::enable_if< fe_degree==-1 >::type *=nullptr)
static void do_mass(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const AlignedVector< Number > &coefficients, const Number *values_in, Number *scratch_data, Number *values_out)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
void load(const Number *ptr)
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:95
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & data
static void do_forward(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const Number *values_in, Number *values_out, const unsigned int basis_size_1_variable=numbers::invalid_unsigned_int, const unsigned int basis_size_2_variable=numbers::invalid_unsigned_int)
static bool run(const unsigned int n_desired_components, const AlignedVector< Number > &inverse_shape, const Number *in_array, Number *out_array, typename std::enable_if< fe_degree !=-1 >::type *=nullptr)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void hermite_grad_vectorized_indexed(T0 &temp_1, T0 &temp_2, const T1 src_ptr_1, const T1 src_ptr_2, const T2 &grad_weight, const T3 &indices_1, const T3 &indices_2)
const std::array< unsigned int, VectorizedArrayType::size()> face_orientations
static void interpolate(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< Number > &data, const Number *input, Number *output, const bool do_gradients, const unsigned int face_no)
const std::vector< ArrayView< const Number > > * sm_ptr
void default_operation(const T0 &temp1, const unsigned int comp)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
const std::array< unsigned int, VectorizedArrayType::size()> face_nos
static bool run(const unsigned int, const AlignedVector< Number > &, const AlignedVector< Number > &, const Number *, Number *, typename std::enable_if< fe_degree==-1 >::type *=nullptr)
size_type size() const
const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & data
static ::ExceptionBase & ExcNotImplemented()
Processor(const unsigned int n_components, const bool integrate, const Number2 *global_vector_ptr, const std::vector< ArrayView< const Number >> *sm_ptr, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, const MatrixFreeFunctions::DoFInfo &dof_info, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool do_values, const bool do_gradients, const unsigned int active_fe_index, const unsigned int first_selected_component, const std::array< unsigned int, VectorizedArrayType::size()> cells, const std::array< unsigned int, VectorizedArrayType::size()> face_nos, const unsigned int subface_index, const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, const std::array< unsigned int, VectorizedArrayType::size()> face_orientations, const Table< 2, unsigned int > &orientation_map)
std::vector< UnivariateShapeData< Number > > data
Definition: shape_info.h:386
const std::array< unsigned int, VectorizedArrayType::size()> cells
void value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
static bool run(const unsigned int n_components, const unsigned int n_face_orientations, const Number2 *src_ptr, const std::vector< ArrayView< const Number >> *sm_ptr, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, const MatrixFreeFunctions::DoFInfo &dof_info, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool evaluate_values, const bool evaluate_gradients, const unsigned int active_fe_index, const unsigned int first_selected_component, const std::array< unsigned int, VectorizedArrayType::size()> cells, const std::array< unsigned int, VectorizedArrayType::size()> face_nos, const unsigned int subface_index, const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, const std::array< unsigned int, VectorizedArrayType::size()> face_orientations, const Table< 2, unsigned int > &orientation_map)
static bool fe_face_evaluation_process_and_io(Processor &proc)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const internal::MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs_actual, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data)
EvaluationFlags
The EvaluationFlags enum.
const std::array< unsigned int, VectorizedArrayType::size()> face_nos
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:135
const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool add_into_values_array)
Processor(VectorizedArrayType *values_array, const unsigned int n_components, const bool integrate, Number2 *global_vector_ptr, const std::vector< ArrayView< const Number >> *sm_ptr, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, const MatrixFreeFunctions::DoFInfo &dof_info, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool do_values, const bool do_gradients, const unsigned int active_fe_index, const unsigned int first_selected_component, const std::array< unsigned int, VectorizedArrayType::size()> cells, const std::array< unsigned int, VectorizedArrayType::size()> face_nos, const unsigned int subface_index, const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, const std::array< unsigned int, VectorizedArrayType::size()> face_orientations, const Table< 2, unsigned int > &orientation_map)
static ::ExceptionBase & ExcInternalError()
void do_vectorized_scatter_add(const VectorizedArrayType src, const unsigned int *indices, Number2 *dst_ptr)
void do_vectorized_gather(const Number2 *src_ptr, const unsigned int *indices, VectorizedArrayType &dst)