Reference documentation for deal.II version Git e87e13460a 2021-10-28 06:28:57 +0200
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
evaluation_kernels.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_evaluation_kernels_h
18 #define dealii_matrix_free_evaluation_kernels_h
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/utilities.h>
24 
25 #include <deal.II/fe/fe_q.h>
26 
33 
34 
36 
37 
38 // forward declaration
39 template <int, typename, bool, typename>
41 
42 
43 
44 namespace internal
45 {
46  // Select evaluator type from element shape function type
47  template <MatrixFreeFunctions::ElementType element, bool is_long>
49  {};
50 
51  template <bool is_long>
52  struct EvaluatorSelector<MatrixFreeFunctions::tensor_general, is_long>
53  {
54  static const EvaluatorVariant variant = evaluate_general;
55  };
56 
57  template <>
58  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, false>
59  {
60  static const EvaluatorVariant variant = evaluate_symmetric;
61  };
62 
63  template <>
64  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, true>
65  {
66  static const EvaluatorVariant variant = evaluate_evenodd;
67  };
68 
69  template <bool is_long>
70  struct EvaluatorSelector<MatrixFreeFunctions::truncated_tensor, is_long>
71  {
72  static const EvaluatorVariant variant = evaluate_general;
73  };
74 
75  template <>
76  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,
77  false>
78  {
79  static const EvaluatorVariant variant = evaluate_general;
80  };
81 
82  template <>
83  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0, true>
84  {
85  static const EvaluatorVariant variant = evaluate_evenodd;
86  };
87 
88  template <bool is_long>
89  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_collocation,
90  is_long>
91  {
92  static const EvaluatorVariant variant = evaluate_evenodd;
93  };
94 
95 
96 
113  template <MatrixFreeFunctions::ElementType type,
114  int dim,
115  int fe_degree,
116  int n_q_points_1d,
117  typename Number>
119  {
120  static void
121  evaluate(const unsigned int n_components,
122  const EvaluationFlags::EvaluationFlags evaluation_flag,
123  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
124  const Number * values_dofs_actual,
125  Number * values_quad,
126  Number * gradients_quad,
127  Number * hessians_quad,
128  Number * scratch_data);
129 
130  static void
131  integrate(const unsigned int n_components,
132  const EvaluationFlags::EvaluationFlags integration_flag,
133  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
134  Number * values_dofs_actual,
135  Number * values_quad,
136  Number * gradients_quad,
137  Number * hessians_quad,
138  Number * scratch_data,
139  const bool add_into_values_array);
140  };
141 
142 
143 
148  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
149  struct FEEvaluationImpl<MatrixFreeFunctions::tensor_none,
150  dim,
151  fe_degree,
152  n_q_points_1d,
153  Number>
154  {
155  static void
156  evaluate(const unsigned int n_components,
157  const EvaluationFlags::EvaluationFlags evaluation_flag,
158  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
159  const Number * values_dofs_actual,
160  Number * values_quad,
161  Number * gradients_quad,
162  Number * hessians_quad,
163  Number * scratch_data);
164 
165  static void
166  integrate(const unsigned int n_components,
167  const EvaluationFlags::EvaluationFlags integration_flag,
168  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
169  Number * values_dofs_actual,
170  Number * values_quad,
171  Number * gradients_quad,
172  Number * hessians_quad,
173  Number * scratch_data,
174  const bool add_into_values_array);
175  };
176 
177 
178 
179  template <MatrixFreeFunctions::ElementType type,
180  int dim,
181  int fe_degree,
182  int n_q_points_1d,
183  typename Number>
184  inline void
186  const unsigned int n_components,
187  const EvaluationFlags::EvaluationFlags evaluation_flag,
188  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
189  const Number * values_dofs_actual,
190  Number * values_quad,
191  Number * gradients_quad,
192  Number * hessians_quad,
193  Number * scratch_data)
194  {
195  if (evaluation_flag == EvaluationFlags::nothing)
196  return;
197 
198  const EvaluatorVariant variant =
200  using Eval = EvaluatorTensorProduct<variant,
201  dim,
202  fe_degree + 1,
203  n_q_points_1d,
204  Number>;
205  Eval eval(variant == evaluate_evenodd ?
206  shape_info.data.front().shape_values_eo :
207  shape_info.data.front().shape_values,
208  variant == evaluate_evenodd ?
209  shape_info.data.front().shape_gradients_eo :
210  shape_info.data.front().shape_gradients,
211  variant == evaluate_evenodd ?
212  shape_info.data.front().shape_hessians_eo :
213  shape_info.data.front().shape_hessians,
214  shape_info.data.front().fe_degree + 1,
215  shape_info.data.front().n_q_points_1d);
216 
217  const unsigned int temp_size =
218  Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
219  0 :
220  (Eval::n_rows_of_product > Eval::n_columns_of_product ?
221  Eval::n_rows_of_product :
222  Eval::n_columns_of_product);
223  Number *temp1 = scratch_data;
224  Number *temp2;
225  if (temp_size == 0)
226  {
227  temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
228  shape_info.data.front().fe_degree + 1),
229  Utilities::fixed_power<dim>(
230  shape_info.data.front().n_q_points_1d));
231  }
232  else
233  {
234  temp2 = temp1 + temp_size;
235  }
236 
237  const unsigned int n_q_points =
238  temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
239  const unsigned int dofs_per_comp =
241  Utilities::fixed_power<dim>(shape_info.data.front().fe_degree + 1) :
242  shape_info.dofs_per_component_on_cell;
243  const Number *values_dofs = values_dofs_actual;
245  {
246  Number *values_dofs_tmp =
247  scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
248  shape_info.n_q_points));
249  const int degree =
250  fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
251  for (unsigned int c = 0; c < n_components; ++c)
252  for (int i = 0, count_p = 0, count_q = 0;
253  i < (dim > 2 ? degree + 1 : 1);
254  ++i)
255  {
256  for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
257  {
258  for (int k = 0; k < degree + 1 - j - i;
259  ++k, ++count_p, ++count_q)
260  values_dofs_tmp[c * dofs_per_comp + count_q] =
261  values_dofs_actual
262  [c * shape_info.dofs_per_component_on_cell + count_p];
263  for (int k = degree + 1 - j - i; k < degree + 1;
264  ++k, ++count_q)
265  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
266  }
267  for (int j = degree + 1 - i; j < degree + 1; ++j)
268  for (int k = 0; k < degree + 1; ++k, ++count_q)
269  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
270  }
271  values_dofs = values_dofs_tmp;
272  }
273 
274  switch (dim)
275  {
276  case 1:
277  for (unsigned int c = 0; c < n_components; ++c)
278  {
279  if (evaluation_flag & EvaluationFlags::values)
280  eval.template values<0, true, false>(values_dofs, values_quad);
281  if (evaluation_flag & EvaluationFlags::gradients)
282  eval.template gradients<0, true, false>(values_dofs,
283  gradients_quad);
284  if (evaluation_flag & EvaluationFlags::hessians)
285  eval.template hessians<0, true, false>(values_dofs,
286  hessians_quad);
287 
288  // advance the next component in 1D array
289  values_dofs += dofs_per_comp;
290  values_quad += n_q_points;
291  gradients_quad += n_q_points;
292  hessians_quad += n_q_points;
293  }
294  break;
295 
296  case 2:
297  for (unsigned int c = 0; c < n_components; ++c)
298  {
299  // grad x
300  if (evaluation_flag & EvaluationFlags::gradients)
301  {
302  eval.template gradients<0, true, false>(values_dofs, temp1);
303  eval.template values<1, true, false>(temp1, gradients_quad);
304  }
305  if (evaluation_flag & EvaluationFlags::hessians)
306  {
307  // grad xy
308  if (!(evaluation_flag & EvaluationFlags::gradients))
309  eval.template gradients<0, true, false>(values_dofs, temp1);
310  eval.template gradients<1, true, false>(temp1,
311  hessians_quad +
312  2 * n_q_points);
313 
314  // grad xx
315  eval.template hessians<0, true, false>(values_dofs, temp1);
316  eval.template values<1, true, false>(temp1, hessians_quad);
317  }
318 
319  // grad y
320  eval.template values<0, true, false>(values_dofs, temp1);
321  if (evaluation_flag & EvaluationFlags::gradients)
322  eval.template gradients<1, true, false>(temp1,
323  gradients_quad +
324  n_q_points);
325 
326  // grad yy
327  if (evaluation_flag & EvaluationFlags::hessians)
328  eval.template hessians<1, true, false>(temp1,
329  hessians_quad +
330  n_q_points);
331 
332  // val: can use values applied in x
333  if (evaluation_flag & EvaluationFlags::values)
334  eval.template values<1, true, false>(temp1, values_quad);
335 
336  // advance to the next component in 1D array
337  values_dofs += dofs_per_comp;
338  values_quad += n_q_points;
339  gradients_quad += 2 * n_q_points;
340  hessians_quad += 3 * n_q_points;
341  }
342  break;
343 
344  case 3:
345  for (unsigned int c = 0; c < n_components; ++c)
346  {
347  if (evaluation_flag & EvaluationFlags::gradients)
348  {
349  // grad x
350  eval.template gradients<0, true, false>(values_dofs, temp1);
351  eval.template values<1, true, false>(temp1, temp2);
352  eval.template values<2, true, false>(temp2, gradients_quad);
353  }
354 
355  if (evaluation_flag & EvaluationFlags::hessians)
356  {
357  // grad xz
358  if (!(evaluation_flag & EvaluationFlags::gradients))
359  {
360  eval.template gradients<0, true, false>(values_dofs,
361  temp1);
362  eval.template values<1, true, false>(temp1, temp2);
363  }
364  eval.template gradients<2, true, false>(temp2,
365  hessians_quad +
366  4 * n_q_points);
367 
368  // grad xy
369  eval.template gradients<1, true, false>(temp1, temp2);
370  eval.template values<2, true, false>(temp2,
371  hessians_quad +
372  3 * n_q_points);
373 
374  // grad xx
375  eval.template hessians<0, true, false>(values_dofs, temp1);
376  eval.template values<1, true, false>(temp1, temp2);
377  eval.template values<2, true, false>(temp2, hessians_quad);
378  }
379 
380  // grad y
381  eval.template values<0, true, false>(values_dofs, temp1);
382  if (evaluation_flag & EvaluationFlags::gradients)
383  {
384  eval.template gradients<1, true, false>(temp1, temp2);
385  eval.template values<2, true, false>(temp2,
386  gradients_quad +
387  n_q_points);
388  }
389 
390  if (evaluation_flag & EvaluationFlags::hessians)
391  {
392  // grad yz
393  if (!(evaluation_flag & EvaluationFlags::gradients))
394  eval.template gradients<1, true, false>(temp1, temp2);
395  eval.template gradients<2, true, false>(temp2,
396  hessians_quad +
397  5 * n_q_points);
398 
399  // grad yy
400  eval.template hessians<1, true, false>(temp1, temp2);
401  eval.template values<2, true, false>(temp2,
402  hessians_quad +
403  n_q_points);
404  }
405 
406  // grad z: can use the values applied in x direction stored in
407  // temp1
408  eval.template values<1, true, false>(temp1, temp2);
409  if (evaluation_flag & EvaluationFlags::gradients)
410  eval.template gradients<2, true, false>(temp2,
411  gradients_quad +
412  2 * n_q_points);
413 
414  // grad zz: can use the values applied in x and y direction stored
415  // in temp2
416  if (evaluation_flag & EvaluationFlags::hessians)
417  eval.template hessians<2, true, false>(temp2,
418  hessians_quad +
419  2 * n_q_points);
420 
421  // val: can use the values applied in x & y direction stored in
422  // temp2
423  if (evaluation_flag & EvaluationFlags::values)
424  eval.template values<2, true, false>(temp2, values_quad);
425 
426  // advance to the next component in 1D array
427  values_dofs += dofs_per_comp;
428  values_quad += n_q_points;
429  gradients_quad += 3 * n_q_points;
430  hessians_quad += 6 * n_q_points;
431  }
432  break;
433 
434  default:
435  AssertThrow(false, ExcNotImplemented());
436  }
437 
438  // case additional dof for FE_Q_DG0: add values; gradients and second
439  // derivatives evaluate to zero
441  (evaluation_flag & EvaluationFlags::values))
442  {
443  values_quad -= n_components * n_q_points;
444  values_dofs -= n_components * dofs_per_comp;
445  for (unsigned int c = 0; c < n_components; ++c)
446  for (unsigned int q = 0; q < shape_info.n_q_points; ++q)
447  values_quad[c * shape_info.n_q_points + q] +=
448  values_dofs[(c + 1) * shape_info.dofs_per_component_on_cell - 1];
449  }
450  }
451 
452 
453 
454  template <MatrixFreeFunctions::ElementType type,
455  int dim,
456  int fe_degree,
457  int n_q_points_1d,
458  typename Number>
459  inline void
461  const unsigned int n_components,
462  const EvaluationFlags::EvaluationFlags integration_flag,
463  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
464  Number * values_dofs_actual,
465  Number * values_quad,
466  Number * gradients_quad,
467  Number * hessians_quad,
468  Number * scratch_data,
469  const bool add_into_values_array)
470  {
471  const EvaluatorVariant variant =
473  using Eval = EvaluatorTensorProduct<variant,
474  dim,
475  fe_degree + 1,
476  n_q_points_1d,
477  Number>;
478  Eval eval(variant == evaluate_evenodd ?
479  shape_info.data.front().shape_values_eo :
480  shape_info.data.front().shape_values,
481  variant == evaluate_evenodd ?
482  shape_info.data.front().shape_gradients_eo :
483  shape_info.data.front().shape_gradients,
484  variant == evaluate_evenodd ?
485  shape_info.data.front().shape_hessians_eo :
486  shape_info.data.front().shape_hessians,
487  shape_info.data.front().fe_degree + 1,
488  shape_info.data.front().n_q_points_1d);
489 
490  const unsigned int temp_size =
491  Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
492  0 :
493  (Eval::n_rows_of_product > Eval::n_columns_of_product ?
494  Eval::n_rows_of_product :
495  Eval::n_columns_of_product);
496  Number *temp1 = scratch_data;
497  Number *temp2;
498  if (temp_size == 0)
499  {
500  temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
501  shape_info.data.front().fe_degree + 1),
502  Utilities::fixed_power<dim>(
503  shape_info.data.front().n_q_points_1d));
504  }
505  else
506  {
507  temp2 = temp1 + temp_size;
508  }
509 
510  const unsigned int n_q_points =
511  temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
512  const unsigned int dofs_per_comp =
514  Utilities::fixed_power<dim>(shape_info.data.front().fe_degree + 1) :
515  shape_info.dofs_per_component_on_cell;
516  // expand dof_values to tensor product for truncated tensor products
517  Number *values_dofs =
519  scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
520  shape_info.n_q_points)) :
521  values_dofs_actual;
522 
523  switch (dim)
524  {
525  case 1:
526  for (unsigned int c = 0; c < n_components; ++c)
527  {
528  if (integration_flag & EvaluationFlags::values)
529  {
530  if (add_into_values_array == false)
531  eval.template values<0, false, false>(values_quad,
532  values_dofs);
533  else
534  eval.template values<0, false, true>(values_quad,
535  values_dofs);
536  }
537  if (integration_flag & EvaluationFlags::gradients)
538  {
539  if (integration_flag & EvaluationFlags::values ||
540  add_into_values_array == true)
541  eval.template gradients<0, false, true>(gradients_quad,
542  values_dofs);
543  else
544  eval.template gradients<0, false, false>(gradients_quad,
545  values_dofs);
546  }
547  if (integration_flag & EvaluationFlags::hessians)
548  {
549  if (integration_flag & EvaluationFlags::values ||
550  integration_flag & EvaluationFlags::gradients ||
551  add_into_values_array == true)
552  eval.template hessians<0, false, true>(hessians_quad,
553  values_dofs);
554  else
555  eval.template hessians<0, false, false>(hessians_quad,
556  values_dofs);
557  }
558 
559  // advance to the next component in 1D array
560  values_dofs += dofs_per_comp;
561  values_quad += n_q_points;
562  gradients_quad += n_q_points;
563  hessians_quad += n_q_points;
564  }
565  break;
566 
567  case 2:
568  for (unsigned int c = 0; c < n_components; ++c)
569  {
570  if ((integration_flag & EvaluationFlags::values) &&
571  !(integration_flag & EvaluationFlags::gradients))
572  {
573  eval.template values<1, false, false>(values_quad, temp1);
574  if (add_into_values_array == false)
575  eval.template values<0, false, false>(temp1, values_dofs);
576  else
577  eval.template values<0, false, true>(temp1, values_dofs);
578  }
579  if (integration_flag & EvaluationFlags::gradients)
580  {
581  eval.template gradients<1, false, false>(gradients_quad +
582  n_q_points,
583  temp1);
584  if (integration_flag & EvaluationFlags::values)
585  eval.template values<1, false, true>(values_quad, temp1);
586  if (add_into_values_array == false)
587  eval.template values<0, false, false>(temp1, values_dofs);
588  else
589  eval.template values<0, false, true>(temp1, values_dofs);
590  eval.template values<1, false, false>(gradients_quad, temp1);
591  eval.template gradients<0, false, true>(temp1, values_dofs);
592  }
593  if (integration_flag & EvaluationFlags::hessians)
594  {
595  // grad xx
596  eval.template values<1, false, false>(hessians_quad, temp1);
597 
598  if (integration_flag & EvaluationFlags::values ||
599  integration_flag & EvaluationFlags::gradients ||
600  add_into_values_array == true)
601  eval.template hessians<0, false, true>(temp1, values_dofs);
602  else
603  eval.template hessians<0, false, false>(temp1, values_dofs);
604 
605  // grad yy
606  eval.template hessians<1, false, false>(hessians_quad +
607  n_q_points,
608  temp1);
609  eval.template values<0, false, true>(temp1, values_dofs);
610 
611  // grad xy
612  eval.template gradients<1, false, false>(hessians_quad +
613  2 * n_q_points,
614  temp1);
615  eval.template gradients<0, false, true>(temp1, values_dofs);
616  }
617 
618  // advance to the next component in 1D array
619  values_dofs += dofs_per_comp;
620  values_quad += n_q_points;
621  gradients_quad += 2 * n_q_points;
622  hessians_quad += 3 * n_q_points;
623  }
624  break;
625 
626  case 3:
627  for (unsigned int c = 0; c < n_components; ++c)
628  {
629  if ((integration_flag & EvaluationFlags::values) &&
630  !(integration_flag & EvaluationFlags::gradients))
631  {
632  eval.template values<2, false, false>(values_quad, temp1);
633  eval.template values<1, false, false>(temp1, temp2);
634  if (add_into_values_array == false)
635  eval.template values<0, false, false>(temp2, values_dofs);
636  else
637  eval.template values<0, false, true>(temp2, values_dofs);
638  }
639  if (integration_flag & EvaluationFlags::gradients)
640  {
641  eval.template gradients<2, false, false>(gradients_quad +
642  2 * n_q_points,
643  temp1);
644  if (integration_flag & EvaluationFlags::values)
645  eval.template values<2, false, true>(values_quad, temp1);
646  eval.template values<1, false, false>(temp1, temp2);
647  eval.template values<2, false, false>(gradients_quad +
648  n_q_points,
649  temp1);
650  eval.template gradients<1, false, true>(temp1, temp2);
651  if (add_into_values_array == false)
652  eval.template values<0, false, false>(temp2, values_dofs);
653  else
654  eval.template values<0, false, true>(temp2, values_dofs);
655  eval.template values<2, false, false>(gradients_quad, temp1);
656  eval.template values<1, false, false>(temp1, temp2);
657  eval.template gradients<0, false, true>(temp2, values_dofs);
658  }
659  if (integration_flag & EvaluationFlags::hessians)
660  {
661  // grad xx
662  eval.template values<2, false, false>(hessians_quad, temp1);
663  eval.template values<1, false, false>(temp1, temp2);
664 
665  if (integration_flag & EvaluationFlags::values ||
666  integration_flag & EvaluationFlags::gradients ||
667  add_into_values_array == true)
668  eval.template hessians<0, false, true>(temp2, values_dofs);
669  else
670  eval.template hessians<0, false, false>(temp2, values_dofs);
671 
672  // grad yy
673  eval.template values<2, false, false>(hessians_quad +
674  n_q_points,
675  temp1);
676  eval.template hessians<1, false, false>(temp1, temp2);
677  eval.template values<0, false, true>(temp2, values_dofs);
678 
679  // grad zz
680  eval.template hessians<2, false, false>(hessians_quad +
681  2 * n_q_points,
682  temp1);
683  eval.template values<1, false, false>(temp1, temp2);
684  eval.template values<0, false, true>(temp2, values_dofs);
685 
686  // grad xy
687  eval.template values<2, false, false>(hessians_quad +
688  3 * n_q_points,
689  temp1);
690  eval.template gradients<1, false, false>(temp1, temp2);
691  eval.template gradients<0, false, true>(temp2, values_dofs);
692 
693  // grad xz
694  eval.template gradients<2, false, false>(hessians_quad +
695  4 * n_q_points,
696  temp1);
697  eval.template values<1, false, false>(temp1, temp2);
698  eval.template gradients<0, false, true>(temp2, values_dofs);
699 
700  // grad yz
701  eval.template gradients<2, false, false>(hessians_quad +
702  5 * n_q_points,
703  temp1);
704  eval.template gradients<1, false, false>(temp1, temp2);
705  eval.template values<0, false, true>(temp2, values_dofs);
706  }
707 
708  // advance to the next component in 1D array
709  values_dofs += dofs_per_comp;
710  values_quad += n_q_points;
711  gradients_quad += 3 * n_q_points;
712  hessians_quad += 6 * n_q_points;
713  }
714  break;
715 
716  default:
717  AssertThrow(false, ExcNotImplemented());
718  }
719 
720  // case FE_Q_DG0: add values, gradients and second derivatives are zero
722  {
723  values_dofs -= n_components * dofs_per_comp -
724  shape_info.dofs_per_component_on_cell + 1;
725  values_quad -= n_components * n_q_points;
726  if (integration_flag & EvaluationFlags::values)
727  for (unsigned int c = 0; c < n_components; ++c)
728  {
729  values_dofs[0] = values_quad[0];
730  for (unsigned int q = 1; q < shape_info.n_q_points; ++q)
731  values_dofs[0] += values_quad[q];
732  values_dofs += dofs_per_comp;
733  values_quad += n_q_points;
734  }
735  else
736  {
737  for (unsigned int c = 0; c < n_components; ++c)
738  values_dofs[c * shape_info.dofs_per_component_on_cell] = Number();
739  values_dofs += n_components * shape_info.dofs_per_component_on_cell;
740  }
741  }
742 
744  {
745  values_dofs -= dofs_per_comp * n_components;
746  const int degree =
747  fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
748  for (unsigned int c = 0; c < n_components; ++c)
749  for (int i = 0, count_p = 0, count_q = 0;
750  i < (dim > 2 ? degree + 1 : 1);
751  ++i)
752  {
753  for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
754  {
755  for (int k = 0; k < degree + 1 - j - i;
756  ++k, ++count_p, ++count_q)
757  values_dofs_actual[c *
758  shape_info.dofs_per_component_on_cell +
759  count_p] =
760  values_dofs[c * dofs_per_comp + count_q];
761  count_q += j + i;
762  }
763  count_q += i * (degree + 1);
764  }
765  }
766  }
767 
768 
769 
770  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
771  inline void
774  dim,
775  fe_degree,
776  n_q_points_1d,
777  Number>::evaluate(const unsigned int n_components,
778  const EvaluationFlags::EvaluationFlags evaluation_flag,
779  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
780  const Number *values_dofs_actual,
781  Number * values_quad,
782  Number * gradients_quad,
783  Number * hessians_quad,
784  Number * scratch_data)
785  {
786  (void)scratch_data;
787 
788  const unsigned int n_dofs = shape_info.dofs_per_component_on_cell;
789  const unsigned int n_q_points = shape_info.n_q_points;
790 
791  using Eval =
793 
794  if (evaluation_flag & EvaluationFlags::values)
795  {
796  const auto shape_values = shape_info.data.front().shape_values.data();
797  auto values_quad_ptr = values_quad;
798  auto values_dofs_actual_ptr = values_dofs_actual;
799 
800  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
801  for (unsigned int c = 0; c < n_components; ++c)
802  {
803  eval.template values<0, true, false>(values_dofs_actual_ptr,
804  values_quad_ptr);
805 
806  values_quad_ptr += n_q_points;
807  values_dofs_actual_ptr += n_dofs;
808  }
809  }
810 
811  if (evaluation_flag & EvaluationFlags::gradients)
812  {
813  const auto shape_gradients =
814  shape_info.data.front().shape_gradients.data();
815  auto gradients_quad_ptr = gradients_quad;
816  auto values_dofs_actual_ptr = values_dofs_actual;
817 
818  for (unsigned int c = 0; c < n_components; ++c)
819  {
820  for (unsigned int d = 0; d < dim; ++d)
821  {
822  Eval eval(nullptr,
823  shape_gradients + n_q_points * n_dofs * d,
824  nullptr,
825  n_dofs,
826  n_q_points);
827 
828  eval.template gradients<0, true, false>(values_dofs_actual_ptr,
829  gradients_quad_ptr);
830 
831  gradients_quad_ptr += n_q_points;
832  }
833  values_dofs_actual_ptr += n_dofs;
834  }
835  }
836 
837  if (evaluation_flag & EvaluationFlags::hessians)
838  {
839  Assert(false, ExcNotImplemented());
840  (void)hessians_quad;
841  }
842  }
843 
844 
845 
846  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
847  inline void
850  dim,
851  fe_degree,
852  n_q_points_1d,
853  Number>::integrate(const unsigned int n_components,
854  const EvaluationFlags::EvaluationFlags integration_flag,
855  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
856  Number * values_dofs_actual,
857  Number * values_quad,
858  Number * gradients_quad,
859  Number * hessians_quad,
860  Number * scratch_data,
861  const bool add_into_values_array)
862  {
863  // TODO: implement hessians
864  (void)hessians_quad;
865  AssertThrow(!(integration_flag & EvaluationFlags::hessians),
867  (void)scratch_data;
868 
869  const unsigned int n_dofs = shape_info.dofs_per_component_on_cell;
870  const unsigned int n_q_points = shape_info.n_q_points;
871 
872  using Eval =
874 
875  if (integration_flag & EvaluationFlags::values)
876  {
877  const auto shape_values = shape_info.data.front().shape_values.data();
878  auto values_quad_ptr = values_quad;
879  auto values_dofs_actual_ptr = values_dofs_actual;
880 
881  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
882  for (unsigned int c = 0; c < n_components; ++c)
883  {
884  if (add_into_values_array == false)
885  eval.template values<0, false, false>(values_quad_ptr,
886  values_dofs_actual_ptr);
887  else
888  eval.template values<0, false, true>(values_quad_ptr,
889  values_dofs_actual_ptr);
890 
891  values_quad_ptr += n_q_points;
892  values_dofs_actual_ptr += n_dofs;
893  }
894  }
895 
896  if (integration_flag & EvaluationFlags::gradients)
897  {
898  const auto shape_gradients =
899  shape_info.data.front().shape_gradients.data();
900  auto gradients_quad_ptr = gradients_quad;
901  auto values_dofs_actual_ptr = values_dofs_actual;
902 
903  for (unsigned int c = 0; c < n_components; ++c)
904  {
905  for (unsigned int d = 0; d < dim; ++d)
906  {
907  Eval eval(nullptr,
908  shape_gradients + n_q_points * n_dofs * d,
909  nullptr,
910  n_dofs,
911  n_q_points);
912 
913  if ((add_into_values_array == false &&
914  (integration_flag & EvaluationFlags::values) == false) &&
915  d == 0)
916  eval.template gradients<0, false, false>(
917  gradients_quad_ptr, values_dofs_actual_ptr);
918  else
919  eval.template gradients<0, false, true>(
920  gradients_quad_ptr, values_dofs_actual_ptr);
921 
922  gradients_quad_ptr += n_q_points;
923  }
924  values_dofs_actual_ptr += n_dofs;
925  }
926  }
927  }
928 
929 
930 
940  template <EvaluatorVariant variant,
941  EvaluatorQuantity quantity,
942  int dim,
943  int basis_size_1,
944  int basis_size_2,
945  typename Number,
946  typename Number2>
948  {
949  static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
950  "The second dimension must not be smaller than the first");
951 
974 #ifndef DEBUG
976 #endif
977  static void
979  const unsigned int n_components,
980  const AlignedVector<Number2> &transformation_matrix,
981  const Number * values_in,
982  Number * values_out,
983  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
984  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
985  {
986  Assert(
987  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
988  ExcMessage("The second dimension must not be smaller than the first"));
989 
991 
992  // we do recursion until dim==1 or dim==2 and we have
993  // basis_size_1==basis_size_2. The latter optimization increases
994  // optimization possibilities for the compiler but does only work for
995  // aliased pointers if the sizes are equal.
996  constexpr int next_dim =
997  (dim > 2 ||
998  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
999  dim - 1 :
1000  dim;
1001 
1002  EvaluatorTensorProduct<variant,
1003  dim,
1004  basis_size_1,
1005  (basis_size_1 == 0 ? 0 : basis_size_2),
1006  Number,
1007  Number2>
1008  eval_val(transformation_matrix,
1011  basis_size_1_variable,
1012  basis_size_2_variable);
1013  const unsigned int np_1 =
1014  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
1015  const unsigned int np_2 =
1016  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
1017  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
1018  ExcMessage("Cannot transform with 0-point basis"));
1019  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
1020  ExcMessage("Cannot transform with 0-point basis"));
1021 
1022  // run loop backwards to ensure correctness if values_in aliases with
1023  // values_out in case with basis_size_1 < basis_size_2
1024  values_in = values_in + n_components * Utilities::fixed_power<dim>(np_1);
1025  values_out =
1026  values_out + n_components * Utilities::fixed_power<dim>(np_2);
1027  for (unsigned int c = n_components; c != 0; --c)
1028  {
1029  values_in -= Utilities::fixed_power<dim>(np_1);
1030  values_out -= Utilities::fixed_power<dim>(np_2);
1031  if (next_dim < dim)
1032  for (unsigned int q = np_1; q != 0; --q)
1034  variant,
1035  quantity,
1036  next_dim,
1037  basis_size_1,
1038  basis_size_2,
1039  Number,
1040  Number2>::do_forward(1,
1041  transformation_matrix,
1042  values_in +
1043  (q - 1) *
1044  Utilities::fixed_power<next_dim>(np_1),
1045  values_out +
1046  (q - 1) *
1047  Utilities::fixed_power<next_dim>(np_2),
1048  basis_size_1_variable,
1049  basis_size_2_variable);
1050 
1051  // the recursion stops if dim==1 or if dim==2 and
1052  // basis_size_1==basis_size_2 (the latter is used because the
1053  // compiler generates nicer code)
1054  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
1055  {
1056  eval_val.template values<0, true, false>(values_in, values_out);
1057  eval_val.template values<1, true, false>(values_out, values_out);
1058  }
1059  else if (dim == 1)
1060  eval_val.template values<dim - 1, true, false>(values_in,
1061  values_out);
1062  else
1063  eval_val.template values<dim - 1, true, false>(values_out,
1064  values_out);
1065  }
1066  }
1067 
1098 #ifndef DEBUG
1100 #endif
1101  static void
1103  const unsigned int n_components,
1104  const AlignedVector<Number2> &transformation_matrix,
1105  const bool add_into_result,
1106  Number * values_in,
1107  Number * values_out,
1108  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
1109  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
1110  {
1111  Assert(
1112  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
1113  ExcMessage("The second dimension must not be smaller than the first"));
1114  Assert(add_into_result == false || values_in != values_out,
1115  ExcMessage(
1116  "Input and output cannot alias with each other when "
1117  "adding the result of the basis change to existing data"));
1118 
1119  Assert(quantity == EvaluatorQuantity::value ||
1120  quantity == EvaluatorQuantity::hessian,
1121  ExcInternalError());
1122 
1123  constexpr int next_dim =
1124  (dim > 2 ||
1125  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
1126  dim - 1 :
1127  dim;
1128  EvaluatorTensorProduct<variant,
1129  dim,
1130  basis_size_1,
1131  (basis_size_1 == 0 ? 0 : basis_size_2),
1132  Number,
1133  Number2>
1134  eval_val(transformation_matrix,
1135  transformation_matrix,
1136  transformation_matrix,
1137  basis_size_1_variable,
1138  basis_size_2_variable);
1139  const unsigned int np_1 =
1140  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
1141  const unsigned int np_2 =
1142  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
1143  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
1144  ExcMessage("Cannot transform with 0-point basis"));
1145  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
1146  ExcMessage("Cannot transform with 0-point basis"));
1147 
1148  for (unsigned int c = 0; c < n_components; ++c)
1149  {
1150  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
1151  {
1152  if (quantity == EvaluatorQuantity::value)
1153  eval_val.template values<1, false, false>(values_in, values_in);
1154  else
1155  eval_val.template hessians<1, false, false>(values_in,
1156  values_in);
1157 
1158  if (add_into_result)
1159  {
1160  if (quantity == EvaluatorQuantity::value)
1161  eval_val.template values<0, false, true>(values_in,
1162  values_out);
1163  else
1164  eval_val.template hessians<0, false, true>(values_in,
1165  values_out);
1166  }
1167  else
1168  {
1169  if (quantity == EvaluatorQuantity::value)
1170  eval_val.template values<0, false, false>(values_in,
1171  values_out);
1172  else
1173  eval_val.template hessians<0, false, false>(values_in,
1174  values_out);
1175  }
1176  }
1177  else
1178  {
1179  if (dim == 1 && add_into_result)
1180  {
1181  if (quantity == EvaluatorQuantity::value)
1182  eval_val.template values<0, false, true>(values_in,
1183  values_out);
1184  else
1185  eval_val.template hessians<0, false, true>(values_in,
1186  values_out);
1187  }
1188  else if (dim == 1)
1189  {
1190  if (quantity == EvaluatorQuantity::value)
1191  eval_val.template values<0, false, false>(values_in,
1192  values_out);
1193  else
1194  eval_val.template hessians<0, false, false>(values_in,
1195  values_out);
1196  }
1197  else
1198  {
1199  if (quantity == EvaluatorQuantity::value)
1200  eval_val.template values<dim - 1, false, false>(values_in,
1201  values_in);
1202  else
1203  eval_val.template hessians<dim - 1, false, false>(
1204  values_in, values_in);
1205  }
1206  }
1207  if (next_dim < dim)
1208  for (unsigned int q = 0; q < np_1; ++q)
1210  quantity,
1211  next_dim,
1212  basis_size_1,
1213  basis_size_2,
1214  Number,
1215  Number2>::
1216  do_backward(1,
1217  transformation_matrix,
1218  add_into_result,
1219  values_in +
1220  q * Utilities::fixed_power<next_dim>(np_2),
1221  values_out +
1222  q * Utilities::fixed_power<next_dim>(np_1),
1223  basis_size_1_variable,
1224  basis_size_2_variable);
1225 
1226  values_in += Utilities::fixed_power<dim>(np_2);
1227  values_out += Utilities::fixed_power<dim>(np_1);
1228  }
1229  }
1230 
1251  static void
1252  do_mass(const unsigned int n_components,
1253  const AlignedVector<Number2> &transformation_matrix,
1254  const AlignedVector<Number> & coefficients,
1255  const Number * values_in,
1256  Number * scratch_data,
1257  Number * values_out)
1258  {
1259  constexpr int next_dim = dim > 1 ? dim - 1 : dim;
1260  Number * my_scratch =
1261  basis_size_1 != basis_size_2 ? scratch_data : values_out;
1262 
1263  const unsigned int size_per_component = Utilities::pow(basis_size_2, dim);
1264  Assert(coefficients.size() == size_per_component ||
1265  coefficients.size() == n_components * size_per_component,
1266  ExcDimensionMismatch(coefficients.size(), size_per_component));
1267  const unsigned int stride =
1268  coefficients.size() == size_per_component ? 0 : 1;
1269 
1270  for (unsigned int q = basis_size_1; q != 0; --q)
1272  variant,
1274  next_dim,
1275  basis_size_1,
1276  basis_size_2,
1277  Number,
1278  Number2>::do_forward(n_components,
1279  transformation_matrix,
1280  values_in +
1281  (q - 1) *
1282  Utilities::pow(basis_size_1, dim - 1),
1283  my_scratch +
1284  (q - 1) *
1285  Utilities::pow(basis_size_2, dim - 1));
1286  EvaluatorTensorProduct<variant,
1287  dim,
1288  basis_size_1,
1289  basis_size_2,
1290  Number,
1291  Number2>
1292  eval_val(transformation_matrix);
1293  const unsigned int n_inner_blocks =
1294  (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
1295  const unsigned int n_blocks = Utilities::pow(basis_size_2, dim - 1);
1296  for (unsigned int ii = 0; ii < n_blocks; ii += n_inner_blocks)
1297  for (unsigned int c = 0; c < n_components; ++c)
1298  {
1299  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1300  eval_val.template values_one_line<dim - 1, true, false>(
1301  my_scratch + i, my_scratch + i);
1302  for (unsigned int q = 0; q < basis_size_2; ++q)
1303  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1304  my_scratch[i + q * n_blocks + c * size_per_component] *=
1305  coefficients[i + q * n_blocks +
1306  c * stride * size_per_component];
1307  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1308  eval_val.template values_one_line<dim - 1, false, false>(
1309  my_scratch + i, my_scratch + i);
1310  }
1311  for (unsigned int q = 0; q < basis_size_1; ++q)
1313  variant,
1315  next_dim,
1316  basis_size_1,
1317  basis_size_2,
1318  Number,
1319  Number2>::do_backward(n_components,
1320  transformation_matrix,
1321  false,
1322  my_scratch +
1323  q * Utilities::pow(basis_size_2, dim - 1),
1324  values_out +
1325  q * Utilities::pow(basis_size_1, dim - 1));
1326  }
1327  };
1328 
1329 
1330 
1343  template <int dim, int fe_degree, typename Number>
1345  {
1346  static void
1347  evaluate(const unsigned int n_components,
1348  const EvaluationFlags::EvaluationFlags evaluation_flag,
1349  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1350  const Number * values_dofs,
1351  Number * values_quad,
1352  Number * gradients_quad,
1353  Number * hessians_quad,
1354  Number * scratch_data);
1355 
1356  static void
1357  integrate(const unsigned int n_components,
1358  const EvaluationFlags::EvaluationFlags integration_flag,
1359  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1360  Number * values_dofs,
1361  Number * values_quad,
1362  Number * gradients_quad,
1363  Number * hessians_quad,
1364  Number * scratch_data,
1365  const bool add_into_values_array);
1366  };
1367 
1368 
1369 
1370  template <int dim, int fe_degree, typename Number>
1371  inline void
1373  const unsigned int n_components,
1374  const EvaluationFlags::EvaluationFlags evaluation_flag,
1375  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1376  const Number * values_dofs,
1377  Number * values_quad,
1378  Number * gradients_quad,
1379  Number * hessians_quad,
1380  Number *)
1381  {
1383  shape_info.data.front().shape_gradients_collocation_eo.size(),
1384  (fe_degree + 2) / 2 * (fe_degree + 1));
1385 
1387  dim,
1388  fe_degree + 1,
1389  fe_degree + 1,
1390  Number>
1391  eval(AlignedVector<Number>(),
1392  shape_info.data.front().shape_gradients_collocation_eo,
1393  shape_info.data.front().shape_hessians_collocation_eo);
1394  constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
1395 
1396  for (unsigned int c = 0; c < n_components; ++c)
1397  {
1398  if (evaluation_flag & EvaluationFlags::values)
1399  for (unsigned int i = 0; i < n_q_points; ++i)
1400  values_quad[i] = values_dofs[i];
1401  if (evaluation_flag &
1403  {
1404  eval.template gradients<0, true, false>(values_dofs,
1405  gradients_quad);
1406  if (dim > 1)
1407  eval.template gradients<1, true, false>(values_dofs,
1408  gradients_quad +
1409  n_q_points);
1410  if (dim > 2)
1411  eval.template gradients<2, true, false>(values_dofs,
1412  gradients_quad +
1413  2 * n_q_points);
1414  }
1415  if (evaluation_flag & EvaluationFlags::hessians)
1416  {
1417  eval.template hessians<0, true, false>(values_dofs, hessians_quad);
1418  if (dim > 1)
1419  {
1420  eval.template gradients<1, true, false>(gradients_quad,
1421  hessians_quad +
1422  dim * n_q_points);
1423  eval.template hessians<1, true, false>(values_dofs,
1424  hessians_quad +
1425  n_q_points);
1426  }
1427  if (dim > 2)
1428  {
1429  eval.template gradients<2, true, false>(gradients_quad,
1430  hessians_quad +
1431  4 * n_q_points);
1432  eval.template gradients<2, true, false>(
1433  gradients_quad + n_q_points, hessians_quad + 5 * n_q_points);
1434  eval.template hessians<2, true, false>(values_dofs,
1435  hessians_quad +
1436  2 * n_q_points);
1437  }
1438  hessians_quad += (dim * (dim + 1)) / 2 * n_q_points;
1439  }
1440  gradients_quad += dim * n_q_points;
1441  values_quad += n_q_points;
1442  values_dofs += n_q_points;
1443  }
1444  }
1445 
1446 
1447 
1448  template <int dim, int fe_degree, typename Number>
1449  inline void
1451  const unsigned int n_components,
1452  const EvaluationFlags::EvaluationFlags integration_flag,
1453  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1454  Number * values_dofs,
1455  Number * values_quad,
1456  Number * gradients_quad,
1457  Number * hessians_quad,
1458  Number * scratch_data,
1459  const bool add_into_values_array)
1460  {
1462  shape_info.data.front().shape_gradients_collocation_eo.size(),
1463  (fe_degree + 2) / 2 * (fe_degree + 1));
1464 
1466  dim,
1467  fe_degree + 1,
1468  fe_degree + 1,
1469  Number>
1470  eval(AlignedVector<Number>(),
1471  shape_info.data.front().shape_gradients_collocation_eo,
1472  shape_info.data.front().shape_hessians_collocation_eo);
1473  constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
1474  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
1475 
1476  for (unsigned int c = 0; c < n_components; ++c)
1477  {
1478  if (integration_flag & EvaluationFlags::values)
1479  {
1480  if (add_into_values_array == false)
1481  for (unsigned int i = 0; i < n_q_points; ++i)
1482  values_dofs[i] = values_quad[i];
1483  else
1484  for (unsigned int i = 0; i < n_q_points; ++i)
1485  values_dofs[i] += values_quad[i];
1486  }
1487  if (integration_flag & EvaluationFlags::gradients)
1488  {
1489  if (integration_flag & EvaluationFlags::values ||
1490  add_into_values_array == true)
1491  eval.template gradients<0, false, true>(gradients_quad,
1492  values_dofs);
1493  else
1494  eval.template gradients<0, false, false>(gradients_quad,
1495  values_dofs);
1496  if (dim > 1)
1497  eval.template gradients<1, false, true>(gradients_quad +
1498  n_q_points,
1499  values_dofs);
1500  if (dim > 2)
1501  eval.template gradients<2, false, true>(gradients_quad +
1502  2 * n_q_points,
1503  values_dofs);
1504  }
1505  if (integration_flag & EvaluationFlags::hessians)
1506  {
1507  // diagonal
1508  // grad xx
1509  if (integration_flag & EvaluationFlags::values ||
1510  integration_flag & EvaluationFlags::gradients ||
1511  add_into_values_array == true)
1512  eval.template hessians<0, false, true>(hessians_quad,
1513  values_dofs);
1514  else
1515  eval.template hessians<0, false, false>(hessians_quad,
1516  values_dofs);
1517  // grad yy
1518  if (dim > 1)
1519  eval.template hessians<1, false, true>(hessians_quad + n_q_points,
1520  values_dofs);
1521  // grad zz
1522  if (dim > 2)
1523  eval.template hessians<2, false, true>(hessians_quad +
1524  2 * n_q_points,
1525  values_dofs);
1526  // off-diagonal
1527  if (dim == 2)
1528  {
1529  // grad xy
1530  eval.template gradients<0, false, false>(hessians_quad +
1531  2 * n_q_points,
1532  scratch_data);
1533  eval.template gradients<1, false, true>(scratch_data,
1534  values_dofs);
1535  }
1536  if (dim == 3)
1537  {
1538  // grad xy
1539  eval.template gradients<0, false, false>(hessians_quad +
1540  3 * n_q_points,
1541  scratch_data);
1542  eval.template gradients<1, false, true>(scratch_data,
1543  values_dofs);
1544  // grad xz
1545  eval.template gradients<0, false, false>(hessians_quad +
1546  4 * n_q_points,
1547  scratch_data);
1548  eval.template gradients<2, false, true>(scratch_data,
1549  values_dofs);
1550  // grad yz
1551  eval.template gradients<1, false, false>(hessians_quad +
1552  5 * n_q_points,
1553  scratch_data);
1554  eval.template gradients<2, false, true>(scratch_data,
1555  values_dofs);
1556  }
1557  hessians_quad += hdim * n_q_points;
1558  }
1559  gradients_quad += dim * n_q_points;
1560  values_quad += n_q_points;
1561  values_dofs += n_q_points;
1562  }
1563  }
1564 
1565 
1566 
1577  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1579  {
1580  static void
1581  evaluate(const unsigned int n_components,
1582  const EvaluationFlags::EvaluationFlags evaluation_flag,
1583  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1584  const Number * values_dofs,
1585  Number * values_quad,
1586  Number * gradients_quad,
1587  Number * hessians_quad,
1588  Number * scratch_data);
1589 
1590  static void
1591  integrate(const unsigned int n_components,
1592  const EvaluationFlags::EvaluationFlags evaluation_flag,
1593  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1594  Number * values_dofs,
1595  Number * values_quad,
1596  Number * gradients_quad,
1597  Number * hessians_quad,
1598  Number * scratch_data,
1599  const bool add_into_values_array);
1600  };
1601 
1602 
1603 
1604  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1605  inline void
1607  dim,
1608  fe_degree,
1609  n_q_points_1d,
1610  Number>::evaluate(const unsigned int n_components,
1611  const EvaluationFlags::EvaluationFlags evaluation_flag,
1612  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1613  const Number * values_dofs,
1614  Number * values_quad,
1615  Number *gradients_quad,
1616  Number *hessians_quad,
1617  Number *)
1618  {
1619  Assert(n_q_points_1d > fe_degree,
1620  ExcMessage("You lose information when going to a collocation space "
1621  "of lower degree, so the evaluation results would be "
1622  "wrong. Thus, this class does not permit the desired "
1623  "operation."));
1624  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
1625 
1626  for (unsigned int c = 0; c < n_components; ++c)
1627  {
1631  dim,
1632  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1633  n_q_points_1d,
1634  Number,
1635  Number>::do_forward(1,
1636  shape_info.data.front().shape_values_eo,
1637  values_dofs,
1638  values_quad);
1639 
1640  // apply derivatives in the collocation space
1641  if (evaluation_flag &
1644  1,
1645  evaluation_flag &
1647  shape_info,
1648  values_quad,
1649  nullptr,
1650  gradients_quad,
1651  hessians_quad,
1652  nullptr);
1653 
1654  values_dofs += shape_info.dofs_per_component_on_cell;
1655  values_quad += n_q_points;
1656  gradients_quad += dim * n_q_points;
1657  hessians_quad += (dim * (dim + 1)) / 2 * n_q_points;
1658  }
1659  }
1660 
1661 
1662 
1663  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1664  inline void
1666  dim,
1667  fe_degree,
1668  n_q_points_1d,
1669  Number>::integrate(const unsigned int n_components,
1670  const EvaluationFlags::EvaluationFlags integration_flag,
1671  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1672  Number * values_dofs,
1673  Number * values_quad,
1674  Number * gradients_quad,
1675  Number * hessians_quad,
1676  Number * scratch_data,
1677  const bool add_into_values_array)
1678  {
1679  Assert(n_q_points_1d > fe_degree,
1680  ExcMessage("You lose information when going to a collocation space "
1681  "of lower degree, so the evaluation results would be "
1682  "wrong. Thus, this class does not permit the desired "
1683  "operation."));
1685  shape_info.data.front().shape_gradients_collocation_eo.size(),
1686  (n_q_points_1d + 1) / 2 * n_q_points_1d);
1687  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
1688  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
1689 
1690  for (unsigned int c = 0; c < n_components; ++c)
1691  {
1692  // apply derivatives in collocation space
1693  if (integration_flag &
1696  integrate(1,
1697  integration_flag & (EvaluationFlags::gradients |
1699  shape_info,
1700  values_quad,
1701  nullptr,
1702  gradients_quad,
1703  hessians_quad,
1704  scratch_data,
1705  /*add_into_values_array=*/integration_flag &
1707 
1708  // transform back to the original space
1712  dim,
1713  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1714  n_q_points_1d,
1715  Number,
1716  Number>::do_backward(1,
1717  shape_info.data.front().shape_values_eo,
1718  add_into_values_array,
1719  values_quad,
1720  values_dofs);
1721  hessians_quad += hdim * n_q_points;
1722  gradients_quad += dim * n_q_points;
1723  values_quad += n_q_points;
1724  values_dofs += shape_info.dofs_per_component_on_cell;
1725  }
1726  }
1727 
1728 
1729 
1745  template <int dim, typename Number>
1747  {
1748  template <int fe_degree, int n_q_points_1d>
1749  static bool
1750  run(const unsigned int n_components,
1751  const EvaluationFlags::EvaluationFlags evaluation_flag,
1753  Number *values_dofs_actual,
1754  Number *values_quad,
1755  Number *gradients_quad,
1756  Number *hessians_quad,
1757  Number *scratch_data)
1758  {
1759  // We enable a transformation to collocation for derivatives if it gives
1760  // correct results (first condition), if it is the most efficient choice
1761  // in terms of operation counts (second condition) and if we were able to
1762  // initialize the fields in shape_info.templates.h from the polynomials
1763  // (third condition).
1764  static constexpr bool use_collocation =
1765  n_q_points_1d > fe_degree && n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
1766  n_q_points_1d < 200;
1767 
1768  if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
1769  shape_info.element_type ==
1771  {
1773  evaluate(n_components,
1774  evaluation_flag,
1775  shape_info,
1776  values_dofs_actual,
1777  values_quad,
1778  gradients_quad,
1779  hessians_quad,
1780  scratch_data);
1781  }
1782  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
1783  // shape_info.h for more details
1784  else if (fe_degree >= 0 && use_collocation &&
1785  shape_info.element_type <=
1787  {
1789  dim,
1790  fe_degree,
1791  n_q_points_1d,
1792  Number>::evaluate(n_components,
1793  evaluation_flag,
1794  shape_info,
1795  values_dofs_actual,
1796  values_quad,
1797  gradients_quad,
1798  hessians_quad,
1799  scratch_data);
1800  }
1801  else if (fe_degree >= 0 &&
1802  shape_info.element_type <=
1804  {
1807  dim,
1808  fe_degree,
1809  n_q_points_1d,
1810  Number>::evaluate(n_components,
1811  evaluation_flag,
1812  shape_info,
1813  values_dofs_actual,
1814  values_quad,
1815  gradients_quad,
1816  hessians_quad,
1817  scratch_data);
1818  }
1819  else if (shape_info.element_type ==
1821  {
1824  dim,
1825  fe_degree,
1826  n_q_points_1d,
1827  Number>::evaluate(n_components,
1828  evaluation_flag,
1829  shape_info,
1830  values_dofs_actual,
1831  values_quad,
1832  gradients_quad,
1833  hessians_quad,
1834  scratch_data);
1835  }
1836  else if (shape_info.element_type ==
1838  {
1841  dim,
1842  fe_degree,
1843  n_q_points_1d,
1844  Number>::evaluate(n_components,
1845  evaluation_flag,
1846  shape_info,
1847  values_dofs_actual,
1848  values_quad,
1849  gradients_quad,
1850  hessians_quad,
1851  scratch_data);
1852  }
1853  else if (shape_info.element_type ==
1855  {
1857  dim,
1858  fe_degree,
1859  n_q_points_1d,
1860  Number>::evaluate(n_components,
1861  evaluation_flag,
1862  shape_info,
1863  values_dofs_actual,
1864  values_quad,
1865  gradients_quad,
1866  hessians_quad,
1867  scratch_data);
1868  }
1869  else
1870  {
1873  dim,
1874  fe_degree,
1875  n_q_points_1d,
1876  Number>::evaluate(n_components,
1877  evaluation_flag,
1878  shape_info,
1879  values_dofs_actual,
1880  values_quad,
1881  gradients_quad,
1882  hessians_quad,
1883  scratch_data);
1884  }
1885 
1886  return false;
1887  }
1888  };
1889 
1890 
1891 
1907  template <int dim, typename Number>
1909  {
1910  template <int fe_degree, int n_q_points_1d>
1911  static bool
1912  run(const unsigned int n_components,
1913  const EvaluationFlags::EvaluationFlags integration_flag,
1915  Number * values_dofs_actual,
1916  Number * values_quad,
1917  Number * gradients_quad,
1918  Number * hessians_quad,
1919  Number * scratch_data,
1920  const bool sum_into_values_array)
1921  {
1922  // We enable a transformation to collocation for derivatives if it gives
1923  // correct results (first condition), if it is the most efficient choice
1924  // in terms of operation counts (second condition) and if we were able to
1925  // initialize the fields in shape_info.templates.h from the polynomials
1926  // (third condition).
1927  constexpr bool use_collocation = n_q_points_1d > fe_degree &&
1928  n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
1929  n_q_points_1d < 200;
1930 
1931  if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
1932  shape_info.element_type ==
1934  {
1936  integrate(n_components,
1937  integration_flag,
1938  shape_info,
1939  values_dofs_actual,
1940  values_quad,
1941  gradients_quad,
1942  hessians_quad,
1943  scratch_data,
1944  sum_into_values_array);
1945  }
1946  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
1947  // shape_info.h for more details
1948  else if (fe_degree >= 0 && use_collocation &&
1949  shape_info.element_type <=
1951  {
1953  dim,
1954  fe_degree,
1955  n_q_points_1d,
1956  Number>::integrate(n_components,
1957  integration_flag,
1958  shape_info,
1959  values_dofs_actual,
1960  values_quad,
1961  gradients_quad,
1962  hessians_quad,
1963  scratch_data,
1964  sum_into_values_array);
1965  }
1966  else if (fe_degree >= 0 &&
1967  shape_info.element_type <=
1969  {
1972  dim,
1973  fe_degree,
1974  n_q_points_1d,
1975  Number>::integrate(n_components,
1976  integration_flag,
1977  shape_info,
1978  values_dofs_actual,
1979  values_quad,
1980  gradients_quad,
1981  hessians_quad,
1982  scratch_data,
1983  sum_into_values_array);
1984  }
1985  else if (shape_info.element_type ==
1987  {
1990  dim,
1991  fe_degree,
1992  n_q_points_1d,
1993  Number>::integrate(n_components,
1994  integration_flag,
1995  shape_info,
1996  values_dofs_actual,
1997  values_quad,
1998  gradients_quad,
1999  hessians_quad,
2000  scratch_data,
2001  sum_into_values_array);
2002  }
2003  else if (shape_info.element_type ==
2005  {
2008  dim,
2009  fe_degree,
2010  n_q_points_1d,
2011  Number>::integrate(n_components,
2012  integration_flag,
2013  shape_info,
2014  values_dofs_actual,
2015  values_quad,
2016  gradients_quad,
2017  hessians_quad,
2018  scratch_data,
2019  sum_into_values_array);
2020  }
2021  else if (shape_info.element_type ==
2023  {
2025  dim,
2026  fe_degree,
2027  n_q_points_1d,
2028  Number>::integrate(n_components,
2029  integration_flag,
2030  shape_info,
2031  values_dofs_actual,
2032  values_quad,
2033  gradients_quad,
2034  hessians_quad,
2035  scratch_data,
2036  sum_into_values_array);
2037  }
2038  else
2039  {
2042  dim,
2043  fe_degree,
2044  n_q_points_1d,
2045  Number>::integrate(n_components,
2046  integration_flag,
2047  shape_info,
2048  values_dofs_actual,
2049  values_quad,
2050  gradients_quad,
2051  hessians_quad,
2052  scratch_data,
2053  sum_into_values_array);
2054  }
2055 
2056  return false;
2057  }
2058  };
2059 
2060 
2061 
2062  template <bool symmetric_evaluate,
2063  int dim,
2064  int fe_degree,
2065  int n_q_points_1d,
2066  typename Number>
2068  {
2069  // We enable a transformation to collocation for derivatives if it gives
2070  // correct results (first two conditions), if it is the most efficient
2071  // choice in terms of operation counts (third condition) and if we were
2072  // able to initialize the fields in shape_info.templates.h from the
2073  // polynomials (fourth condition).
2074  static constexpr bool use_collocation =
2075  symmetric_evaluate &&
2076  n_q_points_1d > fe_degree &&n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
2077  n_q_points_1d < 200;
2078 
2079  static void
2080  evaluate_in_face(const unsigned int n_components,
2082  Number * values_dofs,
2083  Number * values_quad,
2084  Number * gradients_quad,
2085  Number * hessians_quad,
2086  Number * scratch_data,
2087  const bool evaluate_val,
2088  const bool evaluate_grad,
2089  const bool evaluate_hessian,
2090  const unsigned int subface_index)
2091  {
2092  const AlignedVector<Number> &val0 =
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  const AlignedVector<Number> &val1 =
2099  symmetric_evaluate ?
2100  data.data.front().shape_values_eo :
2101  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
2102  data.data.front().shape_values :
2103  data.data.front().values_within_subface[subface_index / 2]);
2104 
2105  const AlignedVector<Number> &grad0 =
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  const AlignedVector<Number> &grad1 =
2112  symmetric_evaluate ?
2113  data.data.front().shape_gradients_eo :
2114  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
2115  data.data.front().shape_gradients :
2116  data.data.front().gradients_within_subface[subface_index / 2]);
2117 
2118  const AlignedVector<Number> &hessian0 =
2119  symmetric_evaluate ?
2120  data.data.front().shape_hessians_eo :
2121  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
2122  data.data.front().shape_hessians :
2123  data.data.front().hessians_within_subface[subface_index % 2]);
2124  const AlignedVector<Number> &hessian1 =
2125  symmetric_evaluate ?
2126  data.data.front().shape_hessians_eo :
2127  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
2128  data.data.front().shape_hessians :
2129  data.data.front().hessians_within_subface[subface_index / 2]);
2130 
2131  using Eval =
2132  internal::EvaluatorTensorProduct<symmetric_evaluate ?
2135  dim - 1,
2136  fe_degree + 1,
2137  n_q_points_1d,
2138  Number>;
2139  Eval eval0(val0,
2140  grad0,
2141  hessian0,
2142  data.data.front().fe_degree + 1,
2143  data.data.front().n_q_points_1d);
2144  Eval eval1(val1,
2145  grad1,
2146  hessian1,
2147  data.data.front().fe_degree + 1,
2148  data.data.front().n_q_points_1d);
2149 
2150  const unsigned int size_deg =
2151  fe_degree > -1 ?
2152  Utilities::pow(fe_degree + 1, dim - 1) :
2153  (dim > 1 ?
2154  Utilities::fixed_power<dim - 1>(data.data.front().fe_degree + 1) :
2155  1);
2156 
2157  const unsigned int n_q_points = fe_degree > -1 ?
2158  Utilities::pow(n_q_points_1d, dim - 1) :
2159  data.n_q_points_face;
2160 
2161  // keep a copy of the original pointer for the case of the Hessians
2162  Number *values_dofs_ptr = values_dofs;
2163 
2164  if (evaluate_val && !evaluate_grad)
2165  for (unsigned int c = 0; c < n_components; ++c)
2166  {
2167  switch (dim)
2168  {
2169  case 3:
2170  eval0.template values<0, true, false>(values_dofs,
2171  values_quad);
2172  eval1.template values<1, true, false>(values_quad,
2173  values_quad);
2174  break;
2175  case 2:
2176  eval0.template values<0, true, false>(values_dofs,
2177  values_quad);
2178  break;
2179  case 1:
2180  values_quad[0] = values_dofs[0];
2181  break;
2182  default:
2183  Assert(false, ExcNotImplemented());
2184  }
2185  values_dofs += 3 * size_deg;
2186  values_quad += n_q_points;
2187  }
2188  else if (evaluate_grad)
2189  for (unsigned int c = 0; c < n_components; ++c)
2190  {
2191  switch (dim)
2192  {
2193  case 3:
2194  if (use_collocation)
2195  {
2196  eval0.template values<0, true, false>(values_dofs,
2197  values_quad);
2198  eval0.template values<1, true, false>(values_quad,
2199  values_quad);
2202  dim - 1,
2203  n_q_points_1d,
2204  n_q_points_1d,
2205  Number>
2206  eval_grad(
2208  data.data.front().shape_gradients_collocation_eo,
2210  eval_grad.template gradients<0, true, false>(
2211  values_quad, gradients_quad);
2212  eval_grad.template gradients<1, true, false>(
2213  values_quad, gradients_quad + n_q_points);
2214  }
2215  else
2216  {
2217  // grad x
2218  eval0.template gradients<0, true, false>(values_dofs,
2219  scratch_data);
2220  eval1.template values<1, true, false>(scratch_data,
2221  gradients_quad);
2222 
2223  // grad y
2224  eval0.template values<0, true, false>(values_dofs,
2225  scratch_data);
2226  eval1.template gradients<1, true, false>(scratch_data,
2227  gradients_quad +
2228  n_q_points);
2229 
2230  if (evaluate_val == true)
2231  eval1.template values<1, true, false>(scratch_data,
2232  values_quad);
2233  }
2234  // grad z
2235  eval0.template values<0, true, false>(values_dofs + size_deg,
2236  scratch_data);
2237  eval1.template values<1, true, false>(
2238  scratch_data, gradients_quad + (dim - 1) * n_q_points);
2239 
2240  break;
2241  case 2:
2242  eval0.template values<0, true, false>(values_dofs + size_deg,
2243  gradients_quad +
2244  (dim - 1) *
2245  n_q_points);
2246  eval0.template gradients<0, true, false>(values_dofs,
2247  gradients_quad);
2248  if (evaluate_val == true)
2249  eval0.template values<0, true, false>(values_dofs,
2250  values_quad);
2251  break;
2252  case 1:
2253  values_quad[0] = values_dofs[0];
2254  gradients_quad[0] = values_dofs[1];
2255  break;
2256  default:
2257  AssertThrow(false, ExcNotImplemented());
2258  }
2259  values_dofs += 3 * size_deg;
2260  values_quad += n_q_points;
2261  gradients_quad += dim * n_q_points;
2262  }
2263 
2264  if (evaluate_hessian)
2265  {
2266  values_dofs = values_dofs_ptr;
2267  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
2268  for (unsigned int c = 0; c < n_components; ++c)
2269  {
2270  switch (dim)
2271  {
2272  case 3:
2273  // grad xx
2274  eval0.template hessians<0, true, false>(values_dofs,
2275  scratch_data);
2276  eval1.template values<1, true, false>(scratch_data,
2277  hessians_quad);
2278 
2279  // grad yy
2280  eval0.template values<0, true, false>(values_dofs,
2281  scratch_data);
2282  eval1.template hessians<1, true, false>(scratch_data,
2283  hessians_quad +
2284  n_q_points);
2285 
2286  // grad zz
2287  eval0.template values<0, true, false>(values_dofs +
2288  2 * size_deg,
2289  scratch_data);
2290  eval1.template values<1, true, false>(scratch_data,
2291  hessians_quad +
2292  2 * n_q_points);
2293 
2294  // grad xy
2295  eval0.template gradients<0, true, false>(values_dofs,
2296  scratch_data);
2297  eval1.template gradients<1, true, false>(scratch_data,
2298  hessians_quad +
2299  3 * n_q_points);
2300 
2301  // grad xz
2302  eval0.template gradients<0, true, false>(values_dofs +
2303  size_deg,
2304  scratch_data);
2305  eval1.template values<1, true, false>(scratch_data,
2306  hessians_quad +
2307  4 * n_q_points);
2308 
2309  // grad yz
2310  eval0.template values<0, true, false>(values_dofs +
2311  size_deg,
2312  scratch_data);
2313  eval1.template gradients<1, true, false>(scratch_data,
2314  hessians_quad +
2315  5 * n_q_points);
2316 
2317  break;
2318  case 2:
2319  // grad xx
2320  eval0.template hessians<0, true, false>(values_dofs,
2321  hessians_quad);
2322  // grad yy
2323  eval0.template values<0, true, false>(
2324  values_dofs + 2 * size_deg, hessians_quad + n_q_points);
2325  // grad xy
2326  eval0.template gradients<0, true, false>(
2327  values_dofs + size_deg, hessians_quad + 2 * n_q_points);
2328  break;
2329  case 1:
2330  hessians_quad[0] = values_dofs[2];
2331  break;
2332  default:
2333  AssertThrow(false, ExcNotImplemented());
2334  }
2335  values_dofs += 3 * size_deg;
2336  hessians_quad += hdim * n_q_points;
2337  }
2338  }
2339  }
2340 
2341  static void
2342  integrate_in_face(const unsigned int n_components,
2344  Number * values_dofs,
2345  Number * values_quad,
2346  Number * gradients_quad,
2347  Number * hessians_quad,
2348  Number * scratch_data,
2349  const bool integrate_val,
2350  const bool integrate_grad,
2351  const bool integrate_hessian,
2352  const unsigned int subface_index)
2353  {
2354  const AlignedVector<Number> &val0 =
2355  symmetric_evaluate ?
2356  data.data.front().shape_values_eo :
2357  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
2358  data.data.front().shape_values :
2359  data.data.front().values_within_subface[subface_index % 2]);
2360  const AlignedVector<Number> &val1 =
2361  symmetric_evaluate ?
2362  data.data.front().shape_values_eo :
2363  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
2364  data.data.front().shape_values :
2365  data.data.front().values_within_subface[subface_index / 2]);
2366 
2367  const AlignedVector<Number> &grad0 =
2368  symmetric_evaluate ?
2369  data.data.front().shape_gradients_eo :
2370  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
2371  data.data.front().shape_gradients :
2372  data.data.front().gradients_within_subface[subface_index % 2]);
2373  const AlignedVector<Number> &grad1 =
2374  symmetric_evaluate ?
2375  data.data.front().shape_gradients_eo :
2376  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
2377  data.data.front().shape_gradients :
2378  data.data.front().gradients_within_subface[subface_index / 2]);
2379 
2380  const AlignedVector<Number> &hessian0 =
2381  symmetric_evaluate ?
2382  data.data.front().shape_hessians_eo :
2383  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
2384  data.data.front().shape_hessians :
2385  data.data.front().hessians_within_subface[subface_index % 2]);
2386  const AlignedVector<Number> &hessian1 =
2387  symmetric_evaluate ?
2388  data.data.front().shape_hessians_eo :
2389  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
2390  data.data.front().shape_hessians :
2391  data.data.front().hessians_within_subface[subface_index / 2]);
2392 
2393  using Eval =
2394  internal::EvaluatorTensorProduct<symmetric_evaluate ?
2397  dim - 1,
2398  fe_degree + 1,
2399  n_q_points_1d,
2400  Number>;
2401  Eval eval0(val0,
2402  grad0,
2403  hessian0,
2404  data.data.front().fe_degree + 1,
2405  data.data.front().n_q_points_1d);
2406  Eval eval1(val1,
2407  grad1,
2408  hessian1,
2409  data.data.front().fe_degree + 1,
2410  data.data.front().n_q_points_1d);
2411 
2412  const unsigned int size_deg =
2413  fe_degree > -1 ?
2414  Utilities::pow(fe_degree + 1, dim - 1) :
2415  (dim > 1 ?
2416  Utilities::fixed_power<dim - 1>(data.data.front().fe_degree + 1) :
2417  1);
2418 
2419  const unsigned int n_q_points = fe_degree > -1 ?
2420  Utilities::pow(n_q_points_1d, dim - 1) :
2421  data.n_q_points_face;
2422 
2423  // keep a copy of the original pointer for the case of the Hessians
2424  Number *values_dofs_ptr = values_dofs;
2425 
2426  if (integrate_val && !integrate_grad)
2427  for (unsigned int c = 0; c < n_components; ++c)
2428  {
2429  switch (dim)
2430  {
2431  case 3:
2432  eval1.template values<1, false, false>(values_quad,
2433  values_quad);
2434  eval0.template values<0, false, false>(values_quad,
2435  values_dofs);
2436  break;
2437  case 2:
2438  eval0.template values<0, false, false>(values_quad,
2439  values_dofs);
2440  break;
2441  case 1:
2442  values_dofs[0] = values_quad[0];
2443  break;
2444  default:
2445  Assert(false, ExcNotImplemented());
2446  }
2447  values_dofs += 3 * size_deg;
2448  values_quad += n_q_points;
2449  }
2450  else if (integrate_grad)
2451  for (unsigned int c = 0; c < n_components; ++c)
2452  {
2453  switch (dim)
2454  {
2455  case 3:
2456  // grad z
2457  eval1.template values<1, false, false>(gradients_quad +
2458  2 * n_q_points,
2459  gradients_quad +
2460  2 * n_q_points);
2461  eval0.template values<0, false, false>(
2462  gradients_quad + 2 * n_q_points, values_dofs + size_deg);
2463  if (use_collocation)
2464  {
2467  dim - 1,
2468  n_q_points_1d,
2469  n_q_points_1d,
2470  Number>
2471  eval_grad(
2473  data.data.front().shape_gradients_collocation_eo,
2475  if (integrate_val)
2476  eval_grad.template gradients<1, false, true>(
2477  gradients_quad + n_q_points, values_quad);
2478  else
2479  eval_grad.template gradients<1, false, false>(
2480  gradients_quad + n_q_points, values_quad);
2481  eval_grad.template gradients<0, false, true>(
2482  gradients_quad, values_quad);
2483  eval0.template values<1, false, false>(values_quad,
2484  values_quad);
2485  eval0.template values<0, false, false>(values_quad,
2486  values_dofs);
2487  }
2488  else
2489  {
2490  if (integrate_val)
2491  {
2492  eval1.template values<1, false, false>(values_quad,
2493  scratch_data);
2494  eval1.template gradients<1, false, true>(
2495  gradients_quad + n_q_points, scratch_data);
2496  }
2497  else
2498  eval1.template gradients<1, false, false>(
2499  gradients_quad + n_q_points, scratch_data);
2500 
2501  // grad y
2502  eval0.template values<0, false, false>(scratch_data,
2503  values_dofs);
2504 
2505  // grad x
2506  eval1.template values<1, false, false>(gradients_quad,
2507  scratch_data);
2508  eval0.template gradients<0, false, true>(scratch_data,
2509  values_dofs);
2510  }
2511  break;
2512  case 2:
2513  eval0.template values<0, false, false>(
2514  gradients_quad + n_q_points, values_dofs + size_deg);
2515  eval0.template gradients<0, false, false>(gradients_quad,
2516  values_dofs);
2517  if (integrate_val == true)
2518  eval0.template values<0, false, true>(values_quad,
2519  values_dofs);
2520  break;
2521  case 1:
2522  values_dofs[0] = values_quad[0];
2523  values_dofs[1] = gradients_quad[0];
2524  break;
2525  default:
2526  AssertThrow(false, ExcNotImplemented());
2527  }
2528  values_dofs += 3 * size_deg;
2529  values_quad += n_q_points;
2530  gradients_quad += dim * n_q_points;
2531  }
2532 
2533  if (integrate_hessian)
2534  {
2535  values_dofs = values_dofs_ptr;
2536  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
2537  for (unsigned int c = 0; c < n_components; ++c)
2538  {
2539  switch (dim)
2540  {
2541  case 3:
2542  // grad xx
2543  eval1.template values<1, false, false>(hessians_quad,
2544  scratch_data);
2545  if (integrate_val || integrate_grad)
2546  eval0.template hessians<0, false, true>(scratch_data,
2547  values_dofs);
2548  else
2549  eval0.template hessians<0, false, false>(scratch_data,
2550  values_dofs);
2551 
2552  // grad yy
2553  eval1.template hessians<1, false, false>(hessians_quad +
2554  n_q_points,
2555  scratch_data);
2556  eval0.template values<0, false, true>(scratch_data,
2557  values_dofs);
2558 
2559  // grad zz
2560  eval1.template values<1, false, false>(hessians_quad +
2561  2 * n_q_points,
2562  scratch_data);
2563  eval0.template values<0, false, false>(scratch_data,
2564  values_dofs +
2565  2 * size_deg);
2566 
2567  // grad xy
2568  eval1.template gradients<1, false, false>(hessians_quad +
2569  3 * n_q_points,
2570  scratch_data);
2571  eval0.template gradients<0, false, true>(scratch_data,
2572  values_dofs);
2573 
2574  // grad xz
2575  eval1.template values<1, false, false>(hessians_quad +
2576  4 * n_q_points,
2577  scratch_data);
2578  if (integrate_grad)
2579  eval0.template gradients<0, false, true>(scratch_data,
2580  values_dofs +
2581  size_deg);
2582  else
2583  eval0.template gradients<0, false, false>(scratch_data,
2584  values_dofs +
2585  size_deg);
2586 
2587  // grad yz
2588  eval1.template gradients<1, false, false>(hessians_quad +
2589  5 * n_q_points,
2590  scratch_data);
2591  eval0.template values<0, false, true>(scratch_data,
2592  values_dofs +
2593  size_deg);
2594 
2595  break;
2596  case 2:
2597  // grad xx
2598  if (integrate_val || integrate_grad)
2599  eval0.template hessians<0, false, true>(hessians_quad,
2600  values_dofs);
2601  else
2602  eval0.template hessians<0, false, false>(hessians_quad,
2603  values_dofs);
2604 
2605  // grad yy
2606  eval0.template values<0, false, false>(
2607  hessians_quad + n_q_points, values_dofs + 2 * size_deg);
2608  // grad xy
2609  if (integrate_grad)
2610  eval0.template gradients<0, false, true>(
2611  hessians_quad + 2 * n_q_points, values_dofs + size_deg);
2612  else
2613  eval0.template gradients<0, false, false>(
2614  hessians_quad + 2 * n_q_points, values_dofs + size_deg);
2615  break;
2616  case 1:
2617  values_dofs[2] = hessians_quad[0];
2618  if (!integrate_val)
2619  values_dofs[0] = 0;
2620  if (!integrate_grad)
2621  values_dofs[1] = 0;
2622  break;
2623  default:
2624  AssertThrow(false, ExcNotImplemented());
2625  }
2626  values_dofs += 3 * size_deg;
2627  hessians_quad += hdim * n_q_points;
2628  }
2629  }
2630  }
2631  };
2632 
2633 
2634 
2635  template <int dim, int fe_degree, typename Number, bool lex_faces = false>
2637  {
2638  template <bool do_evaluate, bool add_into_output>
2639  static void
2640  interpolate(const unsigned int n_components,
2642  const Number * input,
2643  Number * output,
2644  const bool do_gradients,
2645  const bool do_hessians,
2646  const unsigned int face_no)
2647  {
2648  Assert(static_cast<unsigned int>(fe_degree) ==
2649  data.data.front().fe_degree ||
2650  fe_degree == -1,
2651  ExcInternalError());
2652 
2653  interpolate_generic<do_evaluate, add_into_output>(
2654  n_components,
2655  input,
2656  output,
2657  do_gradients,
2658  do_hessians,
2659  face_no,
2660  data.data.front().fe_degree + 1,
2661  data.data.front().shape_data_on_face,
2663  3 * data.dofs_per_component_on_face);
2664  }
2665 
2669  template <bool do_evaluate, bool add_into_output>
2670  static void
2671  interpolate_quadrature(const unsigned int n_components,
2673  const Number * input,
2674  Number * output,
2675  const bool do_gradients,
2676  const bool do_hessians,
2677  const unsigned int face_no)
2678  {
2679  Assert(static_cast<unsigned int>(fe_degree + 1) ==
2680  data.data.front().quadrature.size() ||
2681  fe_degree == -1,
2682  ExcInternalError());
2683 
2684  interpolate_generic<do_evaluate, add_into_output>(
2685  n_components,
2686  input,
2687  output,
2688  do_gradients,
2689  do_hessians,
2690  face_no,
2691  data.data.front().quadrature.size(),
2692  data.data.front().quadrature_data_on_face,
2693  data.n_q_points,
2694  data.n_q_points_face);
2695  }
2696 
2697  private:
2698  template <bool do_evaluate, bool add_into_output, int face_direction = 0>
2699  static void
2700  interpolate_generic(const unsigned int n_components,
2701  const Number * input,
2702  Number * output,
2703  const bool do_gradients,
2704  const bool do_hessians,
2705  const unsigned int face_no,
2706  const unsigned int n_points_1d,
2707  const std::array<AlignedVector<Number>, 2> &shape_data,
2708  const unsigned int dofs_per_component_on_cell,
2709  const unsigned int dofs_per_component_on_face)
2710  {
2711  if (face_direction == face_no / 2)
2712  {
2714  dim,
2715  fe_degree + 1,
2716  0,
2717  Number>
2718  evalf(shape_data[face_no % 2],
2721  n_points_1d,
2722  0);
2723 
2724  const unsigned int in_stride = do_evaluate ?
2725  dofs_per_component_on_cell :
2726  dofs_per_component_on_face;
2727  const unsigned int out_stride = do_evaluate ?
2728  dofs_per_component_on_face :
2729  dofs_per_component_on_cell;
2730 
2731  for (unsigned int c = 0; c < n_components; ++c)
2732  {
2733  if (do_hessians)
2734  evalf.template apply_face<face_direction,
2735  do_evaluate,
2736  add_into_output,
2737  2,
2738  lex_faces>(input, output);
2739  else if (!do_hessians && do_gradients)
2740  evalf.template apply_face<face_direction,
2741  do_evaluate,
2742  add_into_output,
2743  1,
2744  lex_faces>(input, output);
2745  else
2746  evalf.template apply_face<face_direction,
2747  do_evaluate,
2748  add_into_output,
2749  0,
2750  lex_faces>(input, output);
2751  input += in_stride;
2752  output += out_stride;
2753  }
2754  }
2755  else if (face_direction < dim)
2756  {
2757  interpolate_generic<do_evaluate,
2758  add_into_output,
2759  std::min(face_direction + 1, dim - 1)>(
2760  n_components,
2761  input,
2762  output,
2763  do_gradients,
2764  do_hessians,
2765  face_no,
2766  n_points_1d,
2767  shape_data,
2768  dofs_per_component_on_cell,
2769  dofs_per_component_on_face);
2770  }
2771  }
2772  };
2773 
2774 
2775 
2776  // internal helper function for reading data; base version of different types
2777  template <typename VectorizedArrayType, typename Number2>
2778  void
2779  do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
2780  {
2781  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2782  dst[v] = src_ptr[v];
2783  }
2784 
2785 
2786 
2787  // internal helper function for reading data; specialized version where we
2788  // can use a dedicated load function
2789  template <typename Number, unsigned int width>
2790  void
2792  {
2793  dst.load(src_ptr);
2794  }
2795 
2796 
2797 
2798  // internal helper function for reading data; base version of different types
2799  template <typename VectorizedArrayType, typename Number2>
2800  void
2801  do_vectorized_gather(const Number2 * src_ptr,
2802  const unsigned int * indices,
2803  VectorizedArrayType &dst)
2804  {
2805  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2806  dst[v] = src_ptr[indices[v]];
2807  }
2808 
2809 
2810 
2811  // internal helper function for reading data; specialized version where we
2812  // can use a dedicated gather function
2813  template <typename Number, unsigned int width>
2814  void
2815  do_vectorized_gather(const Number * src_ptr,
2816  const unsigned int * indices,
2818  {
2819  dst.gather(src_ptr, indices);
2820  }
2821 
2822 
2823 
2824  // internal helper function for reading data; base version of different types
2825  template <typename VectorizedArrayType, typename Number2>
2826  void
2827  do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
2828  {
2829  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2830  dst_ptr[v] += src[v];
2831  }
2832 
2833 
2834 
2835  // internal helper function for reading data; specialized version where we
2836  // can use a dedicated load function
2837  template <typename Number, unsigned int width>
2838  void
2840  {
2842  tmp.load(dst_ptr);
2843  (tmp + src).store(dst_ptr);
2844  }
2845 
2846 
2847 
2848  // internal helper function for reading data; base version of different types
2849  template <typename VectorizedArrayType, typename Number2>
2850  void
2851  do_vectorized_scatter_add(const VectorizedArrayType src,
2852  const unsigned int * indices,
2853  Number2 * dst_ptr)
2854  {
2855  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2856  dst_ptr[indices[v]] += src[v];
2857  }
2858 
2859 
2860 
2861  // internal helper function for reading data; specialized version where we
2862  // can use a dedicated gather function
2863  template <typename Number, unsigned int width>
2864  void
2866  const unsigned int * indices,
2867  Number * dst_ptr)
2868  {
2869 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
2870  for (unsigned int v = 0; v < width; ++v)
2871  dst_ptr[indices[v]] += src[v];
2872 #else
2874  tmp.gather(dst_ptr, indices);
2875  (tmp + src).scatter(indices, dst_ptr);
2876 #endif
2877  }
2878 
2879 
2880 
2881  template <typename Number>
2882  void
2883  adjust_for_face_orientation(const unsigned int dim,
2884  const unsigned int n_components,
2885  const unsigned int face_orientation,
2886  const Table<2, unsigned int> &orientation_map,
2887  const bool integrate,
2888  const bool values,
2889  const bool gradients,
2890  const bool hessians,
2891  const unsigned int n_q_points,
2892  Number * tmp_values,
2893  Number * values_quad,
2894  Number * gradients_quad,
2895  Number * hessians_quad)
2896  {
2897  Assert(face_orientation, ExcInternalError());
2898  const unsigned int *orientation = &orientation_map[face_orientation][0];
2899  for (unsigned int c = 0; c < n_components; ++c)
2900  {
2901  if (values == true)
2902  {
2903  if (integrate)
2904  for (unsigned int q = 0; q < n_q_points; ++q)
2905  tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
2906  else
2907  for (unsigned int q = 0; q < n_q_points; ++q)
2908  tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
2909  for (unsigned int q = 0; q < n_q_points; ++q)
2910  values_quad[c * n_q_points + q] = tmp_values[q];
2911  }
2912  if (gradients == true)
2913  for (unsigned int d = 0; d < dim; ++d)
2914  {
2915  if (integrate)
2916  for (unsigned int q = 0; q < n_q_points; ++q)
2917  tmp_values[q] =
2918  gradients_quad[(c * dim + d) * n_q_points + orientation[q]];
2919  else
2920  for (unsigned int q = 0; q < n_q_points; ++q)
2921  tmp_values[orientation[q]] =
2922  gradients_quad[(c * dim + d) * n_q_points + q];
2923  for (unsigned int q = 0; q < n_q_points; ++q)
2924  gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q];
2925  }
2926  if (hessians == true)
2927  {
2928  const unsigned int hdim = (dim * (dim + 1)) / 2;
2929  for (unsigned int d = 0; d < hdim; ++d)
2930  {
2931  if (integrate)
2932  for (unsigned int q = 0; q < n_q_points; ++q)
2933  tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
2934  orientation[q]];
2935  else
2936  for (unsigned int q = 0; q < n_q_points; ++q)
2937  tmp_values[orientation[q]] =
2938  hessians_quad[(c * hdim + d) * n_q_points + q];
2939  for (unsigned int q = 0; q < n_q_points; ++q)
2940  hessians_quad[(c * hdim + d) * n_q_points + q] =
2941  tmp_values[q];
2942  }
2943  }
2944  }
2945  }
2946 
2947 
2948 
2949  template <int dim, typename VectorizedArrayType>
2951  {
2952  template <int fe_degree, int n_q_points_1d>
2953  static bool
2954  run(const unsigned int n_components,
2956  const VectorizedArrayType * values_array,
2957  VectorizedArrayType * values_quad,
2958  VectorizedArrayType * gradients_quad,
2959  VectorizedArrayType * hessians_quad,
2960  VectorizedArrayType * scratch_data,
2961  const bool evaluate_values,
2962  const bool evaluate_gradients,
2963  const bool evaluate_hessians,
2964  const unsigned int face_no,
2965  const unsigned int subface_index,
2966  const unsigned int face_orientation,
2967  const Table<2, unsigned int> &orientation_map)
2968  {
2969  if (data.element_type == MatrixFreeFunctions::tensor_none)
2970  {
2971  const unsigned int n_dofs = data.dofs_per_component_on_cell;
2972  const unsigned int n_q_points = data.n_q_points_faces[face_no];
2973  const auto shape_info = data.data.front();
2974 
2976  1,
2977  0,
2978  0,
2979  VectorizedArrayType,
2980  VectorizedArrayType>;
2981 
2982  if (evaluate_values)
2983  {
2984  const auto shape_values =
2985  &shape_info.shape_values_face(face_no, face_orientation, 0);
2986 
2987  auto values_quad_ptr = values_quad;
2988  auto values_dofs_actual_ptr = values_array;
2989 
2990  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
2991  for (unsigned int c = 0; c < n_components; ++c)
2992  {
2993  eval.template values<0, true, false>(values_dofs_actual_ptr,
2994  values_quad_ptr);
2995 
2996  values_quad_ptr += n_q_points;
2997  values_dofs_actual_ptr += n_dofs;
2998  }
2999  }
3000 
3001  if (evaluate_gradients)
3002  {
3003  auto gradients_quad_ptr = gradients_quad;
3004  auto values_dofs_actual_ptr = values_array;
3005 
3006  std::array<const VectorizedArrayType *, dim> shape_gradients;
3007  for (unsigned int d = 0; d < dim; ++d)
3008  shape_gradients[d] = &shape_info.shape_gradients_face(
3009  face_no, face_orientation, d, 0);
3010 
3011  for (unsigned int c = 0; c < n_components; ++c)
3012  {
3013  for (unsigned int d = 0; d < dim; ++d)
3014  {
3015  Eval eval(nullptr,
3016  shape_gradients[d],
3017  nullptr,
3018  n_dofs,
3019  n_q_points);
3020 
3021  eval.template gradients<0, true, false>(
3022  values_dofs_actual_ptr, gradients_quad_ptr);
3023 
3024  gradients_quad_ptr += n_q_points;
3025  }
3026  values_dofs_actual_ptr += n_dofs;
3027  }
3028  }
3029 
3030  if (evaluate_hessians)
3031  {
3032  AssertThrow(false, ExcNotImplemented());
3033  }
3034 
3035  return true;
3036  }
3037 
3038  constexpr unsigned int static_dofs_per_face =
3039  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
3041  const unsigned int dofs_per_face =
3042  fe_degree > -1 ?
3043  static_dofs_per_face :
3044  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
3045 
3046  VectorizedArrayType *temp1 = scratch_data;
3047 
3049  template interpolate<true, false>(n_components,
3050  data,
3051  values_array,
3052  temp1,
3053  evaluate_gradients,
3054  evaluate_hessians,
3055  face_no);
3056 
3057  const unsigned int n_q_points_1d_actual =
3058  fe_degree > -1 ? n_q_points_1d : 0;
3059  if (fe_degree > -1 &&
3060  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
3061  data.element_type <= MatrixFreeFunctions::tensor_symmetric)
3063  true,
3064  dim,
3065  fe_degree,
3066  n_q_points_1d_actual,
3067  VectorizedArrayType>::evaluate_in_face(n_components,
3068  data,
3069  temp1,
3070  values_quad,
3071  gradients_quad,
3072  hessians_quad,
3073  scratch_data + 3 *
3074  n_components *
3075  dofs_per_face,
3076  evaluate_values,
3077  evaluate_gradients,
3078  evaluate_hessians,
3079  subface_index);
3080  else
3082  false,
3083  dim,
3084  fe_degree,
3085  n_q_points_1d_actual,
3086  VectorizedArrayType>::evaluate_in_face(n_components,
3087  data,
3088  temp1,
3089  values_quad,
3090  gradients_quad,
3091  hessians_quad,
3092  scratch_data + 3 *
3093  n_components *
3094  dofs_per_face,
3095  evaluate_values,
3096  evaluate_gradients,
3097  evaluate_hessians,
3098  subface_index);
3099 
3100  if (face_orientation)
3102  n_components,
3103  face_orientation,
3104  orientation_map,
3105  false,
3106  evaluate_values,
3107  evaluate_gradients,
3108  evaluate_hessians,
3109  data.n_q_points_face,
3110  scratch_data,
3111  values_quad,
3112  gradients_quad,
3113  hessians_quad);
3114 
3115  return false;
3116  }
3117  };
3118 
3119 
3120 
3121  template <int dim, typename VectorizedArrayType>
3123  {
3124  template <int fe_degree, int n_q_points_1d>
3125  static bool
3126  run(const unsigned int n_components,
3128  VectorizedArrayType * values_array,
3129  VectorizedArrayType * values_quad,
3130  VectorizedArrayType * gradients_quad,
3131  VectorizedArrayType * hessians_quad,
3132  VectorizedArrayType * scratch_data,
3133  const bool integrate_values,
3134  const bool integrate_gradients,
3135  const bool integrate_hessians,
3136  const unsigned int face_no,
3137  const unsigned int subface_index,
3138  const unsigned int face_orientation,
3139  const Table<2, unsigned int> &orientation_map)
3140  {
3141  if (data.element_type == MatrixFreeFunctions::tensor_none)
3142  {
3143  const unsigned int n_dofs = data.dofs_per_component_on_cell;
3144  const unsigned int n_q_points = data.n_q_points_faces[face_no];
3145  const auto shape_info = data.data.front();
3146 
3148  1,
3149  0,
3150  0,
3151  VectorizedArrayType,
3152  VectorizedArrayType>;
3153 
3154  if (integrate_values)
3155  {
3156  const auto shape_values =
3157  &shape_info.shape_values_face(face_no, face_orientation, 0);
3158 
3159  auto values_quad_ptr = values_quad;
3160  auto values_dofs_actual_ptr = values_array;
3161 
3162  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
3163  for (unsigned int c = 0; c < n_components; ++c)
3164  {
3165  eval.template values<0, false, false>(values_quad_ptr,
3166  values_dofs_actual_ptr);
3167 
3168  values_quad_ptr += n_q_points;
3169  values_dofs_actual_ptr += n_dofs;
3170  }
3171  }
3172 
3173  if (integrate_gradients)
3174  {
3175  auto gradients_quad_ptr = gradients_quad;
3176  auto values_dofs_actual_ptr = values_array;
3177 
3178  std::array<const VectorizedArrayType *, dim> shape_gradients;
3179  for (unsigned int d = 0; d < dim; ++d)
3180  shape_gradients[d] = &shape_info.shape_gradients_face(
3181  face_no, face_orientation, d, 0);
3182 
3183  for (unsigned int c = 0; c < n_components; ++c)
3184  {
3185  for (unsigned int d = 0; d < dim; ++d)
3186  {
3187  Eval eval(nullptr,
3188  shape_gradients[d],
3189  nullptr,
3190  n_dofs,
3191  n_q_points);
3192 
3193  if ((integrate_values == false) && d == 0)
3194  eval.template gradients<0, false, false>(
3195  gradients_quad_ptr, values_dofs_actual_ptr);
3196  else
3197  eval.template gradients<0, false, true>(
3198  gradients_quad_ptr, values_dofs_actual_ptr);
3199 
3200  gradients_quad_ptr += n_q_points;
3201  }
3202  values_dofs_actual_ptr += n_dofs;
3203  }
3204  }
3205 
3206  if (integrate_hessians)
3207  {
3208  AssertThrow(false, ExcNotImplemented());
3209  }
3210 
3211  return true;
3212  }
3213 
3214  if (face_orientation)
3216  n_components,
3217  face_orientation,
3218  orientation_map,
3219  true,
3220  integrate_values,
3221  integrate_gradients,
3222  integrate_hessians,
3223  data.n_q_points_face,
3224  scratch_data,
3225  values_quad,
3226  gradients_quad,
3227  hessians_quad);
3228 
3229  constexpr unsigned int static_dofs_per_face =
3230  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
3232  const unsigned int dofs_per_face =
3233  fe_degree > -1 ?
3234  static_dofs_per_face :
3235  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
3236 
3237  VectorizedArrayType *temp1 = scratch_data;
3238 
3239  const unsigned int n_q_points_1d_actual =
3240  fe_degree > -1 ? n_q_points_1d : 0;
3241  if (fe_degree > -1 &&
3245  true,
3246  dim,
3247  fe_degree,
3248  n_q_points_1d_actual,
3249  VectorizedArrayType>::integrate_in_face(n_components,
3250  data,
3251  temp1,
3252  values_quad,
3253  gradients_quad,
3254  hessians_quad,
3255  scratch_data +
3256  3 * n_components *
3257  dofs_per_face,
3258  integrate_values,
3259  integrate_gradients,
3260  integrate_hessians,
3261  subface_index);
3262  else
3264  false,
3265  dim,
3266  fe_degree,
3267  n_q_points_1d_actual,
3268  VectorizedArrayType>::integrate_in_face(n_components,
3269  data,
3270  temp1,
3271  values_quad,
3272  gradients_quad,
3273  hessians_quad,
3274  scratch_data +
3275  3 * n_components *
3276  dofs_per_face,
3277  integrate_values,
3278  integrate_gradients,
3279  integrate_hessians,
3280  subface_index);
3281 
3283  template interpolate<false, false>(n_components,
3284  data,
3285  temp1,
3286  values_array,
3287  integrate_gradients,
3288  integrate_hessians,
3289  face_no);
3290  return false;
3291  }
3292  };
3293 
3294 
3295 
3296  template <int n_face_orientations, typename Processor>
3297  static bool
3299  {
3300  auto n_components = proc.n_components;
3301  auto integrate = proc.integrate;
3302  auto global_vector_ptr = proc.global_vector_ptr;
3303  auto &sm_ptr = proc.sm_ptr;
3304  auto &data = proc.data;
3305  auto &dof_info = proc.dof_info;
3306  auto values_quad = proc.values_quad;
3307  auto gradients_quad = proc.gradients_quad;
3308  auto hessians_quad = proc.hessians_quad;
3309  auto scratch_data = proc.scratch_data;
3310  auto do_values = proc.do_values;
3311  auto do_gradients = proc.do_gradients;
3312  auto do_hessians = proc.do_hessians;
3313  auto active_fe_index = proc.active_fe_index;
3314  auto first_selected_component = proc.first_selected_component;
3315  auto cells = proc.cells;
3316  auto face_nos = proc.face_nos;
3317  auto subface_index = proc.subface_index;
3318  auto dof_access_index = proc.dof_access_index;
3319  auto face_orientations = proc.face_orientations;
3320  auto &orientation_map = proc.orientation_map;
3321 
3322  static const int dim = Processor::dim_;
3323  static const int fe_degree = Processor::fe_degree_;
3324  using VectorizedArrayType = typename Processor::VectorizedArrayType_;
3325 
3326  using Number = typename Processor::Number_;
3327  using Number2_ = typename Processor::Number2_;
3328 
3329  const unsigned int cell = cells[0];
3330 
3331  // In the case of integration, we do not need to reshuffle the
3332  // data at the quadrature points to adjust for the face
3333  // orientation if the shape functions are nodal at the cell
3334  // boundaries (and we only requested the integration of the
3335  // values) or Hermite shape functions are used. These cases are
3336  // handled later when the values are written back into the
3337  // glrobal vector.
3338  if (integrate &&
3339  (face_orientations[0] > 0 &&
3340  (subface_index < GeometryInfo<dim>::max_children_per_cell ||
3341  !(((do_gradients == false &&
3342  data.data.front().nodal_at_cell_boundaries == true &&
3343  fe_degree > 0) ||
3344  (data.element_type ==
3346  fe_degree > 1)) &&
3347  (dof_info.index_storage_variants[dof_access_index][cell] ==
3349  interleaved_contiguous ||
3350  dof_info.index_storage_variants[dof_access_index][cell] ==
3352  interleaved_contiguous_strided ||
3353  dof_info.index_storage_variants[dof_access_index][cell] ==
3355  interleaved_contiguous_mixed_strides ||
3356  dof_info.index_storage_variants[dof_access_index][cell] ==
3358  contiguous)))))
3359  {
3360  AssertDimension(n_face_orientations, 1);
3362  n_components,
3363  face_orientations[0],
3364  orientation_map,
3365  true,
3366  do_values,
3367  do_gradients,
3368  do_hessians,
3369  data.n_q_points_face,
3370  scratch_data,
3371  values_quad,
3372  gradients_quad,
3373  hessians_quad);
3374  }
3375 
3376  // we know that the gradient weights for the Hermite case on the
3377  // right (side==1) are the negative from the value at the left
3378  // (side==0), so we only read out one of them.
3379  VectorizedArrayType grad_weight =
3380  (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
3381  data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ?
3382  data.data.front()
3383  .shape_data_on_face[0][fe_degree + (integrate ?
3384  (2 - (face_nos[0] % 2)) :
3385  (1 + (face_nos[0] % 2)))] :
3386  VectorizedArrayType(0.0 /*dummy*/);
3387 
3388  constexpr unsigned int static_dofs_per_component =
3389  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) :
3391  constexpr unsigned int static_dofs_per_face =
3392  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
3394  const unsigned int dofs_per_face =
3395  fe_degree > -1 ? static_dofs_per_face :
3396  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
3397 
3398  VectorizedArrayType *temp1 = scratch_data;
3399 
3400  const unsigned int dummy = 0;
3401 
3402  // re-orientation
3403  std::array<const unsigned int *, n_face_orientations> orientation = {};
3404 
3405  if (n_face_orientations == 1)
3406  orientation[0] = (data.data.front().nodal_at_cell_boundaries == true) ?
3407  &data.face_orientations[face_orientations[0]][0] :
3408  &dummy;
3409  else
3410  {
3411  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3412  {
3413  // the loop breaks once an invalid_unsigned_int is hit for
3414  // all cases except the exterior faces in the ECL loop (where
3415  // some faces might be at the boundaries but others not)
3416  if (cells[v] == numbers::invalid_unsigned_int)
3417  continue;
3418 
3419  orientation[v] =
3420  (data.data.front().nodal_at_cell_boundaries == true) ?
3421  &data.face_orientations[face_orientations[v]][0] :
3422  &dummy;
3423  }
3424  }
3425 
3426  // face_to_cell_index_hermite
3427  std::array<const unsigned int *, n_face_orientations> index_array_hermite =
3428  {};
3429 
3430  if (n_face_orientations == 1)
3431  index_array_hermite[0] =
3432  (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
3433  data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ?
3434  &data.face_to_cell_index_hermite(face_nos[0], 0) :
3435  &dummy;
3436 
3437  if (n_face_orientations > 1 &&
3438  data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
3439  data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite)
3440  {
3441  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3442  {
3443  if (cells[v] == numbers::invalid_unsigned_int)
3444  continue;
3445 
3446  grad_weight[v] =
3447  data.data.front().shape_data_on_face
3448  [0][fe_degree + (integrate ? (2 - (face_nos[v] % 2)) :
3449  (1 + (face_nos[v] % 2)))][v];
3450 
3451  index_array_hermite[v] =
3452  &data.face_to_cell_index_hermite(face_nos[v], 0);
3453  }
3454  }
3455 
3456  // face_to_cell_index_nodal
3457  std::array<const unsigned int *, n_face_orientations> index_array_nodal =
3458  {};
3459 
3460  if (n_face_orientations == 1)
3461  index_array_nodal[0] =
3462  (data.data.front().nodal_at_cell_boundaries == true) ?
3463  &data.face_to_cell_index_nodal(face_nos[0], 0) :
3464  &dummy;
3465 
3466  if (n_face_orientations > 1 &&
3467  (data.data.front().nodal_at_cell_boundaries == true))
3468  {
3469  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3470  {
3471  if (cells[v] == numbers::invalid_unsigned_int)
3472  continue;
3473 
3474  index_array_nodal[v] =
3475  &data.face_to_cell_index_nodal(face_nos[v], 0);
3476  }
3477  }
3478 
3479  const auto reorientate = [&](const unsigned int v, const unsigned int i) {
3480  return (dim < 3 ||
3481  face_orientations[n_face_orientations == 1 ? 0 : v] == 0 ||
3482  subface_index < GeometryInfo<dim>::max_children_per_cell) ?
3483  i :
3484  orientation[v][i];
3485  };
3486 
3487  // this variable keeps track of whether we are able to directly write
3488  // the results into the result (function returns true) or not, requiring
3489  // an additional call to another function
3490  bool accesses_global_vector = true;
3491 
3492  for (unsigned int comp = 0; comp < n_components; ++comp)
3493  {
3494  if (integrate)
3495  proc.in_face_operation(temp1, comp);
3496 
3497  // we can only use the fast functions if we know the polynomial degree
3498  // as a template parameter (fe_degree != -1), and it only makes sense
3499  // to use the functions for at least linear functions for values on
3500  // the faces and quadratic functions for gradients on the faces, so
3501  // include the switch here
3502  if (((do_gradients == false &&
3503  data.data.front().nodal_at_cell_boundaries == true &&
3504  fe_degree > 0) ||
3505  (data.element_type ==
3507  fe_degree > 1)) &&
3508  do_hessians == false)
3509  {
3510  // case 1: contiguous and interleaved indices
3511  if (n_face_orientations == 1 &&
3512  dof_info.index_storage_variants[dof_access_index][cell] ==
3514  interleaved_contiguous)
3515  {
3516  AssertDimension(n_face_orientations, 1);
3517 
3519  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
3520  VectorizedArrayType::size());
3521  Number2_ *vector_ptr =
3522  global_vector_ptr +
3523  dof_info.dof_indices_contiguous[dof_access_index]
3524  [cell *
3525  VectorizedArrayType::size()] +
3526  (dof_info
3527  .component_dof_indices_offset[active_fe_index]
3528  [first_selected_component] +
3529  comp * static_dofs_per_component) *
3530  VectorizedArrayType::size();
3531 
3532  if (fe_degree > 1 && do_gradients == true)
3533  {
3534  for (unsigned int i = 0; i < dofs_per_face; ++i)
3535  {
3536  if (n_face_orientations == 1)
3537  {
3538  const unsigned int ind1 =
3539  index_array_hermite[0][2 * i];
3540  const unsigned int ind2 =
3541  index_array_hermite[0][2 * i + 1];
3542  AssertIndexRange(ind1,
3543  data.dofs_per_component_on_cell);
3544  AssertIndexRange(ind2,
3545  data.dofs_per_component_on_cell);
3546  const unsigned int i_ = reorientate(0, i);
3547  proc.hermite_grad_vectorized(
3548  temp1[i_],
3549  temp1[i_ + dofs_per_face],
3550  vector_ptr + ind1 * VectorizedArrayType::size(),
3551  vector_ptr + ind2 * VectorizedArrayType::size(),
3552  grad_weight);
3553  }
3554  else
3555  {
3556  Assert(false, ExcNotImplemented());
3557  }
3558  }
3559  }
3560  else
3561  {
3562  for (unsigned int i = 0; i < dofs_per_face; ++i)
3563  {
3564  if (n_face_orientations == 1)
3565  {
3566  const unsigned int i_ = reorientate(0, i);
3567  const unsigned int ind = index_array_nodal[0][i];
3568  proc.value_vectorized(
3569  temp1[i_],
3570  vector_ptr + ind * VectorizedArrayType::size());
3571  }
3572  else
3573  {
3574  Assert(false, ExcNotImplemented());
3575  }
3576  }
3577  }
3578  }
3579 
3580  // case 2: contiguous and interleaved indices with fixed stride
3581  else if (n_face_orientations == 1 &&
3582  dof_info.index_storage_variants[dof_access_index][cell] ==
3584  interleaved_contiguous_strided)
3585  {
3586  AssertDimension(n_face_orientations, 1);
3587 
3589  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
3590  VectorizedArrayType::size());
3591  const unsigned int *indices =
3592  &dof_info.dof_indices_contiguous[dof_access_index]
3593  [cell *
3594  VectorizedArrayType::size()];
3595  Number2_ *vector_ptr =
3596  global_vector_ptr +
3597  (comp * static_dofs_per_component +
3598  dof_info
3599  .component_dof_indices_offset[active_fe_index]
3600  [first_selected_component]) *
3601  VectorizedArrayType::size();
3602  if (fe_degree > 1 && do_gradients == true)
3603  {
3604  for (unsigned int i = 0; i < dofs_per_face; ++i)
3605  {
3606  if (n_face_orientations == 1)
3607  {
3608  const unsigned int i_ = reorientate(0, i);
3609  const unsigned int ind1 =
3610  index_array_hermite[0][2 * i] *
3611  VectorizedArrayType::size();
3612  const unsigned int ind2 =
3613  index_array_hermite[0][2 * i + 1] *
3614  VectorizedArrayType::size();
3615  proc.hermite_grad_vectorized_indexed(
3616  temp1[i_],
3617  temp1[i_ + dofs_per_face],
3618  vector_ptr + ind1,
3619  vector_ptr + ind2,
3620  grad_weight,
3621  indices,
3622  indices);
3623  }
3624  else
3625  {
3626  Assert(false, ExcNotImplemented());
3627  }
3628  }
3629  }
3630  else
3631  {
3632  for (unsigned int i = 0; i < dofs_per_face; ++i)
3633  {
3634  if (n_face_orientations == 1)
3635  {
3636  const unsigned int i_ = reorientate(0, i);
3637  const unsigned int ind =
3638  index_array_nodal[0][i] *
3639  VectorizedArrayType::size();
3640  proc.value_vectorized_indexed(temp1[i_],
3641  vector_ptr + ind,
3642  indices);
3643  }
3644  else
3645  {
3646  Assert(false, ExcNotImplemented());
3647  }
3648  }
3649  }
3650  }
3651 
3652  // case 3: contiguous and interleaved indices with mixed stride
3653  else if (n_face_orientations == 1 &&
3654  dof_info.index_storage_variants[dof_access_index][cell] ==
3656  interleaved_contiguous_mixed_strides)
3657  {
3658  AssertDimension(n_face_orientations, 1);
3659 
3660  const unsigned int *strides =
3661  &dof_info.dof_indices_interleave_strides
3662  [dof_access_index][cell * VectorizedArrayType::size()];
3663  unsigned int indices[VectorizedArrayType::size()];
3664  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3665  indices[v] =
3666  dof_info.dof_indices_contiguous
3667  [dof_access_index]
3668  [cell * VectorizedArrayType::size() + v] +
3669  (dof_info
3670  .component_dof_indices_offset[active_fe_index]
3671  [first_selected_component] +
3672  comp * static_dofs_per_component) *
3673  strides[v];
3674  const unsigned int n_filled_lanes =
3675  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
3676 
3677  if (fe_degree > 1 && do_gradients == true)
3678  {
3679  if (n_filled_lanes == VectorizedArrayType::size())
3680  for (unsigned int i = 0; i < dofs_per_face; ++i)
3681  {
3682  if (n_face_orientations == 1)
3683  {
3684  const unsigned int i_ = reorientate(0, i);
3685  unsigned int ind1[VectorizedArrayType::size()];
3687  for (unsigned int v = 0;
3688  v < VectorizedArrayType::size();
3689  ++v)
3690  ind1[v] =
3691  indices[v] +
3692  index_array_hermite[0 /*TODO*/][2 * i] *
3693  strides[v];
3694  unsigned int ind2[VectorizedArrayType::size()];
3696  for (unsigned int v = 0;
3697  v < VectorizedArrayType::size();
3698  ++v)
3699  ind2[v] =
3700  indices[v] +
3701  index_array_hermite[0 /*TODO*/][2 * i + 1] *
3702  strides[v];
3703  proc.hermite_grad_vectorized_indexed(
3704  temp1[i_],
3705  temp1[i_ + dofs_per_face],
3706  global_vector_ptr,
3707  global_vector_ptr,
3708  grad_weight,
3709  ind1,
3710  ind2);
3711  }
3712  else
3713  {
3714  Assert(false, ExcNotImplemented());
3715  }
3716  }
3717  else
3718  {
3719  if (integrate == false)
3720  for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
3721  temp1[i] = VectorizedArrayType();
3722 
3723  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3724  for (unsigned int i = 0; i < dofs_per_face; ++i)
3725  {
3726  const unsigned int i_ =
3727  reorientate(n_face_orientations == 1 ? 0 : v,
3728  i);
3729  proc.hermite_grad(
3730  temp1[i_][v],
3731  temp1[i_ + dofs_per_face][v],
3732  global_vector_ptr
3733  [indices[v] +
3734  index_array_hermite
3735  [n_face_orientations == 1 ? 0 : v]
3736  [2 * i] *
3737  strides[v]],
3738  global_vector_ptr
3739  [indices[v] +
3740  index_array_hermite
3741  [n_face_orientations == 1 ? 0 : v]
3742  [2 * i + 1] *
3743  strides[v]],
3744  grad_weight[n_face_orientations == 1 ? 0 : v]);
3745  }
3746  }
3747  }
3748  else
3749  {
3750  if (n_filled_lanes == VectorizedArrayType::size())
3751  for (unsigned int i = 0; i < dofs_per_face; ++i)
3752  {
3753  if (n_face_orientations == 1)
3754  {
3755  unsigned int ind[VectorizedArrayType::size()];
3757  for (unsigned int v = 0;
3758  v < VectorizedArrayType::size();
3759  ++v)
3760  ind[v] = indices[v] +
3761  index_array_nodal[0][i] * strides[v];
3762  const unsigned int i_ = reorientate(0, i);
3763  proc.value_vectorized_indexed(temp1[i_],
3764  global_vector_ptr,
3765  ind);
3766  }
3767  else
3768  {
3769  Assert(false, ExcNotImplemented());
3770  }
3771  }
3772  else
3773  {
3774  if (integrate == false)
3775  for (unsigned int i = 0; i < dofs_per_face; ++i)
3776  temp1[i] = VectorizedArrayType();
3777 
3778  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3779  for (unsigned int i = 0; i < dofs_per_face; ++i)
3780  proc.value(
3781  temp1[reorientate(
3782  n_face_orientations == 1 ? 0 : v, i)][v],
3783  global_vector_ptr
3784  [indices[v] +
3785  index_array_nodal
3786  [n_face_orientations == 1 ? 0 : v][i] *
3787  strides[v]]);
3788  }
3789  }
3790  }
3791 
3792  // case 4: contiguous indices without interleaving
3793  else if ((n_face_orientations > 1 ||
3794  dof_info.index_storage_variants[dof_access_index][cell] ==
3796  contiguous))
3797  {
3798  const unsigned int *indices =
3799  &dof_info.dof_indices_contiguous[dof_access_index]
3800  [cell *
3801  VectorizedArrayType::size()];
3802  Number2_ *vector_ptr =
3803  global_vector_ptr + comp * static_dofs_per_component +
3804  dof_info
3805  .component_dof_indices_offset[active_fe_index]
3806  [first_selected_component];
3807 
3808  const unsigned int n_filled_lanes =
3809  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
3810 
3811  const bool vectorization_possible =
3812  (n_face_orientations == 1) &&
3813  (n_filled_lanes == VectorizedArrayType::size()) &&
3814  (sm_ptr != nullptr);
3815 
3816  std::array<Number2_ *, VectorizedArrayType::size()>
3817  vector_ptrs = {};
3818 
3819  if (vectorization_possible == false)
3820  {
3821  if (n_face_orientations == 1)
3822  {
3823  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3824  if (sm_ptr == nullptr)
3825  {
3826  vector_ptrs[v] = vector_ptr + indices[v];
3827  }
3828  else
3829  {
3830  const auto &temp =
3831  dof_info.dof_indices_contiguous_sm
3832  [dof_access_index]
3833  [cell * VectorizedArrayType::size() + v];
3834  vector_ptrs[v] = const_cast<Number *>(
3835  sm_ptr->operator[](temp.first).data() +
3836  temp.second + comp * static_dofs_per_component +
3837  dof_info.component_dof_indices_offset
3838  [active_fe_index][first_selected_component]);
3839  }
3840  }
3841  else if (n_face_orientations == VectorizedArrayType::size())
3842  {
3843  for (unsigned int v = 0;
3844  v < VectorizedArrayType::size();
3845  ++v)
3846  if (cells[v] != numbers::invalid_unsigned_int)
3847  {
3848  if (sm_ptr == nullptr)
3849  {
3850  vector_ptrs[v] =
3851  vector_ptr +
3852  dof_info
3853  .dof_indices_contiguous[dof_access_index]
3854  [cells[v]];
3855  }
3856  else
3857  {
3858  const auto &temp =
3859  dof_info.dof_indices_contiguous_sm
3860  [dof_access_index][cells[v]];
3861  vector_ptrs[v] = const_cast<Number *>(
3862  sm_ptr->operator[](temp.first).data() +
3863  temp.second +
3864  comp * static_dofs_per_component +
3865  dof_info.component_dof_indices_offset
3866  [active_fe_index]
3867  [first_selected_component]);
3868  }
3869  }
3870  }
3871  else
3872  {
3873  Assert(false, ExcNotImplemented());
3874  }
3875  }
3876 
3877  if (do_gradients == true &&
3878  data.element_type ==
3880  {
3881  if (vectorization_possible)
3882  for (unsigned int i = 0; i < dofs_per_face; ++i)
3883  {
3884  const unsigned int ind1 =
3885  index_array_hermite[0][2 * i];
3886  const unsigned int ind2 =
3887  index_array_hermite[0][2 * i + 1];
3888  const unsigned int i_ = reorientate(0, i);
3889 
3890  proc.hermite_grad_vectorized_indexed(
3891  temp1[i_],
3892  temp1[i_ + dofs_per_face],
3893  vector_ptr + ind1,
3894  vector_ptr + ind2,
3895  grad_weight,
3896  indices,
3897  indices);
3898  }
3899  else if (n_face_orientations == 1)
3900  for (unsigned int i = 0; i < dofs_per_face; ++i)
3901  {
3902  const unsigned int ind1 =
3903  index_array_hermite[0][2 * i];
3904  const unsigned int ind2 =
3905  index_array_hermite[0][2 * i + 1];
3906  const unsigned int i_ = reorientate(0, i);
3907 
3908  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3909  proc.hermite_grad(temp1[i_][v],
3910  temp1[i_ + dofs_per_face][v],
3911  vector_ptrs[v][ind1],
3912  vector_ptrs[v][ind2],
3913  grad_weight[v]);
3914 
3915  if (integrate == false)
3916  for (unsigned int v = n_filled_lanes;
3917  v < VectorizedArrayType::size();
3918  ++v)
3919  {
3920  temp1[i_][v] = 0.0;
3921  temp1[i_ + dofs_per_face][v] = 0.0;
3922  }
3923  }
3924  else
3925  {
3926  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3927  for (unsigned int i = 0; i < dofs_per_face; ++i)
3928  proc.hermite_grad(
3929  temp1[reorientate(v, i)][v],
3930  temp1[reorientate(v, i) + dofs_per_face][v],
3931  vector_ptrs[v][index_array_hermite[v][2 * i]],
3932  vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
3933  grad_weight[v]);
3934  }
3935  }
3936  else
3937  {
3938  if (vectorization_possible)
3939  for (unsigned int i = 0; i < dofs_per_face; ++i)
3940  {
3941  const unsigned int ind = index_array_nodal[0][i];
3942  const unsigned int i_ = reorientate(0, i);
3943 
3944  proc.value_vectorized_indexed(temp1[i_],
3945  vector_ptr + ind,
3946  indices);
3947  }
3948  else if (n_face_orientations == 1)
3949  for (unsigned int i = 0; i < dofs_per_face; ++i)
3950  {
3951  const unsigned int ind = index_array_nodal[0][i];
3952  const unsigned int i_ = reorientate(0, i);
3953 
3954  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3955  proc.value(temp1[i_][v], vector_ptrs[v][ind]);
3956 
3957  if (integrate == false)
3958  for (unsigned int v = n_filled_lanes;
3959  v < VectorizedArrayType::size();
3960  ++v)
3961  temp1[i_][v] = 0.0;
3962  }
3963  else
3964  for (unsigned int i = 0; i < dofs_per_face; ++i)
3965  {
3966  for (unsigned int v = 0;
3967  v < VectorizedArrayType::size();
3968  ++v)
3969  if (cells[v] != numbers::invalid_unsigned_int)
3970  proc.value(
3971  temp1[reorientate(v, i)][v],
3972  vector_ptrs[v][index_array_nodal[v][i]]);
3973  }
3974  }
3975  }
3976  else
3977  {
3978  // case 5: default vector access
3979  // for the integrate_scatter path (integrate == true), we
3980  // need to only prepare the data in this function for all
3981  // components to later call distribute_local_to_global();
3982  // for the gather_evaluate path (integrate == false), we
3983  // instead want to leave early because we need to get the
3984  // vector data from somewhere else
3985  proc.default_operation(temp1, comp);
3986  if (integrate)
3987  accesses_global_vector = false;
3988  else
3989  return false;
3990  }
3991  }
3992  else
3993  {
3994  // case 5: default vector access
3995  proc.default_operation(temp1, comp);
3996  if (integrate)
3997  accesses_global_vector = false;
3998  else
3999  return false;
4000  }
4001 
4002  if (!integrate)
4003  proc.in_face_operation(temp1, comp);
4004  }
4005 
4006  if (!integrate &&
4007  (face_orientations[0] > 0 &&
4009  {
4010  AssertDimension(n_face_orientations, 1);
4012  n_components,
4013  face_orientations[0],
4014  orientation_map,
4015  false,
4016  do_values,
4017  do_gradients,
4018  do_hessians,
4019  data.n_q_points_face,
4020  scratch_data,
4021  values_quad,
4022  gradients_quad,
4023  hessians_quad);
4024  }
4025 
4026  return accesses_global_vector;
4027  }
4028 
4029 
4030 
4031  template <int dim,
4032  typename Number,
4033  typename VectorizedArrayType,
4034  typename Number2 = Number>
4036  {
4037  template <int fe_degree, int n_q_points_1d>
4038  static bool
4039  run(const unsigned int n_components,
4040  const unsigned int n_face_orientations,
4041  const Number2 * src_ptr,
4042  const std::vector<ArrayView<const Number>> *sm_ptr,
4044  const MatrixFreeFunctions::DoFInfo & dof_info,
4045  VectorizedArrayType * values_quad,
4046  VectorizedArrayType *gradients_quad,
4047  VectorizedArrayType *hessians_quad,
4048  VectorizedArrayType *scratch_data,
4049  const bool evaluate_values,
4050  const bool evaluate_gradients,
4051  const bool evaluate_hessians,
4052  const unsigned int active_fe_index,
4053  const unsigned int first_selected_component,
4054  const std::array<unsigned int, VectorizedArrayType::size()> cells,
4055  const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
4056  const unsigned int subface_index,
4057  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
4058  const std::array<unsigned int, VectorizedArrayType::size()>
4059  face_orientations,
4060  const Table<2, unsigned int> &orientation_map)
4061  {
4062  if (src_ptr == nullptr)
4063  return false;
4064 
4066  return false;
4067 
4068  (void)sm_ptr;
4069 
4070  Processor<fe_degree, n_q_points_1d> p(n_components,
4071  false,
4072  src_ptr,
4073  sm_ptr,
4074  data,
4075  dof_info,
4076  values_quad,
4077  gradients_quad,
4078  hessians_quad,
4079  scratch_data,
4080  evaluate_values,
4081  evaluate_gradients,
4082  evaluate_hessians,
4083  active_fe_index,
4084  first_selected_component,
4085  cells,
4086  face_nos,
4087  subface_index,
4088  dof_access_index,
4089  face_orientations,
4090  orientation_map);
4091 
4092  if (n_face_orientations == VectorizedArrayType::size())
4093  return fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
4094  p);
4095  else
4096  return fe_face_evaluation_process_and_io<1>(p);
4097  }
4098 
4099  private:
4100  template <int fe_degree, int n_q_points_1d>
4101  struct Processor
4102  {
4103  static const int dim_ = dim;
4104  static const int fe_degree_ = fe_degree;
4105  static const int n_q_points_1d_ = n_q_points_1d;
4106  using VectorizedArrayType_ = VectorizedArrayType;
4107  using Number_ = Number;
4108  using Number2_ = const Number2;
4109 
4111  const unsigned int n_components,
4112  const bool integrate,
4113  const Number2 * global_vector_ptr,
4114  const std::vector<ArrayView<const Number>> *sm_ptr,
4116  const MatrixFreeFunctions::DoFInfo & dof_info,
4117  VectorizedArrayType * values_quad,
4118  VectorizedArrayType *gradients_quad,
4119  VectorizedArrayType *hessians_quad,
4120  VectorizedArrayType *scratch_data,
4121  const bool do_values,
4122  const bool do_gradients,
4123  const bool do_hessians,
4124  const unsigned int active_fe_index,
4125  const unsigned int first_selected_component,
4126  const std::array<unsigned int, VectorizedArrayType::size()> cells,
4127  const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
4128  const unsigned int subface_index,
4129  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
4130  const std::array<unsigned int, VectorizedArrayType::size()>
4131  face_orientations,
4132  const Table<2, unsigned int> &orientation_map)
4133  : n_components(n_components)
4134  , integrate(integrate)
4135  , global_vector_ptr(global_vector_ptr)
4136  , sm_ptr(sm_ptr)
4137  , data(data)
4138  , dof_info(dof_info)
4139  , values_quad(values_quad)
4140  , gradients_quad(gradients_quad)
4141  , hessians_quad(hessians_quad)
4142  , scratch_data(scratch_data)
4143  , do_values(do_values)
4144  , do_gradients(do_gradients)
4145  , do_hessians(do_hessians)
4146  , active_fe_index(active_fe_index)
4147  , first_selected_component(first_selected_component)
4148  , cells(cells)
4149  , face_nos(face_nos)
4150  , subface_index(subface_index)
4151  , dof_access_index(dof_access_index)
4152  , face_orientations(face_orientations)
4153  , orientation_map(orientation_map)
4154  {}
4155 
4156  template <typename T0, typename T1, typename T2>
4157  void
4159  T0 & temp_2,
4160  const T1 src_ptr_1,
4161  const T1 src_ptr_2,
4162  const T2 &grad_weight)
4163  {
4164  do_vectorized_read(src_ptr_1, temp_1);
4165  do_vectorized_read(src_ptr_2, temp_2);
4166  temp_2 = grad_weight * (temp_1 - temp_2);
4167  }
4168 
4169  template <typename T1, typename T2>
4170  void
4171  value_vectorized(T1 &temp, const T2 src_ptr)
4172  {
4173  do_vectorized_read(src_ptr, temp);
4174  }
4175 
4176  template <typename T0, typename T1, typename T2, typename T3>
4177  void
4179  T0 & temp_2,
4180  const T1 src_ptr_1,
4181  const T1 src_ptr_2,
4182  const T2 &grad_weight,
4183  const T3 &indices_1,
4184  const T3 &indices_2)
4185  {
4186  do_vectorized_gather(src_ptr_1, indices_1, temp_1);
4187  do_vectorized_gather(src_ptr_2, indices_2, temp_2);
4188  temp_2 = grad_weight * (temp_1 - temp_2);
4189  }
4190 
4191  template <typename T0, typename T1, typename T2>
4192  void
4193  value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
4194  {
4195  do_vectorized_gather(src_ptr, indices, temp);
4196  }
4197 
4198  template <typename T0, typename T1, typename T2>
4199  void
4200  hermite_grad(T0 & temp_1,
4201  T0 & temp_2,
4202  const T1 &src_ptr_1,
4203  const T2 &src_ptr_2,
4204  const T2 &grad_weight)
4205  {
4206  // case 3a)
4207  temp_1 = src_ptr_1;
4208  temp_2 = grad_weight * (temp_1 - src_ptr_2);
4209  }
4210 
4211  template <typename T1, typename T2>
4212  void
4213  value(T1 &temp, const T2 &src_ptr)
4214  {
4215  // case 3b)
4216  temp = src_ptr;
4217  }
4218 
4219  template <typename T1>
4220  void
4221  default_operation(const T1 &, const unsigned int)
4222  {
4223  // case 5)
4224  }
4225 
4226  template <typename T1>
4227  void
4228  in_face_operation(T1 &temp1, const unsigned int comp)
4229  {
4230  const unsigned int dofs_per_face =
4231  fe_degree > -1 ?
4232  Utilities::pow(fe_degree + 1, dim - 1) :
4233  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
4234  const unsigned int n_q_points =
4235  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
4236  data.n_q_points_face;
4237  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4238  if (fe_degree > -1 &&
4239  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
4241  FEFaceEvaluationImpl<true,
4242  dim,
4243  fe_degree,
4244  n_q_points_1d,
4245  VectorizedArrayType>::
4246  evaluate_in_face(/* n_components */ 1,
4247  data,
4248  temp1,
4249  values_quad + comp * n_q_points,
4250  gradients_quad + comp * dim * n_q_points,
4251  hessians_quad + comp * hdim * n_q_points,
4252  scratch_data + 3 * dofs_per_face,
4253  do_values,
4254  do_gradients,
4255  do_hessians,
4256  subface_index);
4257  else
4258  FEFaceEvaluationImpl<false,
4259  dim,
4260  fe_degree,
4261  n_q_points_1d,
4262  VectorizedArrayType>::
4263  evaluate_in_face(/* n_components */ 1,
4264  data,
4265  temp1,
4266  values_quad + comp * n_q_points,
4267  gradients_quad + comp * dim * n_q_points,
4268  hessians_quad + comp * hdim * n_q_points,
4269  scratch_data + 3 * dofs_per_face,
4270  do_values,
4271  do_gradients,
4272  do_hessians,
4273  subface_index);
4274  }
4275 
4276  const unsigned int n_components;
4277  const bool integrate;
4278  const Number2 * global_vector_ptr;
4279  const std::vector<ArrayView<const Number>> *sm_ptr;
4282  VectorizedArrayType * values_quad;
4283  VectorizedArrayType * gradients_quad;
4284  VectorizedArrayType * hessians_quad;
4285  VectorizedArrayType * scratch_data;
4286  const bool do_values;
4287  const bool do_gradients;
4288  const bool do_hessians;
4289  const unsigned int active_fe_index;
4290  const unsigned int first_selected_component;
4291  const std::array<unsigned int, VectorizedArrayType::size()> cells;
4292  const std::array<unsigned int, VectorizedArrayType::size()> face_nos;
4293  const unsigned int subface_index;
4295  const std::array<unsigned int, VectorizedArrayType::size()>
4298  };
4299  };
4300 
4301  template <int dim,
4302  typename Number,
4303  typename VectorizedArrayType,
4304  typename Number2 = Number>
4306  {
4307  template <int fe_degree, int n_q_points_1d>
4308  static bool
4309  run(const unsigned int n_components,
4310  const unsigned int n_face_orientations,
4311  Number2 * dst_ptr,
4312  const std::vector<ArrayView<const Number2>> *sm_ptr,
4314  const MatrixFreeFunctions::DoFInfo & dof_info,
4315  VectorizedArrayType * values_array,
4316  VectorizedArrayType * values_quad,
4317  VectorizedArrayType *gradients_quad,
4318  VectorizedArrayType *hessians_quad,
4319  VectorizedArrayType *scratch_data,
4320  const bool integrate_values,
4321  const bool integrate_gradients,
4322  const bool integrate_hessians,
4323  const unsigned int active_fe_index,
4324  const unsigned int first_selected_component,
4325  const std::array<unsigned int, VectorizedArrayType::size()> cells,
4326  const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
4327  const unsigned int subface_index,
4328  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
4329  const std::array<unsigned int, VectorizedArrayType::size()>
4330  face_orientations,
4331  const Table<2, unsigned int> &orientation_map)
4332  {
4333  (void)sm_ptr;
4334 
4335  if (dst_ptr == nullptr ||
4337  {
4338  AssertDimension(n_face_orientations, 1);
4339 
4340  // for block vectors simply integrate
4342  template run<fe_degree, n_q_points_1d>(n_components,
4343  data,
4344  values_array,
4345  values_quad,
4346  gradients_quad,
4347  hessians_quad,
4348  scratch_data,
4349  integrate_values,
4350  integrate_gradients,
4351  integrate_hessians,
4352  face_nos[0],
4353  subface_index,
4354  face_orientations[0],
4355  orientation_map);
4356 
4357  // default vector access
4358  return false;
4359  }
4360 
4361 
4362  Processor<fe_degree, n_q_points_1d> p(values_array,
4363  n_components,
4364  true,
4365  dst_ptr,
4366  sm_ptr,
4367  data,
4368  dof_info,
4369  values_quad,
4370  gradients_quad,
4371  hessians_quad,
4372  scratch_data,
4373  integrate_values,
4374  integrate_gradients,
4375  integrate_hessians,
4376  active_fe_index,
4377  first_selected_component,
4378  cells,
4379  face_nos,
4380  subface_index,
4381  dof_access_index,
4382  face_orientations,
4383  orientation_map);
4384 
4385  if (n_face_orientations == VectorizedArrayType::size())
4386  return fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
4387  p);
4388  else
4389  return fe_face_evaluation_process_and_io<1>(p);
4390  }
4391 
4392  private:
4393  template <int fe_degree, int n_q_points_1d>
4394  struct Processor
4395  {
4396  static const int dim_ = dim;
4397  static const int fe_degree_ = fe_degree;
4398  static const int n_q_points_1d_ = n_q_points_1d;
4399  using VectorizedArrayType_ = VectorizedArrayType;
4400  using Number_ = Number;
4401  using Number2_ = Number2;
4402 
4403 
4405  VectorizedArrayType * values_array,
4406  const unsigned int n_components,
4407  const bool integrate,
4408  Number2 * global_vector_ptr,
4409  const std::vector<ArrayView<const Number>> *sm_ptr,
4411  const MatrixFreeFunctions::DoFInfo & dof_info,
4412  VectorizedArrayType * values_quad,
4413  VectorizedArrayType *gradients_quad,
4414  VectorizedArrayType *hessians_quad,
4415  VectorizedArrayType *scratch_data,
4416  const bool do_values,
4417  const bool do_gradients,
4418  const bool do_hessians,
4419  const unsigned int active_fe_index,
4420  const unsigned int first_selected_component,
4421  const std::array<unsigned int, VectorizedArrayType::size()> cells,
4422  const std::array<unsigned int, VectorizedArrayType::size()> face_nos,
4423  const unsigned int subface_index,
4424  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
4425  const std::array<unsigned int, VectorizedArrayType::size()>
4426  face_orientations,
4427  const Table<2, unsigned int> &orientation_map)
4428  : values_array(values_array)
4429  , n_components(n_components)
4430  , integrate(integrate)
4431  , global_vector_ptr(global_vector_ptr)
4432  , sm_ptr(sm_ptr)
4433  , data(data)
4434  , dof_info(dof_info)
4435  , values_quad(values_quad)
4436  , gradients_quad(gradients_quad)
4437  , hessians_quad(hessians_quad)
4438  , scratch_data(scratch_data)
4439  , do_values(do_values)
4440  , do_gradients(do_gradients)
4441  , do_hessians(do_hessians)
4442  , active_fe_index(active_fe_index)
4443  , first_selected_component(first_selected_component)
4444  , cells(cells)
4445  , face_nos(face_nos)
4446  , subface_index(subface_index)
4447  , dof_access_index(dof_access_index)
4448  , face_orientations(face_orientations)
4449  , orientation_map(orientation_map)
4450  {}
4451 
4452  template <typename T0, typename T1, typename T2, typename T3, typename T4>
4453  void
4454  hermite_grad_vectorized(const T0 &temp_1,
4455  const T1 &temp_2,
4456  T2 dst_ptr_1,
4457  T3 dst_ptr_2,
4458  const T4 &grad_weight)
4459  {
4460  // case 1a)
4461  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
4462  const VectorizedArrayType grad = grad_weight * temp_2;
4463  do_vectorized_add(val, dst_ptr_1);
4464  do_vectorized_add(grad, dst_ptr_2);
4465  }
4466 
4467  template <typename T0, typename T1>
4468  void
4469  value_vectorized(const T0 &temp, T1 dst_ptr)
4470  {
4471  // case 1b)
4472  do_vectorized_add(temp, dst_ptr);
4473  }
4474 
4475  template <typename T0, typename T1, typename T2, typename T3>
4476  void
4478  const T0 &temp_2,
4479  T1 dst_ptr_1,
4480  T1 dst_ptr_2,
4481  const T2 &grad_weight,
4482  const T3 &indices_1,
4483  const T3 &indices_2)
4484  {
4485  // case 2a)
4486  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
4487  const VectorizedArrayType grad = grad_weight * temp_2;
4488  do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
4489  do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
4490  }
4491 
4492  template <typename T0, typename T1, typename T2>
4493  void
4494  value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
4495  {
4496  // case 2b)
4497  do_vectorized_scatter_add(temp, indices, dst_ptr);
4498  }
4499 
4500  template <typename T0, typename T1, typename T2>
4501  void
4502  hermite_grad(const T0 &temp_1,
4503  const T0 &temp_2,
4504  T1 & dst_ptr_1,
4505  T1 & dst_ptr_2,
4506  const T2 &grad_weight)
4507  {
4508  // case 3a)
4509  const Number val = temp_1 - grad_weight * temp_2;
4510  const Number grad = grad_weight * temp_2;
4511  dst_ptr_1 += val;
4512  dst_ptr_2 += grad;
4513  }
4514 
4515  template <typename T0, typename T1>
4516  void
4517  value(const T0 &temp, T1 &dst_ptr)
4518  {
4519  // case 3b)
4520  dst_ptr += temp;
4521  }
4522 
4523  template <typename T0>
4524  void
4525  default_operation(const T0 &temp1, const unsigned int comp)
4526  {
4527  // case 5: default vector access, must be handled separately, just do
4528  // the face-normal interpolation
4529 
4531  template interpolate<false, false>(
4532  /* n_components */ 1,
4533  data,
4534  temp1,
4535  values_array + comp * data.dofs_per_component_on_cell,
4536  do_gradients,
4537  do_hessians,
4538  face_nos[0]);
4539  }
4540 
4541  template <typename T0>
4542  void
4543  in_face_operation(T0 &temp1, const unsigned int comp)
4544  {
4545  const unsigned int dofs_per_face =
4546  fe_degree > -1 ?
4547  Utilities::pow(fe_degree + 1, dim - 1) :
4548  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
4549  const unsigned int n_q_points =
4550  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
4551  data.n_q_points_face;
4552  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4553  if (fe_degree > -1 &&
4554  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
4555  data.element_type <=
4558  dim,
4559  fe_degree,
4560  n_q_points_1d,
4561  VectorizedArrayType>::
4562  integrate_in_face(/* n_components */ 1,
4563  data,
4564  temp1,
4565  values_quad + comp * n_q_points,
4566  gradients_quad + dim * comp * n_q_points,
4567  hessians_quad + hdim * comp * n_q_points,
4568  scratch_data + 3 * dofs_per_face,
4569  do_values,
4570  do_gradients,
4571  do_hessians,
4572  subface_index);
4573  else
4575  dim,
4576  fe_degree,
4577  n_q_points_1d,
4578  VectorizedArrayType>::
4579  integrate_in_face(/* n_components */ 1,
4580  data,
4581  temp1,
4582  values_quad + comp * n_q_points,
4583  gradients_quad + dim * comp * n_q_points,
4584  hessians_quad + hdim * comp * n_q_points,
4585  scratch_data + 3 * dofs_per_face,
4586  do_values,
4587  do_gradients,
4588  do_hessians,
4589  subface_index);
4590  }
4591 
4592  VectorizedArrayType *values_array;
4593 
4594 
4595  const unsigned int n_components;
4596  const bool integrate;
4598  const std::vector<ArrayView<const Number>> *sm_ptr;
4601  VectorizedArrayType * values_quad;
4602  VectorizedArrayType * gradients_quad;
4603  VectorizedArrayType * hessians_quad;
4604  VectorizedArrayType * scratch_data;
4605  const bool do_values;
4606  const bool do_gradients;
4607  const bool do_hessians;
4608  const unsigned int active_fe_index;
4609  const unsigned int first_selected_component;
4610  const std::array<unsigned int, VectorizedArrayType::size()> cells;
4611  const std::array<unsigned int, VectorizedArrayType::size()> face_nos;
4612  const unsigned int subface_index;
4614  const std::array<unsigned int, VectorizedArrayType::size()>
4617  };
4618  };
4619 
4620 
4621 
4626  template <int dim, typename Number>
4628  {
4629  template <int fe_degree, int = 0>
4630  static bool
4631  run(const unsigned int n_components,
4632  const FEEvaluationBaseData<dim,
4633  typename Number::value_type,
4634  false,
4635  Number> &fe_eval,
4636  const Number * in_array,
4637  Number * out_array,
4638  typename std::enable_if<fe_degree != -1>::type * = nullptr)
4639  {
4640  constexpr unsigned int dofs_per_component =
4641  Utilities::pow(fe_degree + 1, dim);
4642 
4643  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
4644  Assert(fe_eval.get_shape_info().element_type <=
4646  ExcNotImplemented());
4647 
4649  dim,
4650  fe_degree + 1,
4651  fe_degree + 1,
4652  Number>
4653  evaluator(
4656  fe_eval.get_shape_info().data.front().inverse_shape_values_eo);
4657 
4658  for (unsigned int d = 0; d < n_components; ++d)
4659  {
4660  const Number *in = in_array + d * dofs_per_component;
4661  Number * out = out_array + d * dofs_per_component;
4662  // Need to select 'apply' method with hessian slot because values
4663  // assume symmetries that do not exist in the inverse shapes
4664  evaluator.template hessians<0, true, false>(in, out);
4665  if (dim > 1)
4666  evaluator.template hessians<1, true, false>(out, out);
4667  if (dim > 2)
4668  evaluator.template hessians<2, true, false>(out, out);
4669  }
4670  for (unsigned int q = 0; q < dofs_per_component; ++q)
4671  {
4672  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
4673  for (unsigned int d = 0; d < n_components; ++d)
4674  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
4675  }
4676  for (unsigned int d = 0; d < n_components; ++d)
4677  {
4678  Number *out = out_array + d * dofs_per_component;
4679  if (dim > 2)
4680  evaluator.template hessians<2, false, false>(out, out);
4681  if (dim > 1)
4682  evaluator.template hessians<1, false, false>(out, out);
4683  evaluator.template hessians<0, false, false>(out, out);
4684  }
4685  return false;
4686  }
4687 
4688  template <int fe_degree, int = 0>
4689  static bool
4690  run(const unsigned int n_components,
4691  const FEEvaluationBaseData<dim,
4692  typename Number::value_type,
4693  false,
4694  Number> &fe_eval,
4695  const Number * in_array,
4696  Number * out_array,
4697  typename std::enable_if<fe_degree == -1>::type * = nullptr)
4698  {
4699  static_assert(fe_degree == -1, "Only usable for degree -1");
4700  const unsigned int dofs_per_component =
4701  fe_eval.get_shape_info().dofs_per_component_on_cell;
4702 
4703  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
4704 
4705  internal::
4706  EvaluatorTensorProduct<internal::evaluate_general, dim, 0, 0, Number>
4707  evaluator(fe_eval.get_shape_info().data.front().inverse_shape_values,
4710  fe_eval.get_shape_info().data.front().fe_degree + 1,
4711  fe_eval.get_shape_info().data.front().fe_degree + 1);
4712 
4713  for (unsigned int d = 0; d < n_components; ++d)
4714  {
4715  const Number *in = in_array + d * dofs_per_component;
4716  Number * out = out_array + d * dofs_per_component;
4717  // Need to select 'apply' method with hessian slot because values
4718  // assume symmetries that do not exist in the inverse shapes
4719  evaluator.template values<0, true, false>(in, out);
4720  if (dim > 1)
4721  evaluator.template values<1, true, false>(out, out);
4722  if (dim > 2)
4723  evaluator.template values<2, true, false>(out, out);
4724  }
4725  for (unsigned int q = 0; q < dofs_per_component; ++q)
4726  {
4727  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
4728  for (unsigned int d = 0; d < n_components; ++d)
4729  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
4730  }
4731  for (unsigned int d = 0; d < n_components; ++d)
4732  {
4733  Number *out = out_array + d * dofs_per_component;
4734  if (dim > 2)
4735  evaluator.template values<2, false, false>(out, out);
4736  if (dim > 1)
4737  evaluator.template values<1, false, false>(out, out);
4738  evaluator.template values<0, false, false>(out, out);
4739  }
4740  return false;
4741  }
4742  };
4743 
4744 
4745 
4750  template <int dim, typename Number>
4752  {
4753  template <int fe_degree, int = 0>
4754  static bool
4755  run(const unsigned int n_desired_components,
4756  const AlignedVector<Number> &inverse_shape,
4757  const AlignedVector<Number> &inverse_coefficients,
4758  const Number * in_array,
4759  Number * out_array,
4760  typename std::enable_if<fe_degree != -1>::type * = nullptr)
4761  {
4762  constexpr unsigned int dofs_per_component =
4763  Utilities::pow(fe_degree + 1, dim);
4764  Assert(inverse_coefficients.size() > 0 &&
4765  inverse_coefficients.size() % dofs_per_component == 0,
4766  ExcMessage(
4767  "Expected diagonal to be a multiple of scalar dof per cells"));
4768  if (inverse_coefficients.size() != dofs_per_component)
4769  AssertDimension(n_desired_components * dofs_per_component,
4770  inverse_coefficients.size());
4771 
4772  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
4773 
4775  dim,
4776  fe_degree + 1,
4777  fe_degree + 1,
4778  Number>
4779  evaluator(AlignedVector<Number>(),
4781  inverse_shape);
4782 
4783  const unsigned int shift_coefficient =
4784  inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
4785  0;
4786  const Number *inv_coefficient = inverse_coefficients.data();
4787  for (unsigned int d = 0; d < n_desired_components; ++d)
4788  {
4789  const Number *in = in_array + d * dofs_per_component;
4790  Number * out = out_array + d * dofs_per_component;
4791  // Need to select 'apply' method with hessian slot because values
4792  // assume symmetries that do not exist in the inverse shapes
4793  evaluator.template hessians<0, true, false>(in, out);
4794  if (dim > 1)
4795  evaluator.template hessians<1, true, false>(out, out);
4796  if (dim > 2)
4797  evaluator.template hessians<2, true, false>(out, out);
4798 
4799  for (unsigned int q = 0; q < dofs_per_component; ++q)
4800  out[q] *= inv_coefficient[q];
4801 
4802  if (dim > 2)
4803  evaluator.template hessians<2, false, false>(out, out);
4804  if (dim > 1)
4805  evaluator.template hessians<1, false, false>(out, out);
4806  evaluator.template hessians<0, false, false>(out, out);
4807 
4808  inv_coefficient += shift_coefficient;
4809  }
4810  return false;
4811  }
4812 
4816  template <int fe_degree, int = 0>
4817  static bool
4818  run(const unsigned int,
4819  const AlignedVector<Number> &,
4820  const AlignedVector<Number> &,
4821  const Number *,
4822  Number *,
4823  typename std::enable_if<fe_degree == -1>::type * = nullptr)
4824  {
4825  static_assert(fe_degree == -1, "Only usable for degree -1");
4826  Assert(false, ExcNotImplemented());
4827  return false;
4828  }
4829  };
4830 
4831 
4832 
4837  template <int dim, typename Number>
4839  {
4840  template <int fe_degree, int n_q_points_1d>
4841  static bool
4842  run(const unsigned int n_desired_components,
4843  const FEEvaluationBaseData<dim,
4844  typename Number::value_type,
4845  false,
4846  Number> &fe_eval,
4847  const Number * in_array,
4848  Number * out_array)
4849  {
4850  static const bool do_inplace =
4851  fe_degree > -1 && (fe_degree + 1 == n_q_points_1d);
4852 
4853  Assert(fe_eval.get_shape_info().element_type !=
4855  ExcNotImplemented());
4856 
4857  const auto &inverse_shape =
4858  do_inplace ?
4859  fe_eval.get_shape_info().data.front().inverse_shape_values_eo :
4860  fe_eval.get_shape_info().data.front().inverse_shape_values;
4861 
4862  const unsigned int dofs_per_component =
4863  do_inplace ? Utilities::pow(fe_degree + 1, dim) :
4864  fe_eval.get_shape_info().dofs_per_component_on_cell;
4865  const unsigned int n_q_points = do_inplace ?
4866  Utilities::pow(fe_degree + 1, dim) :
4867  fe_eval.get_shape_info().n_q_points;
4868 
4871  dim,
4872  fe_degree + 1,
4873  n_q_points_1d,
4874  Number>
4875  evaluator(AlignedVector<Number>(),
4877  inverse_shape,
4878  fe_eval.get_shape_info().data.front().fe_degree + 1,
4879  fe_eval.get_shape_info().data.front().n_q_points_1d);
4880 
4881  for (unsigned int d = 0; d < n_desired_components; ++d)
4882  {
4883  const Number *in = in_array + d * n_q_points;
4884  Number * out = out_array + d * dofs_per_component;
4885 
4886  auto temp_1 = do_inplace ? out : fe_eval.get_scratch_data().begin();
4887  auto temp_2 = do_inplace ?
4888  out :
4889  (temp_1 + std::max(n_q_points, dofs_per_component));
4890 
4891  if (dim == 3)
4892  {
4893  evaluator.template hessians<2, false, false>(in, temp_1);
4894  evaluator.template hessians<1, false, false>(temp_1, temp_2);
4895  evaluator.template hessians<0, false, false>(temp_2, out);
4896  }
4897  if (dim == 2)
4898  {
4899  evaluator.template hessians<1, false, false>(in, temp_1);
4900  evaluator.template hessians<0, false, false>(temp_1, out);
4901  }
4902  if (dim == 1)
4903  evaluator.template hessians<0, false, false>(in, out);
4904  }
4905  return false;
4906  }
4907  };
4908 
4909 } // end of namespace internal
4910 
4911 
4913 
4914 #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)
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 const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
static bool run(const unsigned int n_desired_components, const FEEvaluationBaseData< dim, typename Number::value_type, false, Number > &fe_eval, const Number *in_array, Number *out_array)
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 *hessians_quad, Number *scratch_data, const bool add_into_values_array)
void hermite_grad(const T0 &temp_1, const T0 &temp_2, T1 &dst_ptr_1, T1 &dst_ptr_2, const T2 &grad_weight)
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 *hessians_quad, Number *scratch_data, const bool sum_into_values_array)
void do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
void do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
void hermite_grad(T0 &temp_1, T0 &temp_2, const T1 &src_ptr_1, const T2 &src_ptr_2, const T2 &grad_weight)
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 *hessians_quad, VectorizedArrayType *scratch_data, const bool do_values, const bool do_gradients, const bool do_hessians, 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)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
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()
static void interpolate_quadrature(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< Number > &data, const Number *input, Number *output, const bool do_gradients, const bool do_hessians, const unsigned int face_no)
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 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)
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 bool hessians, const unsigned int n_q_points, Number *tmp_values, Number *values_quad, Number *gradients_quad, Number *hessians_quad)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
std::vector< unsigned int > n_q_points_faces
Definition: shape_info.h:433
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)
static void interpolate(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< Number > &data, const Number *input, Number *output, const bool do_gradients, const bool do_hessians, const unsigned int face_no)
const std::array< unsigned int, VectorizedArrayType::size()> face_orientations
T fixed_power(const T t)
Definition: utilities.h:1082
const std::vector< ArrayView< const Number > > * sm_ptr
const std::array< unsigned int, VectorizedArrayType::size()> cells
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 *hessians_quad, Number *scratch_data, const bool integrate_val, const bool integrate_grad, const bool integrate_hessian, const unsigned int subface_index)
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:1461
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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 *hessians_quad, Number *scratch_data, const bool add_into_values_array)
static bool run(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, const VectorizedArrayType *values_array, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *hessians_quad, VectorizedArrayType *scratch_data, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians, 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:401
void load(const Number *ptr)
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:98
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:50
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)
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
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:400
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
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 *hessians_quad, VectorizedArrayType *scratch_data, const bool integrate_values, const bool integrate_gradients, const bool integrate_hessians, 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 void interpolate_generic(const unsigned int n_components, const Number *input, Number *output, const bool do_gradients, const bool do_hessians, 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)
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 *hessians_quad, Number *scratch_data, const bool evaluate_val, const bool evaluate_grad, const bool evaluate_hessian, const unsigned int subface_index)
const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & data
static ::ExceptionBase & ExcNotImplemented()
std::vector< UnivariateShapeData< Number > > data
Definition: shape_info.h:392
const std::array< unsigned int, VectorizedArrayType::size()> cells
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 *hessians_quad, VectorizedArrayType *scratch_data, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians, 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 value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
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 *hessians_quad, VectorizedArrayType *scratch_data, const bool do_values, const bool do_gradients, const bool do_hessians, 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)
static bool run(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, VectorizedArrayType *values_array, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *hessians_quad, VectorizedArrayType *scratch_data, const bool integrate_values, const bool integrate_gradients, const bool integrate_hessians, const unsigned int face_no, const unsigned int subface_index, const unsigned int face_orientation, const Table< 2, unsigned int > &orientation_map)
EvaluationFlags
The EvaluationFlags enum.
const std::array< unsigned int, VectorizedArrayType::size()> face_nos
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:138
const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index
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)