Reference documentation for deal.II version GIT 0980a66d4b 2023-03-23 20:20:03+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
evaluation_kernels.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_evaluation_kernels_h
18 #define dealii_matrix_free_evaluation_kernels_h
19 
20 #include <deal.II/base/config.h>
21 
23 #include <deal.II/base/ndarray.h>
24 #include <deal.II/base/utilities.h>
26 
32 
33 
35 
36 
37 namespace internal
38 {
39  // Select evaluator type from element shape function type
40  template <MatrixFreeFunctions::ElementType element, bool is_long>
42  {};
43 
44  template <bool is_long>
45  struct EvaluatorSelector<MatrixFreeFunctions::tensor_general, is_long>
46  {
47  static const EvaluatorVariant variant = evaluate_general;
48  };
49 
50  template <>
51  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, false>
52  {
53  static const EvaluatorVariant variant = evaluate_symmetric;
54  };
55 
56  template <>
57  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, true>
58  {
59  static const EvaluatorVariant variant = evaluate_evenodd;
60  };
61 
62  template <bool is_long>
63  struct EvaluatorSelector<MatrixFreeFunctions::truncated_tensor, is_long>
64  {
65  static const EvaluatorVariant variant = evaluate_general;
66  };
67 
68  template <>
69  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,
70  false>
71  {
72  static const EvaluatorVariant variant = evaluate_general;
73  };
74 
75  template <>
76  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0, true>
77  {
78  static const EvaluatorVariant variant = evaluate_evenodd;
79  };
80 
81  template <bool is_long>
82  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_collocation,
83  is_long>
84  {
85  static const EvaluatorVariant variant = evaluate_evenodd;
86  };
87 
88  template <bool is_long>
89  struct EvaluatorSelector<MatrixFreeFunctions::tensor_raviart_thomas, is_long>
90  {
91  static const EvaluatorVariant variant = evaluate_raviart_thomas;
92  };
93 
94 
95 
112  template <MatrixFreeFunctions::ElementType type,
113  int dim,
114  int fe_degree,
115  int n_q_points_1d,
116  typename Number>
118  {
119  static const EvaluatorVariant variant =
120  EvaluatorSelector<type, (fe_degree + n_q_points_1d > 4)>::variant;
121 
123  dim,
124  fe_degree + 1,
125  n_q_points_1d,
126  Number>;
127 
128  static void
129  evaluate(const unsigned int n_components,
130  const EvaluationFlags::EvaluationFlags evaluation_flag,
131  const Number * values_dofs_actual,
133 
134  static void
135  integrate(const unsigned int n_components,
136  const EvaluationFlags::EvaluationFlags integration_flag,
137  Number * values_dofs_actual,
139  const bool add_into_values_array);
140 
141  static Eval
144  *univariate_shape_data)
145  {
146  if (variant == evaluate_evenodd)
147  return Eval(univariate_shape_data->shape_values_eo,
148  univariate_shape_data->shape_gradients_eo,
149  univariate_shape_data->shape_hessians_eo,
150  univariate_shape_data->fe_degree + 1,
151  univariate_shape_data->n_q_points_1d);
152  else
153  return Eval(univariate_shape_data->shape_values,
154  univariate_shape_data->shape_gradients,
155  univariate_shape_data->shape_hessians,
156  univariate_shape_data->fe_degree + 1,
157  univariate_shape_data->n_q_points_1d);
158  }
159  };
160 
161 
162 
167  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
168  struct FEEvaluationImpl<MatrixFreeFunctions::tensor_none,
169  dim,
170  fe_degree,
171  n_q_points_1d,
172  Number>
173  {
174  static void
175  evaluate(const unsigned int n_components,
176  const EvaluationFlags::EvaluationFlags evaluation_flag,
177  const Number * values_dofs_actual,
179 
180  static void
181  integrate(const unsigned int n_components,
182  const EvaluationFlags::EvaluationFlags integration_flag,
183  Number * values_dofs_actual,
185  const bool add_into_values_array);
186  };
187 
192  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
193  struct FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
194  dim,
195  fe_degree,
196  n_q_points_1d,
197  Number>
198  {
199  template <bool integrate>
200  static void
201  evaluate_or_integrate(
202  const EvaluationFlags::EvaluationFlags evaluation_flag,
203  Number * values_dofs_actual,
205  const bool add_into_values_array = false);
206 
207  private:
208  template <typename EvalType>
209  static EvalType
212  {
213  return EvalType(shape_data.shape_values,
214  shape_data.shape_gradients,
215  shape_data.shape_hessians);
216  }
217 
218  template <int normal_dir>
219  static void
220  evaluate_tensor_product_per_component(
221  const EvaluationFlags::EvaluationFlags evaluation_flag,
222  Number * values_dofs_actual,
224  const bool add_into_values_array,
225  std::integral_constant<bool, false>);
226 
227  template <int normal_dir>
228  static void
229  evaluate_tensor_product_per_component(
230  const EvaluationFlags::EvaluationFlags evaluation_flag,
231  Number * values_dofs_actual,
233  const bool add_into_values_array,
234  std::integral_constant<bool, true>);
235  };
236 
237 
238 
239  template <MatrixFreeFunctions::ElementType type,
240  int dim,
241  int fe_degree,
242  int n_q_points_1d,
243  typename Number>
244  inline void
246  const unsigned int n_components,
247  const EvaluationFlags::EvaluationFlags evaluation_flag,
248  const Number * values_dofs_actual,
250  {
251  if (evaluation_flag == EvaluationFlags::nothing)
252  return;
253 
254  std::array<const MatrixFreeFunctions::UnivariateShapeData<Number> *, 3>
255  univariate_shape_data;
256 
257  const auto &shape_data = fe_eval.get_shape_info().data;
258 
259  univariate_shape_data.fill(&shape_data.front());
260 
261  if (shape_data.size() == dim)
262  for (int i = 1; i < dim; ++i)
263  univariate_shape_data[i] = &shape_data[i];
264 
265  Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
266  Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
267  Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);
268 
269  const unsigned int temp_size =
270  Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
271  0 :
272  (Eval::n_rows_of_product > Eval::n_columns_of_product ?
273  Eval::n_rows_of_product :
274  Eval::n_columns_of_product);
275  Number *temp1 = fe_eval.get_scratch_data().begin();
276  Number *temp2;
277  if (temp_size == 0)
278  {
279  temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
280  shape_data.front().fe_degree + 1),
281  Utilities::fixed_power<dim>(
282  shape_data.front().n_q_points_1d));
283  }
284  else
285  {
286  temp2 = temp1 + temp_size;
287  }
288 
289  const std::size_t n_q_points = temp_size == 0 ?
290  fe_eval.get_shape_info().n_q_points :
291  Eval::n_columns_of_product;
292  const std::size_t dofs_per_comp =
294  Utilities::pow(shape_data.front().fe_degree + 1, dim) :
296  const Number *values_dofs = values_dofs_actual;
298  {
299  const std::size_t n_dofs_per_comp =
301  Number *values_dofs_tmp =
302  temp1 + 2 * (std::max(n_dofs_per_comp, n_q_points));
303  const int degree =
304  fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
305  for (unsigned int c = 0; c < n_components; ++c)
306  for (int i = 0, count_p = 0, count_q = 0;
307  i < (dim > 2 ? degree + 1 : 1);
308  ++i)
309  {
310  for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
311  {
312  for (int k = 0; k < degree + 1 - j - i;
313  ++k, ++count_p, ++count_q)
314  values_dofs_tmp[c * dofs_per_comp + count_q] =
315  values_dofs_actual[c * n_dofs_per_comp + count_p];
316  for (int k = degree + 1 - j - i; k < degree + 1;
317  ++k, ++count_q)
318  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
319  }
320  for (int j = degree + 1 - i; j < degree + 1; ++j)
321  for (int k = 0; k < degree + 1; ++k, ++count_q)
322  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
323  }
324  values_dofs = values_dofs_tmp;
325  }
326 
327  Number *values_quad = fe_eval.begin_values();
328  Number *gradients_quad = fe_eval.begin_gradients();
329  Number *hessians_quad = fe_eval.begin_hessians();
330 
331  switch (dim)
332  {
333  case 1:
334  for (unsigned int c = 0; c < n_components; ++c)
335  {
336  if (evaluation_flag & EvaluationFlags::values)
337  eval0.template values<0, true, false>(values_dofs, values_quad);
338  if (evaluation_flag & EvaluationFlags::gradients)
339  eval0.template gradients<0, true, false>(values_dofs,
340  gradients_quad);
341  if (evaluation_flag & EvaluationFlags::hessians)
342  eval0.template hessians<0, true, false>(values_dofs,
343  hessians_quad);
344 
345  // advance the next component in 1d array
346  values_dofs += dofs_per_comp;
347  values_quad += n_q_points;
348  gradients_quad += n_q_points;
349  hessians_quad += n_q_points;
350  }
351  break;
352 
353  case 2:
354  for (unsigned int c = 0; c < n_components; ++c)
355  {
356  // grad x
357  if (evaluation_flag & EvaluationFlags::gradients)
358  {
359  eval0.template gradients<0, true, false>(values_dofs, temp1);
360  eval1.template values<1, true, false>(temp1, gradients_quad);
361  }
362  if (evaluation_flag & EvaluationFlags::hessians)
363  {
364  // grad xy
365  if (!(evaluation_flag & EvaluationFlags::gradients))
366  eval0.template gradients<0, true, false>(values_dofs,
367  temp1);
368  eval1.template gradients<1, true, false>(temp1,
369  hessians_quad +
370  2 * n_q_points);
371 
372  // grad xx
373  eval0.template hessians<0, true, false>(values_dofs, temp1);
374  eval1.template values<1, true, false>(temp1, hessians_quad);
375  }
376 
377  // grad y
378  eval0.template values<0, true, false>(values_dofs, temp1);
379  if (evaluation_flag & EvaluationFlags::gradients)
380  eval1.template gradients<1, true, false>(temp1,
381  gradients_quad +
382  n_q_points);
383 
384  // grad yy
385  if (evaluation_flag & EvaluationFlags::hessians)
386  eval1.template hessians<1, true, false>(temp1,
387  hessians_quad +
388  n_q_points);
389 
390  // val: can use values applied in x
391  if (evaluation_flag & EvaluationFlags::values)
392  eval1.template values<1, true, false>(temp1, values_quad);
393 
394  // advance to the next component in 1d array
395  values_dofs += dofs_per_comp;
396  values_quad += n_q_points;
397  gradients_quad += 2 * n_q_points;
398  hessians_quad += 3 * n_q_points;
399  }
400  break;
401 
402  case 3:
403  for (unsigned int c = 0; c < n_components; ++c)
404  {
405  if (evaluation_flag & EvaluationFlags::gradients)
406  {
407  // grad x
408  eval0.template gradients<0, true, false>(values_dofs, temp1);
409  eval1.template values<1, true, false>(temp1, temp2);
410  eval2.template values<2, true, false>(temp2, gradients_quad);
411  }
412 
413  if (evaluation_flag & EvaluationFlags::hessians)
414  {
415  // grad xz
416  if (!(evaluation_flag & EvaluationFlags::gradients))
417  {
418  eval0.template gradients<0, true, false>(values_dofs,
419  temp1);
420  eval1.template values<1, true, false>(temp1, temp2);
421  }
422  eval2.template gradients<2, true, false>(temp2,
423  hessians_quad +
424  4 * n_q_points);
425 
426  // grad xy
427  eval1.template gradients<1, true, false>(temp1, temp2);
428  eval2.template values<2, true, false>(temp2,
429  hessians_quad +
430  3 * n_q_points);
431 
432  // grad xx
433  eval0.template hessians<0, true, false>(values_dofs, temp1);
434  eval1.template values<1, true, false>(temp1, temp2);
435  eval2.template values<2, true, false>(temp2, hessians_quad);
436  }
437 
438  // grad y
439  eval0.template values<0, true, false>(values_dofs, temp1);
440  if (evaluation_flag & EvaluationFlags::gradients)
441  {
442  eval1.template gradients<1, true, false>(temp1, temp2);
443  eval2.template values<2, true, false>(temp2,
444  gradients_quad +
445  n_q_points);
446  }
447 
448  if (evaluation_flag & EvaluationFlags::hessians)
449  {
450  // grad yz
451  if (!(evaluation_flag & EvaluationFlags::gradients))
452  eval1.template gradients<1, true, false>(temp1, temp2);
453  eval2.template gradients<2, true, false>(temp2,
454  hessians_quad +
455  5 * n_q_points);
456 
457  // grad yy
458  eval1.template hessians<1, true, false>(temp1, temp2);
459  eval2.template values<2, true, false>(temp2,
460  hessians_quad +
461  n_q_points);
462  }
463 
464  // grad z: can use the values applied in x direction stored in
465  // temp1
466  eval1.template values<1, true, false>(temp1, temp2);
467  if (evaluation_flag & EvaluationFlags::gradients)
468  eval2.template gradients<2, true, false>(temp2,
469  gradients_quad +
470  2 * n_q_points);
471 
472  // grad zz: can use the values applied in x and y direction stored
473  // in temp2
474  if (evaluation_flag & EvaluationFlags::hessians)
475  eval2.template hessians<2, true, false>(temp2,
476  hessians_quad +
477  2 * n_q_points);
478 
479  // val: can use the values applied in x & y direction stored in
480  // temp2
481  if (evaluation_flag & EvaluationFlags::values)
482  eval2.template values<2, true, false>(temp2, values_quad);
483 
484  // advance to the next component in 1d array
485  values_dofs += dofs_per_comp;
486  values_quad += n_q_points;
487  gradients_quad += 3 * n_q_points;
488  hessians_quad += 6 * n_q_points;
489  }
490  break;
491 
492  default:
493  AssertThrow(false, ExcNotImplemented());
494  }
495 
496  // case additional dof for FE_Q_DG0: add values; gradients and second
497  // derivatives evaluate to zero
499  (evaluation_flag & EvaluationFlags::values))
500  {
501  values_quad -= n_components * n_q_points;
502  values_dofs -= n_components * dofs_per_comp;
503  for (std::size_t c = 0; c < n_components; ++c)
504  for (std::size_t q = 0; q < n_q_points; ++q)
505  values_quad[c * n_q_points + q] +=
506  values_dofs[(c + 1) * dofs_per_comp - 1];
507  }
508  }
509 
510 
511 
512  template <MatrixFreeFunctions::ElementType type,
513  int dim,
514  int fe_degree,
515  int n_q_points_1d,
516  typename Number>
517  inline void
519  const unsigned int n_components,
520  const EvaluationFlags::EvaluationFlags integration_flag,
521  Number * values_dofs_actual,
523  const bool add_into_values_array)
524  {
525  std::array<const MatrixFreeFunctions::UnivariateShapeData<Number> *, 3>
526  univariate_shape_data;
527 
528  const auto &shape_data = fe_eval.get_shape_info().data;
529  univariate_shape_data.fill(&shape_data.front());
530 
531  if (shape_data.size() == dim)
532  for (int i = 1; i < dim; ++i)
533  univariate_shape_data[i] = &shape_data[i];
534 
535  Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
536  Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
537  Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);
538 
539  const unsigned int temp_size =
540  Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
541  0 :
542  (Eval::n_rows_of_product > Eval::n_columns_of_product ?
543  Eval::n_rows_of_product :
544  Eval::n_columns_of_product);
545  Number *temp1 = fe_eval.get_scratch_data().begin();
546  Number *temp2;
547  if (temp_size == 0)
548  {
549  temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
550  shape_data.front().fe_degree + 1),
551  Utilities::fixed_power<dim>(
552  shape_data.front().n_q_points_1d));
553  }
554  else
555  {
556  temp2 = temp1 + temp_size;
557  }
558 
559  const std::size_t n_q_points = temp_size == 0 ?
560  fe_eval.get_shape_info().n_q_points :
561  Eval::n_columns_of_product;
562  const unsigned int dofs_per_comp =
564  Utilities::fixed_power<dim>(shape_data.front().fe_degree + 1) :
566  // expand dof_values to tensor product for truncated tensor products
567  Number *values_dofs =
569  temp1 + 2 * (std::max<std::size_t>(
571  n_q_points)) :
572  values_dofs_actual;
573 
574  Number *values_quad = fe_eval.begin_values();
575  Number *gradients_quad = fe_eval.begin_gradients();
576  Number *hessians_quad = fe_eval.begin_hessians();
577 
578  switch (dim)
579  {
580  case 1:
581  for (unsigned int c = 0; c < n_components; ++c)
582  {
583  if (integration_flag & EvaluationFlags::values)
584  {
585  if (add_into_values_array == false)
586  eval0.template values<0, false, false>(values_quad,
587  values_dofs);
588  else
589  eval0.template values<0, false, true>(values_quad,
590  values_dofs);
591  }
592  if (integration_flag & EvaluationFlags::gradients)
593  {
594  if (integration_flag & EvaluationFlags::values ||
595  add_into_values_array == true)
596  eval0.template gradients<0, false, true>(gradients_quad,
597  values_dofs);
598  else
599  eval0.template gradients<0, false, false>(gradients_quad,
600  values_dofs);
601  }
602  if ((integration_flag & EvaluationFlags::hessians) != 0u)
603  {
604  if ((integration_flag & EvaluationFlags::values) != 0u ||
605  (integration_flag & EvaluationFlags::gradients) != 0u ||
606  add_into_values_array == true)
607  eval0.template hessians<0, false, true>(hessians_quad,
608  values_dofs);
609  else
610  eval0.template hessians<0, false, false>(hessians_quad,
611  values_dofs);
612  }
613 
614  // advance to the next component in 1d array
615  values_dofs += dofs_per_comp;
616  values_quad += n_q_points;
617  gradients_quad += n_q_points;
618  hessians_quad += n_q_points;
619  }
620  break;
621 
622  case 2:
623  for (unsigned int c = 0; c < n_components; ++c)
624  {
625  if ((integration_flag & EvaluationFlags::values) &&
626  !(integration_flag & EvaluationFlags::gradients))
627  {
628  eval1.template values<1, false, false>(values_quad, temp1);
629  if (add_into_values_array == false)
630  eval0.template values<0, false, false>(temp1, values_dofs);
631  else
632  eval0.template values<0, false, true>(temp1, values_dofs);
633  }
634  if (integration_flag & EvaluationFlags::gradients)
635  {
636  eval1.template gradients<1, false, false>(gradients_quad +
637  n_q_points,
638  temp1);
639  if (integration_flag & EvaluationFlags::values)
640  eval1.template values<1, false, true>(values_quad, temp1);
641  if (add_into_values_array == false)
642  eval0.template values<0, false, false>(temp1, values_dofs);
643  else
644  eval0.template values<0, false, true>(temp1, values_dofs);
645  eval1.template values<1, false, false>(gradients_quad, temp1);
646  eval0.template gradients<0, false, true>(temp1, values_dofs);
647  }
648  if ((integration_flag & EvaluationFlags::hessians) != 0u)
649  {
650  // grad xx
651  eval1.template values<1, false, false>(hessians_quad, temp1);
652 
653  if ((integration_flag & EvaluationFlags::values) != 0u ||
654  (integration_flag & EvaluationFlags::gradients) != 0u ||
655  add_into_values_array == true)
656  eval0.template hessians<0, false, true>(temp1, values_dofs);
657  else
658  eval0.template hessians<0, false, false>(temp1,
659  values_dofs);
660 
661  // grad yy
662  eval1.template hessians<1, false, false>(hessians_quad +
663  n_q_points,
664  temp1);
665  eval0.template values<0, false, true>(temp1, values_dofs);
666 
667  // grad xy
668  eval1.template gradients<1, false, false>(hessians_quad +
669  2 * n_q_points,
670  temp1);
671  eval0.template gradients<0, false, true>(temp1, values_dofs);
672  }
673 
674  // advance to the next component in 1d array
675  values_dofs += dofs_per_comp;
676  values_quad += n_q_points;
677  gradients_quad += 2 * n_q_points;
678  hessians_quad += 3 * n_q_points;
679  }
680  break;
681 
682  case 3:
683  for (unsigned int c = 0; c < n_components; ++c)
684  {
685  if ((integration_flag & EvaluationFlags::values) &&
686  !(integration_flag & EvaluationFlags::gradients))
687  {
688  eval2.template values<2, false, false>(values_quad, temp1);
689  eval1.template values<1, false, false>(temp1, temp2);
690  if (add_into_values_array == false)
691  eval0.template values<0, false, false>(temp2, values_dofs);
692  else
693  eval0.template values<0, false, true>(temp2, values_dofs);
694  }
695  if (integration_flag & EvaluationFlags::gradients)
696  {
697  eval2.template gradients<2, false, false>(gradients_quad +
698  2 * n_q_points,
699  temp1);
700  if (integration_flag & EvaluationFlags::values)
701  eval2.template values<2, false, true>(values_quad, temp1);
702  eval1.template values<1, false, false>(temp1, temp2);
703  eval2.template values<2, false, false>(gradients_quad +
704  n_q_points,
705  temp1);
706  eval1.template gradients<1, false, true>(temp1, temp2);
707  if (add_into_values_array == false)
708  eval0.template values<0, false, false>(temp2, values_dofs);
709  else
710  eval0.template values<0, false, true>(temp2, values_dofs);
711  eval2.template values<2, false, false>(gradients_quad, temp1);
712  eval1.template values<1, false, false>(temp1, temp2);
713  eval0.template gradients<0, false, true>(temp2, values_dofs);
714  }
715  if ((integration_flag & EvaluationFlags::hessians) != 0u)
716  {
717  // grad xx
718  eval2.template values<2, false, false>(hessians_quad, temp1);
719  eval1.template values<1, false, false>(temp1, temp2);
720 
721  if ((integration_flag & EvaluationFlags::values) != 0u ||
722  (integration_flag & EvaluationFlags::gradients) != 0u ||
723  add_into_values_array == true)
724  eval0.template hessians<0, false, true>(temp2, values_dofs);
725  else
726  eval0.template hessians<0, false, false>(temp2,
727  values_dofs);
728 
729  // grad yy
730  eval2.template values<2, false, false>(hessians_quad +
731  n_q_points,
732  temp1);
733  eval1.template hessians<1, false, false>(temp1, temp2);
734  eval0.template values<0, false, true>(temp2, values_dofs);
735 
736  // grad zz
737  eval2.template hessians<2, false, false>(hessians_quad +
738  2 * n_q_points,
739  temp1);
740  eval1.template values<1, false, false>(temp1, temp2);
741  eval0.template values<0, false, true>(temp2, values_dofs);
742 
743  // grad xy
744  eval2.template values<2, false, false>(hessians_quad +
745  3 * n_q_points,
746  temp1);
747  eval1.template gradients<1, false, false>(temp1, temp2);
748  eval0.template gradients<0, false, true>(temp2, values_dofs);
749 
750  // grad xz
751  eval2.template gradients<2, false, false>(hessians_quad +
752  4 * n_q_points,
753  temp1);
754  eval1.template values<1, false, false>(temp1, temp2);
755  eval0.template gradients<0, false, true>(temp2, values_dofs);
756 
757  // grad yz
758  eval2.template gradients<2, false, false>(hessians_quad +
759  5 * n_q_points,
760  temp1);
761  eval1.template gradients<1, false, false>(temp1, temp2);
762  eval0.template values<0, false, true>(temp2, values_dofs);
763  }
764 
765  // advance to the next component in 1d array
766  values_dofs += dofs_per_comp;
767  values_quad += n_q_points;
768  gradients_quad += 3 * n_q_points;
769  hessians_quad += 6 * n_q_points;
770  }
771  break;
772 
773  default:
774  AssertThrow(false, ExcNotImplemented());
775  }
776 
777  // case FE_Q_DG0: add values, gradients and second derivatives are zero
779  {
780  values_dofs -= n_components * dofs_per_comp - dofs_per_comp + 1;
781  values_quad -= n_components * n_q_points;
782  if (integration_flag & EvaluationFlags::values)
783  for (unsigned int c = 0; c < n_components; ++c)
784  {
785  values_dofs[0] = values_quad[0];
786  for (unsigned int q = 1; q < n_q_points; ++q)
787  values_dofs[0] += values_quad[q];
788  values_dofs += dofs_per_comp;
789  values_quad += n_q_points;
790  }
791  else
792  {
793  for (unsigned int c = 0; c < n_components; ++c)
794  values_dofs[c * dofs_per_comp] = Number();
795  values_dofs += n_components * dofs_per_comp;
796  }
797  }
798 
800  {
801  const std::size_t n_dofs_per_comp =
803  values_dofs -= dofs_per_comp * n_components;
804  const int degree =
805  fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
806  for (unsigned int c = 0; c < n_components; ++c)
807  for (int i = 0, count_p = 0, count_q = 0;
808  i < (dim > 2 ? degree + 1 : 1);
809  ++i)
810  {
811  for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
812  {
813  for (int k = 0; k < degree + 1 - j - i;
814  ++k, ++count_p, ++count_q)
815  values_dofs_actual[c * n_dofs_per_comp + count_p] =
816  values_dofs[c * dofs_per_comp + count_q];
817  count_q += j + i;
818  }
819  count_q += i * (degree + 1);
820  }
821  }
822  }
823 
824 
825 
826  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
827  inline void
830  dim,
831  fe_degree,
832  n_q_points_1d,
833  Number>::evaluate(const unsigned int n_components,
834  const EvaluationFlags::EvaluationFlags evaluation_flag,
835  const Number * values_dofs_actual,
837  {
838  Assert(!(evaluation_flag & EvaluationFlags::hessians), ExcNotImplemented());
839 
840  const std::size_t n_dofs =
842  const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
843 
844  const auto &shape_data = fe_eval.get_shape_info().data;
845 
846  using Eval =
848 
849  if (evaluation_flag & EvaluationFlags::values)
850  {
851  const auto shape_values = shape_data.front().shape_values.data();
852  auto values_quad_ptr = fe_eval.begin_values();
853  auto values_dofs_actual_ptr = values_dofs_actual;
854 
855  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
856  for (unsigned int c = 0; c < n_components; ++c)
857  {
858  eval.template values<0, true, false>(values_dofs_actual_ptr,
859  values_quad_ptr);
860 
861  values_quad_ptr += n_q_points;
862  values_dofs_actual_ptr += n_dofs;
863  }
864  }
865 
866  if (evaluation_flag & EvaluationFlags::gradients)
867  {
868  const auto shape_gradients = shape_data.front().shape_gradients.data();
869  auto gradients_quad_ptr = fe_eval.begin_gradients();
870  auto values_dofs_actual_ptr = values_dofs_actual;
871 
872  for (unsigned int c = 0; c < n_components; ++c)
873  {
874  for (unsigned int d = 0; d < dim; ++d)
875  {
876  Eval eval(nullptr,
877  shape_gradients + n_q_points * n_dofs * d,
878  nullptr,
879  n_dofs,
880  n_q_points);
881 
882  eval.template gradients<0, true, false>(values_dofs_actual_ptr,
883  gradients_quad_ptr);
884 
885  gradients_quad_ptr += n_q_points;
886  }
887  values_dofs_actual_ptr += n_dofs;
888  }
889  }
890  }
891 
892 
893 
894  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
895  inline void
898  dim,
899  fe_degree,
900  n_q_points_1d,
901  Number>::integrate(const unsigned int n_components,
902  const EvaluationFlags::EvaluationFlags integration_flag,
903  Number * values_dofs_actual,
905  const bool add_into_values_array)
906  {
907  Assert(!(integration_flag & EvaluationFlags::hessians),
909 
910  const std::size_t n_dofs =
912  const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
913 
914  const auto &shape_data = fe_eval.get_shape_info().data;
915 
916  using Eval =
918 
919  if (integration_flag & EvaluationFlags::values)
920  {
921  const auto shape_values = shape_data.front().shape_values.data();
922  auto values_quad_ptr = fe_eval.begin_values();
923  auto values_dofs_actual_ptr = values_dofs_actual;
924 
925  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
926  for (unsigned int c = 0; c < n_components; ++c)
927  {
928  if (add_into_values_array == false)
929  eval.template values<0, false, false>(values_quad_ptr,
930  values_dofs_actual_ptr);
931  else
932  eval.template values<0, false, true>(values_quad_ptr,
933  values_dofs_actual_ptr);
934 
935  values_quad_ptr += n_q_points;
936  values_dofs_actual_ptr += n_dofs;
937  }
938  }
939 
940  if (integration_flag & EvaluationFlags::gradients)
941  {
942  const auto shape_gradients = shape_data.front().shape_gradients.data();
943  auto gradients_quad_ptr = fe_eval.begin_gradients();
944  auto values_dofs_actual_ptr = values_dofs_actual;
945 
946  for (unsigned int c = 0; c < n_components; ++c)
947  {
948  for (unsigned int d = 0; d < dim; ++d)
949  {
950  Eval eval(nullptr,
951  shape_gradients + n_q_points * n_dofs * d,
952  nullptr,
953  n_dofs,
954  n_q_points);
955 
956  if ((add_into_values_array == false &&
957  !(integration_flag & EvaluationFlags::values)) &&
958  d == 0)
959  eval.template gradients<0, false, false>(
960  gradients_quad_ptr, values_dofs_actual_ptr);
961  else
962  eval.template gradients<0, false, true>(
963  gradients_quad_ptr, values_dofs_actual_ptr);
964 
965  gradients_quad_ptr += n_q_points;
966  }
967  values_dofs_actual_ptr += n_dofs;
968  }
969  }
970  }
971 
972 
973  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
974  template <bool integrate>
975  inline void
977  dim,
978  fe_degree,
979  n_q_points_1d,
980  Number>::
981  evaluate_or_integrate(
982  const EvaluationFlags::EvaluationFlags evaluation_flag,
983  Number * values_dofs_actual,
985  const bool add_into_values_array)
986  {
987  if (evaluation_flag == EvaluationFlags::nothing)
988  return;
989 
990  AssertDimension(fe_eval.get_shape_info().data.size(), 2);
991  // First component:
992  evaluate_tensor_product_per_component<0>(
993  evaluation_flag,
994  values_dofs_actual,
995  fe_eval,
996  add_into_values_array,
997  std::integral_constant<bool, integrate>());
998  // Second component :
999  evaluate_tensor_product_per_component<1>(
1000  evaluation_flag,
1001  values_dofs_actual,
1002  fe_eval,
1003  add_into_values_array,
1004  std::integral_constant<bool, integrate>());
1005  if (dim == 3)
1006  {
1007  // Third component
1008  evaluate_tensor_product_per_component<2>(
1009  evaluation_flag,
1010  values_dofs_actual,
1011  fe_eval,
1012  add_into_values_array,
1013  std::integral_constant<bool, integrate>());
1014  }
1015  }
1016 
1017  // Helper function that applies the 1d evaluation kernels.
1018  // std::integral_constant<bool, false> is the interpolation path, and
1019  // std::integral_constant<bool, true> below is the integration path.
1020  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1021  template <int normal_dir>
1022  inline void
1024  dim,
1025  fe_degree,
1026  n_q_points_1d,
1027  Number>::
1028  evaluate_tensor_product_per_component(
1029  const EvaluationFlags::EvaluationFlags evaluation_flag,
1030  Number * values_dofs_actual,
1032  const bool add_into_values_array,
1033  std::integral_constant<bool, false>)
1034  {
1035  (void)add_into_values_array;
1036 
1037  using EvalNormal =
1039  dim,
1040  (fe_degree == -1) ? 1 : fe_degree + 1,
1041  n_q_points_1d,
1042  Number,
1043  normal_dir>;
1044 
1045  using EvalTangent =
1047  dim,
1048  (fe_degree == -1) ? 1 : fe_degree,
1049  n_q_points_1d,
1050  Number,
1051  normal_dir>;
1052  using Eval0 =
1053  typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
1054  using Eval1 =
1055  typename std::conditional<normal_dir == 1, EvalNormal, EvalTangent>::type;
1056  using Eval2 =
1057  typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
1058 
1059  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info =
1060  fe_eval.get_shape_info();
1061  Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
1062  ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]));
1063  Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
1064  ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]));
1065  Eval2 eval2 = create_evaluator_tensor_product<Eval2>(
1066  ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1]));
1067 
1068  Number *temp1 = fe_eval.get_scratch_data().begin();
1069  Number *temp2;
1070 
1071  temp2 =
1072  temp1 +
1073  std::max(Utilities::fixed_power<dim>(shape_info.data[0].fe_degree + 1),
1074  Utilities::fixed_power<dim>(shape_info.data[0].n_q_points_1d));
1075 
1076  const std::size_t n_q_points = shape_info.n_q_points;
1077  const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell;
1078 
1079  // Initial shift depending on component (normal_dir)
1080  Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
1081  Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
1082  Number *gradients_quad =
1083  fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
1084  Number *hessians_quad =
1085  (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
1086  fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
1087 
1088  switch (dim)
1089  {
1090  case 2:
1091  if (evaluation_flag & EvaluationFlags::gradients)
1092  {
1093  eval0.template gradients<0, true, false>(values_dofs, temp1);
1094  eval1.template values<1, true, false>(temp1, gradients_quad);
1095  }
1096  if (evaluation_flag & EvaluationFlags::hessians)
1097  {
1098  // The evaluation/integration here *should* work, however
1099  // the piola transform is not implemented.
1100  AssertThrow(false, ExcNotImplemented());
1101  // grad xy
1102  if (!(evaluation_flag & EvaluationFlags::gradients))
1103  eval0.template gradients<0, true, false>(values_dofs, temp1);
1104  eval1.template gradients<1, true, false>(temp1,
1105  hessians_quad +
1106  2 * n_q_points);
1107 
1108  // grad xx
1109  eval0.template hessians<0, true, false>(values_dofs, temp1);
1110  eval1.template values<1, true, false>(temp1, hessians_quad);
1111  }
1112 
1113  // grad y
1114  eval0.template values<0, true, false>(values_dofs, temp1);
1115  if (evaluation_flag & EvaluationFlags::gradients)
1116  eval1.template gradients<1, true, false>(temp1,
1117  gradients_quad +
1118  n_q_points);
1119 
1120  // grad yy
1121  if (evaluation_flag & EvaluationFlags::hessians)
1122  eval1.template hessians<1, true, false>(temp1,
1123  hessians_quad + n_q_points);
1124 
1125  // val: can use values applied in x
1126  if (evaluation_flag & EvaluationFlags::values)
1127  eval1.template values<1, true, false>(temp1, values_quad);
1128  break;
1129  case 3:
1130  if (evaluation_flag & EvaluationFlags::gradients)
1131  {
1132  // grad x
1133  eval0.template gradients<0, true, false>(values_dofs, temp1);
1134  eval1.template values<1, true, false>(temp1, temp2);
1135  eval2.template values<2, true, false>(temp2, gradients_quad);
1136  }
1137 
1138  if (evaluation_flag & EvaluationFlags::hessians)
1139  {
1140  // The evaluation/integration here *should* work, however
1141  // the piola transform is not implemented.
1142  AssertThrow(false, ExcNotImplemented());
1143  // grad xz
1144  if (!(evaluation_flag & EvaluationFlags::gradients))
1145  {
1146  eval0.template gradients<0, true, false>(values_dofs, temp1);
1147  eval1.template values<1, true, false>(temp1, temp2);
1148  }
1149  eval2.template gradients<2, true, false>(temp2,
1150  hessians_quad +
1151  4 * n_q_points);
1152 
1153  // grad xy
1154  eval1.template gradients<1, true, false>(temp1, temp2);
1155  eval2.template values<2, true, false>(temp2,
1156  hessians_quad +
1157  3 * n_q_points);
1158 
1159  // grad xx
1160  eval0.template hessians<0, true, false>(values_dofs, temp1);
1161  eval1.template values<1, true, false>(temp1, temp2);
1162  eval2.template values<2, true, false>(temp2, hessians_quad);
1163  }
1164 
1165  // grad y
1166  eval0.template values<0, true, false>(values_dofs, temp1);
1167  if (evaluation_flag & EvaluationFlags::gradients)
1168  {
1169  eval1.template gradients<1, true, false>(temp1, temp2);
1170  eval2.template values<2, true, false>(temp2,
1171  gradients_quad +
1172  n_q_points);
1173  }
1174 
1175  if (evaluation_flag & EvaluationFlags::hessians)
1176  {
1177  // grad yz
1178  if (!(evaluation_flag & EvaluationFlags::gradients))
1179  eval1.template gradients<1, true, false>(temp1, temp2);
1180  eval2.template gradients<2, true, false>(temp2,
1181  hessians_quad +
1182  5 * n_q_points);
1183 
1184  // grad yy
1185  eval1.template hessians<1, true, false>(temp1, temp2);
1186  eval2.template values<2, true, false>(temp2,
1187  hessians_quad + n_q_points);
1188  }
1189 
1190  // grad z: can use the values applied in x direction stored in
1191  // temp1
1192  eval1.template values<1, true, false>(temp1, temp2);
1193  if (evaluation_flag & EvaluationFlags::gradients)
1194  eval2.template gradients<2, true, false>(temp2,
1195  gradients_quad +
1196  2 * n_q_points);
1197 
1198  // grad zz: can use the values applied in x and y direction stored
1199  // in temp2
1200  if (evaluation_flag & EvaluationFlags::hessians)
1201  eval2.template hessians<2, true, false>(temp2,
1202  hessians_quad +
1203  2 * n_q_points);
1204 
1205  // val: can use the values applied in x & y direction stored in
1206  // temp2
1207  if (evaluation_flag & EvaluationFlags::values)
1208  eval2.template values<2, true, false>(temp2, values_quad);
1209  break;
1210  default:
1211  AssertThrow(false, ExcNotImplemented());
1212  }
1213  }
1214 
1215  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1216  template <int normal_dir>
1217  inline void
1219  dim,
1220  fe_degree,
1221  n_q_points_1d,
1222  Number>::
1223  evaluate_tensor_product_per_component(
1224  const EvaluationFlags::EvaluationFlags evaluation_flag,
1225  Number * values_dofs_actual,
1227  const bool add_into_values_array,
1228  std::integral_constant<bool, true>)
1229  {
1230  using EvalNormal =
1232  dim,
1233  (fe_degree == -1) ? 1 : fe_degree + 1,
1234  n_q_points_1d,
1235  Number,
1236  normal_dir>;
1237 
1238  using EvalTangent =
1240  dim,
1241  (fe_degree == -1) ? 1 : fe_degree,
1242  n_q_points_1d,
1243  Number,
1244  normal_dir>;
1245  using Eval0 =
1246  typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
1247  using Eval1 =
1248  typename std::conditional<normal_dir == 1, EvalNormal, EvalTangent>::type;
1249  using Eval2 =
1250  typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
1251 
1252  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info =
1253  fe_eval.get_shape_info();
1254  Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
1255  ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]));
1256  Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
1257  ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]));
1258  Eval2 eval2 = create_evaluator_tensor_product<Eval2>(
1259  ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1]));
1260 
1261  Number *temp1 = fe_eval.get_scratch_data().begin();
1262  Number *temp2;
1263 
1264  temp2 =
1265  temp1 +
1266  std::max(Utilities::fixed_power<dim>(shape_info.data[0].fe_degree + 1),
1267  Utilities::fixed_power<dim>(shape_info.data[0].n_q_points_1d));
1268 
1269  const std::size_t n_q_points = shape_info.n_q_points;
1270  const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell;
1271 
1272  // Initial shift depending on component (normal_dir)
1273  Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
1274  Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
1275  Number *gradients_quad =
1276  fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
1277  Number *hessians_quad =
1278  (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
1279  fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
1280 
1281  // Integrate path
1282  switch (dim)
1283  {
1284  case 2:
1285  if ((evaluation_flag & EvaluationFlags::values) &&
1286  !(evaluation_flag & EvaluationFlags::gradients))
1287  {
1288  eval1.template values<1, false, false>(values_quad, temp1);
1289  if (add_into_values_array == false)
1290  eval0.template values<0, false, false>(temp1, values_dofs);
1291  else
1292  eval0.template values<0, false, true>(temp1, values_dofs);
1293  }
1294  if (evaluation_flag & EvaluationFlags::gradients)
1295  {
1296  eval1.template gradients<1, false, false>(gradients_quad +
1297  n_q_points,
1298  temp1);
1299  if ((evaluation_flag & EvaluationFlags::values))
1300  eval1.template values<1, false, true>(values_quad, temp1);
1301  if (add_into_values_array == false)
1302  eval0.template values<0, false, false>(temp1, values_dofs);
1303  else
1304  eval0.template values<0, false, true>(temp1, values_dofs);
1305  eval1.template values<1, false, false>(gradients_quad, temp1);
1306  eval0.template gradients<0, false, true>(temp1, values_dofs);
1307  }
1308  if (evaluation_flag & EvaluationFlags::hessians)
1309  {
1310  // grad xx
1311  eval1.template values<1, false, false>(hessians_quad, temp1);
1312 
1313  if ((evaluation_flag & EvaluationFlags::values) ||
1314  (evaluation_flag & EvaluationFlags::gradients) ||
1315  add_into_values_array == true)
1316  eval0.template hessians<0, false, true>(temp1, values_dofs);
1317  else
1318  eval0.template hessians<0, false, false>(temp1, values_dofs);
1319 
1320  // grad yy
1321  eval1.template hessians<1, false, false>(hessians_quad +
1322  n_q_points,
1323  temp1);
1324  eval0.template values<0, false, true>(temp1, values_dofs);
1325 
1326  // grad xy
1327  eval1.template gradients<1, false, false>(hessians_quad +
1328  2 * n_q_points,
1329  temp1);
1330  eval0.template gradients<0, false, true>(temp1, values_dofs);
1331  }
1332  break;
1333 
1334  case 3:
1335  if ((evaluation_flag & EvaluationFlags::values) &&
1336  !(evaluation_flag & EvaluationFlags::gradients))
1337  {
1338  eval2.template values<2, false, false>(values_quad, temp1);
1339  eval1.template values<1, false, false>(temp1, temp2);
1340  if (add_into_values_array == false)
1341  eval0.template values<0, false, false>(temp2, values_dofs);
1342  else
1343  eval0.template values<0, false, true>(temp2, values_dofs);
1344  }
1345  if (evaluation_flag & EvaluationFlags::gradients)
1346  {
1347  eval2.template gradients<2, false, false>(gradients_quad +
1348  2 * n_q_points,
1349  temp1);
1350  if ((evaluation_flag & EvaluationFlags::values))
1351  eval2.template values<2, false, true>(values_quad, temp1);
1352  eval1.template values<1, false, false>(temp1, temp2);
1353  eval2.template values<2, false, false>(gradients_quad +
1354  n_q_points,
1355  temp1);
1356  eval1.template gradients<1, false, true>(temp1, temp2);
1357  if (add_into_values_array == false)
1358  eval0.template values<0, false, false>(temp2, values_dofs);
1359  else
1360  eval0.template values<0, false, true>(temp2, values_dofs);
1361  eval2.template values<2, false, false>(gradients_quad, temp1);
1362  eval1.template values<1, false, false>(temp1, temp2);
1363  eval0.template gradients<0, false, true>(temp2, values_dofs);
1364  }
1365  if (evaluation_flag & EvaluationFlags::hessians)
1366  {
1367  // grad xx
1368  eval2.template values<2, false, false>(hessians_quad, temp1);
1369  eval1.template values<1, false, false>(temp1, temp2);
1370 
1371  if ((evaluation_flag & EvaluationFlags::values) ||
1372  (evaluation_flag & EvaluationFlags::gradients) ||
1373  add_into_values_array == true)
1374  eval0.template hessians<0, false, true>(temp2, values_dofs);
1375  else
1376  eval0.template hessians<0, false, false>(temp2, values_dofs);
1377 
1378  // grad yy
1379  eval2.template values<2, false, false>(hessians_quad + n_q_points,
1380  temp1);
1381  eval1.template hessians<1, false, false>(temp1, temp2);
1382  eval0.template values<0, false, true>(temp2, values_dofs);
1383 
1384  // grad zz
1385  eval2.template hessians<2, false, false>(hessians_quad +
1386  2 * n_q_points,
1387  temp1);
1388  eval1.template values<1, false, false>(temp1, temp2);
1389  eval0.template values<0, false, true>(temp2, values_dofs);
1390 
1391  // grad xy
1392  eval2.template values<2, false, false>(hessians_quad +
1393  3 * n_q_points,
1394  temp1);
1395  eval1.template gradients<1, false, false>(temp1, temp2);
1396  eval0.template gradients<0, false, true>(temp2, values_dofs);
1397 
1398  // grad xz
1399  eval2.template gradients<2, false, false>(hessians_quad +
1400  4 * n_q_points,
1401  temp1);
1402  eval1.template values<1, false, false>(temp1, temp2);
1403  eval0.template gradients<0, false, true>(temp2, values_dofs);
1404 
1405  // grad yz
1406  eval2.template gradients<2, false, false>(hessians_quad +
1407  5 * n_q_points,
1408  temp1);
1409  eval1.template gradients<1, false, false>(temp1, temp2);
1410  eval0.template values<0, false, true>(temp2, values_dofs);
1411  }
1412 
1413  break;
1414  default:
1415  AssertThrow(false, ExcNotImplemented());
1416  }
1417  }
1418 
1428  template <EvaluatorVariant variant,
1429  EvaluatorQuantity quantity,
1430  int dim,
1431  int basis_size_1,
1432  int basis_size_2,
1433  typename Number,
1434  typename Number2>
1436  {
1437  static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
1438  "The second dimension must not be smaller than the first");
1439 
1462 #ifndef DEBUG
1464 #endif
1465  static void
1467  const unsigned int n_components,
1468  const AlignedVector<Number2> &transformation_matrix,
1469  const Number * values_in,
1470  Number * values_out,
1471  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
1472  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
1473  {
1474  Assert(
1475  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
1476  ExcMessage("The second dimension must not be smaller than the first"));
1477 
1479 
1480  // we do recursion until dim==1 or dim==2 and we have
1481  // basis_size_1==basis_size_2. The latter optimization increases
1482  // optimization possibilities for the compiler but does only work for
1483  // aliased pointers if the sizes are equal.
1484  constexpr int next_dim =
1485  (dim > 2 ||
1486  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
1487  dim - 1 :
1488  dim;
1489 
1490  EvaluatorTensorProduct<variant,
1491  dim,
1492  basis_size_1,
1493  (basis_size_1 == 0 ? 0 : basis_size_2),
1494  Number,
1495  Number2>
1496  eval_val(transformation_matrix,
1499  basis_size_1_variable,
1500  basis_size_2_variable);
1501  const unsigned int np_1 =
1502  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
1503  const unsigned int np_2 =
1504  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
1505  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
1506  ExcMessage("Cannot transform with 0-point basis"));
1507  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
1508  ExcMessage("Cannot transform with 0-point basis"));
1509 
1510  // run loop backwards to ensure correctness if values_in aliases with
1511  // values_out in case with basis_size_1 < basis_size_2
1512  values_in = values_in + n_components * Utilities::fixed_power<dim>(np_1);
1513  values_out =
1514  values_out + n_components * Utilities::fixed_power<dim>(np_2);
1515  for (unsigned int c = n_components; c != 0; --c)
1516  {
1517  values_in -= Utilities::fixed_power<dim>(np_1);
1518  values_out -= Utilities::fixed_power<dim>(np_2);
1519  if (next_dim < dim)
1520  for (unsigned int q = np_1; q != 0; --q)
1522  variant,
1523  quantity,
1524  next_dim,
1525  basis_size_1,
1526  basis_size_2,
1527  Number,
1528  Number2>::do_forward(1,
1529  transformation_matrix,
1530  values_in +
1531  (q - 1) *
1532  Utilities::fixed_power<next_dim>(np_1),
1533  values_out +
1534  (q - 1) *
1535  Utilities::fixed_power<next_dim>(np_2),
1536  basis_size_1_variable,
1537  basis_size_2_variable);
1538 
1539  // the recursion stops if dim==1 or if dim==2 and
1540  // basis_size_1==basis_size_2 (the latter is used because the
1541  // compiler generates nicer code)
1542  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
1543  {
1544  eval_val.template values<0, true, false>(values_in, values_out);
1545  eval_val.template values<1, true, false>(values_out, values_out);
1546  }
1547  else if (dim == 1)
1548  eval_val.template values<dim - 1, true, false>(values_in,
1549  values_out);
1550  else
1551  eval_val.template values<dim - 1, true, false>(values_out,
1552  values_out);
1553  }
1554  }
1555 
1586 #ifndef DEBUG
1588 #endif
1589  static void
1591  const unsigned int n_components,
1592  const AlignedVector<Number2> &transformation_matrix,
1593  const bool add_into_result,
1594  Number * values_in,
1595  Number * values_out,
1596  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
1597  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
1598  {
1599  Assert(
1600  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
1601  ExcMessage("The second dimension must not be smaller than the first"));
1602  Assert(add_into_result == false || values_in != values_out,
1603  ExcMessage(
1604  "Input and output cannot alias with each other when "
1605  "adding the result of the basis change to existing data"));
1606 
1607  Assert(quantity == EvaluatorQuantity::value ||
1608  quantity == EvaluatorQuantity::hessian,
1609  ExcInternalError());
1610 
1611  constexpr int next_dim =
1612  (dim > 2 ||
1613  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
1614  dim - 1 :
1615  dim;
1616  EvaluatorTensorProduct<variant,
1617  dim,
1618  basis_size_1,
1619  (basis_size_1 == 0 ? 0 : basis_size_2),
1620  Number,
1621  Number2>
1622  eval_val(transformation_matrix,
1623  transformation_matrix,
1624  transformation_matrix,
1625  basis_size_1_variable,
1626  basis_size_2_variable);
1627  const unsigned int np_1 =
1628  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
1629  const unsigned int np_2 =
1630  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
1631  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
1632  ExcMessage("Cannot transform with 0-point basis"));
1633  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
1634  ExcMessage("Cannot transform with 0-point basis"));
1635 
1636  for (unsigned int c = 0; c < n_components; ++c)
1637  {
1638  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
1639  {
1640  if (quantity == EvaluatorQuantity::value)
1641  eval_val.template values<1, false, false>(values_in, values_in);
1642  else
1643  eval_val.template hessians<1, false, false>(values_in,
1644  values_in);
1645 
1646  if (add_into_result)
1647  {
1648  if (quantity == EvaluatorQuantity::value)
1649  eval_val.template values<0, false, true>(values_in,
1650  values_out);
1651  else
1652  eval_val.template hessians<0, false, true>(values_in,
1653  values_out);
1654  }
1655  else
1656  {
1657  if (quantity == EvaluatorQuantity::value)
1658  eval_val.template values<0, false, false>(values_in,
1659  values_out);
1660  else
1661  eval_val.template hessians<0, false, false>(values_in,
1662  values_out);
1663  }
1664  }
1665  else
1666  {
1667  if (dim == 1 && add_into_result)
1668  {
1669  if (quantity == EvaluatorQuantity::value)
1670  eval_val.template values<0, false, true>(values_in,
1671  values_out);
1672  else
1673  eval_val.template hessians<0, false, true>(values_in,
1674  values_out);
1675  }
1676  else if (dim == 1)
1677  {
1678  if (quantity == EvaluatorQuantity::value)
1679  eval_val.template values<0, false, false>(values_in,
1680  values_out);
1681  else
1682  eval_val.template hessians<0, false, false>(values_in,
1683  values_out);
1684  }
1685  else
1686  {
1687  if (quantity == EvaluatorQuantity::value)
1688  eval_val.template values<dim - 1, false, false>(values_in,
1689  values_in);
1690  else
1691  eval_val.template hessians<dim - 1, false, false>(
1692  values_in, values_in);
1693  }
1694  }
1695  if (next_dim < dim)
1696  for (unsigned int q = 0; q < np_1; ++q)
1698  quantity,
1699  next_dim,
1700  basis_size_1,
1701  basis_size_2,
1702  Number,
1703  Number2>::
1704  do_backward(1,
1705  transformation_matrix,
1706  add_into_result,
1707  values_in +
1708  q * Utilities::fixed_power<next_dim>(np_2),
1709  values_out +
1710  q * Utilities::fixed_power<next_dim>(np_1),
1711  basis_size_1_variable,
1712  basis_size_2_variable);
1713 
1714  values_in += Utilities::fixed_power<dim>(np_2);
1715  values_out += Utilities::fixed_power<dim>(np_1);
1716  }
1717  }
1718 
1739  static void
1740  do_mass(const unsigned int n_components,
1741  const AlignedVector<Number2> &transformation_matrix,
1742  const AlignedVector<Number> & coefficients,
1743  const Number * values_in,
1744  Number * scratch_data,
1745  Number * values_out)
1746  {
1747  constexpr int next_dim = dim > 1 ? dim - 1 : dim;
1748  Number * my_scratch =
1749  basis_size_1 != basis_size_2 ? scratch_data : values_out;
1750 
1751  const unsigned int size_per_component = Utilities::pow(basis_size_2, dim);
1752  Assert(coefficients.size() == size_per_component ||
1753  coefficients.size() == n_components * size_per_component,
1754  ExcDimensionMismatch(coefficients.size(), size_per_component));
1755  const unsigned int stride =
1756  coefficients.size() == size_per_component ? 0 : 1;
1757 
1758  for (unsigned int q = basis_size_1; q != 0; --q)
1760  variant,
1762  next_dim,
1763  basis_size_1,
1764  basis_size_2,
1765  Number,
1766  Number2>::do_forward(n_components,
1767  transformation_matrix,
1768  values_in +
1769  (q - 1) *
1770  Utilities::pow(basis_size_1, dim - 1),
1771  my_scratch +
1772  (q - 1) *
1773  Utilities::pow(basis_size_2, dim - 1));
1774  EvaluatorTensorProduct<variant,
1775  dim,
1776  basis_size_1,
1777  basis_size_2,
1778  Number,
1779  Number2>
1780  eval_val(transformation_matrix);
1781  const unsigned int n_inner_blocks =
1782  (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
1783  const unsigned int n_blocks = Utilities::pow(basis_size_2, dim - 1);
1784  for (unsigned int ii = 0; ii < n_blocks; ii += n_inner_blocks)
1785  for (unsigned int c = 0; c < n_components; ++c)
1786  {
1787  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1788  eval_val.template values_one_line<dim - 1, true, false>(
1789  my_scratch + i, my_scratch + i);
1790  for (unsigned int q = 0; q < basis_size_2; ++q)
1791  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1792  my_scratch[i + q * n_blocks + c * size_per_component] *=
1793  coefficients[i + q * n_blocks +
1794  c * stride * size_per_component];
1795  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1796  eval_val.template values_one_line<dim - 1, false, false>(
1797  my_scratch + i, my_scratch + i);
1798  }
1799  for (unsigned int q = 0; q < basis_size_1; ++q)
1801  variant,
1803  next_dim,
1804  basis_size_1,
1805  basis_size_2,
1806  Number,
1807  Number2>::do_backward(n_components,
1808  transformation_matrix,
1809  false,
1810  my_scratch +
1811  q * Utilities::pow(basis_size_2, dim - 1),
1812  values_out +
1813  q * Utilities::pow(basis_size_1, dim - 1));
1814  }
1815  };
1816 
1817 
1818 
1831  template <int dim, int fe_degree, typename Number>
1833  {
1834  static void
1835  evaluate(const unsigned int n_components,
1836  const EvaluationFlags::EvaluationFlags evaluation_flag,
1837  const Number * values_dofs,
1839 
1840  static void
1842  const EvaluationFlags::EvaluationFlags evaluation_flag,
1843  const Number * values_dofs,
1844  Number * gradients_quad,
1845  Number * hessians_quad);
1846 
1847  static void
1848  integrate(const unsigned int n_components,
1849  const EvaluationFlags::EvaluationFlags integration_flag,
1850  Number * values_dofs,
1852  const bool add_into_values_array);
1853 
1854  static void
1856  const EvaluationFlags::EvaluationFlags integration_flag,
1857  Number * values_dofs,
1858  Number * gradients_quad,
1859  const Number * hessians_quad,
1860  const bool add_into_values_array);
1861  };
1862 
1863 
1864 
1865  template <int dim, int fe_degree, typename Number>
1866  inline void
1868  const unsigned int n_components,
1869  const EvaluationFlags::EvaluationFlags evaluation_flag,
1870  const Number * values_dofs,
1872  {
1873  constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1874 
1875  for (unsigned int c = 0; c < n_components; ++c)
1876  {
1877  if ((evaluation_flag & EvaluationFlags::values) != 0u)
1878  for (unsigned int i = 0; i < n_points; ++i)
1879  fe_eval.begin_values()[n_points * c + i] =
1880  values_dofs[n_points * c + i];
1881 
1882  do_evaluate(fe_eval.get_shape_info().data.front(),
1883  evaluation_flag,
1884  values_dofs + c * n_points,
1885  fe_eval.begin_gradients() + c * dim * n_points,
1886  fe_eval.begin_hessians() +
1887  c * dim * (dim + 1) / 2 * n_points);
1888  }
1889  }
1890 
1891 
1892 
1893  template <int dim, int fe_degree, typename Number>
1894  inline void
1897  const EvaluationFlags::EvaluationFlags evaluation_flag,
1898  const Number * values_dofs,
1899  Number * gradients_quad,
1900  Number * hessians_quad)
1901  {
1903  (fe_degree + 2) / 2 * (fe_degree + 1));
1904  constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1905 
1907  dim,
1908  fe_degree + 1,
1909  fe_degree + 1,
1910  Number>
1911  eval(AlignedVector<Number>(),
1914  if ((evaluation_flag &
1916  {
1917  eval.template gradients<0, true, false>(values_dofs, gradients_quad);
1918  if (dim > 1)
1919  eval.template gradients<1, true, false>(values_dofs,
1920  gradients_quad + n_points);
1921  if (dim > 2)
1922  eval.template gradients<2, true, false>(values_dofs,
1923  gradients_quad +
1924  2 * n_points);
1925  }
1926  if (evaluation_flag & EvaluationFlags::hessians)
1927  {
1928  eval.template hessians<0, true, false>(values_dofs, hessians_quad);
1929  if (dim > 1)
1930  {
1931  eval.template gradients<1, true, false>(gradients_quad,
1932  hessians_quad +
1933  dim * n_points);
1934  eval.template hessians<1, true, false>(values_dofs,
1935  hessians_quad + n_points);
1936  }
1937  if (dim > 2)
1938  {
1939  eval.template gradients<2, true, false>(gradients_quad,
1940  hessians_quad +
1941  4 * n_points);
1942  eval.template gradients<2, true, false>(gradients_quad + n_points,
1943  hessians_quad +
1944  5 * n_points);
1945  eval.template hessians<2, true, false>(values_dofs,
1946  hessians_quad +
1947  2 * n_points);
1948  }
1949  }
1950  }
1951 
1952 
1953 
1954  template <int dim, int fe_degree, typename Number>
1955  inline void
1957  const unsigned int n_components,
1958  const EvaluationFlags::EvaluationFlags integration_flag,
1959  Number * values_dofs,
1961  const bool add_into_values_array)
1962  {
1963  constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1964 
1965  for (unsigned int c = 0; c < n_components; ++c)
1966  {
1967  if ((integration_flag & EvaluationFlags::values) != 0u)
1968  {
1969  if (add_into_values_array)
1970  for (unsigned int i = 0; i < n_points; ++i)
1971  values_dofs[n_points * c + i] +=
1972  fe_eval.begin_values()[n_points * c + i];
1973  else
1974  for (unsigned int i = 0; i < n_points; ++i)
1975  values_dofs[n_points * c + i] =
1976  fe_eval.begin_values()[n_points * c + i];
1977  }
1978 
1979  do_integrate(fe_eval.get_shape_info().data.front(),
1980  integration_flag,
1981  values_dofs + c * n_points,
1982  fe_eval.begin_gradients() + c * dim * n_points,
1983  fe_eval.begin_hessians() +
1984  c * dim * (dim + 1) / 2 * n_points,
1985  add_into_values_array ||
1986  ((integration_flag & EvaluationFlags::values) != 0u));
1987  }
1988  }
1989 
1990 
1991 
1992  template <int dim, int fe_degree, typename Number>
1993  inline void
1996  const EvaluationFlags::EvaluationFlags integration_flag,
1997  Number * values_dofs,
1998  Number * gradients_quad,
1999  const Number * hessians_quad,
2000  const bool add_into_values_array)
2001  {
2003  (fe_degree + 2) / 2 * (fe_degree + 1));
2004 
2006  dim,
2007  fe_degree + 1,
2008  fe_degree + 1,
2009  Number>
2010  eval(AlignedVector<Number>(),
2013  constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
2014 
2015  if ((integration_flag & EvaluationFlags::hessians) != 0u)
2016  {
2017  // diagonal
2018  // grad xx
2019  if (add_into_values_array == true)
2020  eval.template hessians<0, false, true>(hessians_quad, values_dofs);
2021  else
2022  eval.template hessians<0, false, false>(hessians_quad, values_dofs);
2023  // grad yy
2024  if (dim > 1)
2025  eval.template hessians<1, false, true>(hessians_quad + n_points,
2026  values_dofs);
2027  // grad zz
2028  if (dim > 2)
2029  eval.template hessians<2, false, true>(hessians_quad + 2 * n_points,
2030  values_dofs);
2031  // off-diagonal
2032  if (dim == 2)
2033  {
2034  // grad xy, queue into gradient
2035  if (integration_flag & EvaluationFlags::gradients)
2036  eval.template gradients<1, false, true>(hessians_quad +
2037  2 * n_points,
2038  gradients_quad);
2039  else
2040  eval.template gradients<1, false, false>(hessians_quad +
2041  2 * n_points,
2042  gradients_quad);
2043  }
2044  if (dim == 3)
2045  {
2046  // grad xy, queue into gradient
2047  if (integration_flag & EvaluationFlags::gradients)
2048  eval.template gradients<1, false, true>(hessians_quad +
2049  3 * n_points,
2050  gradients_quad);
2051  else
2052  eval.template gradients<1, false, false>(hessians_quad +
2053  3 * n_points,
2054  gradients_quad);
2055 
2056  // grad xz
2057  eval.template gradients<2, false, true>(hessians_quad +
2058  4 * n_points,
2059  gradients_quad);
2060 
2061  // grad yz
2062  if (integration_flag & EvaluationFlags::gradients)
2063  eval.template gradients<2, false, true>(
2064  hessians_quad + 5 * n_points, gradients_quad + n_points);
2065  else
2066  eval.template gradients<2, false, false>(
2067  hessians_quad + 5 * n_points, gradients_quad + n_points);
2068  }
2069 
2070  // if we did not integrate gradients, set the last slot to zero
2071  // which was not touched before, in order to avoid the if
2072  // statement in the gradients loop below
2073  if ((integration_flag & EvaluationFlags::gradients) == 0u)
2074  for (unsigned int q = 0; q < n_points; ++q)
2075  gradients_quad[(dim - 1) * n_points + q] = Number();
2076  }
2077 
2078  if ((integration_flag &
2080  {
2081  if (add_into_values_array ||
2082  (integration_flag & EvaluationFlags::hessians) != 0u)
2083  eval.template gradients<0, false, true>(gradients_quad, values_dofs);
2084  else
2085  eval.template gradients<0, false, false>(gradients_quad, values_dofs);
2086  if (dim > 1)
2087  eval.template gradients<1, false, true>(gradients_quad + n_points,
2088  values_dofs);
2089  if (dim > 2)
2090  eval.template gradients<2, false, true>(gradients_quad + 2 * n_points,
2091  values_dofs);
2092  }
2093  }
2094 
2095 
2096 
2107  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2109  {
2110  static void
2111  evaluate(const unsigned int n_components,
2112  const EvaluationFlags::EvaluationFlags evaluation_flag,
2113  const Number * values_dofs,
2115 
2116  static void
2117  integrate(const unsigned int n_components,
2118  const EvaluationFlags::EvaluationFlags evaluation_flag,
2119  Number * values_dofs,
2121  const bool add_into_values_array);
2122  };
2123 
2124 
2125 
2126  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2127  inline void
2129  dim,
2130  fe_degree,
2131  n_q_points_1d,
2132  Number>::evaluate(const unsigned int n_components,
2133  const EvaluationFlags::EvaluationFlags evaluation_flag,
2134  const Number * values_dofs,
2136  {
2137  const auto &shape_data = fe_eval.get_shape_info().data.front();
2138 
2139  Assert(n_q_points_1d > fe_degree,
2140  ExcMessage("You lose information when going to a collocation space "
2141  "of lower degree, so the evaluation results would be "
2142  "wrong. Thus, this class does not permit the desired "
2143  "operation."));
2144  constexpr std::size_t n_dofs = Utilities::pow(fe_degree + 1, dim);
2145  constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
2146 
2147  for (unsigned int c = 0; c < n_components; ++c)
2148  {
2152  dim,
2153  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
2154  n_q_points_1d,
2155  Number,
2156  Number>::do_forward(1,
2157  shape_data.shape_values_eo,
2158  values_dofs + c * n_dofs,
2159  fe_eval.begin_values() + c * n_q_points);
2160 
2161  // apply derivatives in the collocation space
2162  if (evaluation_flag &
2165  do_evaluate(shape_data,
2166  evaluation_flag & (EvaluationFlags::gradients |
2168  fe_eval.begin_values() + c * n_q_points,
2169  fe_eval.begin_gradients() + c * dim * n_q_points,
2170  fe_eval.begin_hessians() +
2171  c * dim * (dim + 1) / 2 * n_q_points);
2172  }
2173  }
2174 
2175 
2176 
2177  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2178  inline void
2180  dim,
2181  fe_degree,
2182  n_q_points_1d,
2183  Number>::integrate(const unsigned int n_components,
2184  const EvaluationFlags::EvaluationFlags integration_flag,
2185  Number * values_dofs,
2187  const bool add_into_values_array)
2188  {
2189  const auto &shape_data = fe_eval.get_shape_info().data.front();
2190 
2191  Assert(n_q_points_1d > fe_degree,
2192  ExcMessage("You lose information when going to a collocation space "
2193  "of lower degree, so the evaluation results would be "
2194  "wrong. Thus, this class does not permit the desired "
2195  "operation."));
2196  constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
2197 
2198  for (unsigned int c = 0; c < n_components; ++c)
2199  {
2200  // apply derivatives in collocation space
2201  if (integration_flag &
2204  do_integrate(shape_data,
2205  integration_flag & (EvaluationFlags::gradients |
2207  fe_eval.begin_values() + c * n_q_points,
2208  fe_eval.begin_gradients() + c * dim * n_q_points,
2209  fe_eval.begin_hessians() +
2210  c * dim * (dim + 1) / 2 * n_q_points,
2211  /*add_into_values_array=*/
2212  integration_flag & EvaluationFlags::values);
2213 
2214  // transform back to the original space
2218  dim,
2219  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
2220  n_q_points_1d,
2221  Number,
2222  Number>::do_backward(1,
2223  shape_data.shape_values_eo,
2224  add_into_values_array,
2225  fe_eval.begin_values() + c * n_q_points,
2226  values_dofs +
2227  c * Utilities::pow(fe_degree + 1, dim));
2228  }
2229  }
2230 
2231 
2232 
2240  constexpr bool
2241  use_collocation_evaluation(const unsigned int fe_degree,
2242  const unsigned int n_q_points_1d)
2243  {
2244  return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
2245  (n_q_points_1d <= 3 * fe_degree / 2 + 1);
2246  }
2247 
2248 
2264  template <int dim, typename Number>
2266  {
2267  template <int fe_degree, int n_q_points_1d>
2268  static bool
2269  run(const unsigned int n_components,
2270  const EvaluationFlags::EvaluationFlags evaluation_flag,
2271  const Number * values_dofs,
2273  {
2274  const auto element_type = fe_eval.get_shape_info().element_type;
2276 
2277  Assert(fe_eval.get_shape_info().data.size() == 1 ||
2278  (fe_eval.get_shape_info().data.size() == dim &&
2279  element_type == ElementType::tensor_general) ||
2280  element_type == ElementType::tensor_raviart_thomas,
2281  ExcNotImplemented());
2282 
2283  if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
2285  {
2287  n_components, evaluation_flag, values_dofs, fe_eval);
2288  }
2289  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
2290  // shape_info.h for more details
2291  else if (fe_degree >= 0 &&
2292  use_collocation_evaluation(fe_degree, n_q_points_1d) &&
2293  element_type <= ElementType::tensor_symmetric)
2294  {
2296  dim,
2297  fe_degree,
2298  n_q_points_1d,
2299  Number>::evaluate(n_components,
2300  evaluation_flag,
2301  values_dofs,
2302  fe_eval);
2303  }
2304  else if (fe_degree >= 0 &&
2306  {
2308  dim,
2309  fe_degree,
2310  n_q_points_1d,
2311  Number>::evaluate(n_components,
2312  evaluation_flag,
2313  values_dofs,
2314  fe_eval);
2315  }
2316  else if (element_type == ElementType::tensor_symmetric_plus_dg0)
2317  {
2319  dim,
2320  fe_degree,
2321  n_q_points_1d,
2322  Number>::evaluate(n_components,
2323  evaluation_flag,
2324  values_dofs,
2325  fe_eval);
2326  }
2327  else if (element_type == ElementType::truncated_tensor)
2328  {
2330  dim,
2331  fe_degree,
2332  n_q_points_1d,
2333  Number>::evaluate(n_components,
2334  evaluation_flag,
2335  values_dofs,
2336  fe_eval);
2337  }
2338  else if (element_type == ElementType::tensor_none)
2339  {
2341  dim,
2342  fe_degree,
2343  n_q_points_1d,
2344  Number>::evaluate(n_components,
2345  evaluation_flag,
2346  values_dofs,
2347  fe_eval);
2348  }
2349  else if (element_type == ElementType::tensor_raviart_thomas)
2350  {
2353  dim,
2354  (fe_degree == -1) ? 1 : fe_degree,
2355  (n_q_points_1d < 1) ? 1 : n_q_points_1d,
2356  Number>::template evaluate_or_integrate<false>(evaluation_flag,
2357  const_cast<Number *>(
2358  values_dofs),
2359  fe_eval);
2360  }
2361  else
2362  {
2364  dim,
2365  fe_degree,
2366  n_q_points_1d,
2367  Number>::evaluate(n_components,
2368  evaluation_flag,
2369  values_dofs,
2370  fe_eval);
2371  }
2372 
2373  return false;
2374  }
2375  };
2376 
2377 
2378 
2394  template <int dim, typename Number>
2396  {
2397  template <int fe_degree, int n_q_points_1d>
2398  static bool
2399  run(const unsigned int n_components,
2400  const EvaluationFlags::EvaluationFlags integration_flag,
2401  Number * values_dofs,
2403  const bool sum_into_values_array)
2404  {
2405  const auto element_type = fe_eval.get_shape_info().element_type;
2407 
2408  Assert(fe_eval.get_shape_info().data.size() == 1 ||
2409  (fe_eval.get_shape_info().data.size() == dim &&
2410  element_type == ElementType::tensor_general) ||
2411  element_type == ElementType::tensor_raviart_thomas,
2412  ExcNotImplemented());
2413 
2414  if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
2416  {
2418  n_components,
2419  integration_flag,
2420  values_dofs,
2421  fe_eval,
2422  sum_into_values_array);
2423  }
2424  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
2425  // shape_info.h for more details
2426  else if (fe_degree >= 0 &&
2427  use_collocation_evaluation(fe_degree, n_q_points_1d) &&
2428  element_type <= ElementType::tensor_symmetric)
2429  {
2431  dim,
2432  fe_degree,
2433  n_q_points_1d,
2434  Number>::integrate(n_components,
2435  integration_flag,
2436  values_dofs,
2437  fe_eval,
2438  sum_into_values_array);
2439  }
2440  else if (fe_degree >= 0 &&
2442  {
2444  dim,
2445  fe_degree,
2446  n_q_points_1d,
2447  Number>::integrate(n_components,
2448  integration_flag,
2449  values_dofs,
2450  fe_eval,
2451  sum_into_values_array);
2452  }
2453  else if (element_type == ElementType::tensor_symmetric_plus_dg0)
2454  {
2456  dim,
2457  fe_degree,
2458  n_q_points_1d,
2459  Number>::integrate(n_components,
2460  integration_flag,
2461  values_dofs,
2462  fe_eval,
2463  sum_into_values_array);
2464  }
2465  else if (element_type == ElementType::truncated_tensor)
2466  {
2468  dim,
2469  fe_degree,
2470  n_q_points_1d,
2471  Number>::integrate(n_components,
2472  integration_flag,
2473  values_dofs,
2474  fe_eval,
2475  sum_into_values_array);
2476  }
2477  else if (element_type == ElementType::tensor_none)
2478  {
2480  dim,
2481  fe_degree,
2482  n_q_points_1d,
2483  Number>::integrate(n_components,
2484  integration_flag,
2485  values_dofs,
2486  fe_eval,
2487  sum_into_values_array);
2488  }
2489  else if (element_type == ElementType::tensor_raviart_thomas)
2490  {
2492  dim,
2493  (fe_degree == -1) ? 1 : fe_degree,
2494  (n_q_points_1d < 1) ? 1 : n_q_points_1d,
2495  Number>::
2496  template evaluate_or_integrate<true>(integration_flag,
2497  const_cast<Number *>(
2498  values_dofs),
2499  fe_eval,
2500  sum_into_values_array);
2501  }
2502  else
2503  {
2505  dim,
2506  fe_degree,
2507  n_q_points_1d,
2508  Number>::integrate(n_components,
2509  integration_flag,
2510  values_dofs,
2511  fe_eval,
2512  sum_into_values_array);
2513  }
2514 
2515  return false;
2516  }
2517  };
2518 
2519 
2520 
2521  template <bool symmetric_evaluate,
2522  int dim,
2523  int fe_degree,
2524  int n_q_points_1d,
2525  typename Number>
2527  {
2528  // We enable a transformation to collocation for derivatives if it gives
2529  // correct results (first two conditions), if it is the most efficient
2530  // choice in terms of operation counts (third condition) and if we were
2531  // able to initialize the fields in shape_info.templates.h from the
2532  // polynomials (fourth condition).
2533  using Eval = EvaluatorTensorProduct<symmetric_evaluate ? evaluate_evenodd :
2535  dim - 1,
2536  fe_degree + 1,
2537  n_q_points_1d,
2538  Number>;
2539 
2540  static Eval
2543  const unsigned int subface_index,
2544  const unsigned int direction)
2545  {
2546  if (symmetric_evaluate)
2547  return Eval(data.shape_values_eo,
2548  data.shape_gradients_eo,
2549  data.shape_hessians_eo,
2550  data.fe_degree + 1,
2551  data.n_q_points_1d);
2552  else if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
2553  return Eval(data.shape_values,
2554  data.shape_gradients,
2555  data.shape_hessians,
2556  data.fe_degree + 1,
2557  data.n_q_points_1d);
2558  else
2559  {
2560  const unsigned int index =
2561  direction == 0 ? subface_index % 2 : subface_index / 2;
2562  return Eval(data.values_within_subface[index],
2565  data.fe_degree + 1,
2566  data.n_q_points_1d);
2567  }
2568  }
2569 
2570  static void
2572  const unsigned int n_components,
2573  const EvaluationFlags::EvaluationFlags evaluation_flag,
2575  Number * values_dofs,
2576  Number * values_quad,
2577  Number * gradients_quad,
2578  Number * hessians_quad,
2579  Number * scratch_data,
2580  const unsigned int subface_index)
2581  {
2582  Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
2583  Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
2584 
2585  const std::size_t n_dofs = fe_degree > -1 ?
2586  Utilities::pow(fe_degree + 1, dim - 1) :
2587  Utilities::pow(data.fe_degree + 1, dim - 1);
2588  const std::size_t n_q_points =
2589  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
2590  Utilities::pow(data.n_q_points_1d, dim - 1);
2591 
2592  // keep a copy of the original pointer for the case of the Hessians
2593  Number *values_dofs_ptr = values_dofs;
2594 
2595  if ((evaluation_flag & EvaluationFlags::values) != 0u &&
2596  ((evaluation_flag & EvaluationFlags::gradients) == 0u))
2597  for (unsigned int c = 0; c < n_components; ++c)
2598  {
2599  switch (dim)
2600  {
2601  case 3:
2602  eval0.template values<0, true, false>(values_dofs,
2603  values_quad);
2604  eval1.template values<1, true, false>(values_quad,
2605  values_quad);
2606  break;
2607  case 2:
2608  eval0.template values<0, true, false>(values_dofs,
2609  values_quad);
2610  break;
2611  case 1:
2612  values_quad[0] = values_dofs[0];
2613  break;
2614  default:
2615  Assert(false, ExcNotImplemented());
2616  }
2617  // Note: we always keep storage of values, 1st and 2nd derivatives
2618  // in an array
2619  values_dofs += 3 * n_dofs;
2620  values_quad += n_q_points;
2621  }
2622  else if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2623  for (unsigned int c = 0; c < n_components; ++c)
2624  {
2625  switch (dim)
2626  {
2627  case 3:
2628  if (symmetric_evaluate &&
2629  use_collocation_evaluation(fe_degree, n_q_points_1d))
2630  {
2631  eval0.template values<0, true, false>(values_dofs,
2632  values_quad);
2633  eval0.template values<1, true, false>(values_quad,
2634  values_quad);
2636  dim - 1,
2637  n_q_points_1d,
2638  n_q_points_1d,
2639  Number>
2640  eval_grad(AlignedVector<Number>(),
2643  eval_grad.template gradients<0, true, false>(
2644  values_quad, gradients_quad);
2645  eval_grad.template gradients<1, true, false>(
2646  values_quad, gradients_quad + n_q_points);
2647  }
2648  else
2649  {
2650  // grad x
2651  eval0.template gradients<0, true, false>(values_dofs,
2652  scratch_data);
2653  eval1.template values<1, true, false>(scratch_data,
2654  gradients_quad);
2655 
2656  // grad y
2657  eval0.template values<0, true, false>(values_dofs,
2658  scratch_data);
2659  eval1.template gradients<1, true, false>(scratch_data,
2660  gradients_quad +
2661  n_q_points);
2662 
2663  if ((evaluation_flag & EvaluationFlags::values) != 0u)
2664  eval1.template values<1, true, false>(scratch_data,
2665  values_quad);
2666  }
2667  // grad z
2668  eval0.template values<0, true, false>(values_dofs + n_dofs,
2669  scratch_data);
2670  eval1.template values<1, true, false>(
2671  scratch_data, gradients_quad + (dim - 1) * n_q_points);
2672 
2673  break;
2674  case 2:
2675  eval0.template values<0, true, false>(values_dofs + n_dofs,
2676  gradients_quad +
2677  n_q_points);
2678  eval0.template gradients<0, true, false>(values_dofs,
2679  gradients_quad);
2680  if ((evaluation_flag & EvaluationFlags::values) != 0u)
2681  eval0.template values<0, true, false>(values_dofs,
2682  values_quad);
2683  break;
2684  case 1:
2685  values_quad[0] = values_dofs[0];
2686  gradients_quad[0] = values_dofs[1];
2687  break;
2688  default:
2689  AssertThrow(false, ExcNotImplemented());
2690  }
2691  values_dofs += 3 * n_dofs;
2692  values_quad += n_q_points;
2693  gradients_quad += dim * n_q_points;
2694  }
2695 
2696  if ((evaluation_flag & EvaluationFlags::hessians) != 0u)
2697  {
2698  values_dofs = values_dofs_ptr;
2699  for (unsigned int c = 0; c < n_components; ++c)
2700  {
2701  switch (dim)
2702  {
2703  case 3:
2704  // grad xx
2705  eval0.template hessians<0, true, false>(values_dofs,
2706  scratch_data);
2707  eval1.template values<1, true, false>(scratch_data,
2708  hessians_quad);
2709 
2710  // grad yy
2711  eval0.template values<0, true, false>(values_dofs,
2712  scratch_data);
2713  eval1.template hessians<1, true, false>(scratch_data,
2714  hessians_quad +
2715  n_q_points);
2716 
2717  // grad zz
2718  eval0.template values<0, true, false>(values_dofs +
2719  2 * n_dofs,
2720  scratch_data);
2721  eval1.template values<1, true, false>(scratch_data,
2722  hessians_quad +
2723  2 * n_q_points);
2724 
2725  // grad xy
2726  eval0.template gradients<0, true, false>(values_dofs,
2727  scratch_data);
2728  eval1.template gradients<1, true, false>(scratch_data,
2729  hessians_quad +
2730  3 * n_q_points);
2731 
2732  // grad xz
2733  eval0.template gradients<0, true, false>(values_dofs +
2734  n_dofs,
2735  scratch_data);
2736  eval1.template values<1, true, false>(scratch_data,
2737  hessians_quad +
2738  4 * n_q_points);
2739 
2740  // grad yz
2741  eval0.template values<0, true, false>(values_dofs + n_dofs,
2742  scratch_data);
2743  eval1.template gradients<1, true, false>(scratch_data,
2744  hessians_quad +
2745  5 * n_q_points);
2746 
2747  break;
2748  case 2:
2749  // grad xx
2750  eval0.template hessians<0, true, false>(values_dofs,
2751  hessians_quad);
2752  // grad yy
2753  eval0.template values<0, true, false>(
2754  values_dofs + 2 * n_dofs, hessians_quad + n_q_points);
2755  // grad xy
2756  eval0.template gradients<0, true, false>(
2757  values_dofs + n_dofs, hessians_quad + 2 * n_q_points);
2758  break;
2759  case 1:
2760  hessians_quad[0] = values_dofs[2];
2761  break;
2762  default:
2763  AssertThrow(false, ExcNotImplemented());
2764  }
2765  values_dofs += 3 * n_dofs;
2766  hessians_quad += dim * (dim + 1) / 2 * n_q_points;
2767  }
2768  }
2769  }
2770 
2771  static void
2773  const unsigned int n_components,
2774  const EvaluationFlags::EvaluationFlags integration_flag,
2776  Number * values_dofs,
2777  Number * values_quad,
2778  Number * gradients_quad,
2779  Number * hessians_quad,
2780  Number * scratch_data,
2781  const unsigned int subface_index)
2782  {
2783  Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
2784  Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
2785 
2786  const std::size_t n_dofs =
2787  fe_degree > -1 ?
2788  Utilities::pow(fe_degree + 1, dim - 1) :
2789  (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
2790  const std::size_t n_q_points =
2791  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
2792  Utilities::pow(data.n_q_points_1d, dim - 1);
2793 
2794  // keep a copy of the original pointer for the case of the Hessians
2795  Number *values_dofs_ptr = values_dofs;
2796 
2797  if ((integration_flag & EvaluationFlags::values) != 0u &&
2798  (integration_flag & EvaluationFlags::gradients) == 0u)
2799  for (unsigned int c = 0; c < n_components; ++c)
2800  {
2801  switch (dim)
2802  {
2803  case 3:
2804  eval1.template values<1, false, false>(values_quad,
2805  values_quad);
2806  eval0.template values<0, false, false>(values_quad,
2807  values_dofs);
2808  break;
2809  case 2:
2810  eval0.template values<0, false, false>(values_quad,
2811  values_dofs);
2812  break;
2813  case 1:
2814  values_dofs[0] = values_quad[0];
2815  break;
2816  default:
2817  Assert(false, ExcNotImplemented());
2818  }
2819  values_dofs += 3 * n_dofs;
2820  values_quad += n_q_points;
2821  }
2822  else if ((integration_flag & EvaluationFlags::gradients) != 0u)
2823  for (unsigned int c = 0; c < n_components; ++c)
2824  {
2825  switch (dim)
2826  {
2827  case 3:
2828  // grad z
2829  eval1.template values<1, false, false>(gradients_quad +
2830  2 * n_q_points,
2831  gradients_quad +
2832  2 * n_q_points);
2833  eval0.template values<0, false, false>(gradients_quad +
2834  2 * n_q_points,
2835  values_dofs + n_dofs);
2836  if (symmetric_evaluate &&
2837  use_collocation_evaluation(fe_degree, n_q_points_1d))
2838  {
2840  dim - 1,
2841  n_q_points_1d,
2842  n_q_points_1d,
2843  Number>
2844  eval_grad(AlignedVector<Number>(),
2847  if ((integration_flag & EvaluationFlags::values) != 0u)
2848  eval_grad.template gradients<1, false, true>(
2849  gradients_quad + n_q_points, values_quad);
2850  else
2851  eval_grad.template gradients<1, false, false>(
2852  gradients_quad + n_q_points, values_quad);
2853  eval_grad.template gradients<0, false, true>(
2854  gradients_quad, values_quad);
2855  eval0.template values<1, false, false>(values_quad,
2856  values_quad);
2857  eval0.template values<0, false, false>(values_quad,
2858  values_dofs);
2859  }
2860  else
2861  {
2862  if ((integration_flag & EvaluationFlags::values) != 0u)
2863  {
2864  eval1.template values<1, false, false>(values_quad,
2865  scratch_data);
2866  eval1.template gradients<1, false, true>(
2867  gradients_quad + n_q_points, scratch_data);
2868  }
2869  else
2870  eval1.template gradients<1, false, false>(
2871  gradients_quad + n_q_points, scratch_data);
2872 
2873  // grad y
2874  eval0.template values<0, false, false>(scratch_data,
2875  values_dofs);
2876 
2877  // grad x
2878  eval1.template values<1, false, false>(gradients_quad,
2879  scratch_data);
2880  eval0.template gradients<0, false, true>(scratch_data,
2881  values_dofs);
2882  }
2883  break;
2884  case 2:
2885  eval0.template values<0, false, false>(gradients_quad +
2886  n_q_points,
2887  values_dofs + n_dofs);
2888  eval0.template gradients<0, false, false>(gradients_quad,
2889  values_dofs);
2890  if ((integration_flag & EvaluationFlags::values) != 0u)
2891  eval0.template values<0, false, true>(values_quad,
2892  values_dofs);
2893  break;
2894  case 1:
2895  values_dofs[0] = values_quad[0];
2896  values_dofs[1] = gradients_quad[0];
2897  break;
2898  default:
2899  AssertThrow(false, ExcNotImplemented());
2900  }
2901  values_dofs += 3 * n_dofs;
2902  values_quad += n_q_points;
2903  gradients_quad += dim * n_q_points;
2904  }
2905 
2906  if ((integration_flag & EvaluationFlags::hessians) != 0u)
2907  {
2908  values_dofs = values_dofs_ptr;
2909  for (unsigned int c = 0; c < n_components; ++c)
2910  {
2911  switch (dim)
2912  {
2913  case 3:
2914  // grad xx
2915  eval1.template values<1, false, false>(hessians_quad,
2916  scratch_data);
2917  if ((integration_flag & (EvaluationFlags::values |
2919  eval0.template hessians<0, false, true>(scratch_data,
2920  values_dofs);
2921  else
2922  eval0.template hessians<0, false, false>(scratch_data,
2923  values_dofs);
2924 
2925  // grad yy
2926  eval1.template hessians<1, false, false>(hessians_quad +
2927  n_q_points,
2928  scratch_data);
2929  eval0.template values<0, false, true>(scratch_data,
2930  values_dofs);
2931 
2932  // grad zz
2933  eval1.template values<1, false, false>(hessians_quad +
2934  2 * n_q_points,
2935  scratch_data);
2936  eval0.template values<0, false, false>(scratch_data,
2937  values_dofs +
2938  2 * n_dofs);
2939 
2940  // grad xy
2941  eval1.template gradients<1, false, false>(hessians_quad +
2942  3 * n_q_points,
2943  scratch_data);
2944  eval0.template gradients<0, false, true>(scratch_data,
2945  values_dofs);
2946 
2947  // grad xz
2948  eval1.template values<1, false, false>(hessians_quad +
2949  4 * n_q_points,
2950  scratch_data);
2951  if ((integration_flag & EvaluationFlags::gradients) != 0u)
2952  eval0.template gradients<0, false, true>(scratch_data,
2953  values_dofs +
2954  n_dofs);
2955  else
2956  eval0.template gradients<0, false, false>(scratch_data,
2957  values_dofs +
2958  n_dofs);
2959 
2960  // grad yz
2961  eval1.template gradients<1, false, false>(hessians_quad +
2962  5 * n_q_points,
2963  scratch_data);
2964  eval0.template values<0, false, true>(scratch_data,
2965  values_dofs + n_dofs);
2966 
2967  break;
2968  case 2:
2969  // grad xx
2970  if ((integration_flag & (EvaluationFlags::values |
2972  eval0.template hessians<0, false, true>(hessians_quad,
2973  values_dofs);
2974  else
2975  eval0.template hessians<0, false, false>(hessians_quad,
2976  values_dofs);
2977 
2978  // grad yy
2979  eval0.template values<0, false, false>(
2980  hessians_quad + n_q_points, values_dofs + 2 * n_dofs);
2981  // grad xy
2982  if ((integration_flag & EvaluationFlags::gradients) != 0u)
2983  eval0.template gradients<0, false, true>(
2984  hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
2985  else
2986  eval0.template gradients<0, false, false>(
2987  hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
2988  break;
2989  case 1:
2990  values_dofs[2] = hessians_quad[0];
2991  if ((integration_flag & EvaluationFlags::values) == 0u)
2992  values_dofs[0] = 0;
2993  if ((integration_flag & EvaluationFlags::gradients) == 0u)
2994  values_dofs[1] = 0;
2995  break;
2996  default:
2997  AssertThrow(false, ExcNotImplemented());
2998  }
2999  values_dofs += 3 * n_dofs;
3000  hessians_quad += dim * (dim + 1) / 2 * n_q_points;
3001  }
3002  }
3003  }
3004  };
3005 
3006  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
3008  {
3010  dim - 1,
3011  fe_degree,
3012  n_q_points_1d,
3013  Number>;
3014  template <typename EvalType>
3015  static EvalType
3018  const unsigned int subface_index,
3019  const unsigned int direction)
3020  {
3021  if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
3022  return EvalType(data.shape_values,
3023  data.shape_gradients,
3024  data.shape_hessians);
3025  else
3026  {
3027  const unsigned int index =
3028  direction == 0 ? subface_index % 2 : subface_index / 2;
3029  return EvalType(data.values_within_subface[index],
3032  }
3033  }
3034 
3035  template <bool integrate>
3036  static void
3038  const EvaluationFlags::EvaluationFlags evaluation_flag,
3039  Number * values_dofs,
3041  Number * scratch_data,
3042  const unsigned int subface_index,
3043  const unsigned int face_no)
3044  {
3045  const unsigned int face_direction = face_no / 2;
3046 
3047  // We first evaluate the anisotropic faces, i.e the faces where
3048  // face_direction != component. Note that the call order here is not
3049  // important, since the pointers are shifted accordingly within the
3050  // function. However, this is the order in which the components will be in
3051  // the quadrature points. Furthermore, the isotropic faces have no "normal
3052  // direction" but we still pass in normal_dir = 2 since this is used for
3053  // the pointers.
3054  // -----------------------------------------------------------------------------------
3055  // | | Anisotropic faces | Isotropic faces|
3056  // | Face dir | comp, coords, normal_dir | comp, coords, normal_dir | comp, coords |
3057  // | --------------------------------------------------------------------------------|
3058  // | 0 | 1, y, 0 | - | 0, y |
3059  // | 1 | 0, x, 0 | - | 1, x |
3060  // | --------------------------------------------------------------------------------|
3061  // | 0 | 1, yz, 0 | 2, yz, 1 | 0, yz |
3062  // | 1 | 2, zx, 0 | 0, zx, 1 | 1, zx |
3063  // | 2 | 0, xy, 0 | 1, xy, 1 | 2, xy |
3064  // -----------------------------------------------------------------------------------
3065  evaluate_in_face_apply<0>(values_dofs,
3066  fe_eval,
3067  scratch_data,
3068  evaluation_flag,
3069  face_direction,
3070  subface_index,
3071  std::integral_constant<bool, integrate>());
3072 
3073  if (dim == 3)
3074  evaluate_in_face_apply<1>(values_dofs,
3075  fe_eval,
3076  scratch_data,
3077  evaluation_flag,
3078  face_direction,
3079  subface_index,
3080  std::integral_constant<bool, integrate>());
3081 
3082  evaluate_in_face_apply<2>(values_dofs,
3083  fe_eval,
3084  scratch_data,
3085  evaluation_flag,
3086  face_direction,
3087  subface_index,
3088  std::integral_constant<bool, integrate>());
3089  }
3090 
3091  /*
3092  * Helper function which applies the 1d kernels for on one
3093  * component in a face. normal_dir indicates the direction of the continuous
3094  * component of the RT space. std::integral_constant<bool, false> is the
3095  * evaluation path, and std::integral_constant<bool, true> below is the
3096  * integration path. These two functions can be fused together since all
3097  * offsets and pointers are the exact same.
3098  */
3099  template <int normal_dir>
3100  static inline void
3102  Number * values_dofs,
3104  Number * scratch_data,
3105  const EvaluationFlags::EvaluationFlags evaluation_flag,
3106  const unsigned int face_direction,
3107  const unsigned int subface_index,
3108  std::integral_constant<bool, false>)
3109  {
3110  using EvalNormal =
3112  dim - 1,
3113  (fe_degree == -1) ? 1 : fe_degree + 1,
3114  n_q_points_1d,
3115  Number,
3116  normal_dir>;
3117  using EvalTangent =
3119  dim - 1,
3120  (fe_degree == -1) ? 1 : fe_degree,
3121  n_q_points_1d,
3122  Number,
3123  normal_dir>;
3124 
3125  using TempEval0 = typename std::
3126  conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
3127  using TempEval1 = typename std::
3128  conditional<normal_dir == 0, EvalTangent, EvalNormal>::type;
3129  using Eval0 = typename std::
3130  conditional<normal_dir == 2, EvalGeneral, TempEval0>::type;
3131  using Eval1 = typename std::
3132  conditional<normal_dir == 2, EvalGeneral, TempEval1>::type;
3133 
3134  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info =
3135  fe_eval.get_shape_info();
3136  Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
3137  ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]),
3138  subface_index,
3139  0);
3140  Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
3141  ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]),
3142  subface_index,
3143  1);
3144 
3145  constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
3146  const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face;
3147  const std::size_t n_dofs_normal =
3148  n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
3149  const std::size_t dofs_stride =
3150  (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
3151  n_dofs_tangent;
3152 
3153  static constexpr ::ndarray<unsigned int, 3, 3> component_table = {
3154  {{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
3155  const unsigned int component =
3156  (dim == 2 && normal_dir == 0 && face_direction == 1) ?
3157  0 :
3158  component_table[face_direction][normal_dir];
3159 
3160  // Initial offsets
3161  values_dofs +=
3162  3 * ((component == 0) ?
3163  0 :
3164  ((component == 1) ?
3165  ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) :
3166  ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent :
3167  n_dofs_normal + n_dofs_tangent)));
3168  const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir;
3169  Number *values_quad = fe_eval.begin_values() + n_q_points * shift;
3170  Number *gradients_quad =
3171  fe_eval.begin_gradients() + dim * n_q_points * shift;
3172  Number *hessians_quad =
3173  fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * shift;
3174 
3175  // Evaluation path
3176  if ((evaluation_flag & EvaluationFlags::values) &&
3177  !(evaluation_flag & EvaluationFlags::gradients))
3178  {
3179  switch (dim)
3180  {
3181  case 3:
3182  eval0.template values<0, true, false>(values_dofs, values_quad);
3183  eval1.template values<1, true, false>(values_quad, values_quad);
3184  break;
3185  case 2:
3186  eval0.template values<0, true, false>(values_dofs, values_quad);
3187  break;
3188  default:
3189  Assert(false, ExcNotImplemented());
3190  }
3191  }
3192  else if (evaluation_flag & EvaluationFlags::gradients)
3193  {
3194  switch (dim)
3195  {
3196  case 3:
3197  // grad x
3198  eval0.template gradients<0, true, false>(values_dofs,
3199  scratch_data);
3200  eval1.template values<1, true, false>(scratch_data,
3201  gradients_quad);
3202 
3203  // grad y
3204  eval0.template values<0, true, false>(values_dofs,
3205  scratch_data);
3206  eval1.template gradients<1, true, false>(scratch_data,
3207  gradients_quad +
3208  n_q_points);
3209 
3210  if (evaluation_flag & EvaluationFlags::values)
3211  eval1.template values<1, true, false>(scratch_data,
3212  values_quad);
3213 
3214  // grad z
3215  eval0.template values<0, true, false>(values_dofs + dofs_stride,
3216  scratch_data);
3217  eval1.template values<1, true, false>(scratch_data,
3218  gradients_quad +
3219  2 * n_q_points);
3220 
3221  break;
3222  case 2:
3223  eval0.template values<0, true, false>(values_dofs + dofs_stride,
3224  gradients_quad +
3225  n_q_points);
3226  eval0.template gradients<0, true, false>(values_dofs,
3227  gradients_quad);
3228  if ((evaluation_flag & EvaluationFlags::values))
3229  eval0.template values<0, true, false>(values_dofs,
3230  values_quad);
3231  break;
3232  default:
3233  AssertThrow(false, ExcNotImplemented());
3234  }
3235  }
3236 
3237  if (evaluation_flag & EvaluationFlags::hessians)
3238  {
3239  switch (dim)
3240  {
3241  case 3:
3242  // grad xx
3243  eval0.template hessians<0, true, false>(values_dofs,
3244  scratch_data);
3245  eval1.template values<1, true, false>(scratch_data,
3246  hessians_quad);
3247 
3248  // grad yy
3249  eval0.template values<0, true, false>(values_dofs,
3250  scratch_data);
3251  eval1.template hessians<1, true, false>(scratch_data,
3252  hessians_quad +
3253  n_q_points);
3254 
3255  // grad zz
3256  eval0.template values<0, true, false>(values_dofs +
3257  2 * dofs_stride,
3258  scratch_data);
3259  eval1.template values<1, true, false>(scratch_data,
3260  hessians_quad +
3261  2 * n_q_points);
3262 
3263  // grad xy
3264  eval0.template gradients<0, true, false>(values_dofs,
3265  scratch_data);
3266  eval1.template gradients<1, true, false>(scratch_data,
3267  hessians_quad +
3268  3 * n_q_points);
3269 
3270  // grad xz
3271  eval0.template gradients<0, true, false>(values_dofs +
3272  dofs_stride,
3273  scratch_data);
3274  eval1.template values<1, true, false>(scratch_data,
3275  hessians_quad +
3276  4 * n_q_points);
3277 
3278  // grad yz
3279  eval0.template values<0, true, false>(values_dofs + dofs_stride,
3280  scratch_data);
3281  eval1.template gradients<1, true, false>(scratch_data,
3282  hessians_quad +
3283  5 * n_q_points);
3284 
3285  break;
3286  case 2:
3287  // grad xx
3288  eval0.template hessians<0, true, false>(values_dofs,
3289  hessians_quad);
3290  // grad yy
3291  eval0.template values<0, true, false>(
3292  values_dofs + 2 * dofs_stride, hessians_quad + n_q_points);
3293  // grad xy
3294  eval0.template gradients<0, true, false>(
3295  values_dofs + dofs_stride, hessians_quad + 2 * n_q_points);
3296  break;
3297  default:
3298  AssertThrow(false, ExcNotImplemented());
3299  }
3300  }
3301  }
3302 
3303  template <int normal_dir>
3304  static inline void
3306  Number * values_dofs,
3308  Number * scratch_data,
3309  const EvaluationFlags::EvaluationFlags evaluation_flag,
3310  const unsigned int face_direction,
3311  const unsigned int subface_index,
3312  std::integral_constant<bool, true>)
3313  {
3314  using EvalNormal =
3316  dim - 1,
3317  (fe_degree == -1) ? 1 : fe_degree + 1,
3318  n_q_points_1d,
3319  Number,
3320  normal_dir>;
3321  using EvalTangent =
3323  dim - 1,
3324  (fe_degree == -1) ? 1 : fe_degree,
3325  n_q_points_1d,
3326  Number,
3327  normal_dir>;
3328 
3329  using TempEval0 = typename std::
3330  conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
3331  using TempEval1 = typename std::
3332  conditional<normal_dir == 0, EvalTangent, EvalNormal>::type;
3333  using Eval0 = typename std::
3334  conditional<normal_dir == 2, EvalGeneral, TempEval0>::type;
3335  using Eval1 = typename std::
3336  conditional<normal_dir == 2, EvalGeneral, TempEval1>::type;
3337 
3338  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info =
3339  fe_eval.get_shape_info();
3340  Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
3341  ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]),
3342  subface_index,
3343  0);
3344  Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
3345  ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]),
3346  subface_index,
3347  1);
3348 
3349  constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
3350  const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face;
3351  const std::size_t n_dofs_normal =
3352  n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
3353  const std::size_t dofs_stride =
3354  (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
3355  n_dofs_tangent;
3356 
3357  static constexpr ::ndarray<unsigned int, 3, 3> component_table = {
3358  {{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
3359  const unsigned int component =
3360  (dim == 2 && normal_dir == 0 && face_direction == 1) ?
3361  0 :
3362  component_table[face_direction][normal_dir];
3363 
3364  // Initial offsets
3365  values_dofs +=
3366  3 * ((component == 0) ?
3367  0 :
3368  ((component == 1) ?
3369  ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) :
3370  ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent :
3371  n_dofs_normal + n_dofs_tangent)));
3372  const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir;
3373  Number *values_quad = fe_eval.begin_values() + n_q_points * shift;
3374  Number *gradients_quad =
3375  fe_eval.begin_gradients() + dim * n_q_points * shift;
3376  Number *hessians_quad =
3377  fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * shift;
3378 
3379  // Integration path
3380  if ((evaluation_flag & EvaluationFlags::values) &&
3381  !(evaluation_flag & EvaluationFlags::gradients))
3382  {
3383  switch (dim)
3384  {
3385  case 3:
3386  eval1.template values<1, false, false>(values_quad,
3387  values_quad);
3388  eval0.template values<0, false, false>(values_quad,
3389  values_dofs);
3390  break;
3391  case 2:
3392  eval0.template values<0, false, false>(values_quad,
3393  values_dofs);
3394  break;
3395  default:
3396  Assert(false, ExcNotImplemented());
3397  }
3398  }
3399  else if (evaluation_flag & EvaluationFlags::gradients)
3400  {
3401  switch (dim)
3402  {
3403  case 3:
3404  // grad z
3405  eval1.template values<1, false, false>(gradients_quad +
3406  2 * n_q_points,
3407  gradients_quad +
3408  2 * n_q_points);
3409  eval0.template values<0, false, false>(
3410  gradients_quad + 2 * n_q_points, values_dofs + dofs_stride);
3411 
3412  if (evaluation_flag & EvaluationFlags::values)
3413  {
3414  eval1.template values<1, false, false>(values_quad,
3415  scratch_data);
3416  eval1.template gradients<1, false, true>(gradients_quad +
3417  n_q_points,
3418  scratch_data);
3419  }
3420  else
3421  eval1.template gradients<1, false, false>(gradients_quad +
3422  n_q_points,
3423  scratch_data);
3424 
3425  // grad y
3426  eval0.template values<0, false, false>(scratch_data,
3427  values_dofs);
3428 
3429  // grad x
3430  eval1.template values<1, false, false>(gradients_quad,
3431  scratch_data);
3432  eval0.template gradients<0, false, true>(scratch_data,
3433  values_dofs);
3434 
3435  break;
3436  case 2:
3437  eval0.template values<0, false, false>(
3438  gradients_quad + n_q_points, values_dofs + dofs_stride);
3439  eval0.template gradients<0, false, false>(gradients_quad,
3440  values_dofs);
3441  if (evaluation_flag & EvaluationFlags::values)
3442  eval0.template values<0, false, true>(values_quad,
3443  values_dofs);
3444  break;
3445  default:
3446  AssertThrow(false, ExcNotImplemented());
3447  }
3448  }
3449 
3450  if (evaluation_flag & EvaluationFlags::hessians)
3451  {
3452  switch (dim)
3453  {
3454  case 3:
3455  // grad xx
3456  eval1.template values<1, false, false>(hessians_quad,
3457  scratch_data);
3458  if ((evaluation_flag &
3460  eval0.template hessians<0, false, true>(scratch_data,
3461  values_dofs);
3462  else
3463  eval0.template hessians<0, false, false>(scratch_data,
3464  values_dofs);
3465 
3466  // grad yy
3467  eval1.template hessians<1, false, false>(hessians_quad +
3468  n_q_points,
3469  scratch_data);
3470  eval0.template values<0, false, true>(scratch_data,
3471  values_dofs);
3472 
3473  // grad zz
3474  eval1.template values<1, false, false>(hessians_quad +
3475  2 * n_q_points,
3476  scratch_data);
3477  eval0.template values<0, false, false>(scratch_data,
3478  values_dofs +
3479  2 * dofs_stride);
3480 
3481  // grad xy
3482  eval1.template gradients<1, false, false>(hessians_quad +
3483  3 * n_q_points,
3484  scratch_data);
3485  eval0.template gradients<0, false, true>(scratch_data,
3486  values_dofs);
3487 
3488  // grad xz
3489  eval1.template values<1, false, false>(hessians_quad +
3490  4 * n_q_points,
3491  scratch_data);
3492  if ((evaluation_flag & EvaluationFlags::gradients))
3493  eval0.template gradients<0, false, true>(scratch_data,
3494  values_dofs +
3495  dofs_stride);
3496  else
3497  eval0.template gradients<0, false, false>(scratch_data,
3498  values_dofs +
3499  dofs_stride);
3500 
3501  // grad yz
3502  eval1.template gradients<1, false, false>(hessians_quad +
3503  5 * n_q_points,
3504  scratch_data);
3505  eval0.template values<0, false, true>(scratch_data,
3506  values_dofs +
3507  dofs_stride);
3508 
3509  break;
3510  case 2:
3511  // grad xx
3512  if (evaluation_flag &
3514  eval0.template hessians<0, false, true>(hessians_quad,
3515  values_dofs);
3516  else
3517  eval0.template hessians<0, false, false>(hessians_quad,
3518  values_dofs);
3519 
3520  // grad yy
3521  eval0.template values<0, false, false>(
3522  hessians_quad + n_q_points, values_dofs + 2 * dofs_stride);
3523  // grad xy
3524  if ((evaluation_flag & EvaluationFlags::gradients))
3525  eval0.template gradients<0, false, true>(
3526  hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
3527  else
3528  eval0.template gradients<0, false, false>(
3529  hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
3530  break;
3531  default:
3532  AssertThrow(false, ExcNotImplemented());
3533  }
3534  }
3535  }
3536  };
3537 
3538 
3539  template <int dim, int fe_degree, typename Number>
3541  {
3542  template <bool do_evaluate, bool add_into_output>
3543  static void
3544  interpolate(const unsigned int n_components,
3546  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
3547  const Number * input,
3548  Number * output,
3549  const unsigned int face_no)
3550  {
3551  Assert(static_cast<unsigned int>(fe_degree) ==
3552  shape_info.data.front().fe_degree ||
3553  fe_degree == -1,
3554  ExcInternalError());
3556  interpolate_generic_raviart_thomas<do_evaluate, add_into_output>(
3557  n_components, input, output, flags, face_no, shape_info);
3558  else
3559  interpolate_generic<do_evaluate, add_into_output>(
3560  n_components,
3561  input,
3562  output,
3563  flags,
3564  face_no,
3565  shape_info.data.front().fe_degree + 1,
3566  shape_info.data.front().shape_data_on_face,
3567  shape_info.dofs_per_component_on_cell,
3568  3 * shape_info.dofs_per_component_on_face);
3569  }
3570 
3574  template <bool do_evaluate, bool add_into_output>
3575  static void
3577  const unsigned int n_components,
3579  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
3580  const Number * input,
3581  Number * output,
3582  const unsigned int face_no)
3583  {
3584  Assert(static_cast<unsigned int>(fe_degree + 1) ==
3585  shape_info.data.front().n_q_points_1d ||
3586  fe_degree == -1,
3587  ExcInternalError());
3588 
3589  interpolate_generic<do_evaluate, add_into_output>(
3590  n_components,
3591  input,
3592  output,
3593  flags,
3594  face_no,
3595  shape_info.data.front().quadrature.size(),
3596  shape_info.data.front().quadrature_data_on_face,
3597  shape_info.n_q_points,
3598  shape_info.n_q_points_face);
3599  }
3600 
3601  private:
3602  template <bool do_evaluate, bool add_into_output, int face_direction = 0>
3603  static void
3604  interpolate_generic(const unsigned int n_components,
3605  const Number * input,
3606  Number * output,
3608  const unsigned int face_no,
3609  const unsigned int n_points_1d,
3610  const std::array<AlignedVector<Number>, 2> &shape_data,
3611  const unsigned int dofs_per_component_on_cell,
3612  const unsigned int dofs_per_component_on_face)
3613  {
3614  if (face_direction == face_no / 2)
3615  {
3617  dim,
3618  fe_degree + 1,
3619  0,
3620  Number>
3621  evalf(shape_data[face_no % 2],
3624  n_points_1d,
3625  0);
3626 
3627  const unsigned int in_stride = do_evaluate ?
3628  dofs_per_component_on_cell :
3629  dofs_per_component_on_face;
3630  const unsigned int out_stride = do_evaluate ?
3631  dofs_per_component_on_face :
3632  dofs_per_component_on_cell;
3633 
3634  for (unsigned int c = 0; c < n_components; ++c)
3635  {
3636  if (flag & EvaluationFlags::hessians)
3637  evalf.template apply_face<face_direction,
3638  do_evaluate,
3639  add_into_output,
3640  2>(input, output);
3641  else if (flag & EvaluationFlags::gradients)
3642  evalf.template apply_face<face_direction,
3643  do_evaluate,
3644  add_into_output,
3645  1>(input, output);
3646  else
3647  evalf.template apply_face<face_direction,
3648  do_evaluate,
3649  add_into_output,
3650  0>(input, output);
3651  input += in_stride;
3652  output += out_stride;
3653  }
3654  }
3655  else if (face_direction < dim)
3656  {
3657  interpolate_generic<do_evaluate,
3658  add_into_output,
3659  std::min(face_direction + 1, dim - 1)>(
3660  n_components,
3661  input,
3662  output,
3663  flag,
3664  face_no,
3665  n_points_1d,
3666  shape_data,
3667  dofs_per_component_on_cell,
3668  dofs_per_component_on_face);
3669  }
3670  }
3671 
3672  template <typename EvalType>
3673  static EvalType
3676  const unsigned int face_no)
3677  {
3678  return EvalType(data.shape_data_on_face[face_no % 2],
3681  }
3682 
3683  template <bool do_evaluate,
3684  bool add_into_output,
3685  int face_direction = 0,
3686  int max_derivative = 0>
3687  static void
3689  const unsigned int n_components,
3690  const Number * input,
3691  Number * output,
3693  const unsigned int face_no,
3694  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info)
3695  {
3696  if (dim == 1)
3697  {
3698  // This should never happen since the FE_RaviartThomasNodal is not
3699  // defined for dim = 1. It prevents compiler warnings of infinite
3700  // recursion.
3701  Assert(false, ExcInternalError());
3702  return;
3703  }
3704 
3705  bool increase_max_der = false;
3706  if ((flag & EvaluationFlags::hessians && max_derivative < 2) ||
3707  (flag & EvaluationFlags::gradients && max_derivative < 1))
3708  increase_max_der = true;
3709 
3710  if (face_direction == face_no / 2 && !increase_max_der)
3711  {
3713  add_into_output,
3714  face_direction,
3715  max_derivative>(
3716  shape_info, face_no, input, output);
3717  }
3718  else if (face_direction == face_no / 2)
3719  {
3720  // Only increase max_derivative
3722  add_into_output,
3723  face_direction,
3724  std::min(max_derivative + 1, 2)>(
3725  n_components, input, output, flag, face_no, shape_info);
3726  }
3727  else if (face_direction < dim)
3728  {
3729  if (increase_max_der)
3730  {
3732  do_evaluate,
3733  add_into_output,
3734  std::min(face_direction + 1, dim - 1),
3735  std::min(max_derivative + 1, 2)>(
3736  n_components, input, output, flag, face_no, shape_info);
3737  }
3738  else
3739  {
3741  add_into_output,
3742  std::min(face_direction + 1,
3743  dim - 1),
3744  max_derivative>(
3745  n_components, input, output, flag, face_no, shape_info);
3746  }
3747  }
3748  }
3749 
3750  /* Help function for interpolate_generic_raviart_thomas */
3751  template <bool do_evaluate,
3752  bool add_into_output,
3753  int face_direction,
3754  int max_derivative>
3755  static inline void
3757  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
3758  const unsigned int face_no,
3759  const Number * input,
3760  Number * output)
3761  {
3762  // These types are evaluators in either normal or tangential direction
3763  // depending on the face direction, with different normal directions for
3764  // the different components.
3765  using Evalf0 = typename std::conditional<
3766  face_direction == 0,
3768  dim,
3769  (fe_degree == -1) ? 1 : fe_degree + 1,
3770  0,
3771  Number,
3772  0>,
3774  dim,
3775  (fe_degree == -1) ? 1 : fe_degree,
3776  0,
3777  Number,
3778  0>>::type;
3779  using Evalf1 = typename std::conditional<
3780  face_direction == 1,
3782  dim,
3783  (fe_degree == -1) ? 1 : fe_degree + 1,
3784  0,
3785  Number,
3786  1>,
3788  dim,
3789  (fe_degree == -1) ? 1 : fe_degree,
3790  0,
3791  Number,
3792  1>>::type;
3793  using Evalf2 = typename std::conditional<
3794  face_direction == 2,
3796  dim,
3797  (fe_degree == -1) ? 1 : fe_degree + 1,
3798  0,
3799  Number,
3800  2>,
3802  dim,
3803  (fe_degree == -1) ? 1 : fe_degree,
3804  0,
3805  Number,
3806  2>>::type;
3807 
3808  Evalf0 evalf0 =
3809  create_evaluator_tensor_product<Evalf0>((face_direction == 0) ?
3810  shape_info.data[0] :
3811  shape_info.data[1],
3812  face_no);
3813  Evalf1 evalf1 =
3814  create_evaluator_tensor_product<Evalf1>((face_direction == 1) ?
3815  shape_info.data[0] :
3816  shape_info.data[1],
3817  face_no);
3818  Evalf2 evalf2 =
3819  create_evaluator_tensor_product<Evalf2>((face_direction == 2) ?
3820  shape_info.data[0] :
3821  shape_info.data[1],
3822  face_no);
3823 
3824  const unsigned int dofs_per_component_on_cell =
3825  shape_info.dofs_per_component_on_cell;
3826  const unsigned int dofs_per_component_on_face =
3827  3 * shape_info.dofs_per_component_on_face;
3828 
3829  // NOTE! dofs_per_component_on_face is in the tangent direction,
3830  // i.e (fe.degree+1)*fe.degree. Normal faces are only
3831  // fe.degree*fe.degree
3832  const unsigned int in_stride =
3833  do_evaluate ? dofs_per_component_on_cell : dofs_per_component_on_face;
3834  const unsigned int out_stride =
3835  do_evaluate ? dofs_per_component_on_face : dofs_per_component_on_cell;
3836 
3837  const unsigned int in_stride_after_normal =
3838  do_evaluate ?
3839  dofs_per_component_on_cell :
3840  dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2);
3841  const unsigned int out_stride_after_normal =
3842  do_evaluate ?
3843  dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2) :
3844  dofs_per_component_on_cell;
3845 
3846  evalf0.template apply_face<face_direction,
3847  do_evaluate,
3848  add_into_output,
3849  max_derivative>(input, output);
3850  // stride to next component
3851  input += (face_direction == 0) ? in_stride_after_normal : in_stride;
3852  output += (face_direction == 0) ? out_stride_after_normal : out_stride;
3853 
3854  evalf1.template apply_face<face_direction,
3855  do_evaluate,
3856  add_into_output,
3857  max_derivative>(input, output);
3858 
3859  if (dim == 3)
3860  {
3861  // stride to next component
3862  input += (face_direction == 1) ? in_stride_after_normal : in_stride;
3863  output +=
3864  (face_direction == 1) ? out_stride_after_normal : out_stride;
3865 
3866  evalf2.template apply_face<face_direction,
3867  do_evaluate,
3868  add_into_output,
3869  max_derivative>(input, output);
3870  }
3871  }
3872  };
3873 
3874 
3875 
3876  // internal helper function for reading data; base version of different types
3877  template <typename VectorizedArrayType, typename Number2>
3878  void
3879  do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
3880  {
3881  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3882  dst[v] = src_ptr[v];
3883  }
3884 
3885 
3886 
3887  // internal helper function for reading data; specialized version where we
3888  // can use a dedicated load function
3889  template <typename Number, std::size_t width>
3890  void
3892  {
3893  dst.load(src_ptr);
3894  }
3895 
3896 
3897 
3898  // internal helper function for reading data; base version of different types
3899  template <typename VectorizedArrayType, typename Number2>
3900  void
3901  do_vectorized_gather(const Number2 * src_ptr,
3902  const unsigned int * indices,
3903  VectorizedArrayType &dst)
3904  {
3905  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3906  dst[v] = src_ptr[indices[v]];
3907  }
3908 
3909 
3910 
3911  // internal helper function for reading data; specialized version where we
3912  // can use a dedicated gather function
3913  template <typename Number, std::size_t width>
3914  void
3915  do_vectorized_gather(const Number * src_ptr,
3916  const unsigned int * indices,
3918  {
3919  dst.gather(src_ptr, indices);
3920  }
3921 
3922 
3923 
3924  // internal helper function for reading data; base version of different types
3925  template <typename VectorizedArrayType, typename Number2>
3926  void
3927  do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
3928  {
3929  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3930  dst_ptr[v] += src[v];
3931  }
3932 
3933 
3934 
3935  // internal helper function for reading data; specialized version where we
3936  // can use a dedicated load function
3937  template <typename Number, std::size_t width>
3938  void
3940  {
3942  tmp.load(dst_ptr);
3943  (tmp + src).store(dst_ptr);
3944  }
3945 
3946 
3947 
3948  // internal helper function for reading data; base version of different types
3949  template <typename VectorizedArrayType, typename Number2>
3950  void
3951  do_vectorized_scatter_add(const VectorizedArrayType src,
3952  const unsigned int * indices,
3953  Number2 * dst_ptr)
3954  {
3955  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3956  dst_ptr[indices[v]] += src[v];
3957  }
3958 
3959 
3960 
3961  // internal helper function for reading data; specialized version where we
3962  // can use a dedicated gather function
3963  template <typename Number, std::size_t width>
3964  void
3966  const unsigned int * indices,
3967  Number * dst_ptr)
3968  {
3969 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
3970  for (unsigned int v = 0; v < width; ++v)
3971  dst_ptr[indices[v]] += src[v];
3972 #else
3974  tmp.gather(dst_ptr, indices);
3975  (tmp + src).scatter(indices, dst_ptr);
3976 #endif
3977  }
3978 
3979 
3980 
3981  template <typename Number>
3982  void
3983  adjust_for_face_orientation(const unsigned int dim,
3984  const unsigned int n_components,
3986  const unsigned int *orientation,
3987  const bool integrate,
3988  const std::size_t n_q_points,
3989  Number * tmp_values,
3990  Number * values_quad,
3991  Number * gradients_quad,
3992  Number * hessians_quad)
3993  {
3994  for (unsigned int c = 0; c < n_components; ++c)
3995  {
3996  if (flag & EvaluationFlags::values)
3997  {
3998  if (integrate)
3999  for (unsigned int q = 0; q < n_q_points; ++q)
4000  tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
4001  else
4002  for (unsigned int q = 0; q < n_q_points; ++q)
4003  tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
4004  for (unsigned int q = 0; q < n_q_points; ++q)
4005  values_quad[c * n_q_points + q] = tmp_values[q];
4006  }
4007  if (flag & EvaluationFlags::gradients)
4008  for (unsigned int d = 0; d < dim; ++d)
4009  {
4010  if (integrate)
4011  for (unsigned int q = 0; q < n_q_points; ++q)
4012  tmp_values[q] =
4013  gradients_quad[(c * dim + d) * n_q_points + orientation[q]];
4014  else
4015  for (unsigned int q = 0; q < n_q_points; ++q)
4016  tmp_values[orientation[q]] =
4017  gradients_quad[(c * dim + d) * n_q_points + q];
4018  for (unsigned int q = 0; q < n_q_points; ++q)
4019  gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q];
4020  }
4021  if (flag & EvaluationFlags::hessians)
4022  {
4023  const unsigned int hdim = (dim * (dim + 1)) / 2;
4024  for (unsigned int d = 0; d < hdim; ++d)
4025  {
4026  if (integrate)
4027  for (unsigned int q = 0; q < n_q_points; ++q)
4028  tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
4029  orientation[q]];
4030  else
4031  for (unsigned int q = 0; q < n_q_points; ++q)
4032  tmp_values[orientation[q]] =
4033  hessians_quad[(c * hdim + d) * n_q_points + q];
4034  for (unsigned int q = 0; q < n_q_points; ++q)
4035  hessians_quad[(c * hdim + d) * n_q_points + q] =
4036  tmp_values[q];
4037  }
4038  }
4039  }
4040  }
4041 
4042 
4043 
4044  template <typename Number, typename VectorizedArrayType>
4045  void
4047  const unsigned int dim,
4048  const unsigned int n_components,
4049  const unsigned int v,
4051  const unsigned int * orientation,
4052  const bool integrate,
4053  const std::size_t n_q_points,
4054  Number * tmp_values,
4055  VectorizedArrayType * values_quad,
4056  VectorizedArrayType * gradients_quad = nullptr,
4057  VectorizedArrayType * hessians_quad = nullptr)
4058  {
4059  for (unsigned int c = 0; c < n_components; ++c)
4060  {
4061  if (flag & EvaluationFlags::values)
4062  {
4063  if (integrate)
4064  for (unsigned int q = 0; q < n_q_points; ++q)
4065  tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v];
4066  else
4067  for (unsigned int q = 0; q < n_q_points; ++q)
4068  tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v];
4069  for (unsigned int q = 0; q < n_q_points; ++q)
4070  values_quad[c * n_q_points + q][v] = tmp_values[q];
4071  }
4072  if (flag & EvaluationFlags::gradients)
4073  for (unsigned int d = 0; d < dim; ++d)
4074  {
4075  Assert(gradients_quad != nullptr, ExcInternalError());
4076  if (integrate)
4077  for (unsigned int q = 0; q < n_q_points; ++q)
4078  tmp_values[q] = gradients_quad[(c * dim + d) * n_q_points +
4079  orientation[q]][v];
4080  else
4081  for (unsigned int q = 0; q < n_q_points; ++q)
4082  tmp_values[orientation[q]] =
4083  gradients_quad[(c * dim + d) * n_q_points + q][v];
4084  for (unsigned int q = 0; q < n_q_points; ++q)
4085  gradients_quad[(c * dim + d) * n_q_points + q][v] =
4086  tmp_values[q];
4087  }
4088  if (flag & EvaluationFlags::hessians)
4089  {
4090  Assert(hessians_quad != nullptr, ExcInternalError());
4091  const unsigned int hdim = (dim * (dim + 1)) / 2;
4092  for (unsigned int d = 0; d < hdim; ++d)
4093  {
4094  if (integrate)
4095  for (unsigned int q = 0; q < n_q_points; ++q)
4096  tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
4097  orientation[q]][v];
4098  else
4099  for (unsigned int q = 0; q < n_q_points; ++q)
4100  tmp_values[orientation[q]] =
4101  hessians_quad[(c * hdim + d) * n_q_points + q][v];
4102  for (unsigned int q = 0; q < n_q_points; ++q)
4103  hessians_quad[(c * hdim + d) * n_q_points + q][v] =
4104  tmp_values[q];
4105  }
4106  }
4107  }
4108  }
4109 
4110 
4111 
4112  template <int dim, typename Number>
4114  {
4115  template <int fe_degree, int n_q_points_1d>
4116  static bool
4117  run(const unsigned int n_components,
4118  const EvaluationFlags::EvaluationFlags evaluation_flag,
4119  const Number * values_dofs,
4121  {
4122  const auto &shape_info = fe_eval.get_shape_info();
4123  const auto &shape_data = shape_info.data.front();
4124 
4125  if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
4126  {
4127  Assert((fe_eval.get_dof_access_index() ==
4129  fe_eval.is_interior_face() == false) == false,
4130  ExcNotImplemented());
4131 
4132  const unsigned int face_no = fe_eval.get_face_no();
4133  const unsigned int face_orientation = fe_eval.get_face_orientation();
4134  const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
4135  const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
4136 
4137  using Eval =
4139 
4140  if (evaluation_flag & EvaluationFlags::values)
4141  {
4142  const auto shape_values =
4143  &shape_data.shape_values_face(face_no, face_orientation, 0);
4144 
4145  auto values_quad_ptr = fe_eval.begin_values();
4146  auto values_dofs_actual_ptr = values_dofs;
4147 
4148  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
4149  for (unsigned int c = 0; c < n_components; ++c)
4150  {
4151  eval.template values<0, true, false>(values_dofs_actual_ptr,
4152  values_quad_ptr);
4153 
4154  values_quad_ptr += n_q_points;
4155  values_dofs_actual_ptr += n_dofs;
4156  }
4157  }
4158 
4159  if (evaluation_flag & EvaluationFlags::gradients)
4160  {
4161  auto gradients_quad_ptr = fe_eval.begin_gradients();
4162  auto values_dofs_actual_ptr = values_dofs;
4163 
4164  std::array<const Number *, dim> shape_gradients;
4165  for (unsigned int d = 0; d < dim; ++d)
4166  shape_gradients[d] = &shape_data.shape_gradients_face(
4167  face_no, face_orientation, d, 0);
4168 
4169  for (unsigned int c = 0; c < n_components; ++c)
4170  {
4171  for (unsigned int d = 0; d < dim; ++d)
4172  {
4173  Eval eval(nullptr,
4174  shape_gradients[d],
4175  nullptr,
4176  n_dofs,
4177  n_q_points);
4178 
4179  eval.template gradients<0, true, false>(
4180  values_dofs_actual_ptr, gradients_quad_ptr);
4181 
4182  gradients_quad_ptr += n_q_points;
4183  }
4184  values_dofs_actual_ptr += n_dofs;
4185  }
4186  }
4187 
4188  Assert(!(evaluation_flag & EvaluationFlags::hessians),
4189  ExcNotImplemented());
4190 
4191  return true;
4192  }
4193 
4194  const unsigned int dofs_per_face =
4195  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
4196  Utilities::pow(shape_data.fe_degree + 1, dim - 1);
4197 
4198  // Note: we always keep storage of values, 1st and 2nd derivatives in an
4199  // array, so reserve space for all three here
4200  Number *temp = fe_eval.get_scratch_data().begin();
4201  Number *scratch_data = temp + 3 * n_components * dofs_per_face;
4202 
4203  bool use_vectorization = true;
4204 
4205  if (fe_eval.get_dof_access_index() ==
4207  fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
4208  for (unsigned int v = 0; v < Number::size(); ++v)
4209  if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
4210  fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
4211  use_vectorization = false;
4212 
4213  if (use_vectorization == false)
4214  {
4215  for (unsigned int v = 0; v < Number::size(); ++v)
4216  {
4217  // the loop breaks once an invalid_unsigned_int is hit for
4218  // all cases except the exterior faces in the ECL loop (where
4219  // some faces might be at the boundaries but others not)
4220  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4221  {
4222  for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4223  ++i)
4224  temp[i][v] = 0;
4225  continue;
4226  }
4227 
4229  template interpolate<true, false>(n_components,
4230  evaluation_flag,
4231  shape_info,
4232  values_dofs,
4233  scratch_data,
4234  fe_eval.get_face_no(v));
4235 
4236  for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4237  ++i)
4238  temp[i][v] = scratch_data[i][v];
4239  }
4240  }
4241  else
4243  template interpolate<true, false>(n_components,
4244  evaluation_flag,
4245  shape_info,
4246  values_dofs,
4247  temp,
4248  fe_eval.get_face_no());
4249 
4250  const unsigned int subface_index = fe_eval.get_subface_index();
4251  constexpr unsigned int n_q_points_1d_actual =
4252  fe_degree > -1 ? n_q_points_1d : 0;
4253 
4254  if (fe_degree >= 1 &&
4255  shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
4256  {
4258  (fe_degree == -1) ? 1 : fe_degree,
4259  (n_q_points_1d < 1) ? 1 :
4260  n_q_points_1d,
4261  Number>::
4262  template evaluate_or_integrate_in_face<false>(
4263  evaluation_flag,
4264  temp,
4265  fe_eval,
4266  scratch_data,
4267  subface_index,
4268  fe_eval.get_face_no());
4269  }
4270  else if (fe_degree > -1 &&
4271  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
4272  shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
4273  FEFaceEvaluationImpl<true,
4274  dim,
4275  fe_degree,
4276  n_q_points_1d_actual,
4277  Number>::evaluate_in_face(n_components,
4278  evaluation_flag,
4279  shape_data,
4280  temp,
4281  fe_eval.begin_values(),
4282  fe_eval
4283  .begin_gradients(),
4284  fe_eval.begin_hessians(),
4285  scratch_data,
4286  subface_index);
4287  else
4288  FEFaceEvaluationImpl<false,
4289  dim,
4290  fe_degree,
4291  n_q_points_1d_actual,
4292  Number>::evaluate_in_face(n_components,
4293  evaluation_flag,
4294  shape_data,
4295  temp,
4296  fe_eval.begin_values(),
4297  fe_eval
4298  .begin_gradients(),
4299  fe_eval.begin_hessians(),
4300  scratch_data,
4301  subface_index);
4302 
4303  if (use_vectorization == false)
4304  {
4305  for (unsigned int v = 0; v < Number::size(); ++v)
4306  {
4307  // the loop breaks once an invalid_unsigned_int is hit for
4308  // all cases except the exterior faces in the ECL loop (where
4309  // some faces might be at the boundaries but others not)
4310  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4311  continue;
4312 
4313  if (fe_eval.get_face_orientation(v) != 0)
4315  dim,
4316  n_components,
4317  v,
4318  evaluation_flag,
4320  fe_eval.get_face_orientation(v), 0),
4321  false,
4322  shape_info.n_q_points_face,
4323  &temp[0][0],
4324  fe_eval.begin_values(),
4325  fe_eval.begin_gradients(),
4326  fe_eval.begin_hessians());
4327  }
4328  }
4329  else if (fe_eval.get_face_orientation() != 0)
4331  dim,
4332  n_components,
4333  evaluation_flag,
4335  fe_eval.get_face_orientation(), 0),
4336  false,
4337  shape_info.n_q_points_face,
4338  temp,
4339  fe_eval.begin_values(),
4340  fe_eval.begin_gradients(),
4341  fe_eval.begin_hessians());
4342 
4343  return false;
4344  }
4345  };
4346 
4347 
4348 
4349  template <int dim, typename Number>
4351  {
4352  template <int fe_degree, int n_q_points_1d>
4353  static bool
4354  run(const unsigned int n_components,
4355  const EvaluationFlags::EvaluationFlags integration_flag,
4356  Number * values_dofs,
4358  {
4359  const auto &shape_info = fe_eval.get_shape_info();
4360  const auto &shape_data = shape_info.data.front();
4361 
4362  if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
4363  {
4364  Assert((fe_eval.get_dof_access_index() ==
4366  fe_eval.is_interior_face() == false) == false,
4367  ExcNotImplemented());
4368 
4369  const unsigned int face_no = fe_eval.get_face_no();
4370  const unsigned int face_orientation = fe_eval.get_face_orientation();
4371  const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
4372  const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
4373 
4374  using Eval =
4376 
4377  if (integration_flag & EvaluationFlags::values)
4378  {
4379  const auto shape_values =
4380  &shape_data.shape_values_face(face_no, face_orientation, 0);
4381 
4382  auto values_quad_ptr = fe_eval.begin_values();
4383  auto values_dofs_actual_ptr = values_dofs;
4384 
4385  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
4386  for (unsigned int c = 0; c < n_components; ++c)
4387  {
4388  eval.template values<0, false, false>(values_quad_ptr,
4389  values_dofs_actual_ptr);
4390 
4391  values_quad_ptr += n_q_points;
4392  values_dofs_actual_ptr += n_dofs;
4393  }
4394  }
4395 
4396  if (integration_flag & EvaluationFlags::gradients)
4397  {
4398  auto gradients_quad_ptr = fe_eval.begin_gradients();
4399  auto values_dofs_actual_ptr = values_dofs;
4400 
4401  std::array<const Number *, dim> shape_gradients;
4402  for (unsigned int d = 0; d < dim; ++d)
4403  shape_gradients[d] = &shape_data.shape_gradients_face(
4404  face_no, face_orientation, d, 0);
4405 
4406  for (unsigned int c = 0; c < n_components; ++c)
4407  {
4408  for (unsigned int d = 0; d < dim; ++d)
4409  {
4410  Eval eval(nullptr,
4411  shape_gradients[d],
4412  nullptr,
4413  n_dofs,
4414  n_q_points);
4415 
4416  if (!(integration_flag & EvaluationFlags::values) &&
4417  d == 0)
4418  eval.template gradients<0, false, false>(
4419  gradients_quad_ptr, values_dofs_actual_ptr);
4420  else
4421  eval.template gradients<0, false, true>(
4422  gradients_quad_ptr, values_dofs_actual_ptr);
4423 
4424  gradients_quad_ptr += n_q_points;
4425  }
4426  values_dofs_actual_ptr += n_dofs;
4427  }
4428  }
4429 
4430  Assert(!(integration_flag & EvaluationFlags::hessians),
4431  ExcNotImplemented());
4432 
4433  return true;
4434  }
4435 
4436  const unsigned int dofs_per_face =
4437  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
4438  Utilities::pow(shape_data.fe_degree + 1, dim - 1);
4439 
4440  Number *temp = fe_eval.get_scratch_data().begin();
4441  Number *scratch_data = temp + 3 * n_components * dofs_per_face;
4442 
4443  bool use_vectorization = true;
4444 
4445  if (fe_eval.get_dof_access_index() ==
4447  fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
4448  use_vectorization =
4449  fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
4450  std::all_of(fe_eval.get_cell_ids().begin() + 1,
4451  fe_eval.get_cell_ids().end(),
4452  [&](const auto &v) {
4453  return v == fe_eval.get_cell_ids()[0] ||
4454  v == numbers::invalid_unsigned_int;
4455  });
4456 
4457  if (use_vectorization == false)
4458  {
4459  for (unsigned int v = 0; v < Number::size(); ++v)
4460  {
4461  // the loop breaks once an invalid_unsigned_int is hit for
4462  // all cases except the exterior faces in the ECL loop (where
4463  // some faces might be at the boundaries but others not)
4464  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4465  continue;
4466 
4467  if (fe_eval.get_face_orientation(v) != 0)
4469  dim,
4470  n_components,
4471  v,
4472  integration_flag,
4474  fe_eval.get_face_orientation(v), 0),
4475  true,
4476  shape_info.n_q_points_face,
4477  &temp[0][0],
4478  fe_eval.begin_values(),
4479  fe_eval.begin_gradients(),
4480  fe_eval.begin_hessians());
4481  }
4482  }
4483  else if (fe_eval.get_face_orientation() != 0)
4485  dim,
4486  n_components,
4487  integration_flag,
4489  fe_eval.get_face_orientation(), 0),
4490  true,
4491  shape_info.n_q_points_face,
4492  temp,
4493  fe_eval.begin_values(),
4494  fe_eval.begin_gradients(),
4495  fe_eval.begin_hessians());
4496 
4497  const unsigned int n_q_points_1d_actual =
4498  fe_degree > -1 ? n_q_points_1d : 0;
4499  const unsigned int subface_index = fe_eval.get_subface_index();
4500 
4501  if (fe_degree >= 1 &&
4502  shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
4503  {
4505  (fe_degree == -1) ? 1 : fe_degree,
4506  (n_q_points_1d < 1) ? 1 :
4507  n_q_points_1d,
4508  Number>::
4509  template evaluate_or_integrate_in_face<true>(integration_flag,
4510  temp,
4511  fe_eval,
4512  scratch_data,
4513  subface_index,
4514  fe_eval.get_face_no());
4515  }
4516  else if (fe_degree > -1 &&
4517  fe_eval.get_subface_index() >=
4519  shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
4521  true,
4522  dim,
4523  fe_degree,
4524  n_q_points_1d_actual,
4525  Number>::integrate_in_face(n_components,
4526  integration_flag,
4527  shape_data,
4528  temp,
4529  fe_eval.begin_values(),
4530  fe_eval.begin_gradients(),
4531  fe_eval.begin_hessians(),
4532  scratch_data,
4533  subface_index);
4534  else
4536  false,
4537  dim,
4538  fe_degree,
4539  n_q_points_1d_actual,
4540  Number>::integrate_in_face(n_components,
4541  integration_flag,
4542  shape_data,
4543  temp,
4544  fe_eval.begin_values(),
4545  fe_eval.begin_gradients(),
4546  fe_eval.begin_hessians(),
4547  scratch_data,
4548  subface_index);
4549 
4550  if (use_vectorization == false)
4551  {
4552  for (unsigned int v = 0; v < Number::size(); ++v)
4553  {
4554  // the loop breaks once an invalid_unsigned_int is hit for
4555  // all cases except the exterior faces in the ECL loop (where
4556  // some faces might be at the boundaries but others not)
4557  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4558  continue;
4559 
4561  template interpolate<false, false>(n_components,
4562  integration_flag,
4563  shape_info,
4564  values_dofs,
4565  scratch_data,
4566  fe_eval.get_face_no(v));
4567 
4568  for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4569  ++i)
4570  temp[i][v] = scratch_data[i][v];
4571  }
4572  }
4573  else
4575  template interpolate<false, false>(n_components,
4576  integration_flag,
4577  shape_info,
4578  temp,
4579  values_dofs,
4580  fe_eval.get_face_no());
4581  return false;
4582  }
4583  };
4584 
4585 
4586 
4587  template <int n_face_orientations,
4588  typename Processor,
4589  typename EvaluationData,
4590  const bool check_face_orientations = false>
4591  void
4593  Processor & proc,
4594  const unsigned int n_components,
4595  const EvaluationFlags::EvaluationFlags evaluation_flag,
4596  typename Processor::Number2_ * global_vector_ptr,
4597  const std::vector<ArrayView<const typename Processor::Number2_>> *sm_ptr,
4598  const EvaluationData & fe_eval,
4599  typename Processor::VectorizedArrayType_ * temp1)
4600  {
4601  constexpr int dim = Processor::dim_;
4602  constexpr int fe_degree = Processor::fe_degree_;
4603  using VectorizedArrayType = typename Processor::VectorizedArrayType_;
4604  constexpr int n_lanes = VectorizedArrayType::size();
4605 
4606  using Number = typename Processor::Number_;
4607  using Number2_ = typename Processor::Number2_;
4608 
4609  const auto & shape_data = fe_eval.get_shape_info().data.front();
4610  constexpr bool integrate = Processor::do_integrate;
4611  const unsigned int face_no = fe_eval.get_face_no();
4612  const auto & dof_info = fe_eval.get_dof_info();
4613  const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
4614  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index =
4615  fe_eval.get_dof_access_index();
4616  AssertIndexRange(cell,
4617  dof_info.index_storage_variants[dof_access_index].size());
4618  constexpr unsigned int dofs_per_face =
4619  Utilities::pow(fe_degree + 1, dim - 1);
4620  const unsigned int subface_index = fe_eval.get_subface_index();
4621 
4622  const unsigned int n_filled_lanes =
4623  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
4624 
4625  bool all_faces_are_same = n_filled_lanes == n_lanes;
4626  if (n_face_orientations == n_lanes)
4627  for (unsigned int v = 1; v < n_lanes; ++v)
4628  if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) ||
4629  fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0))
4630  {
4631  all_faces_are_same = false;
4632  break;
4633  }
4634 
4635  // check for re-orientation ...
4636  std::array<const unsigned int *, n_face_orientations> orientation = {};
4637 
4638  if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same &&
4639  fe_eval.is_interior_face() == 0)
4640  for (unsigned int v = 0; v < n_lanes; ++v)
4641  {
4642  // the loop breaks once an invalid_unsigned_int is hit for
4643  // all cases except the exterior faces in the ECL loop (where
4644  // some faces might be at the boundaries but others not)
4645  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4646  continue;
4647 
4648  if (shape_data.nodal_at_cell_boundaries &&
4649  fe_eval.get_face_orientation(v) != 0)
4650  {
4651  // ... and in case we detect a re-orientation, go to the other
4652  // version of this function that actually allows for this
4653  if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
4654  check_face_orientations == false)
4655  {
4656  fe_face_evaluation_process_and_io<n_face_orientations,
4657  Processor,
4658  EvaluationData,
4659  true>(proc,
4660  n_components,
4661  evaluation_flag,
4662  global_vector_ptr,
4663  sm_ptr,
4664  fe_eval,
4665  temp1);
4666  return;
4667  }
4668  orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
4669  fe_eval.get_face_orientation(v), 0);
4670  }
4671  }
4672  else if (dim == 3 && fe_eval.get_face_orientation() != 0)
4673  {
4674  // go to the other version of this function
4675  if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
4676  check_face_orientations == false)
4677  {
4678  fe_face_evaluation_process_and_io<n_face_orientations,
4679  Processor,
4680  EvaluationData,
4681  true>(proc,
4682  n_components,
4683  evaluation_flag,
4684  global_vector_ptr,
4685  sm_ptr,
4686  fe_eval,
4687  temp1);
4688  return;
4689  }
4690  for (unsigned int v = 0; v < n_face_orientations; ++v)
4691  orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
4692  fe_eval.get_face_orientation(), 0);
4693  }
4694 
4695  // we know that the gradient weights for the Hermite case on the
4696  // right (side==1) are the negative from the value at the left
4697  // (side==0), so we only read out one of them.
4698  VectorizedArrayType grad_weight =
4699  shape_data
4700  .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) :
4701  (1 + face_no % 2))];
4702 
4703  // face_to_cell_index_hermite
4704  std::array<const unsigned int *, n_face_orientations> index_array_hermite =
4705  {};
4706  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4707  {
4708  if (n_face_orientations == 1)
4709  index_array_hermite[0] =
4710  &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0);
4711  else
4712  {
4713  for (unsigned int v = 0; v < n_lanes; ++v)
4714  {
4715  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4716  continue;
4717 
4718  const auto face_no = fe_eval.get_face_no(v);
4719 
4720  grad_weight[v] =
4721  shape_data.shape_data_on_face[0][fe_degree +
4722  (integrate ?
4723  (2 - (face_no % 2)) :
4724  (1 + (face_no % 2)))][0];
4725 
4726  index_array_hermite[v] =
4727  &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no,
4728  0);
4729  }
4730  }
4731  }
4732 
4733  // face_to_cell_index_nodal
4734  std::array<const unsigned int *, n_face_orientations> index_array_nodal =
4735  {};
4736  if (shape_data.nodal_at_cell_boundaries == true)
4737  {
4738  if (n_face_orientations == 1)
4739  index_array_nodal[0] =
4740  &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0);
4741  else
4742  {
4743  for (unsigned int v = 0; v < n_lanes; ++v)
4744  {
4745  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4746  continue;
4747 
4748  const auto face_no = fe_eval.get_face_no(v);
4749 
4750  index_array_nodal[v] =
4751  &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no,
4752  0);
4753  }
4754  }
4755  }
4756 
4757 
4758  const auto reorientate = [&](const unsigned int v, const unsigned int i) {
4759  return (!check_face_orientations || orientation[v] == nullptr) ?
4760  i :
4761  orientation[v][i];
4762  };
4763 
4764  const unsigned int cell_index =
4765  dof_access_index == MatrixFreeFunctions::DoFInfo::dof_access_cell ?
4766  fe_eval.get_cell_ids()[0] :
4767  cell * n_lanes;
4768  const unsigned int *dof_indices =
4769  &dof_info.dof_indices_contiguous[dof_access_index][cell_index];
4770 
4771  for (unsigned int comp = 0; comp < n_components; ++comp)
4772  {
4773  const std::size_t index_offset =
4774  dof_info.component_dof_indices_offset
4775  [fe_eval.get_active_fe_index()]
4776  [fe_eval.get_first_selected_component()] +
4777  comp * Utilities::pow(fe_degree + 1, dim);
4778 
4779  // case 1: contiguous and interleaved indices
4780  if (n_face_orientations == 1 &&
4781  dof_info.index_storage_variants[dof_access_index][cell] ==
4783  interleaved_contiguous)
4784  {
4786  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
4787  n_lanes);
4788  Number2_ *vector_ptr =
4789  global_vector_ptr + dof_indices[0] + index_offset * n_lanes;
4790 
4791  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4792  {
4793  for (unsigned int i = 0; i < dofs_per_face; ++i)
4794  {
4795  Assert(n_face_orientations == 1, ExcNotImplemented());
4796 
4797  const unsigned int ind1 = index_array_hermite[0][2 * i];
4798  const unsigned int ind2 = index_array_hermite[0][2 * i + 1];
4799  const unsigned int i_ = reorientate(0, i);
4800  proc.hermite_grad_vectorized(temp1[i_],
4801  temp1[i_ + dofs_per_face],
4802  vector_ptr + ind1 * n_lanes,
4803  vector_ptr + ind2 * n_lanes,
4804  grad_weight);
4805  }
4806  }
4807  else
4808  {
4809  for (unsigned int i = 0; i < dofs_per_face; ++i)
4810  {
4811  Assert(n_face_orientations == 1, ExcNotImplemented());
4812 
4813  const unsigned int i_ = reorientate(0, i);
4814  const unsigned int ind = index_array_nodal[0][i];
4815  proc.value_vectorized(temp1[i_],
4816  vector_ptr + ind * n_lanes);
4817  }
4818  }
4819  }
4820 
4821  // case 2: contiguous and interleaved indices with fixed stride
4822  else if (n_face_orientations == 1 &&
4823  dof_info.index_storage_variants[dof_access_index][cell] ==
4825  interleaved_contiguous_strided)
4826  {
4828  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
4829  n_lanes);
4830  Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes;
4831  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4832  {
4833  for (unsigned int i = 0; i < dofs_per_face; ++i)
4834  {
4835  Assert(n_face_orientations == 1, ExcNotImplemented());
4836 
4837  const unsigned int i_ = reorientate(0, i);
4838  const unsigned int ind1 =
4839  index_array_hermite[0][2 * i] * n_lanes;
4840  const unsigned int ind2 =
4841  index_array_hermite[0][2 * i + 1] * n_lanes;
4842  proc.hermite_grad_vectorized_indexed(
4843  temp1[i_],
4844  temp1[i_ + dofs_per_face],
4845  vector_ptr + ind1,
4846  vector_ptr + ind2,
4847  grad_weight,
4848  dof_indices,
4849  dof_indices);
4850  }
4851  }
4852  else
4853  {
4854  for (unsigned int i = 0; i < dofs_per_face; ++i)
4855  {
4856  Assert(n_face_orientations == 1, ExcNotImplemented());
4857 
4858  const unsigned int i_ = reorientate(0, i);
4859  const unsigned int ind = index_array_nodal[0][i] * n_lanes;
4860  proc.value_vectorized_indexed(temp1[i_],
4861  vector_ptr + ind,
4862  dof_indices);
4863  }
4864  }
4865  }
4866 
4867  // case 3: contiguous and interleaved indices with mixed stride
4868  else if (n_face_orientations == 1 &&
4869  dof_info.index_storage_variants[dof_access_index][cell] ==
4871  interleaved_contiguous_mixed_strides)
4872  {
4873  const unsigned int *strides =
4874  &dof_info.dof_indices_interleave_strides[dof_access_index]
4875  [cell * n_lanes];
4876  unsigned int indices[n_lanes];
4877  for (unsigned int v = 0; v < n_lanes; ++v)
4878  indices[v] = dof_indices[v] + index_offset * strides[v];
4879  const unsigned int n_filled_lanes =
4880  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
4881 
4882  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4883  {
4884  if (n_filled_lanes == n_lanes)
4885  for (unsigned int i = 0; i < dofs_per_face; ++i)
4886  {
4887  Assert(n_face_orientations == 1, ExcNotImplemented());
4888 
4889  const unsigned int i_ = reorientate(0, i);
4890  unsigned int ind1[n_lanes];
4892  for (unsigned int v = 0; v < n_lanes; ++v)
4893  ind1[v] = indices[v] +
4894  index_array_hermite[0][2 * i] * strides[v];
4895  unsigned int ind2[n_lanes];
4897  for (unsigned int v = 0; v < n_lanes; ++v)
4898  ind2[v] =
4899  indices[v] +
4900  // TODO
4901  index_array_hermite[0][2 * i + 1] * strides[v];
4902  proc.hermite_grad_vectorized_indexed(
4903  temp1[i_],
4904  temp1[i_ + dofs_per_face],
4905  global_vector_ptr,
4906  global_vector_ptr,
4907  grad_weight,
4908  ind1,
4909  ind2);
4910  }
4911  else
4912  {
4913  if (integrate == false)
4914  for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
4915  temp1[i] = VectorizedArrayType();
4916 
4917  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4918  for (unsigned int i = 0; i < dofs_per_face; ++i)
4919  {
4920  const unsigned int i_ =
4921  reorientate(n_face_orientations == 1 ? 0 : v, i);
4922  proc.hermite_grad(
4923  temp1[i_][v],
4924  temp1[i_ + dofs_per_face][v],
4925  global_vector_ptr
4926  [indices[v] +
4927  index_array_hermite
4928  [n_face_orientations == 1 ? 0 : v][2 * i] *
4929  strides[v]],
4930  global_vector_ptr
4931  [indices[v] +
4932  index_array_hermite[n_face_orientations == 1 ?
4933  0 :
4934  v][2 * i + 1] *
4935  strides[v]],
4936  grad_weight[n_face_orientations == 1 ? 0 : v]);
4937  }
4938  }
4939  }
4940  else
4941  {
4942  if (n_filled_lanes == n_lanes)
4943  for (unsigned int i = 0; i < dofs_per_face; ++i)
4944  {
4945  Assert(n_face_orientations == 1, ExcInternalError());
4946  unsigned int ind[n_lanes];
4948  for (unsigned int v = 0; v < n_lanes; ++v)
4949  ind[v] =
4950  indices[v] + index_array_nodal[0][i] * strides[v];
4951  const unsigned int i_ = reorientate(0, i);
4952  proc.value_vectorized_indexed(temp1[i_],
4953  global_vector_ptr,
4954  ind);
4955  }
4956  else
4957  {
4958  if (integrate == false)
4959  for (unsigned int i = 0; i < dofs_per_face; ++i)
4960  temp1[i] = VectorizedArrayType();
4961 
4962  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4963  for (unsigned int i = 0; i < dofs_per_face; ++i)
4964  proc.value(
4965  temp1[reorientate(n_face_orientations == 1 ? 0 : v,
4966  i)][v],
4967  global_vector_ptr
4968  [indices[v] +
4969  index_array_nodal[n_face_orientations == 1 ? 0 : v]
4970  [i] *
4971  strides[v]]);
4972  }
4973  }
4974  }
4975 
4976  // case 4: contiguous indices without interleaving
4977  else if (n_face_orientations > 1 ||
4978  dof_info.index_storage_variants[dof_access_index][cell] ==
4980  contiguous)
4981  {
4982  Number2_ *vector_ptr = global_vector_ptr + index_offset;
4983 
4984  const bool vectorization_possible =
4985  all_faces_are_same && (sm_ptr == nullptr);
4986 
4987  std::array<Number2_ *, n_lanes> vector_ptrs;
4988  std::array<unsigned int, n_lanes> reordered_indices;
4989 
4990  if (vectorization_possible == false)
4991  {
4992  vector_ptrs = {};
4993  if (n_face_orientations == 1)
4994  {
4995  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4996  if (sm_ptr == nullptr)
4997  {
4998  vector_ptrs[v] = vector_ptr + dof_indices[v];
4999  }
5000  else
5001  {
5002  const auto &temp =
5003  dof_info
5004  .dof_indices_contiguous_sm[dof_access_index]
5005  [cell * n_lanes + v];
5006  vector_ptrs[v] = const_cast<Number2_ *>(
5007  sm_ptr->operator[](temp.first).data() +
5008  temp.second + index_offset);
5009  }
5010  }
5011  else if (n_face_orientations == n_lanes)
5012  {
5013  const auto &cells = fe_eval.get_cell_ids();
5014  for (unsigned int v = 0; v < n_lanes; ++v)
5015  if (cells[v] != numbers::invalid_unsigned_int)
5016  {
5017  if (sm_ptr == nullptr)
5018  {
5019  vector_ptrs[v] =
5020  vector_ptr +
5021  dof_info
5022  .dof_indices_contiguous[dof_access_index]
5023  [cells[v]];
5024  }
5025  else
5026  {
5027  const auto &temp =
5028  dof_info
5029  .dof_indices_contiguous_sm[dof_access_index]
5030  [cells[v]];
5031  vector_ptrs[v] = const_cast<Number2_ *>(
5032  sm_ptr->operator[](temp.first).data() +
5033  temp.second + index_offset);
5034  }
5035  }
5036  }
5037  else
5038  {
5039  Assert(false, ExcNotImplemented());
5040  }
5041  }
5042  else if (n_face_orientations == n_lanes)
5043  {
5044  for (unsigned int v = 0; v < n_lanes; ++v)
5045  reordered_indices[v] =
5046  dof_info.dof_indices_contiguous[dof_access_index]
5047  [fe_eval.get_cell_ids()[v]];
5048  dof_indices = reordered_indices.data();
5049  }
5050 
5051  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
5052  {
5053  if (vectorization_possible)
5054  for (unsigned int i = 0; i < dofs_per_face; ++i)
5055  {
5056  const unsigned int ind1 = index_array_hermite[0][2 * i];
5057  const unsigned int ind2 =
5058  index_array_hermite[0][2 * i + 1];
5059  const unsigned int i_ = reorientate(0, i);
5060 
5061  proc.hermite_grad_vectorized_indexed(
5062  temp1[i_],
5063  temp1[i_ + dofs_per_face],
5064  vector_ptr + ind1,
5065  vector_ptr + ind2,
5066  grad_weight,
5067  dof_indices,
5068  dof_indices);
5069  }
5070  else if (n_face_orientations == 1)
5071  for (unsigned int i = 0; i < dofs_per_face; ++i)
5072  {
5073  const unsigned int ind1 = index_array_hermite[0][2 * i];
5074  const unsigned int ind2 =
5075  index_array_hermite[0][2 * i + 1];
5076  const unsigned int i_ = reorientate(0, i);
5077 
5078  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5079  proc.hermite_grad(temp1[i_][v],
5080  temp1[i_ + dofs_per_face][v],
5081  vector_ptrs[v][ind1],
5082  vector_ptrs[v][ind2],
5083  grad_weight[v]);
5084 
5085  if (integrate == false)
5086  for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
5087  {
5088  temp1[i][v] = 0.0;
5089  temp1[i + dofs_per_face][v] = 0.0;
5090  }
5091  }
5092  else
5093  {
5094  if (integrate == false && n_filled_lanes < n_lanes)
5095  for (unsigned int i = 0; i < dofs_per_face; ++i)
5096  temp1[i] = temp1[i + dofs_per_face] = Number();
5097 
5098  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5099  for (unsigned int i = 0; i < dofs_per_face; ++i)
5100  proc.hermite_grad(
5101  temp1[reorientate(v, i)][v],
5102  temp1[reorientate(v, i) + dofs_per_face][v],
5103  vector_ptrs[v][index_array_hermite[v][2 * i]],
5104  vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
5105  grad_weight[v]);
5106  }
5107  }
5108  else
5109  {
5110  if (vectorization_possible)
5111  for (unsigned int i = 0; i < dofs_per_face; ++i)
5112  {
5113  const unsigned int ind = index_array_nodal[0][i];
5114  const unsigned int i_ = reorientate(0, i);
5115 
5116  proc.value_vectorized_indexed(temp1[i_],
5117  vector_ptr + ind,
5118  dof_indices);
5119  }
5120  else if (n_face_orientations == 1)
5121  for (unsigned int i = 0; i < dofs_per_face; ++i)
5122  {
5123  const unsigned int ind = index_array_nodal[0][i];
5124  const unsigned int i_ = reorientate(0, i);
5125 
5126  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5127  proc.value(temp1[i_][v], vector_ptrs[v][ind]);
5128 
5129  if (integrate == false)
5130  for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
5131  temp1[i_][v] = 0.0;
5132  }
5133  else
5134  {
5135  if (integrate == false && n_filled_lanes < n_lanes)
5136  for (unsigned int i = 0; i < dofs_per_face; ++i)
5137  temp1[i] = Number();
5138 
5139  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5140  for (unsigned int i = 0; i < dofs_per_face; ++i)
5141  proc.value(temp1[reorientate(v, i)][v],
5142  vector_ptrs[v][index_array_nodal[v][i]]);
5143  }
5144  }
5145  }
5146  else
5147  {
5148  // We should not end up here, this should be caught by
5149  // FEFaceEvaluationImplGatherEvaluateSelector::supports()
5150  Assert(false, ExcInternalError());
5151  }
5152  temp1 += 3 * dofs_per_face;
5153  }
5154  }
5155 
5156 
5157 
5158  template <int dim, typename Number2, typename VectorizedArrayType>
5160  {
5161  using Number = typename VectorizedArrayType::value_type;
5162 
5163  template <int fe_degree, int n_q_points_1d>
5164  static bool
5165  run(const unsigned int n_components,
5166  const EvaluationFlags::EvaluationFlags evaluation_flag,
5167  const Number2 * src_ptr,
5168  const std::vector<ArrayView<const Number2>> * sm_ptr,
5170  {
5171  Assert(fe_degree > -1, ExcInternalError());
5172  Assert(fe_eval.get_shape_info().element_type <=
5174  ExcInternalError());
5175 
5176  const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
5177 
5178  VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
5179  VectorizedArrayType *scratch_data =
5180  temp + 3 * n_components * dofs_per_face;
5181 
5183 
5184  if (fe_eval.get_dof_access_index() ==
5186  fe_eval.is_interior_face() == false)
5187  fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
5188  p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
5189  else
5190  fe_face_evaluation_process_and_io<1>(
5191  p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
5192 
5193  const unsigned int subface_index = fe_eval.get_subface_index();
5194 
5195  if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
5196  FEFaceEvaluationImpl<true,
5197  dim,
5198  fe_degree,
5199  n_q_points_1d,
5200  VectorizedArrayType>::
5201  evaluate_in_face(n_components,
5202  evaluation_flag,
5203  fe_eval.get_shape_info().data.front(),
5204  temp,
5205  fe_eval.begin_values(),
5206  fe_eval.begin_gradients(),
5207  fe_eval.begin_hessians(),
5208  scratch_data,
5209  subface_index);
5210  else
5211  FEFaceEvaluationImpl<false,
5212  dim,
5213  fe_degree,
5214  n_q_points_1d,
5215  VectorizedArrayType>::
5216  evaluate_in_face(n_components,
5217  evaluation_flag,
5218  fe_eval.get_shape_info().data.front(),
5219  temp,
5220  fe_eval.begin_values(),
5221  fe_eval.begin_gradients(),
5222  fe_eval.begin_hessians(),
5223  scratch_data,
5224  subface_index);
5225 
5226  // re-orientation for cases not possible with above algorithm
5227  if (subface_index < GeometryInfo<dim>::max_children_per_cell)
5228  {
5229  if (fe_eval.get_dof_access_index() ==
5231  fe_eval.is_interior_face() == false)
5232  {
5233  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5234  {
5235  // the loop breaks once an invalid_unsigned_int is hit for
5236  // all cases except the exterior faces in the ECL loop (where
5237  // some faces might be at the boundaries but others not)
5238  if (fe_eval.get_cell_ids()[v] ==
5240  continue;
5241 
5242  if (fe_eval.get_face_orientation(v) != 0)
5244  dim,
5245  n_components,
5246  v,
5247  evaluation_flag,
5249  fe_eval.get_face_orientation(v), 0),
5250  false,
5251  Utilities::pow(n_q_points_1d, dim - 1),
5252  &temp[0][0],
5253  fe_eval.begin_values(),
5254  fe_eval.begin_gradients(),
5255  fe_eval.begin_hessians());
5256  }
5257  }
5258  else if (fe_eval.get_face_orientation() != 0)
5260  dim,
5261  n_components,
5262  evaluation_flag,
5264  fe_eval.get_face_orientation(), 0),
5265  false,
5266  Utilities::pow(n_q_points_1d, dim - 1),
5267  temp,
5268  fe_eval.begin_values(),
5269  fe_eval.begin_gradients(),
5270  fe_eval.begin_hessians());
5271  }
5272 
5273  return false;
5274  }
5275 
5276  static bool
5278  const EvaluationFlags::EvaluationFlags evaluation_flag,
5280  const Number2 * vector_ptr,
5282  {
5283  const unsigned int fe_degree = shape_info.data.front().fe_degree;
5284  if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries ||
5285  (evaluation_flag & EvaluationFlags::gradients &&
5286  (fe_degree < 2 ||
5287  shape_info.data.front().element_type !=
5289  (evaluation_flag & EvaluationFlags::hessians) ||
5290  vector_ptr == nullptr ||
5291  shape_info.data.front().element_type >
5293  storage <
5295  return false;
5296  else
5297  return true;
5298  }
5299 
5300  private:
5301  template <int fe_degree>
5302  struct Processor
5303  {
5304  static const bool do_integrate = false;
5305  static const int dim_ = dim;
5306  static const int fe_degree_ = fe_degree;
5307  using VectorizedArrayType_ = VectorizedArrayType;
5308  using Number_ = Number;
5309  using Number2_ = const Number2;
5310 
5311  template <typename T0, typename T1, typename T2>
5312  void
5314  T0 & temp_2,
5315  const T1 src_ptr_1,
5316  const T1 src_ptr_2,
5317  const T2 &grad_weight)
5318  {
5319  do_vectorized_read(src_ptr_1, temp_1);
5320  do_vectorized_read(src_ptr_2, temp_2);
5321  temp_2 = grad_weight * (temp_1 - temp_2);
5322  }
5323 
5324  template <typename T1, typename T2>
5325  void
5326  value_vectorized(T1 &temp, const T2 src_ptr)
5327  {
5328  do_vectorized_read(src_ptr, temp);
5329  }
5330 
5331  template <typename T0, typename T1, typename T2, typename T3>
5332  void
5334  T0 & temp_2,
5335  const T1 src_ptr_1,
5336  const T1 src_ptr_2,
5337  const T2 &grad_weight,
5338  const T3 &indices_1,
5339  const T3 &indices_2)
5340  {
5341  do_vectorized_gather(src_ptr_1, indices_1, temp_1);
5342  do_vectorized_gather(src_ptr_2, indices_2, temp_2);
5343  temp_2 = grad_weight * (temp_1 - temp_2);
5344  }
5345 
5346  template <typename T0, typename T1, typename T2>
5347  void
5348  value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
5349  {
5350  do_vectorized_gather(src_ptr, indices, temp);
5351  }
5352 
5353  template <typename T0, typename T1, typename T2>
5354  void
5355  hermite_grad(T0 & temp_1,
5356  T0 & temp_2,
5357  const T1 &src_ptr_1,
5358  const T1 &src_ptr_2,
5359  const T2 &grad_weight)
5360  {
5361  // case 3a)
5362  temp_1 = src_ptr_1;
5363  temp_2 = grad_weight * (temp_1 - src_ptr_2);
5364  }
5365 
5366  template <typename T1, typename T2>
5367  void
5368  value(T1 &temp, const T2 &src_ptr)
5369  {
5370  // case 3b)
5371  temp = src_ptr;
5372  }
5373  };
5374  };
5375 
5376 
5377 
5378  template <int dim, typename Number2, typename VectorizedArrayType>
5380  {
5381  using Number = typename VectorizedArrayType::value_type;
5382 
5383  template <int fe_degree, int n_q_points_1d>
5384  static bool
5385  run(const unsigned int n_components,
5386  const EvaluationFlags::EvaluationFlags integration_flag,
5387  Number2 * dst_ptr,
5388  const std::vector<ArrayView<const Number2>> * sm_ptr,
5390  {
5391  Assert(fe_degree > -1, ExcInternalError());
5392  Assert(fe_eval.get_shape_info().element_type <=
5394  ExcInternalError());
5395 
5396  const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
5397 
5398  VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
5399  VectorizedArrayType *scratch_data =
5400  temp + 3 * n_components * dofs_per_face;
5401 
5402  const unsigned int subface_index = fe_eval.get_subface_index();
5403 
5404  // re-orientation for cases not possible with the io function below
5405  if (subface_index < GeometryInfo<dim>::max_children_per_cell)
5406  {
5407  if (fe_eval.get_dof_access_index() ==
5409  fe_eval.is_interior_face() == false)
5410  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5411  {
5412  // the loop breaks once an invalid_unsigned_int is hit for
5413  // all cases except the exterior faces in the ECL loop (where
5414  // some faces might be at the boundaries but others not)
5415  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
5416  continue;
5417 
5418  if (fe_eval.get_face_orientation(v) != 0)
5420  dim,
5421  n_components,
5422  v,
5423  integration_flag,
5425  fe_eval.get_face_orientation(v), 0),
5426  true,
5427  Utilities::pow(n_q_points_1d, dim - 1),
5428  &temp[0][0],
5429  fe_eval.begin_values(),
5430  fe_eval.begin_gradients(),
5431  fe_eval.begin_hessians());
5432  }
5433  else if (fe_eval.get_face_orientation() != 0)
5435  dim,
5436  n_components,
5437  integration_flag,
5439  fe_eval.get_face_orientation(), 0),
5440  true,
5441  Utilities::pow(n_q_points_1d, dim - 1),
5442  temp,
5443  fe_eval.begin_values(),
5444  fe_eval.begin_gradients(),
5445  fe_eval.begin_hessians());
5446  }
5447 
5448  if (fe_degree > -1 && fe_eval.get_subface_index() >=
5449  GeometryInfo<dim - 1>::max_children_per_cell)
5450  FEFaceEvaluationImpl<true,
5451  dim,
5452  fe_degree,
5453  n_q_points_1d,
5454  VectorizedArrayType>::
5455  integrate_in_face(n_components,
5456  integration_flag,
5457  fe_eval.get_shape_info().data.front(),
5458  temp,
5459  fe_eval.begin_values(),
5460  fe_eval.begin_gradients(),
5461  fe_eval.begin_hessians(),
5462  scratch_data,
5463  subface_index);
5464  else
5465  FEFaceEvaluationImpl<false,
5466  dim,
5467  fe_degree,
5468  n_q_points_1d,
5469  VectorizedArrayType>::
5470  integrate_in_face(n_components,
5471  integration_flag,
5472  fe_eval.get_shape_info().data.front(),
5473  temp,
5474  fe_eval.begin_values(),
5475  fe_eval.begin_gradients(),
5476  fe_eval.begin_hessians(),
5477  scratch_data,
5478  subface_index);
5479 
5481 
5482  if (fe_eval.get_dof_access_index() ==
5484  fe_eval.is_interior_face() == false)
5485  fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
5486  p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
5487  else
5488  fe_face_evaluation_process_and_io<1>(
5489  p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
5490 
5491  return false;
5492  }
5493 
5494  private:
5495  template <int fe_degree>
5496  struct Processor
5497  {
5498  static const bool do_integrate = true;
5499  static const int dim_ = dim;
5500  static const int fe_degree_ = fe_degree;
5501  using VectorizedArrayType_ = VectorizedArrayType;
5502  using Number_ = Number;
5503  using Number2_ = Number2;
5504 
5505  template <typename T0, typename T1, typename T2, typename T3, typename T4>
5506  void
5507  hermite_grad_vectorized(const T0 &temp_1,
5508  const T1 &temp_2,
5509  T2 dst_ptr_1,
5510  T3 dst_ptr_2,
5511  const T4 &grad_weight)
5512  {
5513  // case 1a)
5514  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
5515  const VectorizedArrayType grad = grad_weight * temp_2;
5516  do_vectorized_add(val, dst_ptr_1);
5517  do_vectorized_add(grad, dst_ptr_2);
5518  }
5519 
5520  template <typename T0, typename T1>
5521  void
5522  value_vectorized(const T0 &temp, T1 dst_ptr)
5523  {
5524  // case 1b)
5525  do_vectorized_add(temp, dst_ptr);
5526  }
5527 
5528  template <typename T0, typename T1, typename T2, typename T3>
5529  void
5531  const T0 &temp_2,
5532  T1 dst_ptr_1,
5533  T1 dst_ptr_2,
5534  const T2 &grad_weight,
5535  const T3 &indices_1,
5536  const T3 &indices_2)
5537  {
5538  // case 2a)
5539  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
5540  const VectorizedArrayType grad = grad_weight * temp_2;
5541  do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
5542  do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
5543  }
5544 
5545  template <typename T0, typename T1, typename T2>
5546  void
5547  value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
5548  {
5549  // case 2b)
5550  do_vectorized_scatter_add(temp, indices, dst_ptr);
5551  }
5552 
5553  template <typename T0, typename T1, typename T2>
5554  void
5555  hermite_grad(const T0 &temp_1,
5556  const T0 &temp_2,
5557  T1 & dst_ptr_1,
5558  T1 & dst_ptr_2,
5559  const T2 &grad_weight)
5560  {
5561  // case 3a)
5562  const Number val = temp_1 - grad_weight * temp_2;
5563  const Number grad = grad_weight * temp_2;
5564  dst_ptr_1 += val;
5565  dst_ptr_2 += grad;
5566  }
5567 
5568  template <typename T0, typename T1>
5569  void
5570  value(const T0 &temp, T1 &dst_ptr)
5571  {
5572  // case 3b)
5573  dst_ptr += temp;
5574  }
5575  };
5576  };
5577 
5578 
5579 
5584  template <int dim, typename Number>
5586  {
5587  template <int fe_degree, int = 0>
5588  static bool
5589  run(const unsigned int n_components,
5590  const FEEvaluationData<dim, Number, false> &fe_eval,
5591  const Number * in_array,
5592  Number * out_array,
5593  std::enable_if_t<fe_degree != -1> * = nullptr)
5594  {
5595  constexpr unsigned int dofs_per_component =
5596  Utilities::pow(fe_degree + 1, dim);
5597 
5598  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
5599  Assert(fe_eval.get_shape_info().element_type <=
5601  ExcNotImplemented());
5602 
5604  dim,
5605  fe_degree + 1,
5606  fe_degree + 1,
5607  Number>
5608  evaluator(
5611  fe_eval.get_shape_info().data.front().inverse_shape_values_eo);
5612 
5613  for (unsigned int d = 0; d < n_components; ++d)
5614  {
5615  const Number *in = in_array + d * dofs_per_component;
5616  Number * out = out_array + d * dofs_per_component;
5617  // Need to select 'apply' method with hessian slot because values
5618  // assume symmetries that do not exist in the inverse shapes
5619  evaluator.template hessians<0, true, false>(in, out);
5620  if (dim > 1)
5621  evaluator.template hessians<1, true, false>(out, out);
5622  if (dim > 2)
5623  evaluator.template hessians<2, true, false>(out, out);
5624  }
5625  for (unsigned int q = 0; q < dofs_per_component; ++q)
5626  {
5627  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
5628  for (unsigned int d = 0; d < n_components; ++d)
5629  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
5630  }
5631  for (unsigned int d = 0; d < n_components; ++d)
5632  {
5633  Number *out = out_array + d * dofs_per_component;
5634  if (dim > 2)
5635  evaluator.template hessians<2, false, false>(out, out);
5636  if (dim > 1)
5637  evaluator.template hessians<1, false, false>(out, out);
5638  evaluator.template hessians<0, false, false>(out, out);
5639  }
5640  return false;
5641  }
5642 
5643  template <int fe_degree, int = 0>
5644  static bool
5645  run(const unsigned int n_components,
5646  const FEEvaluationData<dim, Number, false> &fe_eval,
5647  const Number * in_array,
5648  Number * out_array,
5649  std::enable_if_t<fe_degree == -1> * = nullptr)
5650  {
5651  static_assert(fe_degree == -1, "Only usable for degree -1");
5652  const unsigned int dofs_per_component =
5654 
5655  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
5656 
5658  fe_eval.get_shape_info().data.front().inverse_shape_values,
5661  fe_eval.get_shape_info().data.front().fe_degree + 1,
5662  fe_eval.get_shape_info().data.front().fe_degree + 1);
5663 
5664  for (unsigned int d = 0; d < n_components; ++d)
5665  {
5666  const Number *in = in_array + d * dofs_per_component;
5667  Number * out = out_array + d * dofs_per_component;
5668  // Need to select 'apply' method with hessian slot because values
5669  // assume symmetries that do not exist in the inverse shapes
5670  evaluator.template values<0, true, false>(in, out);
5671  if (dim > 1)
5672  evaluator.template values<1, true, false>(out, out);
5673  if (dim > 2)
5674  evaluator.template values<2, true, false>(out, out);
5675  }
5676  for (unsigned int q = 0; q < dofs_per_component; ++q)
5677  {
5678  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
5679  for (unsigned int d = 0; d < n_components; ++d)
5680  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
5681  }
5682  for (unsigned int d = 0; d < n_components; ++d)
5683  {
5684  Number *out = out_array + d * dofs_per_component;
5685  if (dim > 2)
5686  evaluator.template values<2, false, false>(out, out);
5687  if (dim > 1)
5688  evaluator.template values<1, false, false>(out, out);
5689  evaluator.template values<0, false, false>(out, out);
5690  }
5691  return false;
5692  }
5693  };
5694 
5695 
5696 
5703  template <int dim, typename Number>
5705  {
5706  template <int fe_degree, int = 0>
5707  static bool
5708  run(const unsigned int n_desired_components,
5709  const FEEvaluationData<dim, Number, false> &fe_eval,
5710  const ArrayView<const Number> & inverse_coefficients,
5711  const bool dyadic_coefficients,
5712  const Number * in_array,
5713  Number * out_array)
5714  {
5715  const unsigned int given_degree =
5716  (fe_degree > -1) ? fe_degree :
5717  fe_eval.get_shape_info().data.front().fe_degree;
5718 
5719  const unsigned int dofs_per_component =
5720  Utilities::pow(given_degree + 1, dim);
5721 
5722  Assert(inverse_coefficients.size() > 0 &&
5723  inverse_coefficients.size() % dofs_per_component == 0,
5724  ExcMessage(
5725  "Expected diagonal to be a multiple of scalar dof per cells"));
5726 
5727  if (!dyadic_coefficients)
5728  {
5729  if (inverse_coefficients.size() != dofs_per_component)
5730  AssertDimension(n_desired_components * dofs_per_component,
5731  inverse_coefficients.size());
5732  }
5733  else
5734  {
5735  AssertDimension(n_desired_components * n_desired_components *
5736  dofs_per_component,
5737  inverse_coefficients.size());
5738  }
5739 
5740  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
5741  Assert(fe_eval.get_shape_info().element_type <=
5743  ExcNotImplemented());
5744 
5746  dim,
5747  fe_degree + 1,
5748  fe_degree + 1,
5749  Number>
5750  evaluator(AlignedVector<Number>(),
5752  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
5753  given_degree + 1,
5754  given_degree + 1);
5755 
5756  const Number *in = in_array;
5757  Number * out = out_array;
5758 
5759  const Number *inv_coefficient = inverse_coefficients.data();
5760 
5761  const unsigned int shift_coefficient =
5762  inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
5763  0;
5764 
5765  const auto n_comp_outer = dyadic_coefficients ? 1 : n_desired_components;
5766  const auto n_comp_inner = dyadic_coefficients ? n_desired_components : 1;
5767 
5768  for (unsigned int d = 0; d < n_comp_outer; ++d)
5769  {
5770  for (unsigned int di = 0; di < n_comp_inner; ++di)
5771  {
5772  const Number *in_ = in + di * dofs_per_component;
5773  Number * out_ = out + di * dofs_per_component;
5774  evaluator.template hessians<0, true, false>(in_, out_);
5775  if (dim > 1)
5776  evaluator.template hessians<1, true, false>(out_, out_);
5777  if (dim > 2)
5778  evaluator.template hessians<2, true, false>(out_, out_);
5779  }
5780  if (dyadic_coefficients)
5781  {
5782  const auto n_coeff_components =
5783  n_desired_components * n_desired_components;
5784  if (n_desired_components == dim)
5785  {
5786  for (unsigned int q = 0; q < dofs_per_component; ++q)
5787  vmult<dim>(&inv_coefficient[q * n_coeff_components],
5788  &in[q],
5789  &out[q],
5790  dofs_per_component);
5791  }
5792  else
5793  {
5794  for (unsigned int q = 0; q < dofs_per_component; ++q)
5795  vmult<-1>(&inv_coefficient[q * n_coeff_components],
5796  &in[q],
5797  &out[q],
5798  dofs_per_component,
5799  n_desired_components);
5800  }
5801  }
5802  else
5803  for (unsigned int q = 0; q < dofs_per_component; ++q)
5804  out[q] *= inv_coefficient[q];
5805 
5806  for (unsigned int di = 0; di < n_comp_inner; ++di)
5807  {
5808  Number *out_ = out + di * dofs_per_component;
5809  if (dim > 2)
5810  evaluator.template hessians<2, false, false>(out_, out_);
5811  if (dim > 1)
5812  evaluator.template hessians<1, false, false>(out_, out_);
5813  evaluator.template hessians<0, false, false>(out_, out_);
5814  }
5815 
5816  in += dofs_per_component;
5817  out += dofs_per_component;
5818  inv_coefficient += shift_coefficient;
5819  }
5820 
5821  return false;
5822  }
5823 
5824  private:
5825  template <int n_components>
5826  static inline void
5827  vmult(const Number * inverse_coefficients,
5828  const Number * src,
5829  Number * dst,
5830  const unsigned int dofs_per_component,
5831  const unsigned int n_given_components = 0)
5832  {
5833  const unsigned int n_desired_components =
5834  (n_components > -1) ? n_components : n_given_components;
5835 
5836  std::array<Number, dim + 2> tmp = {};
5837  Assert(n_desired_components <= dim + 2,
5838  ExcMessage(
5839  "Number of components larger than dim+2 not supported."));
5840 
5841  for (unsigned int d = 0; d < n_desired_components; ++d)
5842  tmp[d] = src[d * dofs_per_component];
5843 
5844  for (unsigned int d1 = 0; d1 < n_desired_components; ++d1)
5845  {
5846  const Number *inv_coeff_row =
5847  &inverse_coefficients[d1 * n_desired_components];
5848  Number sum = inv_coeff_row[0] * tmp[0];
5849  for (unsigned int d2 = 1; d2 < n_desired_components; ++d2)
5850  sum += inv_coeff_row[d2] * tmp[d2];
5851  dst[d1 * dofs_per_component] = sum;
5852  }
5853  }
5854  };
5855 
5856 
5857 
5864  template <int dim, typename Number>
5866  {
5867  template <int fe_degree, int n_q_points_1d>
5868  static bool
5869  run(const unsigned int n_desired_components,
5870  const FEEvaluationData<dim, Number, false> &fe_eval,
5871  const Number * in_array,
5872  Number * out_array)
5873  {
5874  static const bool do_inplace =
5875  fe_degree > -1 && (fe_degree + 1 == n_q_points_1d);
5876 
5877  Assert(fe_eval.get_shape_info().element_type !=
5879  ExcNotImplemented());
5880 
5881  const auto &inverse_shape =
5882  do_inplace ?
5883  fe_eval.get_shape_info().data.front().inverse_shape_values_eo :
5884  fe_eval.get_shape_info().data.front().inverse_shape_values;
5885 
5886  const std::size_t dofs_per_component =
5887  do_inplace ? Utilities::pow(fe_degree + 1, dim) :
5889  const std::size_t n_q_points = do_inplace ?
5890  Utilities::pow(fe_degree + 1, dim) :
5891  fe_eval.get_shape_info().n_q_points;
5892 
5894  dim,
5895  fe_degree + 1,
5896  n_q_points_1d,
5897  Number>
5898  evaluator(AlignedVector<Number>(),
5900  inverse_shape,
5901  fe_eval.get_shape_info().data.front().fe_degree + 1,
5902  fe_eval.get_shape_info().data.front().n_q_points_1d);
5903 
5904  for (unsigned int d = 0; d < n_desired_components; ++d)
5905  {
5906  const Number *in = in_array + d * n_q_points;
5907  Number * out = out_array + d * dofs_per_component;
5908 
5909  auto temp_1 = do_inplace ? out : fe_eval.get_scratch_data().begin();
5910  auto temp_2 = do_inplace ?
5911  out :
5912  (temp_1 + std::max(n_q_points, dofs_per_component));
5913 
5914  if (dim == 3)
5915  {
5916  evaluator.template hessians<2, false, false>(in, temp_1);
5917  evaluator.template hessians<1, false, false>(temp_1, temp_2);
5918  evaluator.template hessians<0, false, false>(temp_2, out);
5919  }
5920  if (dim == 2)
5921  {
5922  evaluator.template hessians<1, false, false>(in, temp_1);
5923  evaluator.template hessians<0, false, false>(temp_1, out);
5924  }
5925  if (dim == 1)
5926  evaluator.template hessians<0, false, false>(in, out);
5927  }
5928  return false;
5929  }
5930  };
5931 
5932 } // end of namespace internal
5933 
5934 
5936 
5937 #endif
size_type size() const
iterator begin() const
Definition: array_view.h:594
value_type * data() const noexcept
Definition: array_view.h:553
std::size_t size() const
Definition: array_view.h:576
std::uint8_t get_face_no(const unsigned int v=0) const
internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex get_dof_access_index() const
Number JxW(const unsigned int q_point) const
const Number * begin_values() const
const Number * begin_hessians() const
const std::array< unsigned int, n_lanes > & get_cell_ids() const
unsigned int get_subface_index() const
const Number * begin_gradients() const
bool is_interior_face() const
ArrayView< Number > get_scratch_data() const
const ShapeInfoType & get_shape_info() const
std::uint8_t get_face_orientation(const unsigned int v=0) const
void gather(const Number *base_ptr, const unsigned int *offsets)
void load(const OtherNumber *ptr)
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:108
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:140
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
unsigned int cell_index
Definition: grid_tools.cc:1191
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
EvaluationFlags
The EvaluationFlags enum.
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2131
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm &mpi_communicator)
T scatter(const MPI_Comm &comm, const std::vector< T > &objects_to_send, const unsigned int root_process=0)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:447
T fixed_power(const T t)
Definition: utilities.h:966
void do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
void adjust_for_face_orientation_per_lane(const unsigned int dim, const unsigned int n_components, const unsigned int v, const EvaluationFlags::EvaluationFlags flag, const unsigned int *orientation, const bool integrate, const std::size_t n_q_points, Number *tmp_values, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad=nullptr, VectorizedArrayType *hessians_quad=nullptr)
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)
void do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
void fe_face_evaluation_process_and_io(Processor &proc, const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, typename Processor::Number2_ *global_vector_ptr, const std::vector< ArrayView< const typename Processor::Number2_ >> *sm_ptr, const EvaluationData &fe_eval, typename Processor::VectorizedArrayType_ *temp1)
void adjust_for_face_orientation(const unsigned int dim, const unsigned int n_components, const EvaluationFlags::EvaluationFlags flag, const unsigned int *orientation, const bool integrate, const std::size_t n_q_points, Number *tmp_values, Number *values_quad, Number *gradients_quad, Number *hessians_quad)
static const unsigned int invalid_unsigned_int
Definition: types.h:213
static bool run(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array, std::enable_if_t< fe_degree==-1 > *=nullptr)
static bool run(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array, std::enable_if_t< fe_degree !=-1 > *=nullptr)
static bool run(const unsigned int n_desired_components, const FEEvaluationData< dim, Number, false > &fe_eval, const ArrayView< const Number > &inverse_coefficients, const bool dyadic_coefficients, const Number *in_array, Number *out_array)
static void vmult(const Number *inverse_coefficients, const Number *src, Number *dst, const unsigned int dofs_per_component, const unsigned int n_given_components=0)
static bool run(const unsigned int n_desired_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array)
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)
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)
static void do_forward(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const Number *values_in, Number *values_out, const unsigned int basis_size_1_variable=numbers::invalid_unsigned_int, const unsigned int basis_size_2_variable=numbers::invalid_unsigned_int)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
static void do_integrate(const MatrixFreeFunctions::UnivariateShapeData< Number > &shape, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, Number *gradients_quad, const Number *hessians_quad, const bool add_into_values_array)
static void do_evaluate(const MatrixFreeFunctions::UnivariateShapeData< Number > &shape, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, Number *gradients_quad, Number *hessians_quad)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static EvalType create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number > &shape_data)