Reference documentation for deal.II version Git 9297d75edf 2020-11-26 18:52:14 +0100
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
evaluation_kernels.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_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  Assert(data.element_type != MatrixFreeFunctions::tensor_none,
2550  ExcNotImplemented());
2551 
2552  constexpr unsigned int static_dofs_per_face =
2553  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
2555  const unsigned int dofs_per_face =
2556  fe_degree > -1 ?
2557  static_dofs_per_face :
2558  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
2559 
2560  VectorizedArrayType *temp1 = scratch_data;
2561 
2563  template interpolate<true, false>(
2564  n_components, data, values_array, temp1, evaluate_gradients, face_no);
2565 
2566  const unsigned int n_q_points_1d_actual =
2567  fe_degree > -1 ? n_q_points_1d : 0;
2568  if (fe_degree > -1 &&
2569  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
2570  data.element_type <= MatrixFreeFunctions::tensor_symmetric)
2572  true,
2573  dim,
2574  fe_degree,
2575  n_q_points_1d_actual,
2576  VectorizedArrayType>::evaluate_in_face(n_components,
2577  data,
2578  temp1,
2579  values_quad,
2580  gradients_quad,
2581  scratch_data + 2 *
2582  n_components *
2583  dofs_per_face,
2584  evaluate_values,
2585  evaluate_gradients,
2586  subface_index);
2587  else
2589  false,
2590  dim,
2591  fe_degree,
2592  n_q_points_1d_actual,
2593  VectorizedArrayType>::evaluate_in_face(n_components,
2594  data,
2595  temp1,
2596  values_quad,
2597  gradients_quad,
2598  scratch_data + 2 *
2599  n_components *
2600  dofs_per_face,
2601  evaluate_values,
2602  evaluate_gradients,
2603  subface_index);
2604 
2605  if (face_orientation)
2607  n_components,
2608  face_orientation,
2609  orientation_map,
2610  false,
2611  evaluate_values,
2612  evaluate_gradients,
2613  data.n_q_points_face,
2614  scratch_data,
2615  values_quad,
2616  gradients_quad);
2617 
2618  return false;
2619  }
2620  };
2621 
2622 
2623 
2624  template <int dim, typename VectorizedArrayType>
2626  {
2627  template <int fe_degree, int n_q_points_1d>
2628  static bool
2629  run(const unsigned int n_components,
2631  VectorizedArrayType * values_array,
2632  VectorizedArrayType * values_quad,
2633  VectorizedArrayType * gradients_quad,
2634  VectorizedArrayType * scratch_data,
2635  const bool integrate_values,
2636  const bool integrate_gradients,
2637  const unsigned int face_no,
2638  const unsigned int subface_index,
2639  const unsigned int face_orientation,
2640  const Table<2, unsigned int> &orientation_map)
2641  {
2642  Assert(data.element_type != MatrixFreeFunctions::tensor_none,
2643  ExcNotImplemented());
2644 
2645  if (face_orientation)
2647  n_components,
2648  face_orientation,
2649  orientation_map,
2650  true,
2651  integrate_values,
2652  integrate_gradients,
2653  data.n_q_points_face,
2654  scratch_data,
2655  values_quad,
2656  gradients_quad);
2657 
2658  constexpr unsigned int static_dofs_per_face =
2659  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
2661  const unsigned int dofs_per_face =
2662  fe_degree > -1 ?
2663  static_dofs_per_face :
2664  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
2665 
2666  VectorizedArrayType *temp1 = scratch_data;
2667 
2668  const unsigned int n_q_points_1d_actual =
2669  fe_degree > -1 ? n_q_points_1d : 0;
2670  if (fe_degree > -1 &&
2674  true,
2675  dim,
2676  fe_degree,
2677  n_q_points_1d_actual,
2678  VectorizedArrayType>::integrate_in_face(n_components,
2679  data,
2680  temp1,
2681  values_quad,
2682  gradients_quad,
2683  scratch_data +
2684  2 * n_components *
2685  dofs_per_face,
2686  integrate_values,
2687  integrate_gradients,
2688  subface_index);
2689  else
2691  false,
2692  dim,
2693  fe_degree,
2694  n_q_points_1d_actual,
2695  VectorizedArrayType>::integrate_in_face(n_components,
2696  data,
2697  temp1,
2698  values_quad,
2699  gradients_quad,
2700  scratch_data +
2701  2 * n_components *
2702  dofs_per_face,
2703  integrate_values,
2704  integrate_gradients,
2705  subface_index);
2706 
2708  template interpolate<false, false>(n_components,
2709  data,
2710  temp1,
2711  values_array,
2712  integrate_gradients,
2713  face_no);
2714  return false;
2715  }
2716  };
2717 
2718 
2719 
2720  template <int n_face_orientations, typename Processor>
2721  static bool
2723  {
2724  auto n_components = proc.n_components;
2725  auto integrate = proc.integrate;
2726  auto global_vector_ptr = proc.global_vector_ptr;
2727  auto &sm_ptr = proc.sm_ptr;
2728  auto &data = proc.data;
2729  auto &dof_info = proc.dof_info;
2730  auto values_quad = proc.values_quad;
2731  auto gradients_quad = proc.gradients_quad;
2732  auto scratch_data = proc.scratch_data;
2733  auto do_values = proc.do_values;
2734  auto do_gradients = proc.do_gradients;
2735  auto active_fe_index = proc.active_fe_index;
2736  auto first_selected_component = proc.first_selected_component;
2737  auto cells = proc.cells;
2738  auto face_nos = proc.face_nos;
2739  auto subface_index = proc.subface_index;
2740  auto dof_access_index = proc.dof_access_index;
2741  auto face_orientations = proc.face_orientations;
2742  auto &orientation_map = proc.orientation_map;
2743 
2744  static const int dim = Processor::dim_;
2745  static const int fe_degree = Processor::fe_degree_;
2746  using VectorizedArrayType = typename Processor::VectorizedArrayType_;
2747 
2748  using Number = typename Processor::Number_;
2749  using Number2_ = typename Processor::Number2_;
2750 
2751  const unsigned int cell = cells[0];
2752 
2753  // In the case of integration, we do not need to reshuffle the
2754  // data at the quadrature points to adjust for the face
2755  // orientation if the shape functions are nodal at the cell
2756  // boundaries (and we only requested the integration of the
2757  // values) or Hermite shape functions are used. These cases are
2758  // handled later when the values are written back into the
2759  // glrobal vector.
2760  if (integrate &&
2761  (face_orientations[0] > 0 &&
2762  (subface_index < GeometryInfo<dim>::max_children_per_cell ||
2763  !(((do_gradients == false &&
2764  data.data.front().nodal_at_cell_boundaries == true &&
2765  fe_degree > 0) ||
2766  (data.element_type ==
2768  fe_degree > 1)) &&
2769  (dof_info.index_storage_variants[dof_access_index][cell] ==
2771  interleaved_contiguous ||
2772  dof_info.index_storage_variants[dof_access_index][cell] ==
2774  interleaved_contiguous_strided ||
2775  dof_info.index_storage_variants[dof_access_index][cell] ==
2777  interleaved_contiguous_mixed_strides ||
2778  dof_info.index_storage_variants[dof_access_index][cell] ==
2780  contiguous)))))
2781  {
2782  AssertDimension(n_face_orientations, 1);
2784  n_components,
2785  face_orientations[0],
2786  orientation_map,
2787  true,
2788  do_values,
2789  do_gradients,
2790  data.n_q_points_face,
2791  scratch_data,
2792  values_quad,
2793  gradients_quad);
2794  }
2795 
2796  // we know that the gradient weights for the Hermite case on the
2797  // right (side==1) are the negative from the value at the left
2798  // (side==0), so we only read out one of them.
2799  VectorizedArrayType grad_weight =
2800  (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
2801  data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ?
2802  data.data.front()
2803  .shape_data_on_face[0][fe_degree + (integrate ?
2804  (2 - (face_nos[0] % 2)) :
2805  (1 + (face_nos[0] % 2)))] :
2806  VectorizedArrayType(0.0 /*dummy*/);
2807 
2808  constexpr unsigned int static_dofs_per_component =
2809  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) :
2811  constexpr unsigned int static_dofs_per_face =
2812  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
2814  const unsigned int dofs_per_face =
2815  fe_degree > -1 ? static_dofs_per_face :
2816  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
2817 
2818  VectorizedArrayType *temp1 = scratch_data;
2819 
2820  const unsigned int dummy = 0;
2821 
2822  // re-orientation
2823  std::array<const unsigned int *, n_face_orientations> orientation = {};
2824 
2825  if (n_face_orientations == 1)
2826  orientation[0] = (data.data.front().nodal_at_cell_boundaries == true) ?
2827  &data.face_orientations[face_orientations[0]][0] :
2828  &dummy;
2829  else
2830  {
2831  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2832  {
2833  // the loop breaks once an invalid_unsigned_int is hit for
2834  // all cases except the exterior faces in the ECL loop (where
2835  // some faces might be at the boundaries but others not)
2836  if (cells[v] == numbers::invalid_unsigned_int)
2837  continue;
2838 
2839  orientation[v] =
2840  (data.data.front().nodal_at_cell_boundaries == true) ?
2841  &data.face_orientations[face_orientations[v]][0] :
2842  &dummy;
2843  }
2844  }
2845 
2846  // face_to_cell_index_hermite
2847  std::array<const unsigned int *, n_face_orientations> index_array_hermite =
2848  {};
2849 
2850  if (n_face_orientations == 1)
2851  index_array_hermite[0] =
2852  (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
2853  data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ?
2854  &data.face_to_cell_index_hermite(face_nos[0], 0) :
2855  &dummy;
2856 
2857  if (n_face_orientations > 1 &&
2858  data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
2859  data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite)
2860  {
2861  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2862  {
2863  if (cells[v] == numbers::invalid_unsigned_int)
2864  continue;
2865 
2866  grad_weight[v] =
2867  data.data.front().shape_data_on_face
2868  [0][fe_degree + (integrate ? (2 - (face_nos[v] % 2)) :
2869  (1 + (face_nos[v] % 2)))][v];
2870 
2871  index_array_hermite[v] =
2872  &data.face_to_cell_index_hermite(face_nos[v], 0);
2873  }
2874  }
2875 
2876  // face_to_cell_index_nodal
2877  std::array<const unsigned int *, n_face_orientations> index_array_nodal =
2878  {};
2879 
2880  if (n_face_orientations == 1)
2881  index_array_nodal[0] =
2882  (data.data.front().nodal_at_cell_boundaries == true) ?
2883  &data.face_to_cell_index_nodal(face_nos[0], 0) :
2884  &dummy;
2885 
2886  if (n_face_orientations > 1 &&
2887  (data.data.front().nodal_at_cell_boundaries == true))
2888  {
2889  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2890  {
2891  if (cells[v] == numbers::invalid_unsigned_int)
2892  continue;
2893 
2894  index_array_nodal[v] =
2895  &data.face_to_cell_index_nodal(face_nos[v], 0);
2896  }
2897  }
2898 
2899  const auto reorientate = [&](const unsigned int v, const unsigned int i) {
2900  return (dim < 3 ||
2901  face_orientations[n_face_orientations == 1 ? 0 : v] == 0 ||
2902  subface_index < GeometryInfo<dim>::max_children_per_cell) ?
2903  i :
2904  orientation[v][i];
2905  };
2906 
2907  // this variable keeps track of whether we are able to directly write
2908  // the results into the result (function returns true) or not, requiring
2909  // an additional call to another function
2910  bool accesses_global_vector = true;
2911 
2912  for (unsigned int comp = 0; comp < n_components; ++comp)
2913  {
2914  if (integrate)
2915  proc.in_face_operation(temp1, comp);
2916 
2917  // we can only use the fast functions if we know the polynomial degree
2918  // as a template parameter (fe_degree != -1), and it only makes sense
2919  // to use the functions for at least linear functions for values on
2920  // the faces and quadratic functions for gradients on the faces, so
2921  // include the switch here
2922  if ((do_gradients == false &&
2923  data.data.front().nodal_at_cell_boundaries == true &&
2924  fe_degree > 0) ||
2925  (data.element_type ==
2927  fe_degree > 1))
2928  {
2929  // case 1: contiguous and interleaved indices
2930  if (n_face_orientations == 1 &&
2931  dof_info.index_storage_variants[dof_access_index][cell] ==
2933  interleaved_contiguous)
2934  {
2935  AssertDimension(n_face_orientations, 1);
2936 
2938  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2939  VectorizedArrayType::size());
2940  Number2_ *vector_ptr =
2941  global_vector_ptr +
2942  dof_info.dof_indices_contiguous[dof_access_index]
2943  [cell *
2944  VectorizedArrayType::size()] +
2945  (dof_info
2946  .component_dof_indices_offset[active_fe_index]
2947  [first_selected_component] +
2948  comp * static_dofs_per_component) *
2949  VectorizedArrayType::size();
2950 
2951  if (fe_degree > 1 && do_gradients == true)
2952  {
2953  for (unsigned int i = 0; i < dofs_per_face; ++i)
2954  {
2955  if (n_face_orientations == 1)
2956  {
2957  const unsigned int ind1 =
2958  index_array_hermite[0][2 * i];
2959  const unsigned int ind2 =
2960  index_array_hermite[0][2 * i + 1];
2961  AssertIndexRange(ind1,
2962  data.dofs_per_component_on_cell);
2963  AssertIndexRange(ind2,
2964  data.dofs_per_component_on_cell);
2965  const unsigned int i_ = reorientate(0, i);
2966  proc.hermite_grad_vectorized(
2967  temp1[i_],
2968  temp1[i_ + dofs_per_face],
2969  vector_ptr + ind1 * VectorizedArrayType::size(),
2970  vector_ptr + ind2 * VectorizedArrayType::size(),
2971  grad_weight);
2972  }
2973  else
2974  {
2975  Assert(false, ExcNotImplemented());
2976  }
2977  }
2978  }
2979  else
2980  {
2981  for (unsigned int i = 0; i < dofs_per_face; ++i)
2982  {
2983  if (n_face_orientations == 1)
2984  {
2985  const unsigned int i_ = reorientate(0, i);
2986  const unsigned int ind = index_array_nodal[0][i];
2987  proc.value_vectorized(
2988  temp1[i_],
2989  vector_ptr + ind * VectorizedArrayType::size());
2990  }
2991  else
2992  {
2993  Assert(false, ExcNotImplemented());
2994  }
2995  }
2996  }
2997  }
2998 
2999  // case 2: contiguous and interleaved indices with fixed stride
3000  else if (n_face_orientations == 1 &&
3001  dof_info.index_storage_variants[dof_access_index][cell] ==
3003  interleaved_contiguous_strided)
3004  {
3005  AssertDimension(n_face_orientations, 1);
3006 
3008  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
3009  VectorizedArrayType::size());
3010  const unsigned int *indices =
3011  &dof_info.dof_indices_contiguous[dof_access_index]
3012  [cell *
3013  VectorizedArrayType::size()];
3014  Number2_ *vector_ptr =
3015  global_vector_ptr +
3016  (comp * static_dofs_per_component +
3017  dof_info
3018  .component_dof_indices_offset[active_fe_index]
3019  [first_selected_component]) *
3020  VectorizedArrayType::size();
3021  if (fe_degree > 1 && do_gradients == true)
3022  {
3023  for (unsigned int i = 0; i < dofs_per_face; ++i)
3024  {
3025  if (n_face_orientations == 1)
3026  {
3027  const unsigned int i_ = reorientate(0, i);
3028  const unsigned int ind1 =
3029  index_array_hermite[0][2 * i] *
3030  VectorizedArrayType::size();
3031  const unsigned int ind2 =
3032  index_array_hermite[0][2 * i + 1] *
3033  VectorizedArrayType::size();
3034  proc.hermite_grad_vectorized_indexed(
3035  temp1[i_],
3036  temp1[i_ + dofs_per_face],
3037  vector_ptr + ind1,
3038  vector_ptr + ind2,
3039  grad_weight,
3040  indices,
3041  indices);
3042  }
3043  else
3044  {
3045  Assert(false, ExcNotImplemented());
3046  }
3047  }
3048  }
3049  else
3050  {
3051  for (unsigned int i = 0; i < dofs_per_face; ++i)
3052  {
3053  if (n_face_orientations == 1)
3054  {
3055  const unsigned int i_ = reorientate(0, i);
3056  const unsigned int ind =
3057  index_array_nodal[0][i] *
3058  VectorizedArrayType::size();
3059  proc.value_vectorized_indexed(temp1[i_],
3060  vector_ptr + ind,
3061  indices);
3062  }
3063  else
3064  {
3065  Assert(false, ExcNotImplemented());
3066  }
3067  }
3068  }
3069  }
3070 
3071  // case 3: contiguous and interleaved indices with mixed stride
3072  else if (n_face_orientations == 1 &&
3073  dof_info.index_storage_variants[dof_access_index][cell] ==
3075  interleaved_contiguous_mixed_strides)
3076  {
3077  AssertDimension(n_face_orientations, 1);
3078 
3079  const unsigned int *strides =
3080  &dof_info.dof_indices_interleave_strides
3081  [dof_access_index][cell * VectorizedArrayType::size()];
3082  unsigned int indices[VectorizedArrayType::size()];
3083  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3084  indices[v] =
3085  dof_info.dof_indices_contiguous
3086  [dof_access_index]
3087  [cell * VectorizedArrayType::size() + v] +
3088  (dof_info
3089  .component_dof_indices_offset[active_fe_index]
3090  [first_selected_component] +
3091  comp * static_dofs_per_component) *
3092  strides[v];
3093  const unsigned int n_filled_lanes =
3094  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
3095 
3096  if (fe_degree > 1 && do_gradients == true)
3097  {
3098  if (n_filled_lanes == VectorizedArrayType::size())
3099  for (unsigned int i = 0; i < dofs_per_face; ++i)
3100  {
3101  if (n_face_orientations == 1)
3102  {
3103  const unsigned int i_ = reorientate(0, i);
3104  unsigned int ind1[VectorizedArrayType::size()];
3106  for (unsigned int v = 0;
3107  v < VectorizedArrayType::size();
3108  ++v)
3109  ind1[v] =
3110  indices[v] +
3111  index_array_hermite[0 /*TODO*/][2 * i] *
3112  strides[v];
3113  unsigned int ind2[VectorizedArrayType::size()];
3115  for (unsigned int v = 0;
3116  v < VectorizedArrayType::size();
3117  ++v)
3118  ind2[v] =
3119  indices[v] +
3120  index_array_hermite[0 /*TODO*/][2 * i + 1] *
3121  strides[v];
3122  proc.hermite_grad_vectorized_indexed(
3123  temp1[i_],
3124  temp1[i_ + dofs_per_face],
3125  global_vector_ptr,
3126  global_vector_ptr,
3127  grad_weight,
3128  ind1,
3129  ind2);
3130  }
3131  else
3132  {
3133  Assert(false, ExcNotImplemented());
3134  }
3135  }
3136  else
3137  {
3138  if (integrate == false)
3139  for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
3140  temp1[i] = VectorizedArrayType();
3141 
3142  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3143  for (unsigned int i = 0; i < dofs_per_face; ++i)
3144  {
3145  const unsigned int i_ =
3146  reorientate(n_face_orientations == 1 ? 0 : v,
3147  i);
3148  proc.hermite_grad(
3149  temp1[i_][v],
3150  temp1[i_ + dofs_per_face][v],
3151  global_vector_ptr
3152  [indices[v] +
3153  index_array_hermite
3154  [n_face_orientations == 1 ? 0 : v]
3155  [2 * i] *
3156  strides[v]],
3157  global_vector_ptr
3158  [indices[v] +
3159  index_array_hermite
3160  [n_face_orientations == 1 ? 0 : v]
3161  [2 * i + 1] *
3162  strides[v]],
3163  grad_weight[n_face_orientations == 1 ? 0 : v]);
3164  }
3165  }
3166  }
3167  else
3168  {
3169  if (n_filled_lanes == VectorizedArrayType::size())
3170  for (unsigned int i = 0; i < dofs_per_face; ++i)
3171  {
3172  if (n_face_orientations == 1)
3173  {
3174  unsigned int ind[VectorizedArrayType::size()];
3176  for (unsigned int v = 0;
3177  v < VectorizedArrayType::size();
3178  ++v)
3179  ind[v] = indices[v] +
3180  index_array_nodal[0][i] * strides[v];
3181  const unsigned int i_ = reorientate(0, i);
3182  proc.value_vectorized_indexed(temp1[i_],
3183  global_vector_ptr,
3184  ind);
3185  }
3186  else
3187  {
3188  Assert(false, ExcNotImplemented());
3189  }
3190  }
3191  else
3192  {
3193  if (integrate == false)
3194  for (unsigned int i = 0; i < dofs_per_face; ++i)
3195  temp1[i] = VectorizedArrayType();
3196 
3197  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3198  for (unsigned int i = 0; i < dofs_per_face; ++i)
3199  proc.value(
3200  temp1[reorientate(
3201  n_face_orientations == 1 ? 0 : v, i)][v],
3202  global_vector_ptr
3203  [indices[v] +
3204  index_array_nodal
3205  [n_face_orientations == 1 ? 0 : v][i] *
3206  strides[v]]);
3207  }
3208  }
3209  }
3210 
3211  // case 4: contiguous indices without interleaving
3212  else if (n_face_orientations > 1 ||
3213  dof_info.index_storage_variants[dof_access_index][cell] ==
3215  contiguous)
3216  {
3217  const unsigned int *indices =
3218  &dof_info.dof_indices_contiguous[dof_access_index]
3219  [cell *
3220  VectorizedArrayType::size()];
3221  Number2_ *vector_ptr =
3222  global_vector_ptr + comp * static_dofs_per_component +
3223  dof_info
3224  .component_dof_indices_offset[active_fe_index]
3225  [first_selected_component];
3226 
3227  const unsigned int n_filled_lanes =
3228  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
3229 
3230  const bool vectorization_possible =
3231  (n_face_orientations == 1) &&
3232  (n_filled_lanes == VectorizedArrayType::size()) &&
3233  (sm_ptr != nullptr);
3234 
3235  std::array<Number2_ *, VectorizedArrayType::size()>
3236  vector_ptrs = {};
3237 
3238  if (vectorization_possible == false)
3239  {
3240  if (n_face_orientations == 1)
3241  {
3242  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3243  if (sm_ptr == nullptr)
3244  {
3245  vector_ptrs[v] = vector_ptr + indices[v];
3246  }
3247  else
3248  {
3249  const auto &temp =
3250  dof_info.dof_indices_contiguous_sm
3251  [dof_access_index]
3252  [cell * VectorizedArrayType::size() + v];
3253  vector_ptrs[v] = const_cast<Number *>(
3254  sm_ptr->operator[](temp.first).data() +
3255  temp.second + comp * static_dofs_per_component +
3256  dof_info.component_dof_indices_offset
3257  [active_fe_index][first_selected_component]);
3258  }
3259  }
3260  else if (n_face_orientations == VectorizedArrayType::size())
3261  {
3262  for (unsigned int v = 0;
3263  v < VectorizedArrayType::size();
3264  ++v)
3265  if (cells[v] != numbers::invalid_unsigned_int)
3266  {
3267  if (sm_ptr == nullptr)
3268  {
3269  vector_ptrs[v] =
3270  vector_ptr +
3271  dof_info
3272  .dof_indices_contiguous[dof_access_index]
3273  [cells[v]];
3274  }
3275  else
3276  {
3277  const auto &temp =
3278  dof_info.dof_indices_contiguous_sm
3279  [dof_access_index][cells[v]];
3280  vector_ptrs[v] = const_cast<Number *>(
3281  sm_ptr->operator[](temp.first).data() +
3282  temp.second +
3283  comp * static_dofs_per_component +
3284  dof_info.component_dof_indices_offset
3285  [active_fe_index]
3286  [first_selected_component]);
3287  }
3288  }
3289  }
3290  else
3291  {
3292  Assert(false, ExcNotImplemented());
3293  }
3294  }
3295 
3296  if (do_gradients == true &&
3297  data.element_type ==
3299  {
3300  if (vectorization_possible)
3301  for (unsigned int i = 0; i < dofs_per_face; ++i)
3302  {
3303  const unsigned int ind1 =
3304  index_array_hermite[0][2 * i];
3305  const unsigned int ind2 =
3306  index_array_hermite[0][2 * i + 1];
3307  const unsigned int i_ = reorientate(0, i);
3308 
3309  proc.hermite_grad_vectorized_indexed(
3310  temp1[i_],
3311  temp1[i_ + dofs_per_face],
3312  vector_ptr + ind1,
3313  vector_ptr + ind2,
3314  grad_weight,
3315  indices,
3316  indices);
3317  }
3318  else if (n_face_orientations == 1)
3319  for (unsigned int i = 0; i < dofs_per_face; ++i)
3320  {
3321  const unsigned int ind1 =
3322  index_array_hermite[0][2 * i];
3323  const unsigned int ind2 =
3324  index_array_hermite[0][2 * i + 1];
3325  const unsigned int i_ = reorientate(0, i);
3326 
3327  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3328  proc.hermite_grad(temp1[i_][v],
3329  temp1[i_ + dofs_per_face][v],
3330  vector_ptrs[v][ind1],
3331  vector_ptrs[v][ind2],
3332  grad_weight[v]);
3333 
3334  if (integrate == false)
3335  for (unsigned int v = n_filled_lanes;
3336  v < VectorizedArrayType::size();
3337  ++v)
3338  {
3339  temp1[i_][v] = 0.0;
3340  temp1[i_ + dofs_per_face][v] = 0.0;
3341  }
3342  }
3343  else
3344  {
3345  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3346  for (unsigned int i = 0; i < dofs_per_face; ++i)
3347  proc.hermite_grad(
3348  temp1[reorientate(v, i)][v],
3349  temp1[reorientate(v, i) + dofs_per_face][v],
3350  vector_ptrs[v][index_array_hermite[v][2 * i]],
3351  vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
3352  grad_weight[v]);
3353  }
3354  }
3355  else
3356  {
3357  if (vectorization_possible)
3358  for (unsigned int i = 0; i < dofs_per_face; ++i)
3359  {
3360  const unsigned int ind = index_array_nodal[0][i];
3361  const unsigned int i_ = reorientate(0, i);
3362 
3363  proc.value_vectorized_indexed(temp1[i_],
3364  vector_ptr + ind,
3365  indices);
3366  }
3367  else if (n_face_orientations == 1)
3368  for (unsigned int i = 0; i < dofs_per_face; ++i)
3369  {
3370  const unsigned int ind = index_array_nodal[0][i];
3371  const unsigned int i_ = reorientate(0, i);
3372 
3373  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3374  proc.value(temp1[i_][v], vector_ptrs[v][ind]);
3375 
3376  if (integrate == false)
3377  for (unsigned int v = n_filled_lanes;
3378  v < VectorizedArrayType::size();
3379  ++v)
3380  temp1[i_][v] = 0.0;
3381  }
3382  else
3383  for (unsigned int i = 0; i < dofs_per_face; ++i)
3384  {
3385  for (unsigned int v = 0;
3386  v < VectorizedArrayType::size();
3387  ++v)
3388  if (cells[v] != numbers::invalid_unsigned_int)
3389  proc.value(
3390  temp1[reorientate(v, i)][v],
3391  vector_ptrs[v][index_array_nodal[v][i]]);
3392  }
3393  }
3394  }
3395  else
3396  {
3397  // case 5: default vector access
3398  AssertDimension(n_face_orientations, 1);
3399 
3400  // for the integrate_scatter path (integrate == true), we
3401  // need to only prepare the data in this function for all
3402  // components to later call distribute_local_to_global();
3403  // for the gather_evaluate path (integrate == false), we
3404  // instead want to leave early because we need to get the
3405  // vector data from somewhere else
3406  proc.default_operation(temp1, comp);
3407  if (integrate)
3408  accesses_global_vector = false;
3409  else
3410  return false;
3411  }
3412  }
3413  else
3414  {
3415  // case 5: default vector access
3416  AssertDimension(n_face_orientations, 1);
3417 
3418  proc.default_operation(temp1, comp);
3419  if (integrate)
3420  accesses_global_vector = false;
3421  else
3422  return false;
3423  }
3424 
3425  if (!integrate)
3426  proc.in_face_operation(temp1, comp);
3427  }
3428 
3429  if (!integrate &&
3430  (face_orientations[0] > 0 &&
3432  {
3433  AssertDimension(n_face_orientations, 1);
3435  n_components,
3436  face_orientations[0],
3437  orientation_map,
3438  false,
3439  do_values,
3440  do_gradients,
3441  data.n_q_points_face,
3442  scratch_data,
3443  values_quad,
3444  gradients_quad);
3445  }
3446 
3447  return accesses_global_vector;
3448  }
3449 
3450 
3451 
3452  template <int dim,
3453  typename Number,
3454  typename VectorizedArrayType,
3455  typename Number2 = Number>
3457  {
3458  template <int fe_degree, int n_q_points_1d>
3459  static bool
3460  run(const unsigned int n_components,
3461  const unsigned int n_face_orientations,
3462  const Number2 * src_ptr,
3463  const std::vector<ArrayView<const Number>> *sm_ptr,
3465  const MatrixFreeFunctions::DoFInfo & dof_info,
3466  VectorizedArrayType * values_quad,
3467  VectorizedArrayType *gradients_quad,
3468  VectorizedArrayType *scratch_data,
3469  const bool evaluate_values,
3470  const bool evaluate_gradients,
3471  const unsigned int active_fe_index,
3472  const unsigned int first_selected_component,
3473  const std::array<unsigned int, VectorizedArrayType::size()> cells,
3474  const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
3475  const unsigned int subface_index,
3476  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
3477  const std::array<unsigned int, VectorizedArrayType::size()>
3478  face_orientations,
3479  const Table<2, unsigned int> &orientation_map)
3480  {
3481  if (src_ptr == nullptr)
3482  {
3483  return false;
3484  }
3485 
3486  (void)sm_ptr;
3487 
3488  Processor<fe_degree, n_q_points_1d> p(n_components,
3489  false,
3490  src_ptr,
3491  sm_ptr,
3492  data,
3493  dof_info,
3494  values_quad,
3495  gradients_quad,
3496  scratch_data,
3497  evaluate_values,
3498  evaluate_gradients,
3499  active_fe_index,
3500  first_selected_component,
3501  cells,
3502  face_nos,
3503  subface_index,
3504  dof_access_index,
3505  face_orientations,
3506  orientation_map);
3507 
3508  if (n_face_orientations == VectorizedArrayType::size())
3509  return fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
3510  p);
3511  else
3512  return fe_face_evaluation_process_and_io<1>(p);
3513  }
3514 
3515  private:
3516  template <int fe_degree, int n_q_points_1d>
3517  struct Processor
3518  {
3519  static const int dim_ = dim;
3520  static const int fe_degree_ = fe_degree;
3521  static const int n_q_points_1d_ = n_q_points_1d;
3522  using VectorizedArrayType_ = VectorizedArrayType;
3523  using Number_ = Number;
3524  using Number2_ = const Number2;
3525 
3527  const unsigned int n_components,
3528  const bool integrate,
3529  const Number2 * global_vector_ptr,
3530  const std::vector<ArrayView<const Number>> *sm_ptr,
3532  const MatrixFreeFunctions::DoFInfo & dof_info,
3533  VectorizedArrayType * values_quad,
3534  VectorizedArrayType *gradients_quad,
3535  VectorizedArrayType *scratch_data,
3536  const bool do_values,
3537  const bool do_gradients,
3538  const unsigned int active_fe_index,
3539  const unsigned int first_selected_component,
3540  const std::array<unsigned int, VectorizedArrayType::size()> cells,
3541  const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
3542  const unsigned int subface_index,
3543  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
3544  const std::array<unsigned int, VectorizedArrayType::size()>
3545  face_orientations,
3546  const Table<2, unsigned int> &orientation_map)
3547  : n_components(n_components)
3548  , integrate(integrate)
3549  , global_vector_ptr(global_vector_ptr)
3550  , sm_ptr(sm_ptr)
3551  , data(data)
3552  , dof_info(dof_info)
3553  , values_quad(values_quad)
3554  , gradients_quad(gradients_quad)
3555  , scratch_data(scratch_data)
3556  , do_values(do_values)
3557  , do_gradients(do_gradients)
3558  , active_fe_index(active_fe_index)
3559  , first_selected_component(first_selected_component)
3560  , cells(cells)
3561  , face_nos(face_nos)
3562  , subface_index(subface_index)
3563  , dof_access_index(dof_access_index)
3564  , face_orientations(face_orientations)
3565  , orientation_map(orientation_map)
3566  {}
3567 
3568  template <typename T0, typename T1, typename T2>
3569  void
3571  T0 & temp_2,
3572  const T1 src_ptr_1,
3573  const T1 src_ptr_2,
3574  const T2 &grad_weight)
3575  {
3576  do_vectorized_read(src_ptr_1, temp_1);
3577  do_vectorized_read(src_ptr_2, temp_2);
3578  temp_2 = grad_weight * (temp_1 - temp_2);
3579  }
3580 
3581  template <typename T1, typename T2>
3582  void
3583  value_vectorized(T1 &temp, const T2 src_ptr)
3584  {
3585  do_vectorized_read(src_ptr, temp);
3586  }
3587 
3588  template <typename T0, typename T1, typename T2, typename T3>
3589  void
3591  T0 & temp_2,
3592  const T1 src_ptr_1,
3593  const T1 src_ptr_2,
3594  const T2 &grad_weight,
3595  const T3 &indices_1,
3596  const T3 &indices_2)
3597  {
3598  do_vectorized_gather(src_ptr_1, indices_1, temp_1);
3599  do_vectorized_gather(src_ptr_2, indices_2, temp_2);
3600  temp_2 = grad_weight * (temp_1 - temp_2);
3601  }
3602 
3603  template <typename T0, typename T1, typename T2>
3604  void
3605  value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
3606  {
3607  do_vectorized_gather(src_ptr, indices, temp);
3608  }
3609 
3610  template <typename T0, typename T1, typename T2>
3611  void
3612  hermite_grad(T0 & temp_1,
3613  T0 & temp_2,
3614  const T1 &src_ptr_1,
3615  const T2 &src_ptr_2,
3616  const T2 &grad_weight)
3617  {
3618  // case 3a)
3619  temp_1 = src_ptr_1;
3620  temp_2 = grad_weight * (temp_1 - src_ptr_2);
3621  }
3622 
3623  template <typename T1, typename T2>
3624  void
3625  value(T1 &temp, const T2 &src_ptr)
3626  {
3627  // case 3b)
3628  temp = src_ptr;
3629  }
3630 
3631  template <typename T1>
3632  void
3633  default_operation(const T1 &, const unsigned int)
3634  {
3635  // case 5)
3636  }
3637 
3638  template <typename T1>
3639  void
3640  in_face_operation(T1 &temp1, const unsigned int comp)
3641  {
3642  const unsigned int dofs_per_face =
3643  fe_degree > -1 ?
3644  Utilities::pow(fe_degree + 1, dim - 1) :
3645  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
3646  const unsigned int n_q_points =
3647  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
3648  data.n_q_points_face;
3649  if (fe_degree > -1 &&
3650  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
3652  FEFaceEvaluationImpl<true,
3653  dim,
3654  fe_degree,
3655  n_q_points_1d,
3656  VectorizedArrayType>::
3657  evaluate_in_face(/* n_components */ 1,
3658  data,
3659  temp1,
3660  values_quad + comp * n_q_points,
3661  gradients_quad + comp * dim * n_q_points,
3662  scratch_data + 2 * dofs_per_face,
3663  do_values,
3664  do_gradients,
3665  subface_index);
3666  else
3667  FEFaceEvaluationImpl<false,
3668  dim,
3669  fe_degree,
3670  n_q_points_1d,
3671  VectorizedArrayType>::
3672  evaluate_in_face(/* n_components */ 1,
3673  data,
3674  temp1,
3675  values_quad + comp * n_q_points,
3676  gradients_quad + comp * dim * n_q_points,
3677  scratch_data + 2 * dofs_per_face,
3678  do_values,
3679  do_gradients,
3680  subface_index);
3681  }
3682 
3683  const unsigned int n_components;
3684  const bool integrate;
3685  const Number2 * global_vector_ptr;
3686  const std::vector<ArrayView<const Number>> *sm_ptr;
3689  VectorizedArrayType * values_quad;
3690  VectorizedArrayType * gradients_quad;
3691  VectorizedArrayType * scratch_data;
3692  const bool do_values;
3693  const bool do_gradients;
3694  const unsigned int active_fe_index;
3695  const unsigned int first_selected_component;
3696  const std::array<unsigned int, VectorizedArrayType::size()> cells;
3697  const std::array<unsigned int, VectorizedArrayType::size()> face_nos;
3698  const unsigned int subface_index;
3700  const std::array<unsigned int, VectorizedArrayType::size()>
3703  };
3704  };
3705 
3706  template <int dim,
3707  typename Number,
3708  typename VectorizedArrayType,
3709  typename Number2 = Number>
3711  {
3712  template <int fe_degree, int n_q_points_1d>
3713  static bool
3714  run(const unsigned int n_components,
3715  const unsigned int n_face_orientations,
3716  Number2 * dst_ptr,
3717  const std::vector<ArrayView<const Number2>> *sm_ptr,
3719  const MatrixFreeFunctions::DoFInfo & dof_info,
3720  VectorizedArrayType * values_array,
3721  VectorizedArrayType * values_quad,
3722  VectorizedArrayType *gradients_quad,
3723  VectorizedArrayType *scratch_data,
3724  const bool integrate_values,
3725  const bool integrate_gradients,
3726  const unsigned int active_fe_index,
3727  const unsigned int first_selected_component,
3728  const std::array<unsigned int, VectorizedArrayType::size()> cells,
3729  const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
3730  const unsigned int subface_index,
3731  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
3732  const std::array<unsigned int, VectorizedArrayType::size()>
3733  face_orientations,
3734  const Table<2, unsigned int> &orientation_map)
3735  {
3736  (void)sm_ptr;
3737 
3738  if (dst_ptr == nullptr)
3739  {
3740  AssertDimension(n_face_orientations, 1);
3741  AssertDimension(n_face_orientations, 1);
3742 
3743  // for block vectors simply integrate
3745  template run<fe_degree, n_q_points_1d>(n_components,
3746  data,
3747  values_array,
3748  values_quad,
3749  gradients_quad,
3750  scratch_data,
3751  integrate_values,
3752  integrate_gradients,
3753  face_nos[0],
3754  subface_index,
3755  face_orientations[0],
3756  orientation_map);
3757 
3758  // default vector access
3759  return false;
3760  }
3761 
3762 
3763  Processor<fe_degree, n_q_points_1d> p(values_array,
3764  n_components,
3765  true,
3766  dst_ptr,
3767  sm_ptr,
3768  data,
3769  dof_info,
3770  values_quad,
3771  gradients_quad,
3772  scratch_data,
3773  integrate_values,
3774  integrate_gradients,
3775  active_fe_index,
3776  first_selected_component,
3777  cells,
3778  face_nos,
3779  subface_index,
3780  dof_access_index,
3781  face_orientations,
3782  orientation_map);
3783 
3784  if (n_face_orientations == VectorizedArrayType::size())
3785  return fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
3786  p);
3787  else
3788  return fe_face_evaluation_process_and_io<1>(p);
3789  }
3790 
3791  private:
3792  template <int fe_degree, int n_q_points_1d>
3793  struct Processor
3794  {
3795  static const int dim_ = dim;
3796  static const int fe_degree_ = fe_degree;
3797  static const int n_q_points_1d_ = n_q_points_1d;
3798  using VectorizedArrayType_ = VectorizedArrayType;
3799  using Number_ = Number;
3800  using Number2_ = Number2;
3801 
3802 
3804  VectorizedArrayType * values_array,
3805  const unsigned int n_components,
3806  const bool integrate,
3807  Number2 * global_vector_ptr,
3808  const std::vector<ArrayView<const Number>> *sm_ptr,
3810  const MatrixFreeFunctions::DoFInfo & dof_info,
3811  VectorizedArrayType * values_quad,
3812  VectorizedArrayType *gradients_quad,
3813  VectorizedArrayType *scratch_data,
3814  const bool do_values,
3815  const bool do_gradients,
3816  const unsigned int active_fe_index,
3817  const unsigned int first_selected_component,
3818  const std::array<unsigned int, VectorizedArrayType::size()> cells,
3819  const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
3820  const unsigned int subface_index,
3821  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
3822  const std::array<unsigned int, VectorizedArrayType::size()>
3823  face_orientations,
3824  const Table<2, unsigned int> &orientation_map)
3825  : values_array(values_array)
3826  , n_components(n_components)
3827  , integrate(integrate)
3828  , global_vector_ptr(global_vector_ptr)
3829  , sm_ptr(sm_ptr)
3830  , data(data)
3831  , dof_info(dof_info)
3832  , values_quad(values_quad)
3833  , gradients_quad(gradients_quad)
3834  , scratch_data(scratch_data)
3835  , do_values(do_values)
3836  , do_gradients(do_gradients)
3837  , active_fe_index(active_fe_index)
3838  , first_selected_component(first_selected_component)
3839  , cells(cells)
3840  , face_nos(face_nos)
3841  , subface_index(subface_index)
3842  , dof_access_index(dof_access_index)
3843  , face_orientations(face_orientations)
3844  , orientation_map(orientation_map)
3845  {}
3846 
3847  template <typename T0, typename T1, typename T2, typename T3, typename T4>
3848  void
3849  hermite_grad_vectorized(const T0 &temp_1,
3850  const T1 &temp_2,
3851  T2 dst_ptr_1,
3852  T3 dst_ptr_2,
3853  const T4 &grad_weight)
3854  {
3855  // case 1a)
3856  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
3857  const VectorizedArrayType grad = grad_weight * temp_2;
3858  do_vectorized_add(val, dst_ptr_1);
3859  do_vectorized_add(grad, dst_ptr_2);
3860  }
3861 
3862  template <typename T0, typename T1>
3863  void
3864  value_vectorized(const T0 &temp, T1 dst_ptr)
3865  {
3866  // case 1b)
3867  do_vectorized_add(temp, dst_ptr);
3868  }
3869 
3870  template <typename T0, typename T1, typename T2, typename T3>
3871  void
3873  const T0 &temp_2,
3874  T1 dst_ptr_1,
3875  T1 dst_ptr_2,
3876  const T2 &grad_weight,
3877  const T3 &indices_1,
3878  const T3 &indices_2)
3879  {
3880  // case 2a)
3881  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
3882  const VectorizedArrayType grad = grad_weight * temp_2;
3883  do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
3884  do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
3885  }
3886 
3887  template <typename T0, typename T1, typename T2>
3888  void
3889  value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
3890  {
3891  // case 2b)
3892  do_vectorized_scatter_add(temp, indices, dst_ptr);
3893  }
3894 
3895  template <typename T0, typename T1, typename T2>
3896  void
3897  hermite_grad(const T0 &temp_1,
3898  const T0 &temp_2,
3899  T1 & dst_ptr_1,
3900  T1 & dst_ptr_2,
3901  const T2 &grad_weight)
3902  {
3903  // case 3a)
3904  const Number val = temp_1 - grad_weight * temp_2;
3905  const Number grad = grad_weight * temp_2;
3906  dst_ptr_1 += val;
3907  dst_ptr_2 += grad;
3908  }
3909 
3910  template <typename T0, typename T1>
3911  void
3912  value(const T0 &temp, T1 &dst_ptr)
3913  {
3914  // case 3b)
3915  dst_ptr += temp;
3916  }
3917 
3918  template <typename T0>
3919  void
3920  default_operation(const T0 &temp1, const unsigned int comp)
3921  {
3922  // case 5: default vector access, must be handled separately, just do
3923  // the face-normal interpolation
3924 
3926  template interpolate<false, false>(
3927  /* n_components */ 1,
3928  data,
3929  temp1,
3930  values_array + comp * data.dofs_per_component_on_cell,
3931  do_gradients,
3932  face_nos[0]);
3933  }
3934 
3935  template <typename T0>
3936  void
3937  in_face_operation(T0 &temp1, const unsigned int comp)
3938  {
3939  const unsigned int dofs_per_face =
3940  fe_degree > -1 ?
3941  Utilities::pow(fe_degree + 1, dim - 1) :
3942  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
3943  const unsigned int n_q_points =
3944  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
3945  data.n_q_points_face;
3946  if (fe_degree > -1 &&
3947  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
3948  data.element_type <=
3951  dim,
3952  fe_degree,
3953  n_q_points_1d,
3954  VectorizedArrayType>::
3955  integrate_in_face(/* n_components */ 1,
3956  data,
3957  temp1,
3958  values_quad + comp * n_q_points,
3959  gradients_quad + dim * comp * n_q_points,
3960  scratch_data + 2 * dofs_per_face,
3961  do_values,
3962  do_gradients,
3963  subface_index);
3964  else
3966  dim,
3967  fe_degree,
3968  n_q_points_1d,
3969  VectorizedArrayType>::
3970  integrate_in_face(/* n_components */ 1,
3971  data,
3972  temp1,
3973  values_quad + comp * n_q_points,
3974  gradients_quad + dim * comp * n_q_points,
3975  scratch_data + 2 * dofs_per_face,
3976  do_values,
3977  do_gradients,
3978  subface_index);
3979  }
3980 
3981  VectorizedArrayType *values_array;
3982 
3983 
3984  const unsigned int n_components;
3985  const bool integrate;
3987  const std::vector<ArrayView<const Number>> *sm_ptr;
3990  VectorizedArrayType * values_quad;
3991  VectorizedArrayType * gradients_quad;
3992  VectorizedArrayType * scratch_data;
3993  const bool do_values;
3994  const bool do_gradients;
3995  const unsigned int active_fe_index;
3996  const unsigned int first_selected_component;
3997  const std::array<unsigned int, VectorizedArrayType::size()> cells;
3998  const std::array<unsigned int, VectorizedArrayType::size()> face_nos;
3999  const unsigned int subface_index;
4001  const std::array<unsigned int, VectorizedArrayType::size()>
4004  };
4005  };
4006 
4007 
4008 
4013  template <int dim, typename Number>
4015  {
4016  template <int fe_degree, int = 0>
4017  static bool
4018  run(const unsigned int n_components,
4019  const FEEvaluationBaseData<dim,
4020  typename Number::value_type,
4021  false,
4022  Number> &fe_eval,
4023  const Number * in_array,
4024  Number * out_array,
4025  typename std::enable_if<fe_degree != -1>::type * = nullptr)
4026  {
4027  constexpr unsigned int dofs_per_component =
4028  Utilities::pow(fe_degree + 1, dim);
4029 
4030  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
4031  Assert(fe_eval.get_shape_info().element_type <=
4033  ExcNotImplemented());
4034 
4036  dim,
4037  fe_degree + 1,
4038  fe_degree + 1,
4039  Number>
4040  evaluator(
4043  fe_eval.get_shape_info().data.front().inverse_shape_values_eo);
4044 
4045  for (unsigned int d = 0; d < n_components; ++d)
4046  {
4047  const Number *in = in_array + d * dofs_per_component;
4048  Number * out = out_array + d * dofs_per_component;
4049  // Need to select 'apply' method with hessian slot because values
4050  // assume symmetries that do not exist in the inverse shapes
4051  evaluator.template hessians<0, true, false>(in, out);
4052  if (dim > 1)
4053  evaluator.template hessians<1, true, false>(out, out);
4054  if (dim > 2)
4055  evaluator.template hessians<2, true, false>(out, out);
4056  }
4057  for (unsigned int q = 0; q < dofs_per_component; ++q)
4058  {
4059  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
4060  for (unsigned int d = 0; d < n_components; ++d)
4061  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
4062  }
4063  for (unsigned int d = 0; d < n_components; ++d)
4064  {
4065  Number *out = out_array + d * dofs_per_component;
4066  if (dim > 2)
4067  evaluator.template hessians<2, false, false>(out, out);
4068  if (dim > 1)
4069  evaluator.template hessians<1, false, false>(out, out);
4070  evaluator.template hessians<0, false, false>(out, out);
4071  }
4072  return false;
4073  }
4074 
4075  template <int fe_degree, int = 0>
4076  static bool
4077  run(const unsigned int n_components,
4078  const FEEvaluationBaseData<dim,
4079  typename Number::value_type,
4080  false,
4081  Number> &fe_eval,
4082  const Number * in_array,
4083  Number * out_array,
4084  typename std::enable_if<fe_degree == -1>::type * = nullptr)
4085  {
4086  static_assert(fe_degree == -1, "Only usable for degree -1");
4087  const unsigned int dofs_per_component =
4088  fe_eval.get_shape_info().dofs_per_component_on_cell;
4089 
4090  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
4091 
4092  internal::
4093  EvaluatorTensorProduct<internal::evaluate_general, dim, 0, 0, Number>
4094  evaluator(fe_eval.get_shape_info().data.front().inverse_shape_values,
4097  fe_eval.get_shape_info().data.front().fe_degree + 1,
4098  fe_eval.get_shape_info().data.front().fe_degree + 1);
4099 
4100  for (unsigned int d = 0; d < n_components; ++d)
4101  {
4102  const Number *in = in_array + d * dofs_per_component;
4103  Number * out = out_array + d * dofs_per_component;
4104  // Need to select 'apply' method with hessian slot because values
4105  // assume symmetries that do not exist in the inverse shapes
4106  evaluator.template values<0, true, false>(in, out);
4107  if (dim > 1)
4108  evaluator.template values<1, true, false>(out, out);
4109  if (dim > 2)
4110  evaluator.template values<2, true, false>(out, out);
4111  }
4112  for (unsigned int q = 0; q < dofs_per_component; ++q)
4113  {
4114  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
4115  for (unsigned int d = 0; d < n_components; ++d)
4116  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
4117  }
4118  for (unsigned int d = 0; d < n_components; ++d)
4119  {
4120  Number *out = out_array + d * dofs_per_component;
4121  if (dim > 2)
4122  evaluator.template values<2, false, false>(out, out);
4123  if (dim > 1)
4124  evaluator.template values<1, false, false>(out, out);
4125  evaluator.template values<0, false, false>(out, out);
4126  }
4127  return false;
4128  }
4129  };
4130 
4131 
4132 
4137  template <int dim, typename Number>
4139  {
4140  template <int fe_degree, int = 0>
4141  static bool
4142  run(const unsigned int n_desired_components,
4143  const AlignedVector<Number> &inverse_shape,
4144  const AlignedVector<Number> &inverse_coefficients,
4145  const Number * in_array,
4146  Number * out_array,
4147  typename std::enable_if<fe_degree != -1>::type * = nullptr)
4148  {
4149  constexpr unsigned int dofs_per_component =
4150  Utilities::pow(fe_degree + 1, dim);
4151  Assert(inverse_coefficients.size() > 0 &&
4152  inverse_coefficients.size() % dofs_per_component == 0,
4153  ExcMessage(
4154  "Expected diagonal to be a multiple of scalar dof per cells"));
4155  if (inverse_coefficients.size() != dofs_per_component)
4156  AssertDimension(n_desired_components * dofs_per_component,
4157  inverse_coefficients.size());
4158 
4159  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
4160 
4162  dim,
4163  fe_degree + 1,
4164  fe_degree + 1,
4165  Number>
4166  evaluator(AlignedVector<Number>(),
4168  inverse_shape);
4169 
4170  const unsigned int shift_coefficient =
4171  inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
4172  0;
4173  const Number *inv_coefficient = inverse_coefficients.data();
4174  for (unsigned int d = 0; d < n_desired_components; ++d)
4175  {
4176  const Number *in = in_array + d * dofs_per_component;
4177  Number * out = out_array + d * dofs_per_component;
4178  // Need to select 'apply' method with hessian slot because values
4179  // assume symmetries that do not exist in the inverse shapes
4180  evaluator.template hessians<0, true, false>(in, out);
4181  if (dim > 1)
4182  evaluator.template hessians<1, true, false>(out, out);
4183  if (dim > 2)
4184  evaluator.template hessians<2, true, false>(out, out);
4185 
4186  for (unsigned int q = 0; q < dofs_per_component; ++q)
4187  out[q] *= inv_coefficient[q];
4188 
4189  if (dim > 2)
4190  evaluator.template hessians<2, false, false>(out, out);
4191  if (dim > 1)
4192  evaluator.template hessians<1, false, false>(out, out);
4193  evaluator.template hessians<0, false, false>(out, out);
4194 
4195  inv_coefficient += shift_coefficient;
4196  }
4197  return false;
4198  }
4199 
4203  template <int fe_degree, int = 0>
4204  static bool
4205  run(const unsigned int,
4206  const AlignedVector<Number> &,
4207  const AlignedVector<Number> &,
4208  const Number *,
4209  Number *,
4210  typename std::enable_if<fe_degree == -1>::type * = nullptr)
4211  {
4212  static_assert(fe_degree == -1, "Only usable for degree -1");
4213  Assert(false, ExcNotImplemented());
4214  return false;
4215  }
4216  };
4217 
4218 
4219 
4224  template <int dim, typename Number>
4226  {
4227  template <int fe_degree, int = 0>
4228  static bool
4229  run(const unsigned int n_desired_components,
4230  const AlignedVector<Number> &inverse_shape,
4231  const Number * in_array,
4232  Number * out_array,
4233  typename std::enable_if<fe_degree != -1>::type * = nullptr)
4234  {
4235  constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim);
4237  dim,
4238  fe_degree + 1,
4239  fe_degree + 1,
4240  Number>
4241  evaluator(AlignedVector<Number>(),
4243  inverse_shape);
4244 
4245  for (unsigned int d = 0; d < n_desired_components; ++d)
4246  {
4247  const Number *in = in_array + d * dofs_per_cell;
4248  Number * out = out_array + d * dofs_per_cell;
4249 
4250  if (dim == 3)
4251  {
4252  evaluator.template hessians<2, false, false>(in, out);
4253  evaluator.template hessians<1, false, false>(out, out);
4254  evaluator.template hessians<0, false, false>(out, out);
4255  }
4256  if (dim == 2)
4257  {
4258  evaluator.template hessians<1, false, false>(in, out);
4259  evaluator.template hessians<0, false, false>(out, out);
4260  }
4261  if (dim == 1)
4262  evaluator.template hessians<0, false, false>(in, out);
4263  }
4264  return false;
4265  }
4266 
4267  template <int fe_degree, int = 0>
4268  static bool
4269  run(const unsigned int,
4270  const AlignedVector<Number> &,
4271  const Number *,
4272  Number *,
4273  typename std::enable_if<fe_degree == -1>::type * = nullptr)
4274  {
4275  static_assert(fe_degree == -1, "Only usable for degree -1");
4276  Assert(false, ExcNotImplemented());
4277  return false;
4278  }
4279  };
4280 
4281 } // end of namespace internal
4282 
4283 
4285 
4286 #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:1623
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:1691
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:1576
static bool run(const unsigned int, const AlignedVector< Number > &, const Number *, Number *, typename std::enable_if< fe_degree==-1 >::type *=nullptr)
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:1045
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:1466
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:372
void load(const Number *ptr)
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:94
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:371
T min(const T &t, const MPI_Comm &mpi_communicator)
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:375
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
T max(const T &t, const MPI_Comm &mpi_communicator)
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:134
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)