Reference documentation for deal.II version GIT 20f059c89a 2022-12-04 14:55:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
evaluation_kernels.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  const std::size_t n_dofs =
840  const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
841 
842  const auto &shape_data = fe_eval.get_shape_info().data;
843 
844  using Eval =
846 
847  if (evaluation_flag & EvaluationFlags::values)
848  {
849  const auto shape_values = shape_data.front().shape_values.data();
850  auto values_quad_ptr = fe_eval.begin_values();
851  auto values_dofs_actual_ptr = values_dofs_actual;
852 
853  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
854  for (unsigned int c = 0; c < n_components; ++c)
855  {
856  eval.template values<0, true, false>(values_dofs_actual_ptr,
857  values_quad_ptr);
858 
859  values_quad_ptr += n_q_points;
860  values_dofs_actual_ptr += n_dofs;
861  }
862  }
863 
864  if (evaluation_flag & EvaluationFlags::gradients)
865  {
866  const auto shape_gradients = shape_data.front().shape_gradients.data();
867  auto gradients_quad_ptr = fe_eval.begin_gradients();
868  auto values_dofs_actual_ptr = values_dofs_actual;
869 
870  for (unsigned int c = 0; c < n_components; ++c)
871  {
872  for (unsigned int d = 0; d < dim; ++d)
873  {
874  Eval eval(nullptr,
875  shape_gradients + n_q_points * n_dofs * d,
876  nullptr,
877  n_dofs,
878  n_q_points);
879 
880  eval.template gradients<0, true, false>(values_dofs_actual_ptr,
881  gradients_quad_ptr);
882 
883  gradients_quad_ptr += n_q_points;
884  }
885  values_dofs_actual_ptr += n_dofs;
886  }
887  }
888 
889  if (evaluation_flag & EvaluationFlags::hessians)
890  Assert(false, ExcNotImplemented());
891  }
892 
893 
894 
895  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
896  inline void
899  dim,
900  fe_degree,
901  n_q_points_1d,
902  Number>::integrate(const unsigned int n_components,
903  const EvaluationFlags::EvaluationFlags integration_flag,
904  Number * values_dofs_actual,
906  const bool add_into_values_array)
907  {
908  // TODO: implement hessians
909  AssertThrow(!(integration_flag & EvaluationFlags::hessians),
911 
912  const std::size_t n_dofs =
914  const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
915 
916  const auto &shape_data = fe_eval.get_shape_info().data;
917 
918  using Eval =
920 
921  if (integration_flag & EvaluationFlags::values)
922  {
923  const auto shape_values = shape_data.front().shape_values.data();
924  auto values_quad_ptr = fe_eval.begin_values();
925  auto values_dofs_actual_ptr = values_dofs_actual;
926 
927  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
928  for (unsigned int c = 0; c < n_components; ++c)
929  {
930  if (add_into_values_array == false)
931  eval.template values<0, false, false>(values_quad_ptr,
932  values_dofs_actual_ptr);
933  else
934  eval.template values<0, false, true>(values_quad_ptr,
935  values_dofs_actual_ptr);
936 
937  values_quad_ptr += n_q_points;
938  values_dofs_actual_ptr += n_dofs;
939  }
940  }
941 
942  if (integration_flag & EvaluationFlags::gradients)
943  {
944  const auto shape_gradients = shape_data.front().shape_gradients.data();
945  auto gradients_quad_ptr = fe_eval.begin_gradients();
946  auto values_dofs_actual_ptr = values_dofs_actual;
947 
948  for (unsigned int c = 0; c < n_components; ++c)
949  {
950  for (unsigned int d = 0; d < dim; ++d)
951  {
952  Eval eval(nullptr,
953  shape_gradients + n_q_points * n_dofs * d,
954  nullptr,
955  n_dofs,
956  n_q_points);
957 
958  if ((add_into_values_array == false &&
959  !(integration_flag & EvaluationFlags::values)) &&
960  d == 0)
961  eval.template gradients<0, false, false>(
962  gradients_quad_ptr, values_dofs_actual_ptr);
963  else
964  eval.template gradients<0, false, true>(
965  gradients_quad_ptr, values_dofs_actual_ptr);
966 
967  gradients_quad_ptr += n_q_points;
968  }
969  values_dofs_actual_ptr += n_dofs;
970  }
971  }
972  }
973 
974 
975  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
976  template <bool integrate>
977  inline void
979  dim,
980  fe_degree,
981  n_q_points_1d,
982  Number>::
983  evaluate_or_integrate(
984  const EvaluationFlags::EvaluationFlags evaluation_flag,
985  Number * values_dofs_actual,
987  const bool add_into_values_array)
988  {
989  if (evaluation_flag == EvaluationFlags::nothing)
990  return;
991 
992  AssertDimension(fe_eval.get_shape_info().data.size(), 2);
993  // First component:
994  evaluate_tensor_product_per_component<0>(
995  evaluation_flag,
996  values_dofs_actual,
997  fe_eval,
998  add_into_values_array,
999  std::integral_constant<bool, integrate>());
1000  // Second component :
1001  evaluate_tensor_product_per_component<1>(
1002  evaluation_flag,
1003  values_dofs_actual,
1004  fe_eval,
1005  add_into_values_array,
1006  std::integral_constant<bool, integrate>());
1007  if (dim == 3)
1008  {
1009  // Third component
1010  evaluate_tensor_product_per_component<2>(
1011  evaluation_flag,
1012  values_dofs_actual,
1013  fe_eval,
1014  add_into_values_array,
1015  std::integral_constant<bool, integrate>());
1016  }
1017  }
1018 
1019  // Helper function that applies the 1d evaluation kernels.
1020  // std::integral_constant<bool, false> is the interpolation path, and
1021  // std::integral_constant<bool, true> below is the integration path.
1022  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1023  template <int normal_dir>
1024  inline void
1026  dim,
1027  fe_degree,
1028  n_q_points_1d,
1029  Number>::
1030  evaluate_tensor_product_per_component(
1031  const EvaluationFlags::EvaluationFlags evaluation_flag,
1032  Number * values_dofs_actual,
1034  const bool add_into_values_array,
1035  std::integral_constant<bool, false>)
1036  {
1037  (void)add_into_values_array;
1038 
1039  using EvalNormal =
1041  dim,
1042  (fe_degree == -1) ? 1 : fe_degree + 1,
1043  n_q_points_1d,
1044  Number,
1045  normal_dir>;
1046 
1047  using EvalTangent =
1049  dim,
1050  (fe_degree == -1) ? 1 : fe_degree,
1051  n_q_points_1d,
1052  Number,
1053  normal_dir>;
1054  using Eval0 =
1055  typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
1056  using Eval1 =
1057  typename std::conditional<normal_dir == 1, EvalNormal, EvalTangent>::type;
1058  using Eval2 =
1059  typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
1060 
1061  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info =
1062  fe_eval.get_shape_info();
1063  Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
1064  ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]));
1065  Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
1066  ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]));
1067  Eval2 eval2 = create_evaluator_tensor_product<Eval2>(
1068  ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1]));
1069 
1070  Number *temp1 = fe_eval.get_scratch_data().begin();
1071  Number *temp2;
1072 
1073  temp2 =
1074  temp1 +
1075  std::max(Utilities::fixed_power<dim>(shape_info.data[0].fe_degree + 1),
1076  Utilities::fixed_power<dim>(shape_info.data[0].n_q_points_1d));
1077 
1078  const std::size_t n_q_points = shape_info.n_q_points;
1079  const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell;
1080 
1081  // Initial shift depending on component (normal_dir)
1082  Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
1083  Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
1084  Number *gradients_quad =
1085  fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
1086  Number *hessians_quad =
1087  (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
1088  fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
1089 
1090  switch (dim)
1091  {
1092  case 2:
1093  if (evaluation_flag & EvaluationFlags::gradients)
1094  {
1095  eval0.template gradients<0, true, false>(values_dofs, temp1);
1096  eval1.template values<1, true, false>(temp1, gradients_quad);
1097  }
1098  if (evaluation_flag & EvaluationFlags::hessians)
1099  {
1100  // The evaluation/integration here *should* work, however
1101  // the piola transform is not implemented.
1102  AssertThrow(false, ExcNotImplemented());
1103  // grad xy
1104  if (!(evaluation_flag & EvaluationFlags::gradients))
1105  eval0.template gradients<0, true, false>(values_dofs, temp1);
1106  eval1.template gradients<1, true, false>(temp1,
1107  hessians_quad +
1108  2 * n_q_points);
1109 
1110  // grad xx
1111  eval0.template hessians<0, true, false>(values_dofs, temp1);
1112  eval1.template values<1, true, false>(temp1, hessians_quad);
1113  }
1114 
1115  // grad y
1116  eval0.template values<0, true, false>(values_dofs, temp1);
1117  if (evaluation_flag & EvaluationFlags::gradients)
1118  eval1.template gradients<1, true, false>(temp1,
1119  gradients_quad +
1120  n_q_points);
1121 
1122  // grad yy
1123  if (evaluation_flag & EvaluationFlags::hessians)
1124  eval1.template hessians<1, true, false>(temp1,
1125  hessians_quad + n_q_points);
1126 
1127  // val: can use values applied in x
1128  if (evaluation_flag & EvaluationFlags::values)
1129  eval1.template values<1, true, false>(temp1, values_quad);
1130  break;
1131  case 3:
1132  if (evaluation_flag & EvaluationFlags::gradients)
1133  {
1134  // grad x
1135  eval0.template gradients<0, true, false>(values_dofs, temp1);
1136  eval1.template values<1, true, false>(temp1, temp2);
1137  eval2.template values<2, true, false>(temp2, gradients_quad);
1138  }
1139 
1140  if (evaluation_flag & EvaluationFlags::hessians)
1141  {
1142  // The evaluation/integration here *should* work, however
1143  // the piola transform is not implemented.
1144  AssertThrow(false, ExcNotImplemented());
1145  // grad xz
1146  if (!(evaluation_flag & EvaluationFlags::gradients))
1147  {
1148  eval0.template gradients<0, true, false>(values_dofs, temp1);
1149  eval1.template values<1, true, false>(temp1, temp2);
1150  }
1151  eval2.template gradients<2, true, false>(temp2,
1152  hessians_quad +
1153  4 * n_q_points);
1154 
1155  // grad xy
1156  eval1.template gradients<1, true, false>(temp1, temp2);
1157  eval2.template values<2, true, false>(temp2,
1158  hessians_quad +
1159  3 * n_q_points);
1160 
1161  // grad xx
1162  eval0.template hessians<0, true, false>(values_dofs, temp1);
1163  eval1.template values<1, true, false>(temp1, temp2);
1164  eval2.template values<2, true, false>(temp2, hessians_quad);
1165  }
1166 
1167  // grad y
1168  eval0.template values<0, true, false>(values_dofs, temp1);
1169  if (evaluation_flag & EvaluationFlags::gradients)
1170  {
1171  eval1.template gradients<1, true, false>(temp1, temp2);
1172  eval2.template values<2, true, false>(temp2,
1173  gradients_quad +
1174  n_q_points);
1175  }
1176 
1177  if (evaluation_flag & EvaluationFlags::hessians)
1178  {
1179  // grad yz
1180  if (!(evaluation_flag & EvaluationFlags::gradients))
1181  eval1.template gradients<1, true, false>(temp1, temp2);
1182  eval2.template gradients<2, true, false>(temp2,
1183  hessians_quad +
1184  5 * n_q_points);
1185 
1186  // grad yy
1187  eval1.template hessians<1, true, false>(temp1, temp2);
1188  eval2.template values<2, true, false>(temp2,
1189  hessians_quad + n_q_points);
1190  }
1191 
1192  // grad z: can use the values applied in x direction stored in
1193  // temp1
1194  eval1.template values<1, true, false>(temp1, temp2);
1195  if (evaluation_flag & EvaluationFlags::gradients)
1196  eval2.template gradients<2, true, false>(temp2,
1197  gradients_quad +
1198  2 * n_q_points);
1199 
1200  // grad zz: can use the values applied in x and y direction stored
1201  // in temp2
1202  if (evaluation_flag & EvaluationFlags::hessians)
1203  eval2.template hessians<2, true, false>(temp2,
1204  hessians_quad +
1205  2 * n_q_points);
1206 
1207  // val: can use the values applied in x & y direction stored in
1208  // temp2
1209  if (evaluation_flag & EvaluationFlags::values)
1210  eval2.template values<2, true, false>(temp2, values_quad);
1211  break;
1212  default:
1213  AssertThrow(false, ExcNotImplemented());
1214  }
1215  }
1216 
1217  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1218  template <int normal_dir>
1219  inline void
1221  dim,
1222  fe_degree,
1223  n_q_points_1d,
1224  Number>::
1225  evaluate_tensor_product_per_component(
1226  const EvaluationFlags::EvaluationFlags evaluation_flag,
1227  Number * values_dofs_actual,
1229  const bool add_into_values_array,
1230  std::integral_constant<bool, true>)
1231  {
1232  using EvalNormal =
1234  dim,
1235  (fe_degree == -1) ? 1 : fe_degree + 1,
1236  n_q_points_1d,
1237  Number,
1238  normal_dir>;
1239 
1240  using EvalTangent =
1242  dim,
1243  (fe_degree == -1) ? 1 : fe_degree,
1244  n_q_points_1d,
1245  Number,
1246  normal_dir>;
1247  using Eval0 =
1248  typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
1249  using Eval1 =
1250  typename std::conditional<normal_dir == 1, EvalNormal, EvalTangent>::type;
1251  using Eval2 =
1252  typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
1253 
1254  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info =
1255  fe_eval.get_shape_info();
1256  Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
1257  ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]));
1258  Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
1259  ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]));
1260  Eval2 eval2 = create_evaluator_tensor_product<Eval2>(
1261  ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1]));
1262 
1263  Number *temp1 = fe_eval.get_scratch_data().begin();
1264  Number *temp2;
1265 
1266  temp2 =
1267  temp1 +
1268  std::max(Utilities::fixed_power<dim>(shape_info.data[0].fe_degree + 1),
1269  Utilities::fixed_power<dim>(shape_info.data[0].n_q_points_1d));
1270 
1271  const std::size_t n_q_points = shape_info.n_q_points;
1272  const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell;
1273 
1274  // Initial shift depending on component (normal_dir)
1275  Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
1276  Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
1277  Number *gradients_quad =
1278  fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
1279  Number *hessians_quad =
1280  (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
1281  fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
1282 
1283  // Integrate path
1284  switch (dim)
1285  {
1286  case 2:
1287  if ((evaluation_flag & EvaluationFlags::values) &&
1288  !(evaluation_flag & EvaluationFlags::gradients))
1289  {
1290  eval1.template values<1, false, false>(values_quad, temp1);
1291  if (add_into_values_array == false)
1292  eval0.template values<0, false, false>(temp1, values_dofs);
1293  else
1294  eval0.template values<0, false, true>(temp1, values_dofs);
1295  }
1296  if (evaluation_flag & EvaluationFlags::gradients)
1297  {
1298  eval1.template gradients<1, false, false>(gradients_quad +
1299  n_q_points,
1300  temp1);
1301  if ((evaluation_flag & EvaluationFlags::values))
1302  eval1.template values<1, false, true>(values_quad, temp1);
1303  if (add_into_values_array == false)
1304  eval0.template values<0, false, false>(temp1, values_dofs);
1305  else
1306  eval0.template values<0, false, true>(temp1, values_dofs);
1307  eval1.template values<1, false, false>(gradients_quad, temp1);
1308  eval0.template gradients<0, false, true>(temp1, values_dofs);
1309  }
1310  if (evaluation_flag & EvaluationFlags::hessians)
1311  {
1312  // grad xx
1313  eval1.template values<1, false, false>(hessians_quad, temp1);
1314 
1315  if ((evaluation_flag & EvaluationFlags::values) ||
1316  (evaluation_flag & EvaluationFlags::gradients) ||
1317  add_into_values_array == true)
1318  eval0.template hessians<0, false, true>(temp1, values_dofs);
1319  else
1320  eval0.template hessians<0, false, false>(temp1, values_dofs);
1321 
1322  // grad yy
1323  eval1.template hessians<1, false, false>(hessians_quad +
1324  n_q_points,
1325  temp1);
1326  eval0.template values<0, false, true>(temp1, values_dofs);
1327 
1328  // grad xy
1329  eval1.template gradients<1, false, false>(hessians_quad +
1330  2 * n_q_points,
1331  temp1);
1332  eval0.template gradients<0, false, true>(temp1, values_dofs);
1333  }
1334  break;
1335 
1336  case 3:
1337  if ((evaluation_flag & EvaluationFlags::values) &&
1338  !(evaluation_flag & EvaluationFlags::gradients))
1339  {
1340  eval2.template values<2, false, false>(values_quad, temp1);
1341  eval1.template values<1, false, false>(temp1, temp2);
1342  if (add_into_values_array == false)
1343  eval0.template values<0, false, false>(temp2, values_dofs);
1344  else
1345  eval0.template values<0, false, true>(temp2, values_dofs);
1346  }
1347  if (evaluation_flag & EvaluationFlags::gradients)
1348  {
1349  eval2.template gradients<2, false, false>(gradients_quad +
1350  2 * n_q_points,
1351  temp1);
1352  if ((evaluation_flag & EvaluationFlags::values))
1353  eval2.template values<2, false, true>(values_quad, temp1);
1354  eval1.template values<1, false, false>(temp1, temp2);
1355  eval2.template values<2, false, false>(gradients_quad +
1356  n_q_points,
1357  temp1);
1358  eval1.template gradients<1, false, true>(temp1, temp2);
1359  if (add_into_values_array == false)
1360  eval0.template values<0, false, false>(temp2, values_dofs);
1361  else
1362  eval0.template values<0, false, true>(temp2, values_dofs);
1363  eval2.template values<2, false, false>(gradients_quad, temp1);
1364  eval1.template values<1, false, false>(temp1, temp2);
1365  eval0.template gradients<0, false, true>(temp2, values_dofs);
1366  }
1367  if (evaluation_flag & EvaluationFlags::hessians)
1368  {
1369  // grad xx
1370  eval2.template values<2, false, false>(hessians_quad, temp1);
1371  eval1.template values<1, false, false>(temp1, temp2);
1372 
1373  if ((evaluation_flag & EvaluationFlags::values) ||
1374  (evaluation_flag & EvaluationFlags::gradients) ||
1375  add_into_values_array == true)
1376  eval0.template hessians<0, false, true>(temp2, values_dofs);
1377  else
1378  eval0.template hessians<0, false, false>(temp2, values_dofs);
1379 
1380  // grad yy
1381  eval2.template values<2, false, false>(hessians_quad + n_q_points,
1382  temp1);
1383  eval1.template hessians<1, false, false>(temp1, temp2);
1384  eval0.template values<0, false, true>(temp2, values_dofs);
1385 
1386  // grad zz
1387  eval2.template hessians<2, false, false>(hessians_quad +
1388  2 * n_q_points,
1389  temp1);
1390  eval1.template values<1, false, false>(temp1, temp2);
1391  eval0.template values<0, false, true>(temp2, values_dofs);
1392 
1393  // grad xy
1394  eval2.template values<2, false, false>(hessians_quad +
1395  3 * n_q_points,
1396  temp1);
1397  eval1.template gradients<1, false, false>(temp1, temp2);
1398  eval0.template gradients<0, false, true>(temp2, values_dofs);
1399 
1400  // grad xz
1401  eval2.template gradients<2, false, false>(hessians_quad +
1402  4 * n_q_points,
1403  temp1);
1404  eval1.template values<1, false, false>(temp1, temp2);
1405  eval0.template gradients<0, false, true>(temp2, values_dofs);
1406 
1407  // grad yz
1408  eval2.template gradients<2, false, false>(hessians_quad +
1409  5 * n_q_points,
1410  temp1);
1411  eval1.template gradients<1, false, false>(temp1, temp2);
1412  eval0.template values<0, false, true>(temp2, values_dofs);
1413  }
1414 
1415  break;
1416  default:
1417  AssertThrow(false, ExcNotImplemented());
1418  }
1419  }
1420 
1430  template <EvaluatorVariant variant,
1431  EvaluatorQuantity quantity,
1432  int dim,
1433  int basis_size_1,
1434  int basis_size_2,
1435  typename Number,
1436  typename Number2>
1438  {
1439  static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
1440  "The second dimension must not be smaller than the first");
1441 
1464 #ifndef DEBUG
1466 #endif
1467  static void
1469  const unsigned int n_components,
1470  const AlignedVector<Number2> &transformation_matrix,
1471  const Number * values_in,
1472  Number * values_out,
1473  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
1474  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
1475  {
1476  Assert(
1477  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
1478  ExcMessage("The second dimension must not be smaller than the first"));
1479 
1481 
1482  // we do recursion until dim==1 or dim==2 and we have
1483  // basis_size_1==basis_size_2. The latter optimization increases
1484  // optimization possibilities for the compiler but does only work for
1485  // aliased pointers if the sizes are equal.
1486  constexpr int next_dim =
1487  (dim > 2 ||
1488  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
1489  dim - 1 :
1490  dim;
1491 
1492  EvaluatorTensorProduct<variant,
1493  dim,
1494  basis_size_1,
1495  (basis_size_1 == 0 ? 0 : basis_size_2),
1496  Number,
1497  Number2>
1498  eval_val(transformation_matrix,
1501  basis_size_1_variable,
1502  basis_size_2_variable);
1503  const unsigned int np_1 =
1504  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
1505  const unsigned int np_2 =
1506  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
1507  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
1508  ExcMessage("Cannot transform with 0-point basis"));
1509  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
1510  ExcMessage("Cannot transform with 0-point basis"));
1511 
1512  // run loop backwards to ensure correctness if values_in aliases with
1513  // values_out in case with basis_size_1 < basis_size_2
1514  values_in = values_in + n_components * Utilities::fixed_power<dim>(np_1);
1515  values_out =
1516  values_out + n_components * Utilities::fixed_power<dim>(np_2);
1517  for (unsigned int c = n_components; c != 0; --c)
1518  {
1519  values_in -= Utilities::fixed_power<dim>(np_1);
1520  values_out -= Utilities::fixed_power<dim>(np_2);
1521  if (next_dim < dim)
1522  for (unsigned int q = np_1; q != 0; --q)
1524  variant,
1525  quantity,
1526  next_dim,
1527  basis_size_1,
1528  basis_size_2,
1529  Number,
1530  Number2>::do_forward(1,
1531  transformation_matrix,
1532  values_in +
1533  (q - 1) *
1534  Utilities::fixed_power<next_dim>(np_1),
1535  values_out +
1536  (q - 1) *
1537  Utilities::fixed_power<next_dim>(np_2),
1538  basis_size_1_variable,
1539  basis_size_2_variable);
1540 
1541  // the recursion stops if dim==1 or if dim==2 and
1542  // basis_size_1==basis_size_2 (the latter is used because the
1543  // compiler generates nicer code)
1544  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
1545  {
1546  eval_val.template values<0, true, false>(values_in, values_out);
1547  eval_val.template values<1, true, false>(values_out, values_out);
1548  }
1549  else if (dim == 1)
1550  eval_val.template values<dim - 1, true, false>(values_in,
1551  values_out);
1552  else
1553  eval_val.template values<dim - 1, true, false>(values_out,
1554  values_out);
1555  }
1556  }
1557 
1588 #ifndef DEBUG
1590 #endif
1591  static void
1593  const unsigned int n_components,
1594  const AlignedVector<Number2> &transformation_matrix,
1595  const bool add_into_result,
1596  Number * values_in,
1597  Number * values_out,
1598  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
1599  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
1600  {
1601  Assert(
1602  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
1603  ExcMessage("The second dimension must not be smaller than the first"));
1604  Assert(add_into_result == false || values_in != values_out,
1605  ExcMessage(
1606  "Input and output cannot alias with each other when "
1607  "adding the result of the basis change to existing data"));
1608 
1609  Assert(quantity == EvaluatorQuantity::value ||
1610  quantity == EvaluatorQuantity::hessian,
1611  ExcInternalError());
1612 
1613  constexpr int next_dim =
1614  (dim > 2 ||
1615  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
1616  dim - 1 :
1617  dim;
1618  EvaluatorTensorProduct<variant,
1619  dim,
1620  basis_size_1,
1621  (basis_size_1 == 0 ? 0 : basis_size_2),
1622  Number,
1623  Number2>
1624  eval_val(transformation_matrix,
1625  transformation_matrix,
1626  transformation_matrix,
1627  basis_size_1_variable,
1628  basis_size_2_variable);
1629  const unsigned int np_1 =
1630  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
1631  const unsigned int np_2 =
1632  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
1633  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
1634  ExcMessage("Cannot transform with 0-point basis"));
1635  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
1636  ExcMessage("Cannot transform with 0-point basis"));
1637 
1638  for (unsigned int c = 0; c < n_components; ++c)
1639  {
1640  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
1641  {
1642  if (quantity == EvaluatorQuantity::value)
1643  eval_val.template values<1, false, false>(values_in, values_in);
1644  else
1645  eval_val.template hessians<1, false, false>(values_in,
1646  values_in);
1647 
1648  if (add_into_result)
1649  {
1650  if (quantity == EvaluatorQuantity::value)
1651  eval_val.template values<0, false, true>(values_in,
1652  values_out);
1653  else
1654  eval_val.template hessians<0, false, true>(values_in,
1655  values_out);
1656  }
1657  else
1658  {
1659  if (quantity == EvaluatorQuantity::value)
1660  eval_val.template values<0, false, false>(values_in,
1661  values_out);
1662  else
1663  eval_val.template hessians<0, false, false>(values_in,
1664  values_out);
1665  }
1666  }
1667  else
1668  {
1669  if (dim == 1 && add_into_result)
1670  {
1671  if (quantity == EvaluatorQuantity::value)
1672  eval_val.template values<0, false, true>(values_in,
1673  values_out);
1674  else
1675  eval_val.template hessians<0, false, true>(values_in,
1676  values_out);
1677  }
1678  else if (dim == 1)
1679  {
1680  if (quantity == EvaluatorQuantity::value)
1681  eval_val.template values<0, false, false>(values_in,
1682  values_out);
1683  else
1684  eval_val.template hessians<0, false, false>(values_in,
1685  values_out);
1686  }
1687  else
1688  {
1689  if (quantity == EvaluatorQuantity::value)
1690  eval_val.template values<dim - 1, false, false>(values_in,
1691  values_in);
1692  else
1693  eval_val.template hessians<dim - 1, false, false>(
1694  values_in, values_in);
1695  }
1696  }
1697  if (next_dim < dim)
1698  for (unsigned int q = 0; q < np_1; ++q)
1700  quantity,
1701  next_dim,
1702  basis_size_1,
1703  basis_size_2,
1704  Number,
1705  Number2>::
1706  do_backward(1,
1707  transformation_matrix,
1708  add_into_result,
1709  values_in +
1710  q * Utilities::fixed_power<next_dim>(np_2),
1711  values_out +
1712  q * Utilities::fixed_power<next_dim>(np_1),
1713  basis_size_1_variable,
1714  basis_size_2_variable);
1715 
1716  values_in += Utilities::fixed_power<dim>(np_2);
1717  values_out += Utilities::fixed_power<dim>(np_1);
1718  }
1719  }
1720 
1741  static void
1742  do_mass(const unsigned int n_components,
1743  const AlignedVector<Number2> &transformation_matrix,
1744  const AlignedVector<Number> & coefficients,
1745  const Number * values_in,
1746  Number * scratch_data,
1747  Number * values_out)
1748  {
1749  constexpr int next_dim = dim > 1 ? dim - 1 : dim;
1750  Number * my_scratch =
1751  basis_size_1 != basis_size_2 ? scratch_data : values_out;
1752 
1753  const unsigned int size_per_component = Utilities::pow(basis_size_2, dim);
1754  Assert(coefficients.size() == size_per_component ||
1755  coefficients.size() == n_components * size_per_component,
1756  ExcDimensionMismatch(coefficients.size(), size_per_component));
1757  const unsigned int stride =
1758  coefficients.size() == size_per_component ? 0 : 1;
1759 
1760  for (unsigned int q = basis_size_1; q != 0; --q)
1762  variant,
1764  next_dim,
1765  basis_size_1,
1766  basis_size_2,
1767  Number,
1768  Number2>::do_forward(n_components,
1769  transformation_matrix,
1770  values_in +
1771  (q - 1) *
1772  Utilities::pow(basis_size_1, dim - 1),
1773  my_scratch +
1774  (q - 1) *
1775  Utilities::pow(basis_size_2, dim - 1));
1776  EvaluatorTensorProduct<variant,
1777  dim,
1778  basis_size_1,
1779  basis_size_2,
1780  Number,
1781  Number2>
1782  eval_val(transformation_matrix);
1783  const unsigned int n_inner_blocks =
1784  (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
1785  const unsigned int n_blocks = Utilities::pow(basis_size_2, dim - 1);
1786  for (unsigned int ii = 0; ii < n_blocks; ii += n_inner_blocks)
1787  for (unsigned int c = 0; c < n_components; ++c)
1788  {
1789  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1790  eval_val.template values_one_line<dim - 1, true, false>(
1791  my_scratch + i, my_scratch + i);
1792  for (unsigned int q = 0; q < basis_size_2; ++q)
1793  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1794  my_scratch[i + q * n_blocks + c * size_per_component] *=
1795  coefficients[i + q * n_blocks +
1796  c * stride * size_per_component];
1797  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1798  eval_val.template values_one_line<dim - 1, false, false>(
1799  my_scratch + i, my_scratch + i);
1800  }
1801  for (unsigned int q = 0; q < basis_size_1; ++q)
1803  variant,
1805  next_dim,
1806  basis_size_1,
1807  basis_size_2,
1808  Number,
1809  Number2>::do_backward(n_components,
1810  transformation_matrix,
1811  false,
1812  my_scratch +
1813  q * Utilities::pow(basis_size_2, dim - 1),
1814  values_out +
1815  q * Utilities::pow(basis_size_1, dim - 1));
1816  }
1817  };
1818 
1819 
1820 
1833  template <int dim, int fe_degree, typename Number>
1835  {
1836  static void
1837  evaluate(const unsigned int n_components,
1838  const EvaluationFlags::EvaluationFlags evaluation_flag,
1839  const Number * values_dofs,
1841 
1842  static void
1844  const EvaluationFlags::EvaluationFlags evaluation_flag,
1845  const Number * values_dofs,
1846  Number * gradients_quad,
1847  Number * hessians_quad);
1848 
1849  static void
1850  integrate(const unsigned int n_components,
1851  const EvaluationFlags::EvaluationFlags integration_flag,
1852  Number * values_dofs,
1854  const bool add_into_values_array);
1855 
1856  static void
1858  const EvaluationFlags::EvaluationFlags integration_flag,
1859  Number * values_dofs,
1860  Number * gradients_quad,
1861  const Number * hessians_quad,
1862  const bool add_into_values_array);
1863  };
1864 
1865 
1866 
1867  template <int dim, int fe_degree, typename Number>
1868  inline void
1870  const unsigned int n_components,
1871  const EvaluationFlags::EvaluationFlags evaluation_flag,
1872  const Number * values_dofs,
1874  {
1875  constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1876 
1877  for (unsigned int c = 0; c < n_components; ++c)
1878  {
1879  if ((evaluation_flag & EvaluationFlags::values) != 0u)
1880  for (unsigned int i = 0; i < n_points; ++i)
1881  fe_eval.begin_values()[n_points * c + i] =
1882  values_dofs[n_points * c + i];
1883 
1884  do_evaluate(fe_eval.get_shape_info().data.front(),
1885  evaluation_flag,
1886  values_dofs + c * n_points,
1887  fe_eval.begin_gradients() + c * dim * n_points,
1888  fe_eval.begin_hessians() +
1889  c * dim * (dim + 1) / 2 * n_points);
1890  }
1891  }
1892 
1893 
1894 
1895  template <int dim, int fe_degree, typename Number>
1896  inline void
1899  const EvaluationFlags::EvaluationFlags evaluation_flag,
1900  const Number * values_dofs,
1901  Number * gradients_quad,
1902  Number * hessians_quad)
1903  {
1905  (fe_degree + 2) / 2 * (fe_degree + 1));
1906  constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1907 
1909  dim,
1910  fe_degree + 1,
1911  fe_degree + 1,
1912  Number>
1913  eval(AlignedVector<Number>(),
1916  if ((evaluation_flag &
1918  {
1919  eval.template gradients<0, true, false>(values_dofs, gradients_quad);
1920  if (dim > 1)
1921  eval.template gradients<1, true, false>(values_dofs,
1922  gradients_quad + n_points);
1923  if (dim > 2)
1924  eval.template gradients<2, true, false>(values_dofs,
1925  gradients_quad +
1926  2 * n_points);
1927  }
1928  if (evaluation_flag & EvaluationFlags::hessians)
1929  {
1930  eval.template hessians<0, true, false>(values_dofs, hessians_quad);
1931  if (dim > 1)
1932  {
1933  eval.template gradients<1, true, false>(gradients_quad,
1934  hessians_quad +
1935  dim * n_points);
1936  eval.template hessians<1, true, false>(values_dofs,
1937  hessians_quad + n_points);
1938  }
1939  if (dim > 2)
1940  {
1941  eval.template gradients<2, true, false>(gradients_quad,
1942  hessians_quad +
1943  4 * n_points);
1944  eval.template gradients<2, true, false>(gradients_quad + n_points,
1945  hessians_quad +
1946  5 * n_points);
1947  eval.template hessians<2, true, false>(values_dofs,
1948  hessians_quad +
1949  2 * n_points);
1950  }
1951  }
1952  }
1953 
1954 
1955 
1956  template <int dim, int fe_degree, typename Number>
1957  inline void
1959  const unsigned int n_components,
1960  const EvaluationFlags::EvaluationFlags integration_flag,
1961  Number * values_dofs,
1963  const bool add_into_values_array)
1964  {
1965  constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1966 
1967  for (unsigned int c = 0; c < n_components; ++c)
1968  {
1969  if ((integration_flag & EvaluationFlags::values) != 0u)
1970  {
1971  if (add_into_values_array)
1972  for (unsigned int i = 0; i < n_points; ++i)
1973  values_dofs[n_points * c + i] +=
1974  fe_eval.begin_values()[n_points * c + i];
1975  else
1976  for (unsigned int i = 0; i < n_points; ++i)
1977  values_dofs[n_points * c + i] =
1978  fe_eval.begin_values()[n_points * c + i];
1979  }
1980 
1981  do_integrate(fe_eval.get_shape_info().data.front(),
1982  integration_flag,
1983  values_dofs + c * n_points,
1984  fe_eval.begin_gradients() + c * dim * n_points,
1985  fe_eval.begin_hessians() +
1986  c * dim * (dim + 1) / 2 * n_points,
1987  add_into_values_array ||
1988  ((integration_flag & EvaluationFlags::values) != 0u));
1989  }
1990  }
1991 
1992 
1993 
1994  template <int dim, int fe_degree, typename Number>
1995  inline void
1998  const EvaluationFlags::EvaluationFlags integration_flag,
1999  Number * values_dofs,
2000  Number * gradients_quad,
2001  const Number * hessians_quad,
2002  const bool add_into_values_array)
2003  {
2005  (fe_degree + 2) / 2 * (fe_degree + 1));
2006 
2008  dim,
2009  fe_degree + 1,
2010  fe_degree + 1,
2011  Number>
2012  eval(AlignedVector<Number>(),
2015  constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
2016 
2017  if ((integration_flag & EvaluationFlags::hessians) != 0u)
2018  {
2019  // diagonal
2020  // grad xx
2021  if (add_into_values_array == true)
2022  eval.template hessians<0, false, true>(hessians_quad, values_dofs);
2023  else
2024  eval.template hessians<0, false, false>(hessians_quad, values_dofs);
2025  // grad yy
2026  if (dim > 1)
2027  eval.template hessians<1, false, true>(hessians_quad + n_points,
2028  values_dofs);
2029  // grad zz
2030  if (dim > 2)
2031  eval.template hessians<2, false, true>(hessians_quad + 2 * n_points,
2032  values_dofs);
2033  // off-diagonal
2034  if (dim == 2)
2035  {
2036  // grad xy, queue into gradient
2037  if (integration_flag & EvaluationFlags::gradients)
2038  eval.template gradients<1, false, true>(hessians_quad +
2039  2 * n_points,
2040  gradients_quad);
2041  else
2042  eval.template gradients<1, false, false>(hessians_quad +
2043  2 * n_points,
2044  gradients_quad);
2045  }
2046  if (dim == 3)
2047  {
2048  // grad xy, queue into gradient
2049  if (integration_flag & EvaluationFlags::gradients)
2050  eval.template gradients<1, false, true>(hessians_quad +
2051  3 * n_points,
2052  gradients_quad);
2053  else
2054  eval.template gradients<1, false, false>(hessians_quad +
2055  3 * n_points,
2056  gradients_quad);
2057 
2058  // grad xz
2059  eval.template gradients<2, false, true>(hessians_quad +
2060  4 * n_points,
2061  gradients_quad);
2062 
2063  // grad yz
2064  if (integration_flag & EvaluationFlags::gradients)
2065  eval.template gradients<2, false, true>(
2066  hessians_quad + 5 * n_points, gradients_quad + n_points);
2067  else
2068  eval.template gradients<2, false, false>(
2069  hessians_quad + 5 * n_points, gradients_quad + n_points);
2070  }
2071 
2072  // if we did not integrate gradients, set the last slot to zero
2073  // which was not touched before, in order to avoid the if
2074  // statement in the gradients loop below
2075  if ((integration_flag & EvaluationFlags::gradients) == 0u)
2076  for (unsigned int q = 0; q < n_points; ++q)
2077  gradients_quad[(dim - 1) * n_points + q] = Number();
2078  }
2079 
2080  if ((integration_flag &
2082  {
2083  if (add_into_values_array ||
2084  (integration_flag & EvaluationFlags::hessians) != 0u)
2085  eval.template gradients<0, false, true>(gradients_quad, values_dofs);
2086  else
2087  eval.template gradients<0, false, false>(gradients_quad, values_dofs);
2088  if (dim > 1)
2089  eval.template gradients<1, false, true>(gradients_quad + n_points,
2090  values_dofs);
2091  if (dim > 2)
2092  eval.template gradients<2, false, true>(gradients_quad + 2 * n_points,
2093  values_dofs);
2094  }
2095  }
2096 
2097 
2098 
2109  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2111  {
2112  static void
2113  evaluate(const unsigned int n_components,
2114  const EvaluationFlags::EvaluationFlags evaluation_flag,
2115  const Number * values_dofs,
2117 
2118  static void
2119  integrate(const unsigned int n_components,
2120  const EvaluationFlags::EvaluationFlags evaluation_flag,
2121  Number * values_dofs,
2123  const bool add_into_values_array);
2124  };
2125 
2126 
2127 
2128  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2129  inline void
2131  dim,
2132  fe_degree,
2133  n_q_points_1d,
2134  Number>::evaluate(const unsigned int n_components,
2135  const EvaluationFlags::EvaluationFlags evaluation_flag,
2136  const Number * values_dofs,
2138  {
2139  const auto &shape_data = fe_eval.get_shape_info().data.front();
2140 
2141  Assert(n_q_points_1d > fe_degree,
2142  ExcMessage("You lose information when going to a collocation space "
2143  "of lower degree, so the evaluation results would be "
2144  "wrong. Thus, this class does not permit the desired "
2145  "operation."));
2146  constexpr std::size_t n_dofs = Utilities::pow(fe_degree + 1, dim);
2147  constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
2148 
2149  for (unsigned int c = 0; c < n_components; ++c)
2150  {
2154  dim,
2155  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
2156  n_q_points_1d,
2157  Number,
2158  Number>::do_forward(1,
2159  shape_data.shape_values_eo,
2160  values_dofs + c * n_dofs,
2161  fe_eval.begin_values() + c * n_q_points);
2162 
2163  // apply derivatives in the collocation space
2164  if (evaluation_flag &
2167  do_evaluate(shape_data,
2168  evaluation_flag & (EvaluationFlags::gradients |
2170  fe_eval.begin_values() + c * n_q_points,
2171  fe_eval.begin_gradients() + c * dim * n_q_points,
2172  fe_eval.begin_hessians() +
2173  c * dim * (dim + 1) / 2 * n_q_points);
2174  }
2175  }
2176 
2177 
2178 
2179  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2180  inline void
2182  dim,
2183  fe_degree,
2184  n_q_points_1d,
2185  Number>::integrate(const unsigned int n_components,
2186  const EvaluationFlags::EvaluationFlags integration_flag,
2187  Number * values_dofs,
2189  const bool add_into_values_array)
2190  {
2191  const auto &shape_data = fe_eval.get_shape_info().data.front();
2192 
2193  Assert(n_q_points_1d > fe_degree,
2194  ExcMessage("You lose information when going to a collocation space "
2195  "of lower degree, so the evaluation results would be "
2196  "wrong. Thus, this class does not permit the desired "
2197  "operation."));
2198  constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
2199 
2200  for (unsigned int c = 0; c < n_components; ++c)
2201  {
2202  // apply derivatives in collocation space
2203  if (integration_flag &
2206  do_integrate(shape_data,
2207  integration_flag & (EvaluationFlags::gradients |
2209  fe_eval.begin_values() + c * n_q_points,
2210  fe_eval.begin_gradients() + c * dim * n_q_points,
2211  fe_eval.begin_hessians() +
2212  c * dim * (dim + 1) / 2 * n_q_points,
2213  /*add_into_values_array=*/
2214  integration_flag & EvaluationFlags::values);
2215 
2216  // transform back to the original space
2220  dim,
2221  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
2222  n_q_points_1d,
2223  Number,
2224  Number>::do_backward(1,
2225  shape_data.shape_values_eo,
2226  add_into_values_array,
2227  fe_eval.begin_values() + c * n_q_points,
2228  values_dofs +
2229  c * Utilities::pow(fe_degree + 1, dim));
2230  }
2231  }
2232 
2233 
2234 
2242  constexpr bool
2243  use_collocation_evaluation(const unsigned int fe_degree,
2244  const unsigned int n_q_points_1d)
2245  {
2246  return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
2247  (n_q_points_1d <= 3 * fe_degree / 2 + 1);
2248  }
2249 
2250 
2266  template <int dim, typename Number>
2268  {
2269  template <int fe_degree, int n_q_points_1d>
2270  static bool
2271  run(const unsigned int n_components,
2272  const EvaluationFlags::EvaluationFlags evaluation_flag,
2273  const Number * values_dofs,
2275  {
2276  const auto element_type = fe_eval.get_shape_info().element_type;
2278 
2279  Assert(fe_eval.get_shape_info().data.size() == 1 ||
2280  (fe_eval.get_shape_info().data.size() == dim &&
2281  element_type == ElementType::tensor_general) ||
2282  element_type == ElementType::tensor_raviart_thomas,
2283  ExcNotImplemented());
2284 
2285  if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
2287  {
2289  n_components, evaluation_flag, values_dofs, fe_eval);
2290  }
2291  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
2292  // shape_info.h for more details
2293  else if (fe_degree >= 0 &&
2294  use_collocation_evaluation(fe_degree, n_q_points_1d) &&
2295  element_type <= ElementType::tensor_symmetric)
2296  {
2298  dim,
2299  fe_degree,
2300  n_q_points_1d,
2301  Number>::evaluate(n_components,
2302  evaluation_flag,
2303  values_dofs,
2304  fe_eval);
2305  }
2306  else if (fe_degree >= 0 &&
2308  {
2310  dim,
2311  fe_degree,
2312  n_q_points_1d,
2313  Number>::evaluate(n_components,
2314  evaluation_flag,
2315  values_dofs,
2316  fe_eval);
2317  }
2318  else if (element_type == ElementType::tensor_symmetric_plus_dg0)
2319  {
2321  dim,
2322  fe_degree,
2323  n_q_points_1d,
2324  Number>::evaluate(n_components,
2325  evaluation_flag,
2326  values_dofs,
2327  fe_eval);
2328  }
2329  else if (element_type == ElementType::truncated_tensor)
2330  {
2332  dim,
2333  fe_degree,
2334  n_q_points_1d,
2335  Number>::evaluate(n_components,
2336  evaluation_flag,
2337  values_dofs,
2338  fe_eval);
2339  }
2340  else if (element_type == ElementType::tensor_none)
2341  {
2343  dim,
2344  fe_degree,
2345  n_q_points_1d,
2346  Number>::evaluate(n_components,
2347  evaluation_flag,
2348  values_dofs,
2349  fe_eval);
2350  }
2351  else if (element_type == ElementType::tensor_raviart_thomas)
2352  {
2355  dim,
2356  (fe_degree == -1) ? 1 : fe_degree,
2357  (n_q_points_1d < 1) ? 1 : n_q_points_1d,
2358  Number>::template evaluate_or_integrate<false>(evaluation_flag,
2359  const_cast<Number *>(
2360  values_dofs),
2361  fe_eval);
2362  }
2363  else
2364  {
2366  dim,
2367  fe_degree,
2368  n_q_points_1d,
2369  Number>::evaluate(n_components,
2370  evaluation_flag,
2371  values_dofs,
2372  fe_eval);
2373  }
2374 
2375  return false;
2376  }
2377  };
2378 
2379 
2380 
2396  template <int dim, typename Number>
2398  {
2399  template <int fe_degree, int n_q_points_1d>
2400  static bool
2401  run(const unsigned int n_components,
2402  const EvaluationFlags::EvaluationFlags integration_flag,
2403  Number * values_dofs,
2405  const bool sum_into_values_array)
2406  {
2407  const auto element_type = fe_eval.get_shape_info().element_type;
2409 
2410  Assert(fe_eval.get_shape_info().data.size() == 1 ||
2411  (fe_eval.get_shape_info().data.size() == dim &&
2412  element_type == ElementType::tensor_general) ||
2413  element_type == ElementType::tensor_raviart_thomas,
2414  ExcNotImplemented());
2415 
2416  if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
2418  {
2420  n_components,
2421  integration_flag,
2422  values_dofs,
2423  fe_eval,
2424  sum_into_values_array);
2425  }
2426  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
2427  // shape_info.h for more details
2428  else if (fe_degree >= 0 &&
2429  use_collocation_evaluation(fe_degree, n_q_points_1d) &&
2430  element_type <= ElementType::tensor_symmetric)
2431  {
2433  dim,
2434  fe_degree,
2435  n_q_points_1d,
2436  Number>::integrate(n_components,
2437  integration_flag,
2438  values_dofs,
2439  fe_eval,
2440  sum_into_values_array);
2441  }
2442  else if (fe_degree >= 0 &&
2444  {
2446  dim,
2447  fe_degree,
2448  n_q_points_1d,
2449  Number>::integrate(n_components,
2450  integration_flag,
2451  values_dofs,
2452  fe_eval,
2453  sum_into_values_array);
2454  }
2455  else if (element_type == ElementType::tensor_symmetric_plus_dg0)
2456  {
2458  dim,
2459  fe_degree,
2460  n_q_points_1d,
2461  Number>::integrate(n_components,
2462  integration_flag,
2463  values_dofs,
2464  fe_eval,
2465  sum_into_values_array);
2466  }
2467  else if (element_type == ElementType::truncated_tensor)
2468  {
2470  dim,
2471  fe_degree,
2472  n_q_points_1d,
2473  Number>::integrate(n_components,
2474  integration_flag,
2475  values_dofs,
2476  fe_eval,
2477  sum_into_values_array);
2478  }
2479  else if (element_type == ElementType::tensor_none)
2480  {
2482  dim,
2483  fe_degree,
2484  n_q_points_1d,
2485  Number>::integrate(n_components,
2486  integration_flag,
2487  values_dofs,
2488  fe_eval,
2489  sum_into_values_array);
2490  }
2491  else if (element_type == ElementType::tensor_raviart_thomas)
2492  {
2494  dim,
2495  (fe_degree == -1) ? 1 : fe_degree,
2496  (n_q_points_1d < 1) ? 1 : n_q_points_1d,
2497  Number>::
2498  template evaluate_or_integrate<true>(integration_flag,
2499  const_cast<Number *>(
2500  values_dofs),
2501  fe_eval,
2502  sum_into_values_array);
2503  }
2504  else
2505  {
2507  dim,
2508  fe_degree,
2509  n_q_points_1d,
2510  Number>::integrate(n_components,
2511  integration_flag,
2512  values_dofs,
2513  fe_eval,
2514  sum_into_values_array);
2515  }
2516 
2517  return false;
2518  }
2519  };
2520 
2521 
2522 
2523  template <bool symmetric_evaluate,
2524  int dim,
2525  int fe_degree,
2526  int n_q_points_1d,
2527  typename Number>
2529  {
2530  // We enable a transformation to collocation for derivatives if it gives
2531  // correct results (first two conditions), if it is the most efficient
2532  // choice in terms of operation counts (third condition) and if we were
2533  // able to initialize the fields in shape_info.templates.h from the
2534  // polynomials (fourth condition).
2535  using Eval = EvaluatorTensorProduct<symmetric_evaluate ? evaluate_evenodd :
2537  dim - 1,
2538  fe_degree + 1,
2539  n_q_points_1d,
2540  Number>;
2541 
2542  static Eval
2545  const unsigned int subface_index,
2546  const unsigned int direction)
2547  {
2548  if (symmetric_evaluate)
2549  return Eval(data.shape_values_eo,
2550  data.shape_gradients_eo,
2551  data.shape_hessians_eo,
2552  data.fe_degree + 1,
2553  data.n_q_points_1d);
2554  else if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
2555  return Eval(data.shape_values,
2556  data.shape_gradients,
2557  data.shape_hessians,
2558  data.fe_degree + 1,
2559  data.n_q_points_1d);
2560  else
2561  {
2562  const unsigned int index =
2563  direction == 0 ? subface_index % 2 : subface_index / 2;
2564  return Eval(data.values_within_subface[index],
2567  data.fe_degree + 1,
2568  data.n_q_points_1d);
2569  }
2570  }
2571 
2572  static void
2574  const unsigned int n_components,
2575  const EvaluationFlags::EvaluationFlags evaluation_flag,
2577  Number * values_dofs,
2578  Number * values_quad,
2579  Number * gradients_quad,
2580  Number * hessians_quad,
2581  Number * scratch_data,
2582  const unsigned int subface_index)
2583  {
2584  Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
2585  Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
2586 
2587  const std::size_t n_dofs = fe_degree > -1 ?
2588  Utilities::pow(fe_degree + 1, dim - 1) :
2589  Utilities::pow(data.fe_degree + 1, dim - 1);
2590  const std::size_t n_q_points =
2591  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
2592  Utilities::pow(data.n_q_points_1d, dim - 1);
2593 
2594  // keep a copy of the original pointer for the case of the Hessians
2595  Number *values_dofs_ptr = values_dofs;
2596 
2597  if ((evaluation_flag & EvaluationFlags::values) != 0u &&
2598  ((evaluation_flag & EvaluationFlags::gradients) == 0u))
2599  for (unsigned int c = 0; c < n_components; ++c)
2600  {
2601  switch (dim)
2602  {
2603  case 3:
2604  eval0.template values<0, true, false>(values_dofs,
2605  values_quad);
2606  eval1.template values<1, true, false>(values_quad,
2607  values_quad);
2608  break;
2609  case 2:
2610  eval0.template values<0, true, false>(values_dofs,
2611  values_quad);
2612  break;
2613  case 1:
2614  values_quad[0] = values_dofs[0];
2615  break;
2616  default:
2617  Assert(false, ExcNotImplemented());
2618  }
2619  // Note: we always keep storage of values, 1st and 2nd derivatives
2620  // in an array
2621  values_dofs += 3 * n_dofs;
2622  values_quad += n_q_points;
2623  }
2624  else if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2625  for (unsigned int c = 0; c < n_components; ++c)
2626  {
2627  switch (dim)
2628  {
2629  case 3:
2630  if (symmetric_evaluate &&
2631  use_collocation_evaluation(fe_degree, n_q_points_1d))
2632  {
2633  eval0.template values<0, true, false>(values_dofs,
2634  values_quad);
2635  eval0.template values<1, true, false>(values_quad,
2636  values_quad);
2638  dim - 1,
2639  n_q_points_1d,
2640  n_q_points_1d,
2641  Number>
2642  eval_grad(AlignedVector<Number>(),
2645  eval_grad.template gradients<0, true, false>(
2646  values_quad, gradients_quad);
2647  eval_grad.template gradients<1, true, false>(
2648  values_quad, gradients_quad + n_q_points);
2649  }
2650  else
2651  {
2652  // grad x
2653  eval0.template gradients<0, true, false>(values_dofs,
2654  scratch_data);
2655  eval1.template values<1, true, false>(scratch_data,
2656  gradients_quad);
2657 
2658  // grad y
2659  eval0.template values<0, true, false>(values_dofs,
2660  scratch_data);
2661  eval1.template gradients<1, true, false>(scratch_data,
2662  gradients_quad +
2663  n_q_points);
2664 
2665  if ((evaluation_flag & EvaluationFlags::values) != 0u)
2666  eval1.template values<1, true, false>(scratch_data,
2667  values_quad);
2668  }
2669  // grad z
2670  eval0.template values<0, true, false>(values_dofs + n_dofs,
2671  scratch_data);
2672  eval1.template values<1, true, false>(
2673  scratch_data, gradients_quad + (dim - 1) * n_q_points);
2674 
2675  break;
2676  case 2:
2677  eval0.template values<0, true, false>(values_dofs + n_dofs,
2678  gradients_quad +
2679  n_q_points);
2680  eval0.template gradients<0, true, false>(values_dofs,
2681  gradients_quad);
2682  if ((evaluation_flag & EvaluationFlags::values) != 0u)
2683  eval0.template values<0, true, false>(values_dofs,
2684  values_quad);
2685  break;
2686  case 1:
2687  values_quad[0] = values_dofs[0];
2688  gradients_quad[0] = values_dofs[1];
2689  break;
2690  default:
2691  AssertThrow(false, ExcNotImplemented());
2692  }
2693  values_dofs += 3 * n_dofs;
2694  values_quad += n_q_points;
2695  gradients_quad += dim * n_q_points;
2696  }
2697 
2698  if ((evaluation_flag & EvaluationFlags::hessians) != 0u)
2699  {
2700  values_dofs = values_dofs_ptr;
2701  for (unsigned int c = 0; c < n_components; ++c)
2702  {
2703  switch (dim)
2704  {
2705  case 3:
2706  // grad xx
2707  eval0.template hessians<0, true, false>(values_dofs,
2708  scratch_data);
2709  eval1.template values<1, true, false>(scratch_data,
2710  hessians_quad);
2711 
2712  // grad yy
2713  eval0.template values<0, true, false>(values_dofs,
2714  scratch_data);
2715  eval1.template hessians<1, true, false>(scratch_data,
2716  hessians_quad +
2717  n_q_points);
2718 
2719  // grad zz
2720  eval0.template values<0, true, false>(values_dofs +
2721  2 * n_dofs,
2722  scratch_data);
2723  eval1.template values<1, true, false>(scratch_data,
2724  hessians_quad +
2725  2 * n_q_points);
2726 
2727  // grad xy
2728  eval0.template gradients<0, true, false>(values_dofs,
2729  scratch_data);
2730  eval1.template gradients<1, true, false>(scratch_data,
2731  hessians_quad +
2732  3 * n_q_points);
2733 
2734  // grad xz
2735  eval0.template gradients<0, true, false>(values_dofs +
2736  n_dofs,
2737  scratch_data);
2738  eval1.template values<1, true, false>(scratch_data,
2739  hessians_quad +
2740  4 * n_q_points);
2741 
2742  // grad yz
2743  eval0.template values<0, true, false>(values_dofs + n_dofs,
2744  scratch_data);
2745  eval1.template gradients<1, true, false>(scratch_data,
2746  hessians_quad +
2747  5 * n_q_points);
2748 
2749  break;
2750  case 2:
2751  // grad xx
2752  eval0.template hessians<0, true, false>(values_dofs,
2753  hessians_quad);
2754  // grad yy
2755  eval0.template values<0, true, false>(
2756  values_dofs + 2 * n_dofs, hessians_quad + n_q_points);
2757  // grad xy
2758  eval0.template gradients<0, true, false>(
2759  values_dofs + n_dofs, hessians_quad + 2 * n_q_points);
2760  break;
2761  case 1:
2762  hessians_quad[0] = values_dofs[2];
2763  break;
2764  default:
2765  AssertThrow(false, ExcNotImplemented());
2766  }
2767  values_dofs += 3 * n_dofs;
2768  hessians_quad += dim * (dim + 1) / 2 * n_q_points;
2769  }
2770  }
2771  }
2772 
2773  static void
2775  const unsigned int n_components,
2776  const EvaluationFlags::EvaluationFlags integration_flag,
2778  Number * values_dofs,
2779  Number * values_quad,
2780  Number * gradients_quad,
2781  Number * hessians_quad,
2782  Number * scratch_data,
2783  const unsigned int subface_index)
2784  {
2785  Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
2786  Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
2787 
2788  const std::size_t n_dofs =
2789  fe_degree > -1 ?
2790  Utilities::pow(fe_degree + 1, dim - 1) :
2791  (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
2792  const std::size_t n_q_points =
2793  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
2794  Utilities::pow(data.n_q_points_1d, dim - 1);
2795 
2796  // keep a copy of the original pointer for the case of the Hessians
2797  Number *values_dofs_ptr = values_dofs;
2798 
2799  if ((integration_flag & EvaluationFlags::values) != 0u &&
2800  (integration_flag & EvaluationFlags::gradients) == 0u)
2801  for (unsigned int c = 0; c < n_components; ++c)
2802  {
2803  switch (dim)
2804  {
2805  case 3:
2806  eval1.template values<1, false, false>(values_quad,
2807  values_quad);
2808  eval0.template values<0, false, false>(values_quad,
2809  values_dofs);
2810  break;
2811  case 2:
2812  eval0.template values<0, false, false>(values_quad,
2813  values_dofs);
2814  break;
2815  case 1:
2816  values_dofs[0] = values_quad[0];
2817  break;
2818  default:
2819  Assert(false, ExcNotImplemented());
2820  }
2821  values_dofs += 3 * n_dofs;
2822  values_quad += n_q_points;
2823  }
2824  else if ((integration_flag & EvaluationFlags::gradients) != 0u)
2825  for (unsigned int c = 0; c < n_components; ++c)
2826  {
2827  switch (dim)
2828  {
2829  case 3:
2830  // grad z
2831  eval1.template values<1, false, false>(gradients_quad +
2832  2 * n_q_points,
2833  gradients_quad +
2834  2 * n_q_points);
2835  eval0.template values<0, false, false>(gradients_quad +
2836  2 * n_q_points,
2837  values_dofs + n_dofs);
2838  if (symmetric_evaluate &&
2839  use_collocation_evaluation(fe_degree, n_q_points_1d))
2840  {
2842  dim - 1,
2843  n_q_points_1d,
2844  n_q_points_1d,
2845  Number>
2846  eval_grad(AlignedVector<Number>(),
2849  if ((integration_flag & EvaluationFlags::values) != 0u)
2850  eval_grad.template gradients<1, false, true>(
2851  gradients_quad + n_q_points, values_quad);
2852  else
2853  eval_grad.template gradients<1, false, false>(
2854  gradients_quad + n_q_points, values_quad);
2855  eval_grad.template gradients<0, false, true>(
2856  gradients_quad, values_quad);
2857  eval0.template values<1, false, false>(values_quad,
2858  values_quad);
2859  eval0.template values<0, false, false>(values_quad,
2860  values_dofs);
2861  }
2862  else
2863  {
2864  if ((integration_flag & EvaluationFlags::values) != 0u)
2865  {
2866  eval1.template values<1, false, false>(values_quad,
2867  scratch_data);
2868  eval1.template gradients<1, false, true>(
2869  gradients_quad + n_q_points, scratch_data);
2870  }
2871  else
2872  eval1.template gradients<1, false, false>(
2873  gradients_quad + n_q_points, scratch_data);
2874 
2875  // grad y
2876  eval0.template values<0, false, false>(scratch_data,
2877  values_dofs);
2878 
2879  // grad x
2880  eval1.template values<1, false, false>(gradients_quad,
2881  scratch_data);
2882  eval0.template gradients<0, false, true>(scratch_data,
2883  values_dofs);
2884  }
2885  break;
2886  case 2:
2887  eval0.template values<0, false, false>(gradients_quad +
2888  n_q_points,
2889  values_dofs + n_dofs);
2890  eval0.template gradients<0, false, false>(gradients_quad,
2891  values_dofs);
2892  if ((integration_flag & EvaluationFlags::values) != 0u)
2893  eval0.template values<0, false, true>(values_quad,
2894  values_dofs);
2895  break;
2896  case 1:
2897  values_dofs[0] = values_quad[0];
2898  values_dofs[1] = gradients_quad[0];
2899  break;
2900  default:
2901  AssertThrow(false, ExcNotImplemented());
2902  }
2903  values_dofs += 3 * n_dofs;
2904  values_quad += n_q_points;
2905  gradients_quad += dim * n_q_points;
2906  }
2907 
2908  if ((integration_flag & EvaluationFlags::hessians) != 0u)
2909  {
2910  values_dofs = values_dofs_ptr;
2911  for (unsigned int c = 0; c < n_components; ++c)
2912  {
2913  switch (dim)
2914  {
2915  case 3:
2916  // grad xx
2917  eval1.template values<1, false, false>(hessians_quad,
2918  scratch_data);
2919  if ((integration_flag & (EvaluationFlags::values |
2921  eval0.template hessians<0, false, true>(scratch_data,
2922  values_dofs);
2923  else
2924  eval0.template hessians<0, false, false>(scratch_data,
2925  values_dofs);
2926 
2927  // grad yy
2928  eval1.template hessians<1, false, false>(hessians_quad +
2929  n_q_points,
2930  scratch_data);
2931  eval0.template values<0, false, true>(scratch_data,
2932  values_dofs);
2933 
2934  // grad zz
2935  eval1.template values<1, false, false>(hessians_quad +
2936  2 * n_q_points,
2937  scratch_data);
2938  eval0.template values<0, false, false>(scratch_data,
2939  values_dofs +
2940  2 * n_dofs);
2941 
2942  // grad xy
2943  eval1.template gradients<1, false, false>(hessians_quad +
2944  3 * n_q_points,
2945  scratch_data);
2946  eval0.template gradients<0, false, true>(scratch_data,
2947  values_dofs);
2948 
2949  // grad xz
2950  eval1.template values<1, false, false>(hessians_quad +
2951  4 * n_q_points,
2952  scratch_data);
2953  if ((integration_flag & EvaluationFlags::gradients) != 0u)
2954  eval0.template gradients<0, false, true>(scratch_data,
2955  values_dofs +
2956  n_dofs);
2957  else
2958  eval0.template gradients<0, false, false>(scratch_data,
2959  values_dofs +
2960  n_dofs);
2961 
2962  // grad yz
2963  eval1.template gradients<1, false, false>(hessians_quad +
2964  5 * n_q_points,
2965  scratch_data);
2966  eval0.template values<0, false, true>(scratch_data,
2967  values_dofs + n_dofs);
2968 
2969  break;
2970  case 2:
2971  // grad xx
2972  if ((integration_flag & (EvaluationFlags::values |
2974  eval0.template hessians<0, false, true>(hessians_quad,
2975  values_dofs);
2976  else
2977  eval0.template hessians<0, false, false>(hessians_quad,
2978  values_dofs);
2979 
2980  // grad yy
2981  eval0.template values<0, false, false>(
2982  hessians_quad + n_q_points, values_dofs + 2 * n_dofs);
2983  // grad xy
2984  if ((integration_flag & EvaluationFlags::gradients) != 0u)
2985  eval0.template gradients<0, false, true>(
2986  hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
2987  else
2988  eval0.template gradients<0, false, false>(
2989  hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
2990  break;
2991  case 1:
2992  values_dofs[2] = hessians_quad[0];
2993  if ((integration_flag & EvaluationFlags::values) == 0u)
2994  values_dofs[0] = 0;
2995  if ((integration_flag & EvaluationFlags::gradients) == 0u)
2996  values_dofs[1] = 0;
2997  break;
2998  default:
2999  AssertThrow(false, ExcNotImplemented());
3000  }
3001  values_dofs += 3 * n_dofs;
3002  hessians_quad += dim * (dim + 1) / 2 * n_q_points;
3003  }
3004  }
3005  }
3006  };
3007 
3008  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
3010  {
3012  dim - 1,
3013  fe_degree,
3014  n_q_points_1d,
3015  Number>;
3016  template <typename EvalType>
3017  static EvalType
3020  const unsigned int subface_index,
3021  const unsigned int direction)
3022  {
3023  if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
3024  return EvalType(data.shape_values,
3025  data.shape_gradients,
3026  data.shape_hessians);
3027  else
3028  {
3029  const unsigned int index =
3030  direction == 0 ? subface_index % 2 : subface_index / 2;
3031  return EvalType(data.values_within_subface[index],
3034  }
3035  }
3036 
3037  template <bool integrate>
3038  static void
3040  const EvaluationFlags::EvaluationFlags evaluation_flag,
3041  Number * values_dofs,
3043  Number * scratch_data,
3044  const unsigned int subface_index,
3045  const unsigned int face_no)
3046  {
3047  const unsigned int face_direction = face_no / 2;
3048 
3049  // We first evaluate the anisotropic faces, i.e the faces where
3050  // face_direction != component. Note that the call order here is not
3051  // important, since the pointers are shifted accordingly within the
3052  // function. However, this is the order in which the components will be in
3053  // the quadrature points. Furthermore, the isotropic faces have no "normal
3054  // direction" but we still pass in normal_dir = 2 since this is used for
3055  // the pointers.
3056  // -----------------------------------------------------------------------------------
3057  // | | Anisotropic faces | Isotropic faces|
3058  // | Face dir | comp, coords, normal_dir | comp, coords, normal_dir | comp, coords |
3059  // | --------------------------------------------------------------------------------|
3060  // | 0 | 1, y, 0 | - | 0, y |
3061  // | 1 | 0, x, 0 | - | 1, x |
3062  // | --------------------------------------------------------------------------------|
3063  // | 0 | 1, yz, 0 | 2, yz, 1 | 0, yz |
3064  // | 1 | 2, zx, 0 | 0, zx, 1 | 1, zx |
3065  // | 2 | 0, xy, 0 | 1, xy, 1 | 2, xy |
3066  // -----------------------------------------------------------------------------------
3067  evaluate_in_face_apply<0>(values_dofs,
3068  fe_eval,
3069  scratch_data,
3070  evaluation_flag,
3071  face_direction,
3072  subface_index,
3073  std::integral_constant<bool, integrate>());
3074 
3075  if (dim == 3)
3076  evaluate_in_face_apply<1>(values_dofs,
3077  fe_eval,
3078  scratch_data,
3079  evaluation_flag,
3080  face_direction,
3081  subface_index,
3082  std::integral_constant<bool, integrate>());
3083 
3084  evaluate_in_face_apply<2>(values_dofs,
3085  fe_eval,
3086  scratch_data,
3087  evaluation_flag,
3088  face_direction,
3089  subface_index,
3090  std::integral_constant<bool, integrate>());
3091  }
3092 
3093  /*
3094  * Helper function which applies the 1D kernels for on one
3095  * component in a face. normal_dir indicates the direction of the continuous
3096  * component of the RT space. std::integral_constant<bool, false> is the
3097  * evaluation path, and std::integral_constant<bool, true> below is the
3098  * integration path. These two functions can be fused together since all
3099  * offsets and pointers are the exact same.
3100  */
3101  template <int normal_dir>
3102  static inline void
3104  Number * values_dofs,
3106  Number * scratch_data,
3107  const EvaluationFlags::EvaluationFlags evaluation_flag,
3108  const unsigned int face_direction,
3109  const unsigned int subface_index,
3110  std::integral_constant<bool, false>)
3111  {
3112  using EvalNormal =
3114  dim - 1,
3115  (fe_degree == -1) ? 1 : fe_degree + 1,
3116  n_q_points_1d,
3117  Number,
3118  normal_dir>;
3119  using EvalTangent =
3121  dim - 1,
3122  (fe_degree == -1) ? 1 : fe_degree,
3123  n_q_points_1d,
3124  Number,
3125  normal_dir>;
3126 
3127  using TempEval0 = typename std::
3128  conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
3129  using TempEval1 = typename std::
3130  conditional<normal_dir == 0, EvalTangent, EvalNormal>::type;
3131  using Eval0 = typename std::
3132  conditional<normal_dir == 2, EvalGeneral, TempEval0>::type;
3133  using Eval1 = typename std::
3134  conditional<normal_dir == 2, EvalGeneral, TempEval1>::type;
3135 
3136  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info =
3137  fe_eval.get_shape_info();
3138  Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
3139  ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]),
3140  subface_index,
3141  0);
3142  Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
3143  ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]),
3144  subface_index,
3145  1);
3146 
3147  constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
3148  const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face;
3149  const std::size_t n_dofs_normal =
3150  n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
3151  const std::size_t dofs_stride =
3152  (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
3153  n_dofs_tangent;
3154 
3155  static constexpr ::ndarray<unsigned int, 3, 3> component_table = {
3156  {{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
3157  const unsigned int component =
3158  (dim == 2 && normal_dir == 0 && face_direction == 1) ?
3159  0 :
3160  component_table[face_direction][normal_dir];
3161 
3162  // Initial offsets
3163  values_dofs +=
3164  3 * ((component == 0) ?
3165  0 :
3166  ((component == 1) ?
3167  ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) :
3168  ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent :
3169  n_dofs_normal + n_dofs_tangent)));
3170  const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir;
3171  Number *values_quad = fe_eval.begin_values() + n_q_points * shift;
3172  Number *gradients_quad =
3173  fe_eval.begin_gradients() + dim * n_q_points * shift;
3174  Number *hessians_quad =
3175  fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * shift;
3176 
3177  // Evaluation path
3178  if ((evaluation_flag & EvaluationFlags::values) &&
3179  !(evaluation_flag & EvaluationFlags::gradients))
3180  {
3181  switch (dim)
3182  {
3183  case 3:
3184  eval0.template values<0, true, false>(values_dofs, values_quad);
3185  eval1.template values<1, true, false>(values_quad, values_quad);
3186  break;
3187  case 2:
3188  eval0.template values<0, true, false>(values_dofs, values_quad);
3189  break;
3190  default:
3191  Assert(false, ExcNotImplemented());
3192  }
3193  }
3194  else if (evaluation_flag & EvaluationFlags::gradients)
3195  {
3196  switch (dim)
3197  {
3198  case 3:
3199  // grad x
3200  eval0.template gradients<0, true, false>(values_dofs,
3201  scratch_data);
3202  eval1.template values<1, true, false>(scratch_data,
3203  gradients_quad);
3204 
3205  // grad y
3206  eval0.template values<0, true, false>(values_dofs,
3207  scratch_data);
3208  eval1.template gradients<1, true, false>(scratch_data,
3209  gradients_quad +
3210  n_q_points);
3211 
3212  if (evaluation_flag & EvaluationFlags::values)
3213  eval1.template values<1, true, false>(scratch_data,
3214  values_quad);
3215 
3216  // grad z
3217  eval0.template values<0, true, false>(values_dofs + dofs_stride,
3218  scratch_data);
3219  eval1.template values<1, true, false>(scratch_data,
3220  gradients_quad +
3221  2 * n_q_points);
3222 
3223  break;
3224  case 2:
3225  eval0.template values<0, true, false>(values_dofs + dofs_stride,
3226  gradients_quad +
3227  n_q_points);
3228  eval0.template gradients<0, true, false>(values_dofs,
3229  gradients_quad);
3230  if ((evaluation_flag & EvaluationFlags::values))
3231  eval0.template values<0, true, false>(values_dofs,
3232  values_quad);
3233  break;
3234  default:
3235  AssertThrow(false, ExcNotImplemented());
3236  }
3237  }
3238 
3239  if (evaluation_flag & EvaluationFlags::hessians)
3240  {
3241  switch (dim)
3242  {
3243  case 3:
3244  // grad xx
3245  eval0.template hessians<0, true, false>(values_dofs,
3246  scratch_data);
3247  eval1.template values<1, true, false>(scratch_data,
3248  hessians_quad);
3249 
3250  // grad yy
3251  eval0.template values<0, true, false>(values_dofs,
3252  scratch_data);
3253  eval1.template hessians<1, true, false>(scratch_data,
3254  hessians_quad +
3255  n_q_points);
3256 
3257  // grad zz
3258  eval0.template values<0, true, false>(values_dofs +
3259  2 * dofs_stride,
3260  scratch_data);
3261  eval1.template values<1, true, false>(scratch_data,
3262  hessians_quad +
3263  2 * n_q_points);
3264 
3265  // grad xy
3266  eval0.template gradients<0, true, false>(values_dofs,
3267  scratch_data);
3268  eval1.template gradients<1, true, false>(scratch_data,
3269  hessians_quad +
3270  3 * n_q_points);
3271 
3272  // grad xz
3273  eval0.template gradients<0, true, false>(values_dofs +
3274  dofs_stride,
3275  scratch_data);
3276  eval1.template values<1, true, false>(scratch_data,
3277  hessians_quad +
3278  4 * n_q_points);
3279 
3280  // grad yz
3281  eval0.template values<0, true, false>(values_dofs + dofs_stride,
3282  scratch_data);
3283  eval1.template gradients<1, true, false>(scratch_data,
3284  hessians_quad +
3285  5 * n_q_points);
3286 
3287  break;
3288  case 2:
3289  // grad xx
3290  eval0.template hessians<0, true, false>(values_dofs,
3291  hessians_quad);
3292  // grad yy
3293  eval0.template values<0, true, false>(
3294  values_dofs + 2 * dofs_stride, hessians_quad + n_q_points);
3295  // grad xy
3296  eval0.template gradients<0, true, false>(
3297  values_dofs + dofs_stride, hessians_quad + 2 * n_q_points);
3298  break;
3299  default:
3300  AssertThrow(false, ExcNotImplemented());
3301  }
3302  }
3303  }
3304 
3305  template <int normal_dir>
3306  static inline void
3308  Number * values_dofs,
3310  Number * scratch_data,
3311  const EvaluationFlags::EvaluationFlags evaluation_flag,
3312  const unsigned int face_direction,
3313  const unsigned int subface_index,
3314  std::integral_constant<bool, true>)
3315  {
3316  using EvalNormal =
3318  dim - 1,
3319  (fe_degree == -1) ? 1 : fe_degree + 1,
3320  n_q_points_1d,
3321  Number,
3322  normal_dir>;
3323  using EvalTangent =
3325  dim - 1,
3326  (fe_degree == -1) ? 1 : fe_degree,
3327  n_q_points_1d,
3328  Number,
3329  normal_dir>;
3330 
3331  using TempEval0 = typename std::
3332  conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
3333  using TempEval1 = typename std::
3334  conditional<normal_dir == 0, EvalTangent, EvalNormal>::type;
3335  using Eval0 = typename std::
3336  conditional<normal_dir == 2, EvalGeneral, TempEval0>::type;
3337  using Eval1 = typename std::
3338  conditional<normal_dir == 2, EvalGeneral, TempEval1>::type;
3339 
3340  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info =
3341  fe_eval.get_shape_info();
3342  Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
3343  ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]),
3344  subface_index,
3345  0);
3346  Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
3347  ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]),
3348  subface_index,
3349  1);
3350 
3351  constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
3352  const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face;
3353  const std::size_t n_dofs_normal =
3354  n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
3355  const std::size_t dofs_stride =
3356  (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
3357  n_dofs_tangent;
3358 
3359  static constexpr ::ndarray<unsigned int, 3, 3> component_table = {
3360  {{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
3361  const unsigned int component =
3362  (dim == 2 && normal_dir == 0 && face_direction == 1) ?
3363  0 :
3364  component_table[face_direction][normal_dir];
3365 
3366  // Initial offsets
3367  values_dofs +=
3368  3 * ((component == 0) ?
3369  0 :
3370  ((component == 1) ?
3371  ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) :
3372  ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent :
3373  n_dofs_normal + n_dofs_tangent)));
3374  const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir;
3375  Number *values_quad = fe_eval.begin_values() + n_q_points * shift;
3376  Number *gradients_quad =
3377  fe_eval.begin_gradients() + dim * n_q_points * shift;
3378  Number *hessians_quad =
3379  fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * shift;
3380 
3381  // Integration path
3382  if ((evaluation_flag & EvaluationFlags::values) &&
3383  !(evaluation_flag & EvaluationFlags::gradients))
3384  {
3385  switch (dim)
3386  {
3387  case 3:
3388  eval1.template values<1, false, false>(values_quad,
3389  values_quad);
3390  eval0.template values<0, false, false>(values_quad,
3391  values_dofs);
3392  break;
3393  case 2:
3394  eval0.template values<0, false, false>(values_quad,
3395  values_dofs);
3396  break;
3397  default:
3398  Assert(false, ExcNotImplemented());
3399  }
3400  }
3401  else if (evaluation_flag & EvaluationFlags::gradients)
3402  {
3403  switch (dim)
3404  {
3405  case 3:
3406  // grad z
3407  eval1.template values<1, false, false>(gradients_quad +
3408  2 * n_q_points,
3409  gradients_quad +
3410  2 * n_q_points);
3411  eval0.template values<0, false, false>(
3412  gradients_quad + 2 * n_q_points, values_dofs + dofs_stride);
3413 
3414  if (evaluation_flag & EvaluationFlags::values)
3415  {
3416  eval1.template values<1, false, false>(values_quad,
3417  scratch_data);
3418  eval1.template gradients<1, false, true>(gradients_quad +
3419  n_q_points,
3420  scratch_data);
3421  }
3422  else
3423  eval1.template gradients<1, false, false>(gradients_quad +
3424  n_q_points,
3425  scratch_data);
3426 
3427  // grad y
3428  eval0.template values<0, false, false>(scratch_data,
3429  values_dofs);
3430 
3431  // grad x
3432  eval1.template values<1, false, false>(gradients_quad,
3433  scratch_data);
3434  eval0.template gradients<0, false, true>(scratch_data,
3435  values_dofs);
3436 
3437  break;
3438  case 2:
3439  eval0.template values<0, false, false>(
3440  gradients_quad + n_q_points, values_dofs + dofs_stride);
3441  eval0.template gradients<0, false, false>(gradients_quad,
3442  values_dofs);
3443  if (evaluation_flag & EvaluationFlags::values)
3444  eval0.template values<0, false, true>(values_quad,
3445  values_dofs);
3446  break;
3447  default:
3448  AssertThrow(false, ExcNotImplemented());
3449  }
3450  }
3451 
3452  if (evaluation_flag & EvaluationFlags::hessians)
3453  {
3454  switch (dim)
3455  {
3456  case 3:
3457  // grad xx
3458  eval1.template values<1, false, false>(hessians_quad,
3459  scratch_data);
3460  if ((evaluation_flag &
3462  eval0.template hessians<0, false, true>(scratch_data,
3463  values_dofs);
3464  else
3465  eval0.template hessians<0, false, false>(scratch_data,
3466  values_dofs);
3467 
3468  // grad yy
3469  eval1.template hessians<1, false, false>(hessians_quad +
3470  n_q_points,
3471  scratch_data);
3472  eval0.template values<0, false, true>(scratch_data,
3473  values_dofs);
3474 
3475  // grad zz
3476  eval1.template values<1, false, false>(hessians_quad +
3477  2 * n_q_points,
3478  scratch_data);
3479  eval0.template values<0, false, false>(scratch_data,
3480  values_dofs +
3481  2 * dofs_stride);
3482 
3483  // grad xy
3484  eval1.template gradients<1, false, false>(hessians_quad +
3485  3 * n_q_points,
3486  scratch_data);
3487  eval0.template gradients<0, false, true>(scratch_data,
3488  values_dofs);
3489 
3490  // grad xz
3491  eval1.template values<1, false, false>(hessians_quad +
3492  4 * n_q_points,
3493  scratch_data);
3494  if ((evaluation_flag & EvaluationFlags::gradients))
3495  eval0.template gradients<0, false, true>(scratch_data,
3496  values_dofs +
3497  dofs_stride);
3498  else
3499  eval0.template gradients<0, false, false>(scratch_data,
3500  values_dofs +
3501  dofs_stride);
3502 
3503  // grad yz
3504  eval1.template gradients<1, false, false>(hessians_quad +
3505  5 * n_q_points,
3506  scratch_data);
3507  eval0.template values<0, false, true>(scratch_data,
3508  values_dofs +
3509  dofs_stride);
3510 
3511  break;
3512  case 2:
3513  // grad xx
3514  if (evaluation_flag &
3516  eval0.template hessians<0, false, true>(hessians_quad,
3517  values_dofs);
3518  else
3519  eval0.template hessians<0, false, false>(hessians_quad,
3520  values_dofs);
3521 
3522  // grad yy
3523  eval0.template values<0, false, false>(
3524  hessians_quad + n_q_points, values_dofs + 2 * dofs_stride);
3525  // grad xy
3526  if ((evaluation_flag & EvaluationFlags::gradients))
3527  eval0.template gradients<0, false, true>(
3528  hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
3529  else
3530  eval0.template gradients<0, false, false>(
3531  hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
3532  break;
3533  default:
3534  AssertThrow(false, ExcNotImplemented());
3535  }
3536  }
3537  }
3538  };
3539 
3540 
3541  template <int dim, int fe_degree, typename Number>
3543  {
3544  template <bool do_evaluate, bool add_into_output>
3545  static void
3546  interpolate(const unsigned int n_components,
3548  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
3549  const Number * input,
3550  Number * output,
3551  const unsigned int face_no)
3552  {
3553  Assert(static_cast<unsigned int>(fe_degree) ==
3554  shape_info.data.front().fe_degree ||
3555  fe_degree == -1,
3556  ExcInternalError());
3558  interpolate_generic_raviart_thomas<do_evaluate, add_into_output>(
3559  n_components, input, output, flags, face_no, shape_info);
3560  else
3561  interpolate_generic<do_evaluate, add_into_output>(
3562  n_components,
3563  input,
3564  output,
3565  flags,
3566  face_no,
3567  shape_info.data.front().fe_degree + 1,
3568  shape_info.data.front().shape_data_on_face,
3569  shape_info.dofs_per_component_on_cell,
3570  3 * shape_info.dofs_per_component_on_face);
3571  }
3572 
3576  template <bool do_evaluate, bool add_into_output>
3577  static void
3579  const unsigned int n_components,
3581  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
3582  const Number * input,
3583  Number * output,
3584  const unsigned int face_no)
3585  {
3586  Assert(static_cast<unsigned int>(fe_degree + 1) ==
3587  shape_info.data.front().n_q_points_1d ||
3588  fe_degree == -1,
3589  ExcInternalError());
3590 
3591  interpolate_generic<do_evaluate, add_into_output>(
3592  n_components,
3593  input,
3594  output,
3595  flags,
3596  face_no,
3597  shape_info.data.front().quadrature.size(),
3598  shape_info.data.front().quadrature_data_on_face,
3599  shape_info.n_q_points,
3600  shape_info.n_q_points_face);
3601  }
3602 
3603  private:
3604  template <bool do_evaluate, bool add_into_output, int face_direction = 0>
3605  static void
3606  interpolate_generic(const unsigned int n_components,
3607  const Number * input,
3608  Number * output,
3610  const unsigned int face_no,
3611  const unsigned int n_points_1d,
3612  const std::array<AlignedVector<Number>, 2> &shape_data,
3613  const unsigned int dofs_per_component_on_cell,
3614  const unsigned int dofs_per_component_on_face)
3615  {
3616  if (face_direction == face_no / 2)
3617  {
3619  dim,
3620  fe_degree + 1,
3621  0,
3622  Number>
3623  evalf(shape_data[face_no % 2],
3626  n_points_1d,
3627  0);
3628 
3629  const unsigned int in_stride = do_evaluate ?
3630  dofs_per_component_on_cell :
3631  dofs_per_component_on_face;
3632  const unsigned int out_stride = do_evaluate ?
3633  dofs_per_component_on_face :
3634  dofs_per_component_on_cell;
3635 
3636  for (unsigned int c = 0; c < n_components; ++c)
3637  {
3638  if (flag & EvaluationFlags::hessians)
3639  evalf.template apply_face<face_direction,
3640  do_evaluate,
3641  add_into_output,
3642  2>(input, output);
3643  else if (flag & EvaluationFlags::gradients)
3644  evalf.template apply_face<face_direction,
3645  do_evaluate,
3646  add_into_output,
3647  1>(input, output);
3648  else
3649  evalf.template apply_face<face_direction,
3650  do_evaluate,
3651  add_into_output,
3652  0>(input, output);
3653  input += in_stride;
3654  output += out_stride;
3655  }
3656  }
3657  else if (face_direction < dim)
3658  {
3659  interpolate_generic<do_evaluate,
3660  add_into_output,
3661  std::min(face_direction + 1, dim - 1)>(
3662  n_components,
3663  input,
3664  output,
3665  flag,
3666  face_no,
3667  n_points_1d,
3668  shape_data,
3669  dofs_per_component_on_cell,
3670  dofs_per_component_on_face);
3671  }
3672  }
3673 
3674  template <typename EvalType>
3675  static EvalType
3678  const unsigned int face_no)
3679  {
3680  return EvalType(data.shape_data_on_face[face_no % 2],
3683  }
3684 
3685  template <bool do_evaluate,
3686  bool add_into_output,
3687  int face_direction = 0,
3688  int max_derivative = 0>
3689  static void
3691  const unsigned int n_components,
3692  const Number * input,
3693  Number * output,
3695  const unsigned int face_no,
3696  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info)
3697  {
3698  if (dim == 1)
3699  {
3700  // This should never happen since the FE_RaviartThomasNodal is not
3701  // defined for dim = 1. It prevents compiler warnings of infinite
3702  // recursion.
3703  Assert(false, ExcInternalError());
3704  return;
3705  }
3706 
3707  bool increase_max_der = false;
3708  if ((flag & EvaluationFlags::hessians && max_derivative < 2) ||
3709  (flag & EvaluationFlags::gradients && max_derivative < 1))
3710  increase_max_der = true;
3711 
3712  if (face_direction == face_no / 2 && !increase_max_der)
3713  {
3715  add_into_output,
3716  face_direction,
3717  max_derivative>(
3718  shape_info, face_no, input, output);
3719  }
3720  else if (face_direction == face_no / 2)
3721  {
3722  // Only increase max_derivative
3724  add_into_output,
3725  face_direction,
3726  std::min(max_derivative + 1, 2)>(
3727  n_components, input, output, flag, face_no, shape_info);
3728  }
3729  else if (face_direction < dim)
3730  {
3731  if (increase_max_der)
3732  {
3734  do_evaluate,
3735  add_into_output,
3736  std::min(face_direction + 1, dim - 1),
3737  std::min(max_derivative + 1, 2)>(
3738  n_components, input, output, flag, face_no, shape_info);
3739  }
3740  else
3741  {
3743  add_into_output,
3744  std::min(face_direction + 1,
3745  dim - 1),
3746  max_derivative>(
3747  n_components, input, output, flag, face_no, shape_info);
3748  }
3749  }
3750  }
3751 
3752  /* Help function for interpolate_generic_raviart_thomas */
3753  template <bool do_evaluate,
3754  bool add_into_output,
3755  int face_direction,
3756  int max_derivative>
3757  static inline void
3759  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
3760  const unsigned int face_no,
3761  const Number * input,
3762  Number * output)
3763  {
3764  // These types are evaluators in either normal or tangential direction
3765  // depending on the face direction, with different normal directions for
3766  // the different components.
3767  using Evalf0 = typename std::conditional<
3768  face_direction == 0,
3770  dim,
3771  (fe_degree == -1) ? 1 : fe_degree + 1,
3772  0,
3773  Number,
3774  0>,
3776  dim,
3777  (fe_degree == -1) ? 1 : fe_degree,
3778  0,
3779  Number,
3780  0>>::type;
3781  using Evalf1 = typename std::conditional<
3782  face_direction == 1,
3784  dim,
3785  (fe_degree == -1) ? 1 : fe_degree + 1,
3786  0,
3787  Number,
3788  1>,
3790  dim,
3791  (fe_degree == -1) ? 1 : fe_degree,
3792  0,
3793  Number,
3794  1>>::type;
3795  using Evalf2 = typename std::conditional<
3796  face_direction == 2,
3798  dim,
3799  (fe_degree == -1) ? 1 : fe_degree + 1,
3800  0,
3801  Number,
3802  2>,
3804  dim,
3805  (fe_degree == -1) ? 1 : fe_degree,
3806  0,
3807  Number,
3808  2>>::type;
3809 
3810  Evalf0 evalf0 =
3811  create_evaluator_tensor_product<Evalf0>((face_direction == 0) ?
3812  shape_info.data[0] :
3813  shape_info.data[1],
3814  face_no);
3815  Evalf1 evalf1 =
3816  create_evaluator_tensor_product<Evalf1>((face_direction == 1) ?
3817  shape_info.data[0] :
3818  shape_info.data[1],
3819  face_no);
3820  Evalf2 evalf2 =
3821  create_evaluator_tensor_product<Evalf2>((face_direction == 2) ?
3822  shape_info.data[0] :
3823  shape_info.data[1],
3824  face_no);
3825 
3826  const unsigned int dofs_per_component_on_cell =
3827  shape_info.dofs_per_component_on_cell;
3828  const unsigned int dofs_per_component_on_face =
3829  3 * shape_info.dofs_per_component_on_face;
3830 
3831  // NOTE! dofs_per_component_on_face is in the tangent direction,
3832  // i.e (fe.degree+1)*fe.degree. Normal faces are only
3833  // fe.degree*fe.degree
3834  const unsigned int in_stride =
3835  do_evaluate ? dofs_per_component_on_cell : dofs_per_component_on_face;
3836  const unsigned int out_stride =
3837  do_evaluate ? dofs_per_component_on_face : dofs_per_component_on_cell;
3838 
3839  const unsigned int in_stride_after_normal =
3840  do_evaluate ?
3841  dofs_per_component_on_cell :
3842  dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2);
3843  const unsigned int out_stride_after_normal =
3844  do_evaluate ?
3845  dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2) :
3846  dofs_per_component_on_cell;
3847 
3848  evalf0.template apply_face<face_direction,
3849  do_evaluate,
3850  add_into_output,
3851  max_derivative>(input, output);
3852  // stride to next component
3853  input += (face_direction == 0) ? in_stride_after_normal : in_stride;
3854  output += (face_direction == 0) ? out_stride_after_normal : out_stride;
3855 
3856  evalf1.template apply_face<face_direction,
3857  do_evaluate,
3858  add_into_output,
3859  max_derivative>(input, output);
3860 
3861  if (dim == 3)
3862  {
3863  // stride to next component
3864  input += (face_direction == 1) ? in_stride_after_normal : in_stride;
3865  output +=
3866  (face_direction == 1) ? out_stride_after_normal : out_stride;
3867 
3868  evalf2.template apply_face<face_direction,
3869  do_evaluate,
3870  add_into_output,
3871  max_derivative>(input, output);
3872  }
3873  }
3874  };
3875 
3876 
3877 
3878  // internal helper function for reading data; base version of different types
3879  template <typename VectorizedArrayType, typename Number2>
3880  void
3881  do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
3882  {
3883  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3884  dst[v] = src_ptr[v];
3885  }
3886 
3887 
3888 
3889  // internal helper function for reading data; specialized version where we
3890  // can use a dedicated load function
3891  template <typename Number, std::size_t width>
3892  void
3894  {
3895  dst.load(src_ptr);
3896  }
3897 
3898 
3899 
3900  // internal helper function for reading data; base version of different types
3901  template <typename VectorizedArrayType, typename Number2>
3902  void
3903  do_vectorized_gather(const Number2 * src_ptr,
3904  const unsigned int * indices,
3905  VectorizedArrayType &dst)
3906  {
3907  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3908  dst[v] = src_ptr[indices[v]];
3909  }
3910 
3911 
3912 
3913  // internal helper function for reading data; specialized version where we
3914  // can use a dedicated gather function
3915  template <typename Number, std::size_t width>
3916  void
3917  do_vectorized_gather(const Number * src_ptr,
3918  const unsigned int * indices,
3920  {
3921  dst.gather(src_ptr, indices);
3922  }
3923 
3924 
3925 
3926  // internal helper function for reading data; base version of different types
3927  template <typename VectorizedArrayType, typename Number2>
3928  void
3929  do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
3930  {
3931  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3932  dst_ptr[v] += src[v];
3933  }
3934 
3935 
3936 
3937  // internal helper function for reading data; specialized version where we
3938  // can use a dedicated load function
3939  template <typename Number, std::size_t width>
3940  void
3942  {
3944  tmp.load(dst_ptr);
3945  (tmp + src).store(dst_ptr);
3946  }
3947 
3948 
3949 
3950  // internal helper function for reading data; base version of different types
3951  template <typename VectorizedArrayType, typename Number2>
3952  void
3953  do_vectorized_scatter_add(const VectorizedArrayType src,
3954  const unsigned int * indices,
3955  Number2 * dst_ptr)
3956  {
3957  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3958  dst_ptr[indices[v]] += src[v];
3959  }
3960 
3961 
3962 
3963  // internal helper function for reading data; specialized version where we
3964  // can use a dedicated gather function
3965  template <typename Number, std::size_t width>
3966  void
3968  const unsigned int * indices,
3969  Number * dst_ptr)
3970  {
3971 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
3972  for (unsigned int v = 0; v < width; ++v)
3973  dst_ptr[indices[v]] += src[v];
3974 #else
3976  tmp.gather(dst_ptr, indices);
3977  (tmp + src).scatter(indices, dst_ptr);
3978 #endif
3979  }
3980 
3981 
3982 
3983  template <typename Number>
3984  void
3985  adjust_for_face_orientation(const unsigned int dim,
3986  const unsigned int n_components,
3988  const unsigned int *orientation,
3989  const bool integrate,
3990  const std::size_t n_q_points,
3991  Number * tmp_values,
3992  Number * values_quad,
3993  Number * gradients_quad,
3994  Number * hessians_quad)
3995  {
3996  for (unsigned int c = 0; c < n_components; ++c)
3997  {
3998  if (flag & EvaluationFlags::values)
3999  {
4000  if (integrate)
4001  for (unsigned int q = 0; q < n_q_points; ++q)
4002  tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
4003  else
4004  for (unsigned int q = 0; q < n_q_points; ++q)
4005  tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
4006  for (unsigned int q = 0; q < n_q_points; ++q)
4007  values_quad[c * n_q_points + q] = tmp_values[q];
4008  }
4009  if (flag & EvaluationFlags::gradients)
4010  for (unsigned int d = 0; d < dim; ++d)
4011  {
4012  if (integrate)
4013  for (unsigned int q = 0; q < n_q_points; ++q)
4014  tmp_values[q] =
4015  gradients_quad[(c * dim + d) * n_q_points + orientation[q]];
4016  else
4017  for (unsigned int q = 0; q < n_q_points; ++q)
4018  tmp_values[orientation[q]] =
4019  gradients_quad[(c * dim + d) * n_q_points + q];
4020  for (unsigned int q = 0; q < n_q_points; ++q)
4021  gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q];
4022  }
4023  if (flag & EvaluationFlags::hessians)
4024  {
4025  const unsigned int hdim = (dim * (dim + 1)) / 2;
4026  for (unsigned int d = 0; d < hdim; ++d)
4027  {
4028  if (integrate)
4029  for (unsigned int q = 0; q < n_q_points; ++q)
4030  tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
4031  orientation[q]];
4032  else
4033  for (unsigned int q = 0; q < n_q_points; ++q)
4034  tmp_values[orientation[q]] =
4035  hessians_quad[(c * hdim + d) * n_q_points + q];
4036  for (unsigned int q = 0; q < n_q_points; ++q)
4037  hessians_quad[(c * hdim + d) * n_q_points + q] =
4038  tmp_values[q];
4039  }
4040  }
4041  }
4042  }
4043 
4044 
4045 
4046  template <typename Number, typename VectorizedArrayType>
4047  void
4049  const unsigned int dim,
4050  const unsigned int n_components,
4051  const unsigned int v,
4053  const unsigned int * orientation,
4054  const bool integrate,
4055  const std::size_t n_q_points,
4056  Number * tmp_values,
4057  VectorizedArrayType * values_quad,
4058  VectorizedArrayType * gradients_quad = nullptr,
4059  VectorizedArrayType * hessians_quad = nullptr)
4060  {
4061  for (unsigned int c = 0; c < n_components; ++c)
4062  {
4063  if (flag & EvaluationFlags::values)
4064  {
4065  if (integrate)
4066  for (unsigned int q = 0; q < n_q_points; ++q)
4067  tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v];
4068  else
4069  for (unsigned int q = 0; q < n_q_points; ++q)
4070  tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v];
4071  for (unsigned int q = 0; q < n_q_points; ++q)
4072  values_quad[c * n_q_points + q][v] = tmp_values[q];
4073  }
4074  if (flag & EvaluationFlags::gradients)
4075  for (unsigned int d = 0; d < dim; ++d)
4076  {
4077  Assert(gradients_quad != nullptr, ExcInternalError());
4078  if (integrate)
4079  for (unsigned int q = 0; q < n_q_points; ++q)
4080  tmp_values[q] = gradients_quad[(c * dim + d) * n_q_points +
4081  orientation[q]][v];
4082  else
4083  for (unsigned int q = 0; q < n_q_points; ++q)
4084  tmp_values[orientation[q]] =
4085  gradients_quad[(c * dim + d) * n_q_points + q][v];
4086  for (unsigned int q = 0; q < n_q_points; ++q)
4087  gradients_quad[(c * dim + d) * n_q_points + q][v] =
4088  tmp_values[q];
4089  }
4090  if (flag & EvaluationFlags::hessians)
4091  {
4092  Assert(hessians_quad != nullptr, ExcInternalError());
4093  const unsigned int hdim = (dim * (dim + 1)) / 2;
4094  for (unsigned int d = 0; d < hdim; ++d)
4095  {
4096  if (integrate)
4097  for (unsigned int q = 0; q < n_q_points; ++q)
4098  tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
4099  orientation[q]][v];
4100  else
4101  for (unsigned int q = 0; q < n_q_points; ++q)
4102  tmp_values[orientation[q]] =
4103  hessians_quad[(c * hdim + d) * n_q_points + q][v];
4104  for (unsigned int q = 0; q < n_q_points; ++q)
4105  hessians_quad[(c * hdim + d) * n_q_points + q][v] =
4106  tmp_values[q];
4107  }
4108  }
4109  }
4110  }
4111 
4112 
4113 
4114  template <int dim, typename Number>
4116  {
4117  template <int fe_degree, int n_q_points_1d>
4118  static bool
4119  run(const unsigned int n_components,
4120  const EvaluationFlags::EvaluationFlags evaluation_flag,
4121  const Number * values_dofs,
4123  {
4124  const auto &shape_info = fe_eval.get_shape_info();
4125  const auto &shape_data = shape_info.data.front();
4126 
4127  if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
4128  {
4129  Assert((fe_eval.get_dof_access_index() ==
4131  fe_eval.is_interior_face() == false) == false,
4132  ExcNotImplemented());
4133 
4134  const unsigned int face_no = fe_eval.get_face_no();
4135  const unsigned int face_orientation = fe_eval.get_face_orientation();
4136  const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
4137  const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
4138 
4139  using Eval =
4141 
4142  if (evaluation_flag & EvaluationFlags::values)
4143  {
4144  const auto shape_values =
4145  &shape_data.shape_values_face(face_no, face_orientation, 0);
4146 
4147  auto values_quad_ptr = fe_eval.begin_values();
4148  auto values_dofs_actual_ptr = values_dofs;
4149 
4150  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
4151  for (unsigned int c = 0; c < n_components; ++c)
4152  {
4153  eval.template values<0, true, false>(values_dofs_actual_ptr,
4154  values_quad_ptr);
4155 
4156  values_quad_ptr += n_q_points;
4157  values_dofs_actual_ptr += n_dofs;
4158  }
4159  }
4160 
4161  if (evaluation_flag & EvaluationFlags::gradients)
4162  {
4163  auto gradients_quad_ptr = fe_eval.begin_gradients();
4164  auto values_dofs_actual_ptr = values_dofs;
4165 
4166  std::array<const Number *, dim> shape_gradients;
4167  for (unsigned int d = 0; d < dim; ++d)
4168  shape_gradients[d] = &shape_data.shape_gradients_face(
4169  face_no, face_orientation, d, 0);
4170 
4171  for (unsigned int c = 0; c < n_components; ++c)
4172  {
4173  for (unsigned int d = 0; d < dim; ++d)
4174  {
4175  Eval eval(nullptr,
4176  shape_gradients[d],
4177  nullptr,
4178  n_dofs,
4179  n_q_points);
4180 
4181  eval.template gradients<0, true, false>(
4182  values_dofs_actual_ptr, gradients_quad_ptr);
4183 
4184  gradients_quad_ptr += n_q_points;
4185  }
4186  values_dofs_actual_ptr += n_dofs;
4187  }
4188  }
4189 
4190  Assert(!(evaluation_flag & EvaluationFlags::hessians),
4191  ExcNotImplemented());
4192 
4193  return true;
4194  }
4195 
4196  const unsigned int dofs_per_face =
4197  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
4198  Utilities::pow(shape_data.fe_degree + 1, dim - 1);
4199 
4200  // Note: we always keep storage of values, 1st and 2nd derivatives in an
4201  // array, so reserve space for all three here
4202  Number *temp = fe_eval.get_scratch_data().begin();
4203  Number *scratch_data = temp + 3 * n_components * dofs_per_face;
4204 
4205  bool use_vectorization = true;
4206 
4207  if (fe_eval.get_dof_access_index() ==
4209  fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
4210  for (unsigned int v = 0; v < Number::size(); ++v)
4211  if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
4212  fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
4213  use_vectorization = false;
4214 
4215  if (use_vectorization == false)
4216  {
4217  for (unsigned int v = 0; v < Number::size(); ++v)
4218  {
4219  // the loop breaks once an invalid_unsigned_int is hit for
4220  // all cases except the exterior faces in the ECL loop (where
4221  // some faces might be at the boundaries but others not)
4222  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4223  {
4224  for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4225  ++i)
4226  temp[i][v] = 0;
4227  continue;
4228  }
4229 
4231  template interpolate<true, false>(n_components,
4232  evaluation_flag,
4233  shape_info,
4234  values_dofs,
4235  scratch_data,
4236  fe_eval.get_face_no(v));
4237 
4238  for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4239  ++i)
4240  temp[i][v] = scratch_data[i][v];
4241  }
4242  }
4243  else
4245  template interpolate<true, false>(n_components,
4246  evaluation_flag,
4247  shape_info,
4248  values_dofs,
4249  temp,
4250  fe_eval.get_face_no());
4251 
4252  const unsigned int subface_index = fe_eval.get_subface_index();
4253  constexpr unsigned int n_q_points_1d_actual =
4254  fe_degree > -1 ? n_q_points_1d : 0;
4255 
4256  if (fe_degree >= 1 &&
4257  shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
4258  {
4260  (fe_degree == -1) ? 1 : fe_degree,
4261  (n_q_points_1d < 1) ? 1 :
4262  n_q_points_1d,
4263  Number>::
4264  template evaluate_or_integrate_in_face<false>(
4265  evaluation_flag,
4266  temp,
4267  fe_eval,
4268  scratch_data,
4269  subface_index,
4270  fe_eval.get_face_no());
4271  }
4272  else if (fe_degree > -1 &&
4273  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
4274  shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
4275  FEFaceEvaluationImpl<true,
4276  dim,
4277  fe_degree,
4278  n_q_points_1d_actual,
4279  Number>::evaluate_in_face(n_components,
4280  evaluation_flag,
4281  shape_data,
4282  temp,
4283  fe_eval.begin_values(),
4284  fe_eval
4285  .begin_gradients(),
4286  fe_eval.begin_hessians(),
4287  scratch_data,
4288  subface_index);
4289  else
4290  FEFaceEvaluationImpl<false,
4291  dim,
4292  fe_degree,
4293  n_q_points_1d_actual,
4294  Number>::evaluate_in_face(n_components,
4295  evaluation_flag,
4296  shape_data,
4297  temp,
4298  fe_eval.begin_values(),
4299  fe_eval
4300  .begin_gradients(),
4301  fe_eval.begin_hessians(),
4302  scratch_data,
4303  subface_index);
4304 
4305  if (use_vectorization == false)
4306  {
4307  for (unsigned int v = 0; v < Number::size(); ++v)
4308  {
4309  // the loop breaks once an invalid_unsigned_int is hit for
4310  // all cases except the exterior faces in the ECL loop (where
4311  // some faces might be at the boundaries but others not)
4312  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4313  continue;
4314 
4315  if (fe_eval.get_face_orientation(v) != 0)
4317  dim,
4318  n_components,
4319  v,
4320  evaluation_flag,
4322  fe_eval.get_face_orientation(v), 0),
4323  false,
4324  shape_info.n_q_points_face,
4325  &temp[0][0],
4326  fe_eval.begin_values(),
4327  fe_eval.begin_gradients(),
4328  fe_eval.begin_hessians());
4329  }
4330  }
4331  else if (fe_eval.get_face_orientation() != 0)
4333  dim,
4334  n_components,
4335  evaluation_flag,
4337  fe_eval.get_face_orientation(), 0),
4338  false,
4339  shape_info.n_q_points_face,
4340  temp,
4341  fe_eval.begin_values(),
4342  fe_eval.begin_gradients(),
4343  fe_eval.begin_hessians());
4344 
4345  return false;
4346  }
4347  };
4348 
4349 
4350 
4351  template <int dim, typename Number>
4353  {
4354  template <int fe_degree, int n_q_points_1d>
4355  static bool
4356  run(const unsigned int n_components,
4357  const EvaluationFlags::EvaluationFlags integration_flag,
4358  Number * values_dofs,
4360  {
4361  const auto &shape_info = fe_eval.get_shape_info();
4362  const auto &shape_data = shape_info.data.front();
4363 
4364  if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
4365  {
4366  Assert((fe_eval.get_dof_access_index() ==
4368  fe_eval.is_interior_face() == false) == false,
4369  ExcNotImplemented());
4370 
4371  const unsigned int face_no = fe_eval.get_face_no();
4372  const unsigned int face_orientation = fe_eval.get_face_orientation();
4373  const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
4374  const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
4375 
4376  using Eval =
4378 
4379  if (integration_flag & EvaluationFlags::values)
4380  {
4381  const auto shape_values =
4382  &shape_data.shape_values_face(face_no, face_orientation, 0);
4383 
4384  auto values_quad_ptr = fe_eval.begin_values();
4385  auto values_dofs_actual_ptr = values_dofs;
4386 
4387  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
4388  for (unsigned int c = 0; c < n_components; ++c)
4389  {
4390  eval.template values<0, false, false>(values_quad_ptr,
4391  values_dofs_actual_ptr);
4392 
4393  values_quad_ptr += n_q_points;
4394  values_dofs_actual_ptr += n_dofs;
4395  }
4396  }
4397 
4398  if (integration_flag & EvaluationFlags::gradients)
4399  {
4400  auto gradients_quad_ptr = fe_eval.begin_gradients();
4401  auto values_dofs_actual_ptr = values_dofs;
4402 
4403  std::array<const Number *, dim> shape_gradients;
4404  for (unsigned int d = 0; d < dim; ++d)
4405  shape_gradients[d] = &shape_data.shape_gradients_face(
4406  face_no, face_orientation, d, 0);
4407 
4408  for (unsigned int c = 0; c < n_components; ++c)
4409  {
4410  for (unsigned int d = 0; d < dim; ++d)
4411  {
4412  Eval eval(nullptr,
4413  shape_gradients[d],
4414  nullptr,
4415  n_dofs,
4416  n_q_points);
4417 
4418  if (!(integration_flag & EvaluationFlags::values) &&
4419  d == 0)
4420  eval.template gradients<0, false, false>(
4421  gradients_quad_ptr, values_dofs_actual_ptr);
4422  else
4423  eval.template gradients<0, false, true>(
4424  gradients_quad_ptr, values_dofs_actual_ptr);
4425 
4426  gradients_quad_ptr += n_q_points;
4427  }
4428  values_dofs_actual_ptr += n_dofs;
4429  }
4430  }
4431 
4432  Assert(!(integration_flag & EvaluationFlags::hessians),
4433  ExcNotImplemented());
4434 
4435  return true;
4436  }
4437 
4438  const unsigned int dofs_per_face =
4439  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
4440  Utilities::pow(shape_data.fe_degree + 1, dim - 1);
4441 
4442  Number *temp = fe_eval.get_scratch_data().begin();
4443  Number *scratch_data = temp + 3 * n_components * dofs_per_face;
4444 
4445  bool use_vectorization = true;
4446 
4447  if (fe_eval.get_dof_access_index() ==
4449  fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
4450  use_vectorization =
4451  fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
4452  std::all_of(fe_eval.get_cell_ids().begin() + 1,
4453  fe_eval.get_cell_ids().end(),
4454  [&](const auto &v) {
4455  return v == fe_eval.get_cell_ids()[0] ||
4456  v == numbers::invalid_unsigned_int;
4457  });
4458 
4459  if (use_vectorization == false)
4460  {
4461  for (unsigned int v = 0; v < Number::size(); ++v)
4462  {
4463  // the loop breaks once an invalid_unsigned_int is hit for
4464  // all cases except the exterior faces in the ECL loop (where
4465  // some faces might be at the boundaries but others not)
4466  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4467  continue;
4468 
4469  if (fe_eval.get_face_orientation(v) != 0)
4471  dim,
4472  n_components,
4473  v,
4474  integration_flag,
4476  fe_eval.get_face_orientation(v), 0),
4477  true,
4478  shape_info.n_q_points_face,
4479  &temp[0][0],
4480  fe_eval.begin_values(),
4481  fe_eval.begin_gradients(),
4482  fe_eval.begin_hessians());
4483  }
4484  }
4485  else if (fe_eval.get_face_orientation() != 0)
4487  dim,
4488  n_components,
4489  integration_flag,
4491  fe_eval.get_face_orientation(), 0),
4492  true,
4493  shape_info.n_q_points_face,
4494  temp,
4495  fe_eval.begin_values(),
4496  fe_eval.begin_gradients(),
4497  fe_eval.begin_hessians());
4498 
4499  const unsigned int n_q_points_1d_actual =
4500  fe_degree > -1 ? n_q_points_1d : 0;
4501  const unsigned int subface_index = fe_eval.get_subface_index();
4502 
4503  if (fe_degree >= 1 &&
4504  shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
4505  {
4507  (fe_degree == -1) ? 1 : fe_degree,
4508  (n_q_points_1d < 1) ? 1 :
4509  n_q_points_1d,
4510  Number>::
4511  template evaluate_or_integrate_in_face<true>(integration_flag,
4512  temp,
4513  fe_eval,
4514  scratch_data,
4515  subface_index,
4516  fe_eval.get_face_no());
4517  }
4518  else if (fe_degree > -1 &&
4519  fe_eval.get_subface_index() >=
4521  shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
4523  true,
4524  dim,
4525  fe_degree,
4526  n_q_points_1d_actual,
4527  Number>::integrate_in_face(n_components,
4528  integration_flag,
4529  shape_data,
4530  temp,
4531  fe_eval.begin_values(),
4532  fe_eval.begin_gradients(),
4533  fe_eval.begin_hessians(),
4534  scratch_data,
4535  subface_index);
4536  else
4538  false,
4539  dim,
4540  fe_degree,
4541  n_q_points_1d_actual,
4542  Number>::integrate_in_face(n_components,
4543  integration_flag,
4544  shape_data,
4545  temp,
4546  fe_eval.begin_values(),
4547  fe_eval.begin_gradients(),
4548  fe_eval.begin_hessians(),
4549  scratch_data,
4550  subface_index);
4551 
4552  if (use_vectorization == false)
4553  {
4554  for (unsigned int v = 0; v < Number::size(); ++v)
4555  {
4556  // the loop breaks once an invalid_unsigned_int is hit for
4557  // all cases except the exterior faces in the ECL loop (where
4558  // some faces might be at the boundaries but others not)
4559  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4560  continue;
4561 
4563  template interpolate<false, false>(n_components,
4564  integration_flag,
4565  shape_info,
4566  values_dofs,
4567  scratch_data,
4568  fe_eval.get_face_no(v));
4569 
4570  for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4571  ++i)
4572  temp[i][v] = scratch_data[i][v];
4573  }
4574  }
4575  else
4577  template interpolate<false, false>(n_components,
4578  integration_flag,
4579  shape_info,
4580  temp,
4581  values_dofs,
4582  fe_eval.get_face_no());
4583  return false;
4584  }
4585  };
4586 
4587 
4588 
4589  template <int n_face_orientations,
4590  typename Processor,
4591  typename EvaluationData,
4592  const bool check_face_orientations = false>
4593  void
4595  Processor & proc,
4596  const unsigned int n_components,
4597  const EvaluationFlags::EvaluationFlags evaluation_flag,
4598  typename Processor::Number2_ * global_vector_ptr,
4599  const std::vector<ArrayView<const typename Processor::Number2_>> *sm_ptr,
4600  const EvaluationData & fe_eval,
4601  typename Processor::VectorizedArrayType_ * temp1)
4602  {
4603  constexpr int dim = Processor::dim_;
4604  constexpr int fe_degree = Processor::fe_degree_;
4605  using VectorizedArrayType = typename Processor::VectorizedArrayType_;
4606  constexpr int n_lanes = VectorizedArrayType::size();
4607 
4608  using Number = typename Processor::Number_;
4609  using Number2_ = typename Processor::Number2_;
4610 
4611  const auto & shape_data = fe_eval.get_shape_info().data.front();
4612  constexpr bool integrate = Processor::do_integrate;
4613  const unsigned int face_no = fe_eval.get_face_no();
4614  const auto & dof_info = fe_eval.get_dof_info();
4615  const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
4616  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index =
4617  fe_eval.get_dof_access_index();
4618  AssertIndexRange(cell,
4619  dof_info.index_storage_variants[dof_access_index].size());
4620  constexpr unsigned int dofs_per_face =
4621  Utilities::pow(fe_degree + 1, dim - 1);
4622  const unsigned int subface_index = fe_eval.get_subface_index();
4623 
4624  const unsigned int n_filled_lanes =
4625  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
4626 
4627  bool all_faces_are_same = n_filled_lanes == n_lanes;
4628  if (n_face_orientations == n_lanes)
4629  for (unsigned int v = 1; v < n_lanes; ++v)
4630  if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) ||
4631  fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0))
4632  {
4633  all_faces_are_same = false;
4634  break;
4635  }
4636 
4637  // check for re-orientation ...
4638  std::array<const unsigned int *, n_face_orientations> orientation = {};
4639 
4640  if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same &&
4641  fe_eval.is_interior_face() == 0)
4642  for (unsigned int v = 0; v < n_lanes; ++v)
4643  {
4644  // the loop breaks once an invalid_unsigned_int is hit for
4645  // all cases except the exterior faces in the ECL loop (where
4646  // some faces might be at the boundaries but others not)
4647  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4648  continue;
4649 
4650  if (shape_data.nodal_at_cell_boundaries &&
4651  fe_eval.get_face_orientation(v) != 0)
4652  {
4653  // ... and in case we detect a re-orientation, go to the other
4654  // version of this function that actually allows for this
4655  if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
4656  check_face_orientations == false)
4657  {
4658  fe_face_evaluation_process_and_io<n_face_orientations,
4659  Processor,
4660  EvaluationData,
4661  true>(proc,
4662  n_components,
4663  evaluation_flag,
4664  global_vector_ptr,
4665  sm_ptr,
4666  fe_eval,
4667  temp1);
4668  return;
4669  }
4670  orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
4671  fe_eval.get_face_orientation(v), 0);
4672  }
4673  }
4674  else if (dim == 3 && fe_eval.get_face_orientation() != 0)
4675  {
4676  // go to the other version of this function
4677  if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
4678  check_face_orientations == false)
4679  {
4680  fe_face_evaluation_process_and_io<n_face_orientations,
4681  Processor,
4682  EvaluationData,
4683  true>(proc,
4684  n_components,
4685  evaluation_flag,
4686  global_vector_ptr,
4687  sm_ptr,
4688  fe_eval,
4689  temp1);
4690  return;
4691  }
4692  for (unsigned int v = 0; v < n_face_orientations; ++v)
4693  orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
4694  fe_eval.get_face_orientation(), 0);
4695  }
4696 
4697  // we know that the gradient weights for the Hermite case on the
4698  // right (side==1) are the negative from the value at the left
4699  // (side==0), so we only read out one of them.
4700  VectorizedArrayType grad_weight =
4701  shape_data
4702  .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) :
4703  (1 + face_no % 2))];
4704 
4705  // face_to_cell_index_hermite
4706  std::array<const unsigned int *, n_face_orientations> index_array_hermite =
4707  {};
4708  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4709  {
4710  if (n_face_orientations == 1)
4711  index_array_hermite[0] =
4712  &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0);
4713  else
4714  {
4715  for (unsigned int v = 0; v < n_lanes; ++v)
4716  {
4717  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4718  continue;
4719 
4720  const auto face_no = fe_eval.get_face_no(v);
4721 
4722  grad_weight[v] =
4723  shape_data.shape_data_on_face[0][fe_degree +
4724  (integrate ?
4725  (2 - (face_no % 2)) :
4726  (1 + (face_no % 2)))][0];
4727 
4728  index_array_hermite[v] =
4729  &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no,
4730  0);
4731  }
4732  }
4733  }
4734 
4735  // face_to_cell_index_nodal
4736  std::array<const unsigned int *, n_face_orientations> index_array_nodal =
4737  {};
4738  if (shape_data.nodal_at_cell_boundaries == true)
4739  {
4740  if (n_face_orientations == 1)
4741  index_array_nodal[0] =
4742  &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0);
4743  else
4744  {
4745  for (unsigned int v = 0; v < n_lanes; ++v)
4746  {
4747  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4748  continue;
4749 
4750  const auto face_no = fe_eval.get_face_no(v);
4751 
4752  index_array_nodal[v] =
4753  &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no,
4754  0);
4755  }
4756  }
4757  }
4758 
4759 
4760  const auto reorientate = [&](const unsigned int v, const unsigned int i) {
4761  return (!check_face_orientations || orientation[v] == nullptr) ?
4762  i :
4763  orientation[v][i];
4764  };
4765 
4766  const unsigned int cell_index =
4767  dof_access_index == MatrixFreeFunctions::DoFInfo::dof_access_cell ?
4768  fe_eval.get_cell_ids()[0] :
4769  cell * n_lanes;
4770  const unsigned int *dof_indices =
4771  &dof_info.dof_indices_contiguous[dof_access_index][cell_index];
4772 
4773  for (unsigned int comp = 0; comp < n_components; ++comp)
4774  {
4775  const std::size_t index_offset =
4776  dof_info.component_dof_indices_offset
4777  [fe_eval.get_active_fe_index()]
4778  [fe_eval.get_first_selected_component()] +
4779  comp * Utilities::pow(fe_degree + 1, dim);
4780 
4781  // case 1: contiguous and interleaved indices
4782  if (n_face_orientations == 1 &&
4783  dof_info.index_storage_variants[dof_access_index][cell] ==
4785  interleaved_contiguous)
4786  {
4788  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
4789  n_lanes);
4790  Number2_ *vector_ptr =
4791  global_vector_ptr + dof_indices[0] + index_offset * n_lanes;
4792 
4793  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4794  {
4795  for (unsigned int i = 0; i < dofs_per_face; ++i)
4796  {
4797  Assert(n_face_orientations == 1, ExcNotImplemented());
4798 
4799  const unsigned int ind1 = index_array_hermite[0][2 * i];
4800  const unsigned int ind2 = index_array_hermite[0][2 * i + 1];
4801  const unsigned int i_ = reorientate(0, i);
4802  proc.hermite_grad_vectorized(temp1[i_],
4803  temp1[i_ + dofs_per_face],
4804  vector_ptr + ind1 * n_lanes,
4805  vector_ptr + ind2 * n_lanes,
4806  grad_weight);
4807  }
4808  }
4809  else
4810  {
4811  for (unsigned int i = 0; i < dofs_per_face; ++i)
4812  {
4813  Assert(n_face_orientations == 1, ExcNotImplemented());
4814 
4815  const unsigned int i_ = reorientate(0, i);
4816  const unsigned int ind = index_array_nodal[0][i];
4817  proc.value_vectorized(temp1[i_],
4818  vector_ptr + ind * n_lanes);
4819  }
4820  }
4821  }
4822 
4823  // case 2: contiguous and interleaved indices with fixed stride
4824  else if (n_face_orientations == 1 &&
4825  dof_info.index_storage_variants[dof_access_index][cell] ==
4827  interleaved_contiguous_strided)
4828  {
4830  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
4831  n_lanes);
4832  Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes;
4833  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4834  {
4835  for (unsigned int i = 0; i < dofs_per_face; ++i)
4836  {
4837  Assert(n_face_orientations == 1, ExcNotImplemented());
4838 
4839  const unsigned int i_ = reorientate(0, i);
4840  const unsigned int ind1 =
4841  index_array_hermite[0][2 * i] * n_lanes;
4842  const unsigned int ind2 =
4843  index_array_hermite[0][2 * i + 1] * n_lanes;
4844  proc.hermite_grad_vectorized_indexed(
4845  temp1[i_],
4846  temp1[i_ + dofs_per_face],
4847  vector_ptr + ind1,
4848  vector_ptr + ind2,
4849  grad_weight,
4850  dof_indices,
4851  dof_indices);
4852  }
4853  }
4854  else
4855  {
4856  for (unsigned int i = 0; i < dofs_per_face; ++i)
4857  {
4858  Assert(n_face_orientations == 1, ExcNotImplemented());
4859 
4860  const unsigned int i_ = reorientate(0, i);
4861  const unsigned int ind = index_array_nodal[0][i] * n_lanes;
4862  proc.value_vectorized_indexed(temp1[i_],
4863  vector_ptr + ind,
4864  dof_indices);
4865  }
4866  }
4867  }
4868 
4869  // case 3: contiguous and interleaved indices with mixed stride
4870  else if (n_face_orientations == 1 &&
4871  dof_info.index_storage_variants[dof_access_index][cell] ==
4873  interleaved_contiguous_mixed_strides)
4874  {
4875  const unsigned int *strides =
4876  &dof_info.dof_indices_interleave_strides[dof_access_index]
4877  [cell * n_lanes];
4878  unsigned int indices[n_lanes];
4879  for (unsigned int v = 0; v < n_lanes; ++v)
4880  indices[v] = dof_indices[v] + index_offset * strides[v];
4881  const unsigned int n_filled_lanes =
4882  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
4883 
4884  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4885  {
4886  if (n_filled_lanes == n_lanes)
4887  for (unsigned int i = 0; i < dofs_per_face; ++i)
4888  {
4889  Assert(n_face_orientations == 1, ExcNotImplemented());
4890 
4891  const unsigned int i_ = reorientate(0, i);
4892  unsigned int ind1[n_lanes];
4894  for (unsigned int v = 0; v < n_lanes; ++v)
4895  ind1[v] = indices[v] +
4896  index_array_hermite[0][2 * i] * strides[v];
4897  unsigned int ind2[n_lanes];
4899  for (unsigned int v = 0; v < n_lanes; ++v)
4900  ind2[v] =
4901  indices[v] +
4902  // TODO
4903  index_array_hermite[0][2 * i + 1] * strides[v];
4904  proc.hermite_grad_vectorized_indexed(
4905  temp1[i_],
4906  temp1[i_ + dofs_per_face],
4907  global_vector_ptr,
4908  global_vector_ptr,
4909  grad_weight,
4910  ind1,
4911  ind2);
4912  }
4913  else
4914  {
4915  if (integrate == false)
4916  for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
4917  temp1[i] = VectorizedArrayType();
4918 
4919  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4920  for (unsigned int i = 0; i < dofs_per_face; ++i)
4921  {
4922  const unsigned int i_ =
4923  reorientate(n_face_orientations == 1 ? 0 : v, i);
4924  proc.hermite_grad(
4925  temp1[i_][v],
4926  temp1[i_ + dofs_per_face][v],
4927  global_vector_ptr
4928  [indices[v] +
4929  index_array_hermite
4930  [n_face_orientations == 1 ? 0 : v][2 * i] *
4931  strides[v]],
4932  global_vector_ptr
4933  [indices[v] +
4934  index_array_hermite[n_face_orientations == 1 ?
4935  0 :
4936  v][2 * i + 1] *
4937  strides[v]],
4938  grad_weight[n_face_orientations == 1 ? 0 : v]);
4939  }
4940  }
4941  }
4942  else
4943  {
4944  if (n_filled_lanes == n_lanes)
4945  for (unsigned int i = 0; i < dofs_per_face; ++i)
4946  {
4947  Assert(n_face_orientations == 1, ExcInternalError());
4948  unsigned int ind[n_lanes];
4950  for (unsigned int v = 0; v < n_lanes; ++v)
4951  ind[v] =
4952  indices[v] + index_array_nodal[0][i] * strides[v];
4953  const unsigned int i_ = reorientate(0, i);
4954  proc.value_vectorized_indexed(temp1[i_],
4955  global_vector_ptr,
4956  ind);
4957  }
4958  else
4959  {
4960  if (integrate == false)
4961  for (unsigned int i = 0; i < dofs_per_face; ++i)
4962  temp1[i] = VectorizedArrayType();
4963 
4964  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4965  for (unsigned int i = 0; i < dofs_per_face; ++i)
4966  proc.value(
4967  temp1[reorientate(n_face_orientations == 1 ? 0 : v,
4968  i)][v],
4969  global_vector_ptr
4970  [indices[v] +
4971  index_array_nodal[n_face_orientations == 1 ? 0 : v]
4972  [i] *
4973  strides[v]]);
4974  }
4975  }
4976  }
4977 
4978  // case 4: contiguous indices without interleaving
4979  else if (n_face_orientations > 1 ||
4980  dof_info.index_storage_variants[dof_access_index][cell] ==
4982  contiguous)
4983  {
4984  Number2_ *vector_ptr = global_vector_ptr + index_offset;
4985 
4986  const bool vectorization_possible =
4987  all_faces_are_same && (sm_ptr == nullptr);
4988 
4989  std::array<Number2_ *, n_lanes> vector_ptrs;
4990  std::array<unsigned int, n_lanes> reordered_indices;
4991 
4992  if (vectorization_possible == false)
4993  {
4994  vector_ptrs = {};
4995  if (n_face_orientations == 1)
4996  {
4997  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4998  if (sm_ptr == nullptr)
4999  {
5000  vector_ptrs[v] = vector_ptr + dof_indices[v];
5001  }
5002  else
5003  {
5004  const auto &temp =
5005  dof_info
5006  .dof_indices_contiguous_sm[dof_access_index]
5007  [cell * n_lanes + v];
5008  vector_ptrs[v] = const_cast<Number2_ *>(
5009  sm_ptr->operator[](temp.first).data() +
5010  temp.second + index_offset);
5011  }
5012  }
5013  else if (n_face_orientations == n_lanes)
5014  {
5015  const auto &cells = fe_eval.get_cell_ids();
5016  for (unsigned int v = 0; v < n_lanes; ++v)
5017  if (cells[v] != numbers::invalid_unsigned_int)
5018  {
5019  if (sm_ptr == nullptr)
5020  {
5021  vector_ptrs[v] =
5022  vector_ptr +
5023  dof_info
5024  .dof_indices_contiguous[dof_access_index]
5025  [cells[v]];
5026  }
5027  else
5028  {
5029  const auto &temp =
5030  dof_info
5031  .dof_indices_contiguous_sm[dof_access_index]
5032  [cells[v]];
5033  vector_ptrs[v] = const_cast<Number2_ *>(
5034  sm_ptr->operator[](temp.first).data() +
5035  temp.second + index_offset);
5036  }
5037  }
5038  }
5039  else
5040  {
5041  Assert(false, ExcNotImplemented());
5042  }
5043  }
5044  else if (n_face_orientations == n_lanes)
5045  {
5046  for (unsigned int v = 0; v < n_lanes; ++v)
5047  reordered_indices[v] =
5048  dof_info.dof_indices_contiguous[dof_access_index]
5049  [fe_eval.get_cell_ids()[v]];
5050  dof_indices = reordered_indices.data();
5051  }
5052 
5053  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
5054  {
5055  if (vectorization_possible)
5056  for (unsigned int i = 0; i < dofs_per_face; ++i)
5057  {
5058  const unsigned int ind1 = index_array_hermite[0][2 * i];
5059  const unsigned int ind2 =
5060  index_array_hermite[0][2 * i + 1];
5061  const unsigned int i_ = reorientate(0, i);
5062 
5063  proc.hermite_grad_vectorized_indexed(
5064  temp1[i_],
5065  temp1[i_ + dofs_per_face],
5066  vector_ptr + ind1,
5067  vector_ptr + ind2,
5068  grad_weight,
5069  dof_indices,
5070  dof_indices);
5071  }
5072  else if (n_face_orientations == 1)
5073  for (unsigned int i = 0; i < dofs_per_face; ++i)
5074  {
5075  const unsigned int ind1 = index_array_hermite[0][2 * i];
5076  const unsigned int ind2 =
5077  index_array_hermite[0][2 * i + 1];
5078  const unsigned int i_ = reorientate(0, i);
5079 
5080  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5081  proc.hermite_grad(temp1[i_][v],
5082  temp1[i_ + dofs_per_face][v],
5083  vector_ptrs[v][ind1],
5084  vector_ptrs[v][ind2],
5085  grad_weight[v]);
5086 
5087  if (integrate == false)
5088  for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
5089  {
5090  temp1[i][v] = 0.0;
5091  temp1[i + dofs_per_face][v] = 0.0;
5092  }
5093  }
5094  else
5095  {
5096  if (integrate == false && n_filled_lanes < n_lanes)
5097  for (unsigned int i = 0; i < dofs_per_face; ++i)
5098  temp1[i] = temp1[i + dofs_per_face] = Number();
5099 
5100  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5101  for (unsigned int i = 0; i < dofs_per_face; ++i)
5102  proc.hermite_grad(
5103  temp1[reorientate(v, i)][v],
5104  temp1[reorientate(v, i) + dofs_per_face][v],
5105  vector_ptrs[v][index_array_hermite[v][2 * i]],
5106  vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
5107  grad_weight[v]);
5108  }
5109  }
5110  else
5111  {
5112  if (vectorization_possible)
5113  for (unsigned int i = 0; i < dofs_per_face; ++i)
5114  {
5115  const unsigned int ind = index_array_nodal[0][i];
5116  const unsigned int i_ = reorientate(0, i);
5117 
5118  proc.value_vectorized_indexed(temp1[i_],
5119  vector_ptr + ind,
5120  dof_indices);
5121  }
5122  else if (n_face_orientations == 1)
5123  for (unsigned int i = 0; i < dofs_per_face; ++i)
5124  {
5125  const unsigned int ind = index_array_nodal[0][i];
5126  const unsigned int i_ = reorientate(0, i);
5127 
5128  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5129  proc.value(temp1[i_][v], vector_ptrs[v][ind]);
5130 
5131  if (integrate == false)
5132  for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
5133  temp1[i_][v] = 0.0;
5134  }
5135  else
5136  {
5137  if (integrate == false && n_filled_lanes < n_lanes)
5138  for (unsigned int i = 0; i < dofs_per_face; ++i)
5139  temp1[i] = Number();
5140 
5141  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5142  for (unsigned int i = 0; i < dofs_per_face; ++i)
5143  proc.value(temp1[reorientate(v, i)][v],
5144  vector_ptrs[v][index_array_nodal[v][i]]);
5145  }
5146  }
5147  }
5148  else
5149  {
5150  // We should not end up here, this should be caught by
5151  // FEFaceEvaluationImplGatherEvaluateSelector::supports()
5152  Assert(false, ExcInternalError());
5153  }
5154  temp1 += 3 * dofs_per_face;
5155  }
5156  }
5157 
5158 
5159 
5160  template <int dim, typename Number2, typename VectorizedArrayType>
5162  {
5163  using Number = typename VectorizedArrayType::value_type;
5164 
5165  template <int fe_degree, int n_q_points_1d>
5166  static bool
5167  run(const unsigned int n_components,
5168  const EvaluationFlags::EvaluationFlags evaluation_flag,
5169  const Number2 * src_ptr,
5170  const std::vector<ArrayView<const Number2>> * sm_ptr,
5172  {
5173  Assert(fe_degree > -1, ExcInternalError());
5174  Assert(fe_eval.get_shape_info().element_type <=
5176  ExcInternalError());
5177 
5178  const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
5179 
5180  VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
5181  VectorizedArrayType *scratch_data =
5182  temp + 3 * n_components * dofs_per_face;
5183 
5185 
5186  if (fe_eval.get_dof_access_index() ==
5188  fe_eval.is_interior_face() == false)
5189  fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
5190  p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
5191  else
5192  fe_face_evaluation_process_and_io<1>(
5193  p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
5194 
5195  const unsigned int subface_index = fe_eval.get_subface_index();
5196 
5197  if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
5198  FEFaceEvaluationImpl<true,
5199  dim,
5200  fe_degree,
5201  n_q_points_1d,
5202  VectorizedArrayType>::
5203  evaluate_in_face(n_components,
5204  evaluation_flag,
5205  fe_eval.get_shape_info().data.front(),
5206  temp,
5207  fe_eval.begin_values(),
5208  fe_eval.begin_gradients(),
5209  fe_eval.begin_hessians(),
5210  scratch_data,
5211  subface_index);
5212  else
5213  FEFaceEvaluationImpl<false,
5214  dim,
5215  fe_degree,
5216  n_q_points_1d,
5217  VectorizedArrayType>::
5218  evaluate_in_face(n_components,
5219  evaluation_flag,
5220  fe_eval.get_shape_info().data.front(),
5221  temp,
5222  fe_eval.begin_values(),
5223  fe_eval.begin_gradients(),
5224  fe_eval.begin_hessians(),
5225  scratch_data,
5226  subface_index);
5227 
5228  // re-orientation for cases not possible with above algorithm
5229  if (subface_index < GeometryInfo<dim>::max_children_per_cell)
5230  {
5231  if (fe_eval.get_dof_access_index() ==
5233  fe_eval.is_interior_face() == false)
5234  {
5235  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5236  {
5237  // the loop breaks once an invalid_unsigned_int is hit for
5238  // all cases except the exterior faces in the ECL loop (where
5239  // some faces might be at the boundaries but others not)
5240  if (fe_eval.get_cell_ids()[v] ==
5242  continue;
5243 
5244  if (fe_eval.get_face_orientation(v) != 0)
5246  dim,
5247  n_components,
5248  v,
5249  evaluation_flag,
5251  fe_eval.get_face_orientation(v), 0),
5252  false,
5253  Utilities::pow(n_q_points_1d, dim - 1),
5254  &temp[0][0],
5255  fe_eval.begin_values(),
5256  fe_eval.begin_gradients(),
5257  fe_eval.begin_hessians());
5258  }
5259  }
5260  else if (fe_eval.get_face_orientation() != 0)
5262  dim,
5263  n_components,
5264  evaluation_flag,
5266  fe_eval.get_face_orientation(), 0),
5267  false,
5268  Utilities::pow(n_q_points_1d, dim - 1),
5269  temp,
5270  fe_eval.begin_values(),
5271  fe_eval.begin_gradients(),
5272  fe_eval.begin_hessians());
5273  }
5274 
5275  return false;
5276  }
5277 
5278  static bool
5280  const EvaluationFlags::EvaluationFlags evaluation_flag,
5282  const Number2 * vector_ptr,
5284  {
5285  const unsigned int fe_degree = shape_info.data.front().fe_degree;
5286  if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries ||
5287  (evaluation_flag & EvaluationFlags::gradients &&
5288  (fe_degree < 2 ||
5289  shape_info.data.front().element_type !=
5291  (evaluation_flag & EvaluationFlags::hessians) ||
5292  vector_ptr == nullptr ||
5293  shape_info.data.front().element_type >
5295  storage <
5297  return false;
5298  else
5299  return true;
5300  }
5301 
5302  private:
5303  template <int fe_degree>
5304  struct Processor
5305  {
5306  static const bool do_integrate = false;
5307  static const int dim_ = dim;
5308  static const int fe_degree_ = fe_degree;
5309  using VectorizedArrayType_ = VectorizedArrayType;
5310  using Number_ = Number;
5311  using Number2_ = const Number2;
5312 
5313  template <typename T0, typename T1, typename T2>
5314  void
5316  T0 & temp_2,
5317  const T1 src_ptr_1,
5318  const T1 src_ptr_2,
5319  const T2 &grad_weight)
5320  {
5321  do_vectorized_read(src_ptr_1, temp_1);
5322  do_vectorized_read(src_ptr_2, temp_2);
5323  temp_2 = grad_weight * (temp_1 - temp_2);
5324  }
5325 
5326  template <typename T1, typename T2>
5327  void
5328  value_vectorized(T1 &temp, const T2 src_ptr)
5329  {
5330  do_vectorized_read(src_ptr, temp);
5331  }
5332 
5333  template <typename T0, typename T1, typename T2, typename T3>
5334  void
5336  T0 & temp_2,
5337  const T1 src_ptr_1,
5338  const T1 src_ptr_2,
5339  const T2 &grad_weight,
5340  const T3 &indices_1,
5341  const T3 &indices_2)
5342  {
5343  do_vectorized_gather(src_ptr_1, indices_1, temp_1);
5344  do_vectorized_gather(src_ptr_2, indices_2, temp_2);
5345  temp_2 = grad_weight * (temp_1 - temp_2);
5346  }
5347 
5348  template <typename T0, typename T1, typename T2>
5349  void
5350  value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
5351  {
5352  do_vectorized_gather(src_ptr, indices, temp);
5353  }
5354 
5355  template <typename T0, typename T1, typename T2>
5356  void
5357  hermite_grad(T0 & temp_1,
5358  T0 & temp_2,
5359  const T1 &src_ptr_1,
5360  const T1 &src_ptr_2,
5361  const T2 &grad_weight)
5362  {
5363  // case 3a)
5364  temp_1 = src_ptr_1;
5365  temp_2 = grad_weight * (temp_1 - src_ptr_2);
5366  }
5367 
5368  template <typename T1, typename T2>
5369  void
5370  value(T1 &temp, const T2 &src_ptr)
5371  {
5372  // case 3b)
5373  temp = src_ptr;
5374  }
5375  };
5376  };
5377 
5378 
5379 
5380  template <int dim, typename Number2, typename VectorizedArrayType>
5382  {
5383  using Number = typename VectorizedArrayType::value_type;
5384 
5385  template <int fe_degree, int n_q_points_1d>
5386  static bool
5387  run(const unsigned int n_components,
5388  const EvaluationFlags::EvaluationFlags integration_flag,
5389  Number2 * dst_ptr,
5390  const std::vector<ArrayView<const Number2>> * sm_ptr,
5392  {
5393  Assert(fe_degree > -1, ExcInternalError());
5394  Assert(fe_eval.get_shape_info().element_type <=
5396  ExcInternalError());
5397 
5398  const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
5399 
5400  VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
5401  VectorizedArrayType *scratch_data =
5402  temp + 3 * n_components * dofs_per_face;
5403 
5404  const unsigned int subface_index = fe_eval.get_subface_index();
5405 
5406  // re-orientation for cases not possible with the io function below
5407  if (subface_index < GeometryInfo<dim>::max_children_per_cell)
5408  {
5409  if (fe_eval.get_dof_access_index() ==
5411  fe_eval.is_interior_face() == false)
5412  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5413  {
5414  // the loop breaks once an invalid_unsigned_int is hit for
5415  // all cases except the exterior faces in the ECL loop (where
5416  // some faces might be at the boundaries but others not)
5417  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
5418  continue;
5419 
5420  if (fe_eval.get_face_orientation(v) != 0)
5422  dim,
5423  n_components,
5424  v,
5425  integration_flag,
5427  fe_eval.get_face_orientation(v), 0),
5428  true,
5429  Utilities::pow(n_q_points_1d, dim - 1),
5430  &temp[0][0],
5431  fe_eval.begin_values(),
5432  fe_eval.begin_gradients(),
5433  fe_eval.begin_hessians());
5434  }
5435  else if (fe_eval.get_face_orientation() != 0)
5437  dim,
5438  n_components,
5439  integration_flag,
5441  fe_eval.get_face_orientation(), 0),
5442  true,
5443  Utilities::pow(n_q_points_1d, dim - 1),
5444  temp,
5445  fe_eval.begin_values(),
5446  fe_eval.begin_gradients(),
5447  fe_eval.begin_hessians());
5448  }
5449 
5450  if (fe_degree > -1 && fe_eval.get_subface_index() >=
5451  GeometryInfo<dim - 1>::max_children_per_cell)
5452  FEFaceEvaluationImpl<true,
5453  dim,
5454  fe_degree,
5455  n_q_points_1d,
5456  VectorizedArrayType>::
5457  integrate_in_face(n_components,
5458  integration_flag,
5459  fe_eval.get_shape_info().data.front(),
5460  temp,
5461  fe_eval.begin_values(),
5462  fe_eval.begin_gradients(),
5463  fe_eval.begin_hessians(),
5464  scratch_data,
5465  subface_index);
5466  else
5467  FEFaceEvaluationImpl<false,
5468  dim,
5469  fe_degree,
5470  n_q_points_1d,
5471  VectorizedArrayType>::
5472  integrate_in_face(n_components,
5473  integration_flag,
5474  fe_eval.get_shape_info().data.front(),
5475  temp,
5476  fe_eval.begin_values(),
5477  fe_eval.begin_gradients(),
5478  fe_eval.begin_hessians(),
5479  scratch_data,
5480  subface_index);
5481 
5483 
5484  if (fe_eval.get_dof_access_index() ==
5486  fe_eval.is_interior_face() == false)
5487  fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
5488  p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
5489  else
5490  fe_face_evaluation_process_and_io<1>(
5491  p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
5492 
5493  return false;
5494  }
5495 
5496  private:
5497  template <int fe_degree>
5498  struct Processor
5499  {
5500  static const bool do_integrate = true;
5501  static const int dim_ = dim;
5502  static const int fe_degree_ = fe_degree;
5503  using VectorizedArrayType_ = VectorizedArrayType;
5504  using Number_ = Number;
5505  using Number2_ = Number2;
5506 
5507  template <typename T0, typename T1, typename T2, typename T3, typename T4>
5508  void
5509  hermite_grad_vectorized(const T0 &temp_1,
5510  const T1 &temp_2,
5511  T2 dst_ptr_1,
5512  T3 dst_ptr_2,
5513  const T4 &grad_weight)
5514  {
5515  // case 1a)
5516  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
5517  const VectorizedArrayType grad = grad_weight * temp_2;
5518  do_vectorized_add(val, dst_ptr_1);
5519  do_vectorized_add(grad, dst_ptr_2);
5520  }
5521 
5522  template <typename T0, typename T1>
5523  void
5524  value_vectorized(const T0 &temp, T1 dst_ptr)
5525  {
5526  // case 1b)
5527  do_vectorized_add(temp, dst_ptr);
5528  }
5529 
5530  template <typename T0, typename T1, typename T2, typename T3>
5531  void
5533  const T0 &temp_2,
5534  T1 dst_ptr_1,
5535  T1 dst_ptr_2,
5536  const T2 &grad_weight,
5537  const T3 &indices_1,
5538  const T3 &indices_2)
5539  {
5540  // case 2a)
5541  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
5542  const VectorizedArrayType grad = grad_weight * temp_2;
5543  do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
5544  do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
5545  }
5546 
5547  template <typename T0, typename T1, typename T2>
5548  void
5549  value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
5550  {
5551  // case 2b)
5552  do_vectorized_scatter_add(temp, indices, dst_ptr);
5553  }
5554 
5555  template <typename T0, typename T1, typename T2>
5556  void
5557  hermite_grad(const T0 &temp_1,
5558  const T0 &temp_2,
5559  T1 & dst_ptr_1,
5560  T1 & dst_ptr_2,
5561  const T2 &grad_weight)
5562  {
5563  // case 3a)
5564  const Number val = temp_1 - grad_weight * temp_2;
5565  const Number grad = grad_weight * temp_2;
5566  dst_ptr_1 += val;
5567  dst_ptr_2 += grad;
5568  }
5569 
5570  template <typename T0, typename T1>
5571  void
5572  value(const T0 &temp, T1 &dst_ptr)
5573  {
5574  // case 3b)
5575  dst_ptr += temp;
5576  }
5577  };
5578  };
5579 
5580 
5581 
5586  template <int dim, typename Number>
5588  {
5589  template <int fe_degree, int = 0>
5590  static bool
5591  run(const unsigned int n_components,
5592  const FEEvaluationData<dim, Number, false> &fe_eval,
5593  const Number * in_array,
5594  Number * out_array,
5595  std::enable_if_t<fe_degree != -1> * = nullptr)
5596  {
5597  constexpr unsigned int dofs_per_component =
5598  Utilities::pow(fe_degree + 1, dim);
5599 
5600  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
5601  Assert(fe_eval.get_shape_info().element_type <=
5603  ExcNotImplemented());
5604 
5606  dim,
5607  fe_degree + 1,
5608  fe_degree + 1,
5609  Number>
5610  evaluator(
5613  fe_eval.get_shape_info().data.front().inverse_shape_values_eo);
5614 
5615  for (unsigned int d = 0; d < n_components; ++d)
5616  {
5617  const Number *in = in_array + d * dofs_per_component;
5618  Number * out = out_array + d * dofs_per_component;
5619  // Need to select 'apply' method with hessian slot because values
5620  // assume symmetries that do not exist in the inverse shapes
5621  evaluator.template hessians<0, true, false>(in, out);
5622  if (dim > 1)
5623  evaluator.template hessians<1, true, false>(out, out);
5624  if (dim > 2)
5625  evaluator.template hessians<2, true, false>(out, out);
5626  }
5627  for (unsigned int q = 0; q < dofs_per_component; ++q)
5628  {
5629  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
5630  for (unsigned int d = 0; d < n_components; ++d)
5631  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
5632  }
5633  for (unsigned int d = 0; d < n_components; ++d)
5634  {
5635  Number *out = out_array + d * dofs_per_component;
5636  if (dim > 2)
5637  evaluator.template hessians<2, false, false>(out, out);
5638  if (dim > 1)
5639  evaluator.template hessians<1, false, false>(out, out);
5640  evaluator.template hessians<0, false, false>(out, out);
5641  }
5642  return false;
5643  }
5644 
5645  template <int fe_degree, int = 0>
5646  static bool
5647  run(const unsigned int n_components,
5648  const FEEvaluationData<dim, Number, false> &fe_eval,
5649  const Number * in_array,
5650  Number * out_array,
5651  std::enable_if_t<fe_degree == -1> * = nullptr)
5652  {
5653  static_assert(fe_degree == -1, "Only usable for degree -1");
5654  const unsigned int dofs_per_component =
5656 
5657  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
5658 
5660  fe_eval.get_shape_info().data.front().inverse_shape_values,
5663  fe_eval.get_shape_info().data.front().fe_degree + 1,
5664  fe_eval.get_shape_info().data.front().fe_degree + 1);
5665 
5666  for (unsigned int d = 0; d < n_components; ++d)
5667  {
5668  const Number *in = in_array + d * dofs_per_component;
5669  Number * out = out_array + d * dofs_per_component;
5670  // Need to select 'apply' method with hessian slot because values
5671  // assume symmetries that do not exist in the inverse shapes
5672  evaluator.template values<0, true, false>(in, out);
5673  if (dim > 1)
5674  evaluator.template values<1, true, false>(out, out);
5675  if (dim > 2)
5676  evaluator.template values<2, true, false>(out, out);
5677  }
5678  for (unsigned int q = 0; q < dofs_per_component; ++q)
5679  {
5680  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
5681  for (unsigned int d = 0; d < n_components; ++d)
5682  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
5683  }
5684  for (unsigned int d = 0; d < n_components; ++d)
5685  {
5686  Number *out = out_array + d * dofs_per_component;
5687  if (dim > 2)
5688  evaluator.template values<2, false, false>(out, out);
5689  if (dim > 1)
5690  evaluator.template values<1, false, false>(out, out);
5691  evaluator.template values<0, false, false>(out, out);
5692  }
5693  return false;
5694  }
5695  };
5696 
5697 
5698 
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 AlignedVector<Number> &inverse_shape,
5710  const AlignedVector<Number> &inverse_coefficients,
5711  const Number * in_array,
5712  Number * out_array,
5713  std::enable_if_t<fe_degree != -1> * = nullptr)
5714  {
5715  constexpr unsigned int dofs_per_component =
5716  Utilities::pow(fe_degree + 1, dim);
5717  Assert(inverse_coefficients.size() > 0 &&
5718  inverse_coefficients.size() % dofs_per_component == 0,
5719  ExcMessage(
5720  "Expected diagonal to be a multiple of scalar dof per cells"));
5721  if (inverse_coefficients.size() != dofs_per_component)
5722  AssertDimension(n_desired_components * dofs_per_component,
5723  inverse_coefficients.size());
5724 
5725  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
5726 
5728  dim,
5729  fe_degree + 1,
5730  fe_degree + 1,
5731  Number>
5732  evaluator(AlignedVector<Number>(),
5734  inverse_shape);
5735 
5736  const unsigned int shift_coefficient =
5737  inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
5738  0;
5739  const Number *inv_coefficient = inverse_coefficients.data();
5740  for (unsigned int d = 0; d < n_desired_components; ++d)
5741  {
5742  const Number *in = in_array + d * dofs_per_component;
5743  Number * out = out_array + d * dofs_per_component;
5744  // Need to select 'apply' method with hessian slot because values
5745  // assume symmetries that do not exist in the inverse shapes
5746  evaluator.template hessians<0, true, false>(in, out);
5747  if (dim > 1)
5748  evaluator.template hessians<1, true, false>(out, out);
5749  if (dim > 2)
5750  evaluator.template hessians<2, true, false>(out, out);
5751 
5752  for (unsigned int q = 0; q < dofs_per_component; ++q)
5753  out[q] *= inv_coefficient[q];
5754 
5755  if (dim > 2)
5756  evaluator.template hessians<2, false, false>(out, out);
5757  if (dim > 1)
5758  evaluator.template hessians<1, false, false>(out, out);
5759  evaluator.template hessians<0, false, false>(out, out);
5760 
5761  inv_coefficient += shift_coefficient;
5762  }
5763  return false;
5764  }
5765 
5769  template <int fe_degree, int = 0>
5770  static bool
5771  run(const unsigned int,
5772  const AlignedVector<Number> &,
5773  const AlignedVector<Number> &,
5774  const Number *,
5775  Number *,
5776  std::enable_if_t<fe_degree == -1> * = nullptr)
5777  {
5778  static_assert(fe_degree == -1, "Only usable for degree -1");
5779  Assert(false, ExcNotImplemented());
5780  return false;
5781  }
5782  };
5783 
5784 
5785 
5790  template <int dim, typename Number>
5792  {
5793  template <int fe_degree, int n_q_points_1d>
5794  static bool
5795  run(const unsigned int n_desired_components,
5796  const FEEvaluationData<dim, Number, false> &fe_eval,
5797  const Number * in_array,
5798  Number * out_array)
5799  {
5800  static const bool do_inplace =
5801  fe_degree > -1 && (fe_degree + 1 == n_q_points_1d);
5802 
5803  Assert(fe_eval.get_shape_info().element_type !=
5805  ExcNotImplemented());
5806 
5807  const auto &inverse_shape =
5808  do_inplace ?
5809  fe_eval.get_shape_info().data.front().inverse_shape_values_eo :
5810  fe_eval.get_shape_info().data.front().inverse_shape_values;
5811 
5812  const std::size_t dofs_per_component =
5813  do_inplace ? Utilities::pow(fe_degree + 1, dim) :
5815  const std::size_t n_q_points = do_inplace ?
5816  Utilities::pow(fe_degree + 1, dim) :
5817  fe_eval.get_shape_info().n_q_points;
5818 
5820  dim,
5821  fe_degree + 1,
5822  n_q_points_1d,
5823  Number>
5824  evaluator(AlignedVector<Number>(),
5826  inverse_shape,
5827  fe_eval.get_shape_info().data.front().fe_degree + 1,
5828  fe_eval.get_shape_info().data.front().n_q_points_1d);
5829 
5830  for (unsigned int d = 0; d < n_desired_components; ++d)
5831  {
5832  const Number *in = in_array + d * n_q_points;
5833  Number * out = out_array + d * dofs_per_component;
5834 
5835  auto temp_1 = do_inplace ? out : fe_eval.get_scratch_data().begin();
5836  auto temp_2 = do_inplace ?
5837  out :
5838  (temp_1 + std::max(n_q_points, dofs_per_component));
5839 
5840  if (dim == 3)
5841  {
5842  evaluator.template hessians<2, false, false>(in, temp_1);
5843  evaluator.template hessians<1, false, false>(temp_1, temp_2);
5844  evaluator.template hessians<0, false, false>(temp_2, out);
5845  }
5846  if (dim == 2)
5847  {
5848  evaluator.template hessians<1, false, false>(in, temp_1);
5849  evaluator.template hessians<0, false, false>(temp_1, out);
5850  }
5851  if (dim == 1)
5852  evaluator.template hessians<0, false, false>(in, out);
5853  }
5854  return false;
5855  }
5856  };
5857 
5858 } // end of namespace internal
5859 
5860 
5862 
5863 #endif
pointer data()
size_type size() const
iterator begin() const
Definition: array_view.h:585
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:102
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:142
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
unsigned int cell_index
Definition: grid_tools.cc:1177
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
EvaluationFlags
The EvaluationFlags enum.
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2098
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 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:450
T fixed_power(const T t)
Definition: utilities.h:969
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:206
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, const AlignedVector< Number > &, const AlignedVector< Number > &, const Number *, Number *, std::enable_if_t< fe_degree==-1 > *=nullptr)
static bool run(const unsigned int n_desired_components, const AlignedVector< Number > &inverse_shape, const AlignedVector< Number > &inverse_coefficients, const Number *in_array, Number *out_array, 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 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)
static const EvaluatorVariant variant
static Eval create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number > *univariate_shape_data)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs_actual, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
EvaluatorTensorProduct< variant, dim, fe_degree+1, n_q_points_1d, Number > Eval
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs_actual, 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, true > &fe_eval)
void hermite_grad(T0 &temp_1, T0 &temp_2, const T1 &src_ptr_1, const T1 &src_ptr_2, const T2 &grad_weight)
void value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
void hermite_grad_vectorized(T0 &temp_1, T0 &temp_2, const T1 src_ptr_1, const T1 src_ptr_2, const T2 &grad_weight)
void hermite_grad_vectorized_indexed(T0 &temp_1, T0 &temp_2, const T1 src_ptr_1, const T1 src_ptr_2, const T2 &grad_weight, const T3 &indices_1, const T3 &indices_2)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number2 *src_ptr, const std::vector< ArrayView< const Number2 >> *sm_ptr, FEEvaluationData< dim, VectorizedArrayType, true > &fe_eval)
Definition: