Reference documentation for deal.II version GIT 6da2e5d553 2022-07-01 18: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 && element_type <= ElementType::tensor_symmetric)
2307  {
2309  dim,
2310  fe_degree,
2311  n_q_points_1d,
2312  Number>::evaluate(n_components,
2313  evaluation_flag,
2314  values_dofs,
2315  fe_eval);
2316  }
2317  else if (element_type == ElementType::tensor_symmetric_plus_dg0)
2318  {
2320  dim,
2321  fe_degree,
2322  n_q_points_1d,
2323  Number>::evaluate(n_components,
2324  evaluation_flag,
2325  values_dofs,
2326  fe_eval);
2327  }
2328  else if (element_type == ElementType::truncated_tensor)
2329  {
2331  dim,
2332  fe_degree,
2333  n_q_points_1d,
2334  Number>::evaluate(n_components,
2335  evaluation_flag,
2336  values_dofs,
2337  fe_eval);
2338  }
2339  else if (element_type == ElementType::tensor_none)
2340  {
2342  dim,
2343  fe_degree,
2344  n_q_points_1d,
2345  Number>::evaluate(n_components,
2346  evaluation_flag,
2347  values_dofs,
2348  fe_eval);
2349  }
2350  else if (element_type == ElementType::tensor_raviart_thomas)
2351  {
2354  dim,
2355  (fe_degree == -1) ? 1 : fe_degree,
2356  (n_q_points_1d < 1) ? 1 : n_q_points_1d,
2357  Number>::template evaluate_or_integrate<false>(evaluation_flag,
2358  const_cast<Number *>(
2359  values_dofs),
2360  fe_eval);
2361  }
2362  else
2363  {
2365  dim,
2366  fe_degree,
2367  n_q_points_1d,
2368  Number>::evaluate(n_components,
2369  evaluation_flag,
2370  values_dofs,
2371  fe_eval);
2372  }
2373 
2374  return false;
2375  }
2376  };
2377 
2378 
2379 
2395  template <int dim, typename Number>
2397  {
2398  template <int fe_degree, int n_q_points_1d>
2399  static bool
2400  run(const unsigned int n_components,
2401  const EvaluationFlags::EvaluationFlags integration_flag,
2402  Number * values_dofs,
2404  const bool sum_into_values_array)
2405  {
2406  const auto element_type = fe_eval.get_shape_info().element_type;
2408 
2409  Assert(fe_eval.get_shape_info().data.size() == 1 ||
2410  (fe_eval.get_shape_info().data.size() == dim &&
2411  element_type == ElementType::tensor_general) ||
2412  element_type == ElementType::tensor_raviart_thomas,
2413  ExcNotImplemented());
2414 
2415  if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
2417  {
2419  n_components,
2420  integration_flag,
2421  values_dofs,
2422  fe_eval,
2423  sum_into_values_array);
2424  }
2425  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
2426  // shape_info.h for more details
2427  else if (fe_degree >= 0 &&
2428  use_collocation_evaluation(fe_degree, n_q_points_1d) &&
2429  element_type <= ElementType::tensor_symmetric)
2430  {
2432  dim,
2433  fe_degree,
2434  n_q_points_1d,
2435  Number>::integrate(n_components,
2436  integration_flag,
2437  values_dofs,
2438  fe_eval,
2439  sum_into_values_array);
2440  }
2441  else if (fe_degree >= 0 && element_type <= ElementType::tensor_symmetric)
2442  {
2444  dim,
2445  fe_degree,
2446  n_q_points_1d,
2447  Number>::integrate(n_components,
2448  integration_flag,
2449  values_dofs,
2450  fe_eval,
2451  sum_into_values_array);
2452  }
2453  else if (element_type == ElementType::tensor_symmetric_plus_dg0)
2454  {
2456  dim,
2457  fe_degree,
2458  n_q_points_1d,
2459  Number>::integrate(n_components,
2460  integration_flag,
2461  values_dofs,
2462  fe_eval,
2463  sum_into_values_array);
2464  }
2465  else if (element_type == ElementType::truncated_tensor)
2466  {
2468  dim,
2469  fe_degree,
2470  n_q_points_1d,
2471  Number>::integrate(n_components,
2472  integration_flag,
2473  values_dofs,
2474  fe_eval,
2475  sum_into_values_array);
2476  }
2477  else if (element_type == ElementType::tensor_none)
2478  {
2480  dim,
2481  fe_degree,
2482  n_q_points_1d,
2483  Number>::integrate(n_components,
2484  integration_flag,
2485  values_dofs,
2486  fe_eval,
2487  sum_into_values_array);
2488  }
2489  else if (element_type == ElementType::tensor_raviart_thomas)
2490  {
2492  dim,
2493  (fe_degree == -1) ? 1 : fe_degree,
2494  (n_q_points_1d < 1) ? 1 : n_q_points_1d,
2495  Number>::
2496  template evaluate_or_integrate<true>(integration_flag,
2497  const_cast<Number *>(
2498  values_dofs),
2499  fe_eval,
2500  sum_into_values_array);
2501  }
2502  else
2503  {
2505  dim,
2506  fe_degree,
2507  n_q_points_1d,
2508  Number>::integrate(n_components,
2509  integration_flag,
2510  values_dofs,
2511  fe_eval,
2512  sum_into_values_array);
2513  }
2514 
2515  return false;
2516  }
2517  };
2518 
2519 
2520 
2521  template <bool symmetric_evaluate,
2522  int dim,
2523  int fe_degree,
2524  int n_q_points_1d,
2525  typename Number>
2527  {
2528  // We enable a transformation to collocation for derivatives if it gives
2529  // correct results (first two conditions), if it is the most efficient
2530  // choice in terms of operation counts (third condition) and if we were
2531  // able to initialize the fields in shape_info.templates.h from the
2532  // polynomials (fourth condition).
2533  using Eval = EvaluatorTensorProduct<symmetric_evaluate ? evaluate_evenodd :
2535  dim - 1,
2536  fe_degree + 1,
2537  n_q_points_1d,
2538  Number>;
2539 
2540  static Eval
2543  const unsigned int subface_index,
2544  const unsigned int direction)
2545  {
2546  if (symmetric_evaluate)
2547  return Eval(data.shape_values_eo,
2548  data.shape_gradients_eo,
2549  data.shape_hessians_eo,
2550  data.fe_degree + 1,
2551  data.n_q_points_1d);
2552  else if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
2553  return Eval(data.shape_values,
2554  data.shape_gradients,
2555  data.shape_hessians,
2556  data.fe_degree + 1,
2557  data.n_q_points_1d);
2558  else
2559  {
2560  const unsigned int index =
2561  direction == 0 ? subface_index % 2 : subface_index / 2;
2562  return Eval(data.values_within_subface[index],
2565  data.fe_degree + 1,
2566  data.n_q_points_1d);
2567  }
2568  }
2569 
2570  static void
2572  const unsigned int n_components,
2573  const EvaluationFlags::EvaluationFlags evaluation_flag,
2575  Number * values_dofs,
2576  Number * values_quad,
2577  Number * gradients_quad,
2578  Number * hessians_quad,
2579  Number * scratch_data,
2580  const unsigned int subface_index)
2581  {
2582  Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
2583  Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
2584 
2585  const std::size_t n_dofs = fe_degree > -1 ?
2586  Utilities::pow(fe_degree + 1, dim - 1) :
2587  Utilities::pow(data.fe_degree + 1, dim - 1);
2588  const std::size_t n_q_points =
2589  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
2590  Utilities::pow(data.n_q_points_1d, dim - 1);
2591 
2592  // keep a copy of the original pointer for the case of the Hessians
2593  Number *values_dofs_ptr = values_dofs;
2594 
2595  if ((evaluation_flag & EvaluationFlags::values) != 0u &&
2596  ((evaluation_flag & EvaluationFlags::gradients) == 0u))
2597  for (unsigned int c = 0; c < n_components; ++c)
2598  {
2599  switch (dim)
2600  {
2601  case 3:
2602  eval0.template values<0, true, false>(values_dofs,
2603  values_quad);
2604  eval1.template values<1, true, false>(values_quad,
2605  values_quad);
2606  break;
2607  case 2:
2608  eval0.template values<0, true, false>(values_dofs,
2609  values_quad);
2610  break;
2611  case 1:
2612  values_quad[0] = values_dofs[0];
2613  break;
2614  default:
2615  Assert(false, ExcNotImplemented());
2616  }
2617  // Note: we always keep storage of values, 1st and 2nd derivatives
2618  // in an array
2619  values_dofs += 3 * n_dofs;
2620  values_quad += n_q_points;
2621  }
2622  else if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2623  for (unsigned int c = 0; c < n_components; ++c)
2624  {
2625  switch (dim)
2626  {
2627  case 3:
2628  if (symmetric_evaluate &&
2629  use_collocation_evaluation(fe_degree, n_q_points_1d))
2630  {
2631  eval0.template values<0, true, false>(values_dofs,
2632  values_quad);
2633  eval0.template values<1, true, false>(values_quad,
2634  values_quad);
2636  dim - 1,
2637  n_q_points_1d,
2638  n_q_points_1d,
2639  Number>
2640  eval_grad(AlignedVector<Number>(),
2643  eval_grad.template gradients<0, true, false>(
2644  values_quad, gradients_quad);
2645  eval_grad.template gradients<1, true, false>(
2646  values_quad, gradients_quad + n_q_points);
2647  }
2648  else
2649  {
2650  // grad x
2651  eval0.template gradients<0, true, false>(values_dofs,
2652  scratch_data);
2653  eval1.template values<1, true, false>(scratch_data,
2654  gradients_quad);
2655 
2656  // grad y
2657  eval0.template values<0, true, false>(values_dofs,
2658  scratch_data);
2659  eval1.template gradients<1, true, false>(scratch_data,
2660  gradients_quad +
2661  n_q_points);
2662 
2663  if ((evaluation_flag & EvaluationFlags::values) != 0u)
2664  eval1.template values<1, true, false>(scratch_data,
2665  values_quad);
2666  }
2667  // grad z
2668  eval0.template values<0, true, false>(values_dofs + n_dofs,
2669  scratch_data);
2670  eval1.template values<1, true, false>(
2671  scratch_data, gradients_quad + (dim - 1) * n_q_points);
2672 
2673  break;
2674  case 2:
2675  eval0.template values<0, true, false>(values_dofs + n_dofs,
2676  gradients_quad +
2677  n_q_points);
2678  eval0.template gradients<0, true, false>(values_dofs,
2679  gradients_quad);
2680  if ((evaluation_flag & EvaluationFlags::values) != 0u)
2681  eval0.template values<0, true, false>(values_dofs,
2682  values_quad);
2683  break;
2684  case 1:
2685  values_quad[0] = values_dofs[0];
2686  gradients_quad[0] = values_dofs[1];
2687  break;
2688  default:
2689  AssertThrow(false, ExcNotImplemented());
2690  }
2691  values_dofs += 3 * n_dofs;
2692  values_quad += n_q_points;
2693  gradients_quad += dim * n_q_points;
2694  }
2695 
2696  if ((evaluation_flag & EvaluationFlags::hessians) != 0u)
2697  {
2698  values_dofs = values_dofs_ptr;
2699  for (unsigned int c = 0; c < n_components; ++c)
2700  {
2701  switch (dim)
2702  {
2703  case 3:
2704  // grad xx
2705  eval0.template hessians<0, true, false>(values_dofs,
2706  scratch_data);
2707  eval1.template values<1, true, false>(scratch_data,
2708  hessians_quad);
2709 
2710  // grad yy
2711  eval0.template values<0, true, false>(values_dofs,
2712  scratch_data);
2713  eval1.template hessians<1, true, false>(scratch_data,
2714  hessians_quad +
2715  n_q_points);
2716 
2717  // grad zz
2718  eval0.template values<0, true, false>(values_dofs +
2719  2 * n_dofs,
2720  scratch_data);
2721  eval1.template values<1, true, false>(scratch_data,
2722  hessians_quad +
2723  2 * n_q_points);
2724 
2725  // grad xy
2726  eval0.template gradients<0, true, false>(values_dofs,
2727  scratch_data);
2728  eval1.template gradients<1, true, false>(scratch_data,
2729  hessians_quad +
2730  3 * n_q_points);
2731 
2732  // grad xz
2733  eval0.template gradients<0, true, false>(values_dofs +
2734  n_dofs,
2735  scratch_data);
2736  eval1.template values<1, true, false>(scratch_data,
2737  hessians_quad +
2738  4 * n_q_points);
2739 
2740  // grad yz
2741  eval0.template values<0, true, false>(values_dofs + n_dofs,
2742  scratch_data);
2743  eval1.template gradients<1, true, false>(scratch_data,
2744  hessians_quad +
2745  5 * n_q_points);
2746 
2747  break;
2748  case 2:
2749  // grad xx
2750  eval0.template hessians<0, true, false>(values_dofs,
2751  hessians_quad);
2752  // grad yy
2753  eval0.template values<0, true, false>(
2754  values_dofs + 2 * n_dofs, hessians_quad + n_q_points);
2755  // grad xy
2756  eval0.template gradients<0, true, false>(
2757  values_dofs + n_dofs, hessians_quad + 2 * n_q_points);
2758  break;
2759  case 1:
2760  hessians_quad[0] = values_dofs[2];
2761  break;
2762  default:
2763  AssertThrow(false, ExcNotImplemented());
2764  }
2765  values_dofs += 3 * n_dofs;
2766  hessians_quad += dim * (dim + 1) / 2 * n_q_points;
2767  }
2768  }
2769  }
2770 
2771  static void
2773  const unsigned int n_components,
2774  const EvaluationFlags::EvaluationFlags integration_flag,
2776  Number * values_dofs,
2777  Number * values_quad,
2778  Number * gradients_quad,
2779  Number * hessians_quad,
2780  Number * scratch_data,
2781  const unsigned int subface_index)
2782  {
2783  Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
2784  Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
2785 
2786  const std::size_t n_dofs =
2787  fe_degree > -1 ?
2788  Utilities::pow(fe_degree + 1, dim - 1) :
2789  (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
2790  const std::size_t n_q_points =
2791  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
2792  Utilities::pow(data.n_q_points_1d, dim - 1);
2793 
2794  // keep a copy of the original pointer for the case of the Hessians
2795  Number *values_dofs_ptr = values_dofs;
2796 
2797  if ((integration_flag & EvaluationFlags::values) != 0u &&
2798  (integration_flag & EvaluationFlags::gradients) == 0u)
2799  for (unsigned int c = 0; c < n_components; ++c)
2800  {
2801  switch (dim)
2802  {
2803  case 3:
2804  eval1.template values<1, false, false>(values_quad,
2805  values_quad);
2806  eval0.template values<0, false, false>(values_quad,
2807  values_dofs);
2808  break;
2809  case 2:
2810  eval0.template values<0, false, false>(values_quad,
2811  values_dofs);
2812  break;
2813  case 1:
2814  values_dofs[0] = values_quad[0];
2815  break;
2816  default:
2817  Assert(false, ExcNotImplemented());
2818  }
2819  values_dofs += 3 * n_dofs;
2820  values_quad += n_q_points;
2821  }
2822  else if ((integration_flag & EvaluationFlags::gradients) != 0u)
2823  for (unsigned int c = 0; c < n_components; ++c)
2824  {
2825  switch (dim)
2826  {
2827  case 3:
2828  // grad z
2829  eval1.template values<1, false, false>(gradients_quad +
2830  2 * n_q_points,
2831  gradients_quad +
2832  2 * n_q_points);
2833  eval0.template values<0, false, false>(gradients_quad +
2834  2 * n_q_points,
2835  values_dofs + n_dofs);
2836  if (symmetric_evaluate &&
2837  use_collocation_evaluation(fe_degree, n_q_points_1d))
2838  {
2840  dim - 1,
2841  n_q_points_1d,
2842  n_q_points_1d,
2843  Number>
2844  eval_grad(AlignedVector<Number>(),
2847  if ((integration_flag & EvaluationFlags::values) != 0u)
2848  eval_grad.template gradients<1, false, true>(
2849  gradients_quad + n_q_points, values_quad);
2850  else
2851  eval_grad.template gradients<1, false, false>(
2852  gradients_quad + n_q_points, values_quad);
2853  eval_grad.template gradients<0, false, true>(
2854  gradients_quad, values_quad);
2855  eval0.template values<1, false, false>(values_quad,
2856  values_quad);
2857  eval0.template values<0, false, false>(values_quad,
2858  values_dofs);
2859  }
2860  else
2861  {
2862  if ((integration_flag & EvaluationFlags::values) != 0u)
2863  {
2864  eval1.template values<1, false, false>(values_quad,
2865  scratch_data);
2866  eval1.template gradients<1, false, true>(
2867  gradients_quad + n_q_points, scratch_data);
2868  }
2869  else
2870  eval1.template gradients<1, false, false>(
2871  gradients_quad + n_q_points, scratch_data);
2872 
2873  // grad y
2874  eval0.template values<0, false, false>(scratch_data,
2875  values_dofs);
2876 
2877  // grad x
2878  eval1.template values<1, false, false>(gradients_quad,
2879  scratch_data);
2880  eval0.template gradients<0, false, true>(scratch_data,
2881  values_dofs);
2882  }
2883  break;
2884  case 2:
2885  eval0.template values<0, false, false>(gradients_quad +
2886  n_q_points,
2887  values_dofs + n_dofs);
2888  eval0.template gradients<0, false, false>(gradients_quad,
2889  values_dofs);
2890  if ((integration_flag & EvaluationFlags::values) != 0u)
2891  eval0.template values<0, false, true>(values_quad,
2892  values_dofs);
2893  break;
2894  case 1:
2895  values_dofs[0] = values_quad[0];
2896  values_dofs[1] = gradients_quad[0];
2897  break;
2898  default:
2899  AssertThrow(false, ExcNotImplemented());
2900  }
2901  values_dofs += 3 * n_dofs;
2902  values_quad += n_q_points;
2903  gradients_quad += dim * n_q_points;
2904  }
2905 
2906  if ((integration_flag & EvaluationFlags::hessians) != 0u)
2907  {
2908  values_dofs = values_dofs_ptr;
2909  for (unsigned int c = 0; c < n_components; ++c)
2910  {
2911  switch (dim)
2912  {
2913  case 3:
2914  // grad xx
2915  eval1.template values<1, false, false>(hessians_quad,
2916  scratch_data);
2917  if ((integration_flag & (EvaluationFlags::values |
2919  eval0.template hessians<0, false, true>(scratch_data,
2920  values_dofs);
2921  else
2922  eval0.template hessians<0, false, false>(scratch_data,
2923  values_dofs);
2924 
2925  // grad yy
2926  eval1.template hessians<1, false, false>(hessians_quad +
2927  n_q_points,
2928  scratch_data);
2929  eval0.template values<0, false, true>(scratch_data,
2930  values_dofs);
2931 
2932  // grad zz
2933  eval1.template values<1, false, false>(hessians_quad +
2934  2 * n_q_points,
2935  scratch_data);
2936  eval0.template values<0, false, false>(scratch_data,
2937  values_dofs +
2938  2 * n_dofs);
2939 
2940  // grad xy
2941  eval1.template gradients<1, false, false>(hessians_quad +
2942  3 * n_q_points,
2943  scratch_data);
2944  eval0.template gradients<0, false, true>(scratch_data,
2945  values_dofs);
2946 
2947  // grad xz
2948  eval1.template values<1, false, false>(hessians_quad +
2949  4 * n_q_points,
2950  scratch_data);
2951  if ((integration_flag & EvaluationFlags::gradients) != 0u)
2952  eval0.template gradients<0, false, true>(scratch_data,
2953  values_dofs +
2954  n_dofs);
2955  else
2956  eval0.template gradients<0, false, false>(scratch_data,
2957  values_dofs +
2958  n_dofs);
2959 
2960  // grad yz
2961  eval1.template gradients<1, false, false>(hessians_quad +
2962  5 * n_q_points,
2963  scratch_data);
2964  eval0.template values<0, false, true>(scratch_data,
2965  values_dofs + n_dofs);
2966 
2967  break;
2968  case 2:
2969  // grad xx
2970  if ((integration_flag & (EvaluationFlags::values |
2972  eval0.template hessians<0, false, true>(hessians_quad,
2973  values_dofs);
2974  else
2975  eval0.template hessians<0, false, false>(hessians_quad,
2976  values_dofs);
2977 
2978  // grad yy
2979  eval0.template values<0, false, false>(
2980  hessians_quad + n_q_points, values_dofs + 2 * n_dofs);
2981  // grad xy
2982  if ((integration_flag & EvaluationFlags::gradients) != 0u)
2983  eval0.template gradients<0, false, true>(
2984  hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
2985  else
2986  eval0.template gradients<0, false, false>(
2987  hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
2988  break;
2989  case 1:
2990  values_dofs[2] = hessians_quad[0];
2991  if ((integration_flag & EvaluationFlags::values) == 0u)
2992  values_dofs[0] = 0;
2993  if ((integration_flag & EvaluationFlags::gradients) == 0u)
2994  values_dofs[1] = 0;
2995  break;
2996  default:
2997  AssertThrow(false, ExcNotImplemented());
2998  }
2999  values_dofs += 3 * n_dofs;
3000  hessians_quad += dim * (dim + 1) / 2 * n_q_points;
3001  }
3002  }
3003  }
3004  };
3005 
3006  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
3008  {
3010  dim - 1,
3011  fe_degree,
3012  n_q_points_1d,
3013  Number>;
3014  template <typename EvalType>
3015  static EvalType
3018  const unsigned int subface_index,
3019  const unsigned int direction)
3020  {
3021  if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
3022  return EvalType(data.shape_values,
3023  data.shape_gradients,
3024  data.shape_hessians);
3025  else
3026  {
3027  const unsigned int index =
3028  direction == 0 ? subface_index % 2 : subface_index / 2;
3029  return EvalType(data.values_within_subface[index],
3032  }
3033  }
3034 
3035  template <bool integrate>
3036  static void
3038  const EvaluationFlags::EvaluationFlags evaluation_flag,
3039  Number * values_dofs,
3041  Number * scratch_data,
3042  const unsigned int subface_index,
3043  const unsigned int face_no)
3044  {
3045  const unsigned int face_direction = face_no / 2;
3046 
3047  // We first evaluate the anisotropic faces, i.e the faces where
3048  // face_direction != component. Note that the call order here is not
3049  // important, since the pointers are shifted accordingly within the
3050  // function. However, this is the order in which the components will be in
3051  // the quadrature points. Furthermore, the isotropic faces have no "normal
3052  // direction" but we still pass in normal_dir = 2 since this is used for
3053  // the pointers.
3054  // -----------------------------------------------------------------------------------
3055  // | | Anisotropic faces | Isotropic faces|
3056  // | Face dir | comp, coords, normal_dir | comp, coords, normal_dir | comp, coords |
3057  // | --------------------------------------------------------------------------------|
3058  // | 0 | 1, y, 0 | - | 0, y |
3059  // | 1 | 0, x, 0 | - | 1, x |
3060  // | --------------------------------------------------------------------------------|
3061  // | 0 | 1, yz, 0 | 2, yz, 1 | 0, yz |
3062  // | 1 | 2, zx, 0 | 0, zx, 1 | 1, zx |
3063  // | 2 | 0, xy, 0 | 1, xy, 1 | 2, xy |
3064  // -----------------------------------------------------------------------------------
3065  evaluate_in_face_apply<0>(values_dofs,
3066  fe_eval,
3067  scratch_data,
3068  evaluation_flag,
3069  face_direction,
3070  subface_index,
3071  std::integral_constant<bool, integrate>());
3072 
3073  if (dim == 3)
3074  evaluate_in_face_apply<1>(values_dofs,
3075  fe_eval,
3076  scratch_data,
3077  evaluation_flag,
3078  face_direction,
3079  subface_index,
3080  std::integral_constant<bool, integrate>());
3081 
3082  evaluate_in_face_apply<2>(values_dofs,
3083  fe_eval,
3084  scratch_data,
3085  evaluation_flag,
3086  face_direction,
3087  subface_index,
3088  std::integral_constant<bool, integrate>());
3089  }
3090 
3091  /*
3092  * Helper function which applies the 1D kernels for on one
3093  * component in a face. normal_dir indicates the direction of the continuous
3094  * component of the RT space. std::integral_constant<bool, false> is the
3095  * evaluation path, and std::integral_constant<bool, true> below is the
3096  * integration path. These two functions can be fused together since all
3097  * offsets and pointers are the exact same.
3098  */
3099  template <int normal_dir>
3100  static inline void
3102  Number * values_dofs,
3104  Number * scratch_data,
3105  const EvaluationFlags::EvaluationFlags evaluation_flag,
3106  const unsigned int face_direction,
3107  const unsigned int subface_index,
3108  std::integral_constant<bool, false>)
3109  {
3110  using EvalNormal =
3112  dim - 1,
3113  (fe_degree == -1) ? 1 : fe_degree + 1,
3114  n_q_points_1d,
3115  Number,
3116  normal_dir>;
3117  using EvalTangent =
3119  dim - 1,
3120  (fe_degree == -1) ? 1 : fe_degree,
3121  n_q_points_1d,
3122  Number,
3123  normal_dir>;
3124 
3125  using TempEval0 = typename std::
3126  conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
3127  using TempEval1 = typename std::
3128  conditional<normal_dir == 0, EvalTangent, EvalNormal>::type;
3129  using Eval0 = typename std::
3130  conditional<normal_dir == 2, EvalGeneral, TempEval0>::type;
3131  using Eval1 = typename std::
3132  conditional<normal_dir == 2, EvalGeneral, TempEval1>::type;
3133 
3134  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info =
3135  fe_eval.get_shape_info();
3136  Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
3137  ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]),
3138  subface_index,
3139  0);
3140  Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
3141  ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]),
3142  subface_index,
3143  1);
3144 
3145  constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
3146  const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face;
3147  const std::size_t n_dofs_normal =
3148  n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
3149  const std::size_t dofs_stride =
3150  (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
3151  n_dofs_tangent;
3152 
3153  static constexpr ::ndarray<unsigned int, 3, 3> component_table = {
3154  {{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
3155  const unsigned int component =
3156  (dim == 2 && normal_dir == 0 && face_direction == 1) ?
3157  0 :
3158  component_table[face_direction][normal_dir];
3159 
3160  // Initial offsets
3161  values_dofs +=
3162  3 * ((component == 0) ?
3163  0 :
3164  ((component == 1) ?
3165  ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) :
3166  ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent :
3167  n_dofs_normal + n_dofs_tangent)));
3168  const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir;
3169  Number *values_quad = fe_eval.begin_values() + n_q_points * shift;
3170  Number *gradients_quad =
3171  fe_eval.begin_gradients() + dim * n_q_points * shift;
3172  Number *hessians_quad =
3173  fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * shift;
3174 
3175  // Evaluation path
3176  if ((evaluation_flag & EvaluationFlags::values) &&
3177  !(evaluation_flag & EvaluationFlags::gradients))
3178  {
3179  switch (dim)
3180  {
3181  case 3:
3182  eval0.template values<0, true, false>(values_dofs, values_quad);
3183  eval1.template values<1, true, false>(values_quad, values_quad);
3184  break;
3185  case 2:
3186  eval0.template values<0, true, false>(values_dofs, values_quad);
3187  break;
3188  default:
3189  Assert(false, ExcNotImplemented());
3190  }
3191  }
3192  else if (evaluation_flag & EvaluationFlags::gradients)
3193  {
3194  switch (dim)
3195  {
3196  case 3:
3197  // grad x
3198  eval0.template gradients<0, true, false>(values_dofs,
3199  scratch_data);
3200  eval1.template values<1, true, false>(scratch_data,
3201  gradients_quad);
3202 
3203  // grad y
3204  eval0.template values<0, true, false>(values_dofs,
3205  scratch_data);
3206  eval1.template gradients<1, true, false>(scratch_data,
3207  gradients_quad +
3208  n_q_points);
3209 
3210  if (evaluation_flag & EvaluationFlags::values)
3211  eval1.template values<1, true, false>(scratch_data,
3212  values_quad);
3213 
3214  // grad z
3215  eval0.template values<0, true, false>(values_dofs + dofs_stride,
3216  scratch_data);
3217  eval1.template values<1, true, false>(scratch_data,
3218  gradients_quad +
3219  2 * n_q_points);
3220 
3221  break;
3222  case 2:
3223  eval0.template values<0, true, false>(values_dofs + dofs_stride,
3224  gradients_quad +
3225  n_q_points);
3226  eval0.template gradients<0, true, false>(values_dofs,
3227  gradients_quad);
3228  if ((evaluation_flag & EvaluationFlags::values))
3229  eval0.template values<0, true, false>(values_dofs,
3230  values_quad);
3231  break;
3232  default:
3233  AssertThrow(false, ExcNotImplemented());
3234  }
3235  }
3236 
3237  if (evaluation_flag & EvaluationFlags::hessians)
3238  {
3239  switch (dim)
3240  {
3241  case 3:
3242  // grad xx
3243  eval0.template hessians<0, true, false>(values_dofs,
3244  scratch_data);
3245  eval1.template values<1, true, false>(scratch_data,
3246  hessians_quad);
3247 
3248  // grad yy
3249  eval0.template values<0, true, false>(values_dofs,
3250  scratch_data);
3251  eval1.template hessians<1, true, false>(scratch_data,
3252  hessians_quad +
3253  n_q_points);
3254 
3255  // grad zz
3256  eval0.template values<0, true, false>(values_dofs +
3257  2 * dofs_stride,
3258  scratch_data);
3259  eval1.template values<1, true, false>(scratch_data,
3260  hessians_quad +
3261  2 * n_q_points);
3262 
3263  // grad xy
3264  eval0.template gradients<0, true, false>(values_dofs,
3265  scratch_data);
3266  eval1.template gradients<1, true, false>(scratch_data,
3267  hessians_quad +
3268  3 * n_q_points);
3269 
3270  // grad xz
3271  eval0.template gradients<0, true, false>(values_dofs +
3272  dofs_stride,
3273  scratch_data);
3274  eval1.template values<1, true, false>(scratch_data,
3275  hessians_quad +
3276  4 * n_q_points);
3277 
3278  // grad yz
3279  eval0.template values<0, true, false>(values_dofs + dofs_stride,
3280  scratch_data);
3281  eval1.template gradients<1, true, false>(scratch_data,
3282  hessians_quad +
3283  5 * n_q_points);
3284 
3285  break;
3286  case 2:
3287  // grad xx
3288  eval0.template hessians<0, true, false>(values_dofs,
3289  hessians_quad);
3290  // grad yy
3291  eval0.template values<0, true, false>(
3292  values_dofs + 2 * dofs_stride, hessians_quad + n_q_points);
3293  // grad xy
3294  eval0.template gradients<0, true, false>(
3295  values_dofs + dofs_stride, hessians_quad + 2 * n_q_points);
3296  break;
3297  default:
3298  AssertThrow(false, ExcNotImplemented());
3299  }
3300  }
3301  }
3302 
3303  template <int normal_dir>
3304  static inline void
3306  Number * values_dofs,
3308  Number * scratch_data,
3309  const EvaluationFlags::EvaluationFlags evaluation_flag,
3310  const unsigned int face_direction,
3311  const unsigned int subface_index,
3312  std::integral_constant<bool, true>)
3313  {
3314  using EvalNormal =
3316  dim - 1,
3317  (fe_degree == -1) ? 1 : fe_degree + 1,
3318  n_q_points_1d,
3319  Number,
3320  normal_dir>;
3321  using EvalTangent =
3323  dim - 1,
3324  (fe_degree == -1) ? 1 : fe_degree,
3325  n_q_points_1d,
3326  Number,
3327  normal_dir>;
3328 
3329  using TempEval0 = typename std::
3330  conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
3331  using TempEval1 = typename std::
3332  conditional<normal_dir == 0, EvalTangent, EvalNormal>::type;
3333  using Eval0 = typename std::
3334  conditional<normal_dir == 2, EvalGeneral, TempEval0>::type;
3335  using Eval1 = typename std::
3336  conditional<normal_dir == 2, EvalGeneral, TempEval1>::type;
3337 
3338  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info =
3339  fe_eval.get_shape_info();
3340  Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
3341  ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]),
3342  subface_index,
3343  0);
3344  Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
3345  ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]),
3346  subface_index,
3347  1);
3348 
3349  constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
3350  const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face;
3351  const std::size_t n_dofs_normal =
3352  n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
3353  const std::size_t dofs_stride =
3354  (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
3355  n_dofs_tangent;
3356 
3357  static constexpr ::ndarray<unsigned int, 3, 3> component_table = {
3358  {{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
3359  const unsigned int component =
3360  (dim == 2 && normal_dir == 0 && face_direction == 1) ?
3361  0 :
3362  component_table[face_direction][normal_dir];
3363 
3364  // Initial offsets
3365  values_dofs +=
3366  3 * ((component == 0) ?
3367  0 :
3368  ((component == 1) ?
3369  ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) :
3370  ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent :
3371  n_dofs_normal + n_dofs_tangent)));
3372  const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir;
3373  Number *values_quad = fe_eval.begin_values() + n_q_points * shift;
3374  Number *gradients_quad =
3375  fe_eval.begin_gradients() + dim * n_q_points * shift;
3376  Number *hessians_quad =
3377  fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * shift;
3378 
3379  // Integration path
3380  if ((evaluation_flag & EvaluationFlags::values) &&
3381  !(evaluation_flag & EvaluationFlags::gradients))
3382  {
3383  switch (dim)
3384  {
3385  case 3:
3386  eval1.template values<1, false, false>(values_quad,
3387  values_quad);
3388  eval0.template values<0, false, false>(values_quad,
3389  values_dofs);
3390  break;
3391  case 2:
3392  eval0.template values<0, false, false>(values_quad,
3393  values_dofs);
3394  break;
3395  default:
3396  Assert(false, ExcNotImplemented());
3397  }
3398  }
3399  else if (evaluation_flag & EvaluationFlags::gradients)
3400  {
3401  switch (dim)
3402  {
3403  case 3:
3404  // grad z
3405  eval1.template values<1, false, false>(gradients_quad +
3406  2 * n_q_points,
3407  gradients_quad +
3408  2 * n_q_points);
3409  eval0.template values<0, false, false>(
3410  gradients_quad + 2 * n_q_points, values_dofs + dofs_stride);
3411 
3412  if (evaluation_flag & EvaluationFlags::values)
3413  {
3414  eval1.template values<1, false, false>(values_quad,
3415  scratch_data);
3416  eval1.template gradients<1, false, true>(gradients_quad +
3417  n_q_points,
3418  scratch_data);
3419  }
3420  else
3421  eval1.template gradients<1, false, false>(gradients_quad +
3422  n_q_points,
3423  scratch_data);
3424 
3425  // grad y
3426  eval0.template values<0, false, false>(scratch_data,
3427  values_dofs);
3428 
3429  // grad x
3430  eval1.template values<1, false, false>(gradients_quad,
3431  scratch_data);
3432  eval0.template gradients<0, false, true>(scratch_data,
3433  values_dofs);
3434 
3435  break;
3436  case 2:
3437  eval0.template values<0, false, false>(
3438  gradients_quad + n_q_points, values_dofs + dofs_stride);
3439  eval0.template gradients<0, false, false>(gradients_quad,
3440  values_dofs);
3441  if (evaluation_flag & EvaluationFlags::values)
3442  eval0.template values<0, false, true>(values_quad,
3443  values_dofs);
3444  break;
3445  default:
3446  AssertThrow(false, ExcNotImplemented());
3447  }
3448  }
3449 
3450  if (evaluation_flag & EvaluationFlags::hessians)
3451  {
3452  switch (dim)
3453  {
3454  case 3:
3455  // grad xx
3456  eval1.template values<1, false, false>(hessians_quad,
3457  scratch_data);
3458  if ((evaluation_flag &
3460  eval0.template hessians<0, false, true>(scratch_data,
3461  values_dofs);
3462  else
3463  eval0.template hessians<0, false, false>(scratch_data,
3464  values_dofs);
3465 
3466  // grad yy
3467  eval1.template hessians<1, false, false>(hessians_quad +
3468  n_q_points,
3469  scratch_data);
3470  eval0.template values<0, false, true>(scratch_data,
3471  values_dofs);
3472 
3473  // grad zz
3474  eval1.template values<1, false, false>(hessians_quad +
3475  2 * n_q_points,
3476  scratch_data);
3477  eval0.template values<0, false, false>(scratch_data,
3478  values_dofs +
3479  2 * dofs_stride);
3480 
3481  // grad xy
3482  eval1.template gradients<1, false, false>(hessians_quad +
3483  3 * n_q_points,
3484  scratch_data);
3485  eval0.template gradients<0, false, true>(scratch_data,
3486  values_dofs);
3487 
3488  // grad xz
3489  eval1.template values<1, false, false>(hessians_quad +
3490  4 * n_q_points,
3491  scratch_data);
3492  if ((evaluation_flag & EvaluationFlags::gradients))
3493  eval0.template gradients<0, false, true>(scratch_data,
3494  values_dofs +
3495  dofs_stride);
3496  else
3497  eval0.template gradients<0, false, false>(scratch_data,
3498  values_dofs +
3499  dofs_stride);
3500 
3501  // grad yz
3502  eval1.template gradients<1, false, false>(hessians_quad +
3503  5 * n_q_points,
3504  scratch_data);
3505  eval0.template values<0, false, true>(scratch_data,
3506  values_dofs +
3507  dofs_stride);
3508 
3509  break;
3510  case 2:
3511  // grad xx
3512  if (evaluation_flag &
3514  eval0.template hessians<0, false, true>(hessians_quad,
3515  values_dofs);
3516  else
3517  eval0.template hessians<0, false, false>(hessians_quad,
3518  values_dofs);
3519 
3520  // grad yy
3521  eval0.template values<0, false, false>(
3522  hessians_quad + n_q_points, values_dofs + 2 * dofs_stride);
3523  // grad xy
3524  if ((evaluation_flag & EvaluationFlags::gradients))
3525  eval0.template gradients<0, false, true>(
3526  hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
3527  else
3528  eval0.template gradients<0, false, false>(
3529  hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
3530  break;
3531  default:
3532  AssertThrow(false, ExcNotImplemented());
3533  }
3534  }
3535  }
3536  };
3537 
3538 
3539  template <int dim, int fe_degree, typename Number>
3541  {
3542  template <bool do_evaluate, bool add_into_output>
3543  static void
3544  interpolate(const unsigned int n_components,
3546  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
3547  const Number * input,
3548  Number * output,
3549  const unsigned int face_no)
3550  {
3551  Assert(static_cast<unsigned int>(fe_degree) ==
3552  shape_info.data.front().fe_degree ||
3553  fe_degree == -1,
3554  ExcInternalError());
3556  interpolate_generic_raviart_thomas<do_evaluate, add_into_output>(
3557  n_components, input, output, flags, face_no, shape_info);
3558  else
3559  interpolate_generic<do_evaluate, add_into_output>(
3560  n_components,
3561  input,
3562  output,
3563  flags,
3564  face_no,
3565  shape_info.data.front().fe_degree + 1,
3566  shape_info.data.front().shape_data_on_face,
3567  shape_info.dofs_per_component_on_cell,
3568  3 * shape_info.dofs_per_component_on_face);
3569  }
3570 
3574  template <bool do_evaluate, bool add_into_output>
3575  static void
3577  const unsigned int n_components,
3579  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
3580  const Number * input,
3581  Number * output,
3582  const unsigned int face_no)
3583  {
3584  Assert(static_cast<unsigned int>(fe_degree + 1) ==
3585  shape_info.data.front().n_q_points_1d ||
3586  fe_degree == -1,
3587  ExcInternalError());
3588 
3589  interpolate_generic<do_evaluate, add_into_output>(
3590  n_components,
3591  input,
3592  output,
3593  flags,
3594  face_no,
3595  shape_info.data.front().quadrature.size(),
3596  shape_info.data.front().quadrature_data_on_face,
3597  shape_info.n_q_points,
3598  shape_info.n_q_points_face);
3599  }
3600 
3601  private:
3602  template <bool do_evaluate, bool add_into_output, int face_direction = 0>
3603  static void
3604  interpolate_generic(const unsigned int n_components,
3605  const Number * input,
3606  Number * output,
3608  const unsigned int face_no,
3609  const unsigned int n_points_1d,
3610  const std::array<AlignedVector<Number>, 2> &shape_data,
3611  const unsigned int dofs_per_component_on_cell,
3612  const unsigned int dofs_per_component_on_face)
3613  {
3614  if (face_direction == face_no / 2)
3615  {
3617  dim,
3618  fe_degree + 1,
3619  0,
3620  Number>
3621  evalf(shape_data[face_no % 2],
3624  n_points_1d,
3625  0);
3626 
3627  const unsigned int in_stride = do_evaluate ?
3628  dofs_per_component_on_cell :
3629  dofs_per_component_on_face;
3630  const unsigned int out_stride = do_evaluate ?
3631  dofs_per_component_on_face :
3632  dofs_per_component_on_cell;
3633 
3634  for (unsigned int c = 0; c < n_components; ++c)
3635  {
3636  if (flag & EvaluationFlags::hessians)
3637  evalf.template apply_face<face_direction,
3638  do_evaluate,
3639  add_into_output,
3640  2>(input, output);
3641  else if (flag & EvaluationFlags::gradients)
3642  evalf.template apply_face<face_direction,
3643  do_evaluate,
3644  add_into_output,
3645  1>(input, output);
3646  else
3647  evalf.template apply_face<face_direction,
3648  do_evaluate,
3649  add_into_output,
3650  0>(input, output);
3651  input += in_stride;
3652  output += out_stride;
3653  }
3654  }
3655  else if (face_direction < dim)
3656  {
3657  interpolate_generic<do_evaluate,
3658  add_into_output,
3659  std::min(face_direction + 1, dim - 1)>(
3660  n_components,
3661  input,
3662  output,
3663  flag,
3664  face_no,
3665  n_points_1d,
3666  shape_data,
3667  dofs_per_component_on_cell,
3668  dofs_per_component_on_face);
3669  }
3670  }
3671 
3672  template <typename EvalType>
3673  static EvalType
3676  const unsigned int face_no)
3677  {
3678  return EvalType(data.shape_data_on_face[face_no % 2],
3681  }
3682 
3683  template <bool do_evaluate,
3684  bool add_into_output,
3685  int face_direction = 0,
3686  int max_derivative = 0>
3687  static void
3689  const unsigned int n_components,
3690  const Number * input,
3691  Number * output,
3693  const unsigned int face_no,
3694  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info)
3695  {
3696  if (dim == 1)
3697  {
3698  // This should never happen since the FE_RaviartThomasNodal is not
3699  // defined for dim = 1. It prevents compiler warnings of infinite
3700  // recursion.
3701  Assert(false, ExcInternalError());
3702  return;
3703  }
3704 
3705  bool increase_max_der = false;
3706  if ((flag & EvaluationFlags::hessians && max_derivative < 2) ||
3707  (flag & EvaluationFlags::gradients && max_derivative < 1))
3708  increase_max_der = true;
3709 
3710  if (face_direction == face_no / 2 && !increase_max_der)
3711  {
3713  add_into_output,
3714  face_direction,
3715  max_derivative>(
3716  shape_info, face_no, input, output);
3717  }
3718  else if (face_direction == face_no / 2)
3719  {
3720  // Only increase max_derivative
3722  add_into_output,
3723  face_direction,
3724  std::min(max_derivative + 1, 2)>(
3725  n_components, input, output, flag, face_no, shape_info);
3726  }
3727  else if (face_direction < dim)
3728  {
3729  if (increase_max_der)
3730  {
3732  do_evaluate,
3733  add_into_output,
3734  std::min(face_direction + 1, dim - 1),
3735  std::min(max_derivative + 1, 2)>(
3736  n_components, input, output, flag, face_no, shape_info);
3737  }
3738  else
3739  {
3741  add_into_output,
3742  std::min(face_direction + 1,
3743  dim - 1),
3744  max_derivative>(
3745  n_components, input, output, flag, face_no, shape_info);
3746  }
3747  }
3748  }
3749 
3750  /* Help function for interpolate_generic_raviart_thomas */
3751  template <bool do_evaluate,
3752  bool add_into_output,
3753  int face_direction,
3754  int max_derivative>
3755  static inline void
3757  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
3758  const unsigned int face_no,
3759  const Number * input,
3760  Number * output)
3761  {
3762  // These types are evaluators in either normal or tangential direction
3763  // depending on the face direction, with different normal directions for
3764  // the different components.
3765  using Evalf0 = typename std::conditional<
3766  face_direction == 0,
3768  dim,
3769  (fe_degree == -1) ? 1 : fe_degree + 1,
3770  0,
3771  Number,
3772  0>,
3774  dim,
3775  (fe_degree == -1) ? 1 : fe_degree,
3776  0,
3777  Number,
3778  0>>::type;
3779  using Evalf1 = typename std::conditional<
3780  face_direction == 1,
3782  dim,
3783  (fe_degree == -1) ? 1 : fe_degree + 1,
3784  0,
3785  Number,
3786  1>,
3788  dim,
3789  (fe_degree == -1) ? 1 : fe_degree,
3790  0,
3791  Number,
3792  1>>::type;
3793  using Evalf2 = typename std::conditional<
3794  face_direction == 2,
3796  dim,
3797  (fe_degree == -1) ? 1 : fe_degree + 1,
3798  0,
3799  Number,
3800  2>,
3802  dim,
3803  (fe_degree == -1) ? 1 : fe_degree,
3804  0,
3805  Number,
3806  2>>::type;
3807 
3808  Evalf0 evalf0 =
3809  create_evaluator_tensor_product<Evalf0>((face_direction == 0) ?
3810  shape_info.data[0] :
3811  shape_info.data[1],
3812  face_no);
3813  Evalf1 evalf1 =
3814  create_evaluator_tensor_product<Evalf1>((face_direction == 1) ?
3815  shape_info.data[0] :
3816  shape_info.data[1],
3817  face_no);
3818  Evalf2 evalf2 =
3819  create_evaluator_tensor_product<Evalf2>((face_direction == 2) ?
3820  shape_info.data[0] :
3821  shape_info.data[1],
3822  face_no);
3823 
3824  const unsigned int dofs_per_component_on_cell =
3825  shape_info.dofs_per_component_on_cell;
3826  const unsigned int dofs_per_component_on_face =
3827  3 * shape_info.dofs_per_component_on_face;
3828 
3829  // NOTE! dofs_per_component_on_face is in the tangent direction,
3830  // i.e (fe.degree+1)*fe.degree. Normal faces are only
3831  // fe.degree*fe.degree
3832  const unsigned int in_stride =
3833  do_evaluate ? dofs_per_component_on_cell : dofs_per_component_on_face;
3834  const unsigned int out_stride =
3835  do_evaluate ? dofs_per_component_on_face : dofs_per_component_on_cell;
3836 
3837  const unsigned int in_stride_after_normal =
3838  do_evaluate ?
3839  dofs_per_component_on_cell :
3840  dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2);
3841  const unsigned int out_stride_after_normal =
3842  do_evaluate ?
3843  dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2) :
3844  dofs_per_component_on_cell;
3845 
3846  evalf0.template apply_face<face_direction,
3847  do_evaluate,
3848  add_into_output,
3849  max_derivative>(input, output);
3850  // stride to next component
3851  input += (face_direction == 0) ? in_stride_after_normal : in_stride;
3852  output += (face_direction == 0) ? out_stride_after_normal : out_stride;
3853 
3854  evalf1.template apply_face<face_direction,
3855  do_evaluate,
3856  add_into_output,
3857  max_derivative>(input, output);
3858 
3859  if (dim == 3)
3860  {
3861  // stride to next component
3862  input += (face_direction == 1) ? in_stride_after_normal : in_stride;
3863  output +=
3864  (face_direction == 1) ? out_stride_after_normal : out_stride;
3865 
3866  evalf2.template apply_face<face_direction,
3867  do_evaluate,
3868  add_into_output,
3869  max_derivative>(input, output);
3870  }
3871  }
3872  };
3873 
3874 
3875 
3876  // internal helper function for reading data; base version of different types
3877  template <typename VectorizedArrayType, typename Number2>
3878  void
3879  do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
3880  {
3881  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3882  dst[v] = src_ptr[v];
3883  }
3884 
3885 
3886 
3887  // internal helper function for reading data; specialized version where we
3888  // can use a dedicated load function
3889  template <typename Number, std::size_t width>
3890  void
3892  {
3893  dst.load(src_ptr);
3894  }
3895 
3896 
3897 
3898  // internal helper function for reading data; base version of different types
3899  template <typename VectorizedArrayType, typename Number2>
3900  void
3901  do_vectorized_gather(const Number2 * src_ptr,
3902  const unsigned int * indices,
3903  VectorizedArrayType &dst)
3904  {
3905  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3906  dst[v] = src_ptr[indices[v]];
3907  }
3908 
3909 
3910 
3911  // internal helper function for reading data; specialized version where we
3912  // can use a dedicated gather function
3913  template <typename Number, std::size_t width>
3914  void
3915  do_vectorized_gather(const Number * src_ptr,
3916  const unsigned int * indices,
3918  {
3919  dst.gather(src_ptr, indices);
3920  }
3921 
3922 
3923 
3924  // internal helper function for reading data; base version of different types
3925  template <typename VectorizedArrayType, typename Number2>
3926  void
3927  do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
3928  {
3929  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3930  dst_ptr[v] += src[v];
3931  }
3932 
3933 
3934 
3935  // internal helper function for reading data; specialized version where we
3936  // can use a dedicated load function
3937  template <typename Number, std::size_t width>
3938  void
3940  {
3942  tmp.load(dst_ptr);
3943  (tmp + src).store(dst_ptr);
3944  }
3945 
3946 
3947 
3948  // internal helper function for reading data; base version of different types
3949  template <typename VectorizedArrayType, typename Number2>
3950  void
3951  do_vectorized_scatter_add(const VectorizedArrayType src,
3952  const unsigned int * indices,
3953  Number2 * dst_ptr)
3954  {
3955  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3956  dst_ptr[indices[v]] += src[v];
3957  }
3958 
3959 
3960 
3961  // internal helper function for reading data; specialized version where we
3962  // can use a dedicated gather function
3963  template <typename Number, std::size_t width>
3964  void
3966  const unsigned int * indices,
3967  Number * dst_ptr)
3968  {
3969 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
3970  for (unsigned int v = 0; v < width; ++v)
3971  dst_ptr[indices[v]] += src[v];
3972 #else
3974  tmp.gather(dst_ptr, indices);
3975  (tmp + src).scatter(indices, dst_ptr);
3976 #endif
3977  }
3978 
3979 
3980 
3981  template <typename Number>
3982  void
3983  adjust_for_face_orientation(const unsigned int dim,
3984  const unsigned int n_components,
3986  const unsigned int *orientation,
3987  const bool integrate,
3988  const std::size_t n_q_points,
3989  Number * tmp_values,
3990  Number * values_quad,
3991  Number * gradients_quad,
3992  Number * hessians_quad)
3993  {
3994  for (unsigned int c = 0; c < n_components; ++c)
3995  {
3996  if (flag & EvaluationFlags::values)
3997  {
3998  if (integrate)
3999  for (unsigned int q = 0; q < n_q_points; ++q)
4000  tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
4001  else
4002  for (unsigned int q = 0; q < n_q_points; ++q)
4003  tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
4004  for (unsigned int q = 0; q < n_q_points; ++q)
4005  values_quad[c * n_q_points + q] = tmp_values[q];
4006  }
4007  if (flag & EvaluationFlags::gradients)
4008  for (unsigned int d = 0; d < dim; ++d)
4009  {
4010  if (integrate)
4011  for (unsigned int q = 0; q < n_q_points; ++q)
4012  tmp_values[q] =
4013  gradients_quad[(c * dim + d) * n_q_points + orientation[q]];
4014  else
4015  for (unsigned int q = 0; q < n_q_points; ++q)
4016  tmp_values[orientation[q]] =
4017  gradients_quad[(c * dim + d) * n_q_points + q];
4018  for (unsigned int q = 0; q < n_q_points; ++q)
4019  gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q];
4020  }
4021  if (flag & EvaluationFlags::hessians)
4022  {
4023  const unsigned int hdim = (dim * (dim + 1)) / 2;
4024  for (unsigned int d = 0; d < hdim; ++d)
4025  {
4026  if (integrate)
4027  for (unsigned int q = 0; q < n_q_points; ++q)
4028  tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
4029  orientation[q]];
4030  else
4031  for (unsigned int q = 0; q < n_q_points; ++q)
4032  tmp_values[orientation[q]] =
4033  hessians_quad[(c * hdim + d) * n_q_points + q];
4034  for (unsigned int q = 0; q < n_q_points; ++q)
4035  hessians_quad[(c * hdim + d) * n_q_points + q] =
4036  tmp_values[q];
4037  }
4038  }
4039  }
4040  }
4041 
4042 
4043 
4044  template <typename Number, typename VectorizedArrayType>
4045  void
4047  const unsigned int dim,
4048  const unsigned int n_components,
4049  const unsigned int v,
4051  const unsigned int * orientation,
4052  const bool integrate,
4053  const std::size_t n_q_points,
4054  Number * tmp_values,
4055  VectorizedArrayType * values_quad,
4056  VectorizedArrayType * gradients_quad = nullptr,
4057  VectorizedArrayType * hessians_quad = nullptr)
4058  {
4059  for (unsigned int c = 0; c < n_components; ++c)
4060  {
4061  if (flag & EvaluationFlags::values)
4062  {
4063  if (integrate)
4064  for (unsigned int q = 0; q < n_q_points; ++q)
4065  tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v];
4066  else
4067  for (unsigned int q = 0; q < n_q_points; ++q)
4068  tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v];
4069  for (unsigned int q = 0; q < n_q_points; ++q)
4070  values_quad[c * n_q_points + q][v] = tmp_values[q];
4071  }
4072  if (flag & EvaluationFlags::gradients)
4073  for (unsigned int d = 0; d < dim; ++d)
4074  {
4075  Assert(gradients_quad != nullptr, ExcInternalError());
4076  if (integrate)
4077  for (unsigned int q = 0; q < n_q_points; ++q)
4078  tmp_values[q] = gradients_quad[(c * dim + d) * n_q_points +
4079  orientation[q]][v];
4080  else
4081  for (unsigned int q = 0; q < n_q_points; ++q)
4082  tmp_values[orientation[q]] =
4083  gradients_quad[(c * dim + d) * n_q_points + q][v];
4084  for (unsigned int q = 0; q < n_q_points; ++q)
4085  gradients_quad[(c * dim + d) * n_q_points + q][v] =
4086  tmp_values[q];
4087  }
4088  if (flag & EvaluationFlags::hessians)
4089  {
4090  Assert(hessians_quad != nullptr, ExcInternalError());
4091  const unsigned int hdim = (dim * (dim + 1)) / 2;
4092  for (unsigned int d = 0; d < hdim; ++d)
4093  {
4094  if (integrate)
4095  for (unsigned int q = 0; q < n_q_points; ++q)
4096  tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
4097  orientation[q]][v];
4098  else
4099  for (unsigned int q = 0; q < n_q_points; ++q)
4100  tmp_values[orientation[q]] =
4101  hessians_quad[(c * hdim + d) * n_q_points + q][v];
4102  for (unsigned int q = 0; q < n_q_points; ++q)
4103  hessians_quad[(c * hdim + d) * n_q_points + q][v] =
4104  tmp_values[q];
4105  }
4106  }
4107  }
4108  }
4109 
4110 
4111 
4112  template <int dim, typename Number>
4114  {
4115  template <int fe_degree, int n_q_points_1d>
4116  static bool
4117  run(const unsigned int n_components,
4118  const EvaluationFlags::EvaluationFlags evaluation_flag,
4119  const Number * values_dofs,
4121  {
4122  const auto &shape_info = fe_eval.get_shape_info();
4123  const auto &shape_data = shape_info.data.front();
4124 
4125  if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
4126  {
4127  Assert((fe_eval.get_dof_access_index() ==
4129  fe_eval.is_interior_face() == false) == false,
4130  ExcNotImplemented());
4131 
4132  const unsigned int face_no = fe_eval.get_face_no();
4133  const unsigned int face_orientation = fe_eval.get_face_orientation();
4134  const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
4135  const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
4136 
4137  using Eval =
4139 
4140  if (evaluation_flag & EvaluationFlags::values)
4141  {
4142  const auto shape_values =
4143  &shape_data.shape_values_face(face_no, face_orientation, 0);
4144 
4145  auto values_quad_ptr = fe_eval.begin_values();
4146  auto values_dofs_actual_ptr = values_dofs;
4147 
4148  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
4149  for (unsigned int c = 0; c < n_components; ++c)
4150  {
4151  eval.template values<0, true, false>(values_dofs_actual_ptr,
4152  values_quad_ptr);
4153 
4154  values_quad_ptr += n_q_points;
4155  values_dofs_actual_ptr += n_dofs;
4156  }
4157  }
4158 
4159  if (evaluation_flag & EvaluationFlags::gradients)
4160  {
4161  auto gradients_quad_ptr = fe_eval.begin_gradients();
4162  auto values_dofs_actual_ptr = values_dofs;
4163 
4164  std::array<const Number *, dim> shape_gradients;
4165  for (unsigned int d = 0; d < dim; ++d)
4166  shape_gradients[d] = &shape_data.shape_gradients_face(
4167  face_no, face_orientation, d, 0);
4168 
4169  for (unsigned int c = 0; c < n_components; ++c)
4170  {
4171  for (unsigned int d = 0; d < dim; ++d)
4172  {
4173  Eval eval(nullptr,
4174  shape_gradients[d],
4175  nullptr,
4176  n_dofs,
4177  n_q_points);
4178 
4179  eval.template gradients<0, true, false>(
4180  values_dofs_actual_ptr, gradients_quad_ptr);
4181 
4182  gradients_quad_ptr += n_q_points;
4183  }
4184  values_dofs_actual_ptr += n_dofs;
4185  }
4186  }
4187 
4188  Assert(!(evaluation_flag & EvaluationFlags::hessians),
4189  ExcNotImplemented());
4190 
4191  return true;
4192  }
4193 
4194  const unsigned int dofs_per_face =
4195  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
4196  Utilities::pow(shape_data.fe_degree + 1, dim - 1);
4197 
4198  // Note: we always keep storage of values, 1st and 2nd derivatives in an
4199  // array, so reserve space for all three here
4200  Number *temp = fe_eval.get_scratch_data().begin();
4201  Number *scratch_data = temp + 3 * n_components * dofs_per_face;
4202 
4203  bool use_vectorization = true;
4204 
4205  if (fe_eval.get_dof_access_index() ==
4207  fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
4208  use_vectorization =
4209  fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
4210  std::all_of(fe_eval.get_cell_ids().begin() + 1,
4211  fe_eval.get_cell_ids().end(),
4212  [&](const auto &v) {
4213  return v == fe_eval.get_cell_ids()[0] ||
4214  v == numbers::invalid_unsigned_int;
4215  });
4216 
4217  if (use_vectorization == false)
4218  {
4219  for (unsigned int v = 0; v < Number::size(); ++v)
4220  {
4221  // the loop breaks once an invalid_unsigned_int is hit for
4222  // all cases except the exterior faces in the ECL loop (where
4223  // some faces might be at the boundaries but others not)
4224  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4225  {
4226  for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4227  ++i)
4228  temp[i][v] = 0;
4229  continue;
4230  }
4231 
4233  template interpolate<true, false>(n_components,
4234  evaluation_flag,
4235  shape_info,
4236  values_dofs,
4237  scratch_data,
4238  fe_eval.get_face_no(v));
4239 
4240  for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4241  ++i)
4242  temp[i][v] = scratch_data[i][v];
4243  }
4244  }
4245  else
4247  template interpolate<true, false>(n_components,
4248  evaluation_flag,
4249  shape_info,
4250  values_dofs,
4251  temp,
4252  fe_eval.get_face_no());
4253 
4254  const unsigned int subface_index = fe_eval.get_subface_index();
4255  constexpr unsigned int n_q_points_1d_actual =
4256  fe_degree > -1 ? n_q_points_1d : 0;
4257 
4258  if (fe_degree >= 1 &&
4259  shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
4260  {
4262  (fe_degree == -1) ? 1 : fe_degree,
4263  (n_q_points_1d < 1) ? 1 :
4264  n_q_points_1d,
4265  Number>::
4266  template evaluate_or_integrate_in_face<false>(
4267  evaluation_flag,
4268  temp,
4269  fe_eval,
4270  scratch_data,
4271  subface_index,
4272  fe_eval.get_face_no());
4273  }
4274  else if (fe_degree > -1 &&
4275  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
4276  shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
4277  FEFaceEvaluationImpl<true,
4278  dim,
4279  fe_degree,
4280  n_q_points_1d_actual,
4281  Number>::evaluate_in_face(n_components,
4282  evaluation_flag,
4283  shape_data,
4284  temp,
4285  fe_eval.begin_values(),
4286  fe_eval
4287  .begin_gradients(),
4288  fe_eval.begin_hessians(),
4289  scratch_data,
4290  subface_index);
4291  else
4292  FEFaceEvaluationImpl<false,
4293  dim,
4294  fe_degree,
4295  n_q_points_1d_actual,
4296  Number>::evaluate_in_face(n_components,
4297  evaluation_flag,
4298  shape_data,
4299  temp,
4300  fe_eval.begin_values(),
4301  fe_eval
4302  .begin_gradients(),
4303  fe_eval.begin_hessians(),
4304  scratch_data,
4305  subface_index);
4306 
4307  if (use_vectorization == false)
4308  {
4309  for (unsigned int v = 0; v < Number::size(); ++v)
4310  {
4311  // the loop breaks once an invalid_unsigned_int is hit for
4312  // all cases except the exterior faces in the ECL loop (where
4313  // some faces might be at the boundaries but others not)
4314  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4315  continue;
4316 
4317  if (fe_eval.get_face_orientation(v) != 0)
4319  dim,
4320  n_components,
4321  v,
4322  evaluation_flag,
4324  fe_eval.get_face_orientation(v), 0),
4325  false,
4326  shape_info.n_q_points_face,
4327  &temp[0][0],
4328  fe_eval.begin_values(),
4329  fe_eval.begin_gradients(),
4330  fe_eval.begin_hessians());
4331  }
4332  }
4333  else if (fe_eval.get_face_orientation() != 0)
4335  dim,
4336  n_components,
4337  evaluation_flag,
4339  fe_eval.get_face_orientation(), 0),
4340  false,
4341  shape_info.n_q_points_face,
4342  temp,
4343  fe_eval.begin_values(),
4344  fe_eval.begin_gradients(),
4345  fe_eval.begin_hessians());
4346 
4347  return false;
4348  }
4349  };
4350 
4351 
4352 
4353  template <int dim, typename Number>
4355  {
4356  template <int fe_degree, int n_q_points_1d>
4357  static bool
4358  run(const unsigned int n_components,
4359  const EvaluationFlags::EvaluationFlags integration_flag,
4360  Number * values_dofs,
4362  {
4363  const auto &shape_info = fe_eval.get_shape_info();
4364  const auto &shape_data = shape_info.data.front();
4365 
4366  if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
4367  {
4368  Assert((fe_eval.get_dof_access_index() ==
4370  fe_eval.is_interior_face() == false) == false,
4371  ExcNotImplemented());
4372 
4373  const unsigned int face_no = fe_eval.get_face_no();
4374  const unsigned int face_orientation = fe_eval.get_face_orientation();
4375  const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
4376  const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
4377 
4378  using Eval =
4380 
4381  if (integration_flag & EvaluationFlags::values)
4382  {
4383  const auto shape_values =
4384  &shape_data.shape_values_face(face_no, face_orientation, 0);
4385 
4386  auto values_quad_ptr = fe_eval.begin_values();
4387  auto values_dofs_actual_ptr = values_dofs;
4388 
4389  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
4390  for (unsigned int c = 0; c < n_components; ++c)
4391  {
4392  eval.template values<0, false, false>(values_quad_ptr,
4393  values_dofs_actual_ptr);
4394 
4395  values_quad_ptr += n_q_points;
4396  values_dofs_actual_ptr += n_dofs;
4397  }
4398  }
4399 
4400  if (integration_flag & EvaluationFlags::gradients)
4401  {
4402  auto gradients_quad_ptr = fe_eval.begin_gradients();
4403  auto values_dofs_actual_ptr = values_dofs;
4404 
4405  std::array<const Number *, dim> shape_gradients;
4406  for (unsigned int d = 0; d < dim; ++d)
4407  shape_gradients[d] = &shape_data.shape_gradients_face(
4408  face_no, face_orientation, d, 0);
4409 
4410  for (unsigned int c = 0; c < n_components; ++c)
4411  {
4412  for (unsigned int d = 0; d < dim; ++d)
4413  {
4414  Eval eval(nullptr,
4415  shape_gradients[d],
4416  nullptr,
4417  n_dofs,
4418  n_q_points);
4419 
4420  if (!(integration_flag & EvaluationFlags::values) &&
4421  d == 0)
4422  eval.template gradients<0, false, false>(
4423  gradients_quad_ptr, values_dofs_actual_ptr);
4424  else
4425  eval.template gradients<0, false, true>(
4426  gradients_quad_ptr, values_dofs_actual_ptr);
4427 
4428  gradients_quad_ptr += n_q_points;
4429  }
4430  values_dofs_actual_ptr += n_dofs;
4431  }
4432  }
4433 
4434  Assert(!(integration_flag & EvaluationFlags::hessians),
4435  ExcNotImplemented());
4436 
4437  return true;
4438  }
4439 
4440  const unsigned int dofs_per_face =
4441  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
4442  Utilities::pow(shape_data.fe_degree + 1, dim - 1);
4443 
4444  Number *temp = fe_eval.get_scratch_data().begin();
4445  Number *scratch_data = temp + 3 * n_components * dofs_per_face;
4446 
4447  bool use_vectorization = true;
4448 
4449  if (fe_eval.get_dof_access_index() ==
4451  fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
4452  use_vectorization =
4453  fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
4454  std::all_of(fe_eval.get_cell_ids().begin() + 1,
4455  fe_eval.get_cell_ids().end(),
4456  [&](const auto &v) {
4457  return v == fe_eval.get_cell_ids()[0] ||
4458  v == numbers::invalid_unsigned_int;
4459  });
4460 
4461  if (use_vectorization == false)
4462  {
4463  for (unsigned int v = 0; v < Number::size(); ++v)
4464  {
4465  // the loop breaks once an invalid_unsigned_int is hit for
4466  // all cases except the exterior faces in the ECL loop (where
4467  // some faces might be at the boundaries but others not)
4468  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4469  continue;
4470 
4471  if (fe_eval.get_face_orientation(v) != 0)
4473  dim,
4474  n_components,
4475  v,
4476  integration_flag,
4478  fe_eval.get_face_orientation(v), 0),
4479  true,
4480  shape_info.n_q_points_face,
4481  &temp[0][0],
4482  fe_eval.begin_values(),
4483  fe_eval.begin_gradients(),
4484  fe_eval.begin_hessians());
4485  }
4486  }
4487  else if (fe_eval.get_face_orientation() != 0)
4489  dim,
4490  n_components,
4491  integration_flag,
4493  fe_eval.get_face_orientation(), 0),
4494  true,
4495  shape_info.n_q_points_face,
4496  temp,
4497  fe_eval.begin_values(),
4498  fe_eval.begin_gradients(),
4499  fe_eval.begin_hessians());
4500 
4501  const unsigned int n_q_points_1d_actual =
4502  fe_degree > -1 ? n_q_points_1d : 0;
4503  const unsigned int subface_index = fe_eval.get_subface_index();
4504 
4505  if (fe_degree >= 1 &&
4506  shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
4507  {
4509  (fe_degree == -1) ? 1 : fe_degree,
4510  (n_q_points_1d < 1) ? 1 :
4511  n_q_points_1d,
4512  Number>::
4513  template evaluate_or_integrate_in_face<true>(integration_flag,
4514  temp,
4515  fe_eval,
4516  scratch_data,
4517  subface_index,
4518  fe_eval.get_face_no());
4519  }
4520  else if (fe_degree > -1 &&
4521  fe_eval.get_subface_index() >=
4523  shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
4525  true,
4526  dim,
4527  fe_degree,
4528  n_q_points_1d_actual,
4529  Number>::integrate_in_face(n_components,
4530  integration_flag,
4531  shape_data,
4532  temp,
4533  fe_eval.begin_values(),
4534  fe_eval.begin_gradients(),
4535  fe_eval.begin_hessians(),
4536  scratch_data,
4537  subface_index);
4538  else
4540  false,
4541  dim,
4542  fe_degree,
4543  n_q_points_1d_actual,
4544  Number>::integrate_in_face(n_components,
4545  integration_flag,
4546  shape_data,
4547  temp,
4548  fe_eval.begin_values(),
4549  fe_eval.begin_gradients(),
4550  fe_eval.begin_hessians(),
4551  scratch_data,
4552  subface_index);
4553 
4554  if (use_vectorization == false)
4555  {
4556  for (unsigned int v = 0; v < Number::size(); ++v)
4557  {
4558  // the loop breaks once an invalid_unsigned_int is hit for
4559  // all cases except the exterior faces in the ECL loop (where
4560  // some faces might be at the boundaries but others not)
4561  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4562  continue;
4563 
4565  template interpolate<false, false>(n_components,
4566  integration_flag,
4567  shape_info,
4568  values_dofs,
4569  scratch_data,
4570  fe_eval.get_face_no(v));
4571 
4572  for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4573  ++i)
4574  temp[i][v] = scratch_data[i][v];
4575  }
4576  }
4577  else
4579  template interpolate<false, false>(n_components,
4580  integration_flag,
4581  shape_info,
4582  temp,
4583  values_dofs,
4584  fe_eval.get_face_no());
4585  return false;
4586  }
4587  };
4588 
4589 
4590 
4591  template <int n_face_orientations,
4592  typename Processor,
4593  typename EvaluationData,
4594  const bool check_face_orientations = false>
4595  void
4597  Processor & proc,
4598  const unsigned int n_components,
4599  const EvaluationFlags::EvaluationFlags evaluation_flag,
4600  typename Processor::Number2_ * global_vector_ptr,
4601  const std::vector<ArrayView<const typename Processor::Number2_>> *sm_ptr,
4602  const EvaluationData & fe_eval,
4603  typename Processor::VectorizedArrayType_ * temp1)
4604  {
4605  constexpr int dim = Processor::dim_;
4606  constexpr int fe_degree = Processor::fe_degree_;
4607  using VectorizedArrayType = typename Processor::VectorizedArrayType_;
4608  constexpr int n_lanes = VectorizedArrayType::size();
4609 
4610  using Number = typename Processor::Number_;
4611  using Number2_ = typename Processor::Number2_;
4612 
4613  const auto & shape_data = fe_eval.get_shape_info().data.front();
4614  constexpr bool integrate = Processor::do_integrate;
4615  const unsigned int face_no = fe_eval.get_face_no();
4616  const auto & dof_info = fe_eval.get_dof_info();
4617  const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
4618  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index =
4619  fe_eval.get_dof_access_index();
4620  AssertIndexRange(cell,
4621  dof_info.index_storage_variants[dof_access_index].size());
4622  constexpr unsigned int dofs_per_face =
4623  Utilities::pow(fe_degree + 1, dim - 1);
4624  const unsigned int subface_index = fe_eval.get_subface_index();
4625 
4626  const unsigned int n_filled_lanes =
4627  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
4628 
4629  bool all_faces_are_same = n_filled_lanes == n_lanes;
4630  if (n_face_orientations == n_lanes)
4631  for (unsigned int v = 1; v < n_lanes; ++v)
4632  if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) ||
4633  fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0))
4634  {
4635  all_faces_are_same = false;
4636  break;
4637  }
4638 
4639  // check for re-orientation ...
4640  std::array<const unsigned int *, n_face_orientations> orientation = {};
4641 
4642  if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same &&
4643  fe_eval.is_interior_face() == 0)
4644  for (unsigned int v = 0; v < n_lanes; ++v)
4645  {
4646  // the loop breaks once an invalid_unsigned_int is hit for
4647  // all cases except the exterior faces in the ECL loop (where
4648  // some faces might be at the boundaries but others not)
4649  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4650  continue;
4651 
4652  if (shape_data.nodal_at_cell_boundaries &&
4653  fe_eval.get_face_orientation(v) != 0)
4654  {
4655  // ... and in case we detect a re-orientation, go to the other
4656  // version of this function that actually allows for this
4657  if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
4658  check_face_orientations == false)
4659  {
4660  fe_face_evaluation_process_and_io<n_face_orientations,
4661  Processor,
4662  EvaluationData,
4663  true>(proc,
4664  n_components,
4665  evaluation_flag,
4666  global_vector_ptr,
4667  sm_ptr,
4668  fe_eval,
4669  temp1);
4670  return;
4671  }
4672  orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
4673  fe_eval.get_face_orientation(v), 0);
4674  }
4675  }
4676  else if (dim == 3 && fe_eval.get_face_orientation() != 0)
4677  {
4678  // go to the other version of this function
4679  if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
4680  check_face_orientations == false)
4681  {
4682  fe_face_evaluation_process_and_io<n_face_orientations,
4683  Processor,
4684  EvaluationData,
4685  true>(proc,
4686  n_components,
4687  evaluation_flag,
4688  global_vector_ptr,
4689  sm_ptr,
4690  fe_eval,
4691  temp1);
4692  return;
4693  }
4694  for (unsigned int v = 0; v < n_face_orientations; ++v)
4695  orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
4696  fe_eval.get_face_orientation(), 0);
4697  }
4698 
4699  // we know that the gradient weights for the Hermite case on the
4700  // right (side==1) are the negative from the value at the left
4701  // (side==0), so we only read out one of them.
4702  VectorizedArrayType grad_weight =
4703  shape_data
4704  .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) :
4705  (1 + face_no % 2))];
4706 
4707  // face_to_cell_index_hermite
4708  std::array<const unsigned int *, n_face_orientations> index_array_hermite =
4709  {};
4710  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4711  {
4712  if (n_face_orientations == 1)
4713  index_array_hermite[0] =
4714  &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0);
4715  else
4716  {
4717  for (unsigned int v = 0; v < n_lanes; ++v)
4718  {
4719  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4720  continue;
4721 
4722  const auto face_no = fe_eval.get_face_no(v);
4723 
4724  grad_weight[v] =
4725  shape_data.shape_data_on_face[0][fe_degree +
4726  (integrate ?
4727  (2 - (face_no % 2)) :
4728  (1 + (face_no % 2)))][0];
4729 
4730  index_array_hermite[v] =
4731  &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no,
4732  0);
4733  }
4734  }
4735  }
4736 
4737  // face_to_cell_index_nodal
4738  std::array<const unsigned int *, n_face_orientations> index_array_nodal =
4739  {};
4740  if (shape_data.nodal_at_cell_boundaries == true)
4741  {
4742  if (n_face_orientations == 1)
4743  index_array_nodal[0] =
4744  &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0);
4745  else
4746  {
4747  for (unsigned int v = 0; v < n_lanes; ++v)
4748  {
4749  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4750  continue;
4751 
4752  const auto face_no = fe_eval.get_face_no(v);
4753 
4754  index_array_nodal[v] =
4755  &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no,
4756  0);
4757  }
4758  }
4759  }
4760 
4761 
4762  const auto reorientate = [&](const unsigned int v, const unsigned int i) {
4763  return (!check_face_orientations || orientation[v] == nullptr) ?
4764  i :
4765  orientation[v][i];
4766  };
4767 
4768  const unsigned int cell_index =
4769  dof_access_index == MatrixFreeFunctions::DoFInfo::dof_access_cell ?
4770  fe_eval.get_cell_ids()[0] :
4771  cell * n_lanes;
4772  const unsigned int *dof_indices =
4773  &dof_info.dof_indices_contiguous[dof_access_index][cell_index];
4774 
4775  for (unsigned int comp = 0; comp < n_components; ++comp)
4776  {
4777  const std::size_t index_offset =
4778  dof_info.component_dof_indices_offset
4779  [fe_eval.get_active_fe_index()]
4780  [fe_eval.get_first_selected_component()] +
4781  comp * Utilities::pow(fe_degree + 1, dim);
4782 
4783  // case 1: contiguous and interleaved indices
4784  if (n_face_orientations == 1 &&
4785  dof_info.index_storage_variants[dof_access_index][cell] ==
4787  interleaved_contiguous)
4788  {
4790  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
4791  n_lanes);
4792  Number2_ *vector_ptr =
4793  global_vector_ptr + dof_indices[0] + index_offset * n_lanes;
4794 
4795  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4796  {
4797  for (unsigned int i = 0; i < dofs_per_face; ++i)
4798  {
4799  Assert(n_face_orientations == 1, ExcNotImplemented());
4800 
4801  const unsigned int ind1 = index_array_hermite[0][2 * i];
4802  const unsigned int ind2 = index_array_hermite[0][2 * i + 1];
4803  const unsigned int i_ = reorientate(0, i);
4804  proc.hermite_grad_vectorized(temp1[i_],
4805  temp1[i_ + dofs_per_face],
4806  vector_ptr + ind1 * n_lanes,
4807  vector_ptr + ind2 * n_lanes,
4808  grad_weight);
4809  }
4810  }
4811  else
4812  {
4813  for (unsigned int i = 0; i < dofs_per_face; ++i)
4814  {
4815  Assert(n_face_orientations == 1, ExcNotImplemented());
4816 
4817  const unsigned int i_ = reorientate(0, i);
4818  const unsigned int ind = index_array_nodal[0][i];
4819  proc.value_vectorized(temp1[i_],
4820  vector_ptr + ind * n_lanes);
4821  }
4822  }
4823  }
4824 
4825  // case 2: contiguous and interleaved indices with fixed stride
4826  else if (n_face_orientations == 1 &&
4827  dof_info.index_storage_variants[dof_access_index][cell] ==
4829  interleaved_contiguous_strided)
4830  {
4832  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
4833  n_lanes);
4834  Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes;
4835  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4836  {
4837  for (unsigned int i = 0; i < dofs_per_face; ++i)
4838  {
4839  Assert(n_face_orientations == 1, ExcNotImplemented());
4840 
4841  const unsigned int i_ = reorientate(0, i);
4842  const unsigned int ind1 =
4843  index_array_hermite[0][2 * i] * n_lanes;
4844  const unsigned int ind2 =
4845  index_array_hermite[0][2 * i + 1] * n_lanes;
4846  proc.hermite_grad_vectorized_indexed(
4847  temp1[i_],
4848  temp1[i_ + dofs_per_face],
4849  vector_ptr + ind1,
4850  vector_ptr + ind2,
4851  grad_weight,
4852  dof_indices,
4853  dof_indices);
4854  }
4855  }
4856  else
4857  {
4858  for (unsigned int i = 0; i < dofs_per_face; ++i)
4859  {
4860  Assert(n_face_orientations == 1, ExcNotImplemented());
4861 
4862  const unsigned int i_ = reorientate(0, i);
4863  const unsigned int ind = index_array_nodal[0][i] * n_lanes;
4864  proc.value_vectorized_indexed(temp1[i_],
4865  vector_ptr + ind,
4866  dof_indices);
4867  }
4868  }
4869  }
4870 
4871  // case 3: contiguous and interleaved indices with mixed stride
4872  else if (n_face_orientations == 1 &&
4873  dof_info.index_storage_variants[dof_access_index][cell] ==
4875  interleaved_contiguous_mixed_strides)
4876  {
4877  const unsigned int *strides =
4878  &dof_info.dof_indices_interleave_strides[dof_access_index]
4879  [cell * n_lanes];
4880  unsigned int indices[n_lanes];
4881  for (unsigned int v = 0; v < n_lanes; ++v)
4882  indices[v] = dof_indices[v] + index_offset * strides[v];
4883  const unsigned int n_filled_lanes =
4884  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
4885 
4886  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4887  {
4888  if (n_filled_lanes == n_lanes)
4889  for (unsigned int i = 0; i < dofs_per_face; ++i)
4890  {
4891  Assert(n_face_orientations == 1, ExcNotImplemented());
4892 
4893  const unsigned int i_ = reorientate(0, i);
4894  unsigned int ind1[n_lanes];
4896  for (unsigned int v = 0; v < n_lanes; ++v)
4897  ind1[v] = indices[v] +
4898  index_array_hermite[0][2 * i] * strides[v];
4899  unsigned int ind2[n_lanes];
4901  for (unsigned int v = 0; v < n_lanes; ++v)
4902  ind2[v] =
4903  indices[v] +
4904  // TODO
4905  index_array_hermite[0][2 * i + 1] * strides[v];
4906  proc.hermite_grad_vectorized_indexed(
4907  temp1[i_],
4908  temp1[i_ + dofs_per_face],
4909  global_vector_ptr,
4910  global_vector_ptr,
4911  grad_weight,
4912  ind1,
4913  ind2);
4914  }
4915  else
4916  {
4917  if (integrate == false)
4918  for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
4919  temp1[i] = VectorizedArrayType();
4920 
4921  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4922  for (unsigned int i = 0; i < dofs_per_face; ++i)
4923  {
4924  const unsigned int i_ =
4925  reorientate(n_face_orientations == 1 ? 0 : v, i);
4926  proc.hermite_grad(
4927  temp1[i_][v],
4928  temp1[i_ + dofs_per_face][v],
4929  global_vector_ptr
4930  [indices[v] +
4931  index_array_hermite
4932  [n_face_orientations == 1 ? 0 : v][2 * i] *
4933  strides[v]],
4934  global_vector_ptr
4935  [indices[v] +
4936  index_array_hermite[n_face_orientations == 1 ?
4937  0 :
4938  v][2 * i + 1] *
4939  strides[v]],
4940  grad_weight[n_face_orientations == 1 ? 0 : v]);
4941  }
4942  }
4943  }
4944  else
4945  {
4946  if (n_filled_lanes == n_lanes)
4947  for (unsigned int i = 0; i < dofs_per_face; ++i)
4948  {
4949  Assert(n_face_orientations == 1, ExcInternalError());
4950  unsigned int ind[n_lanes];
4952  for (unsigned int v = 0; v < n_lanes; ++v)
4953  ind[v] =
4954  indices[v] + index_array_nodal[0][i] * strides[v];
4955  const unsigned int i_ = reorientate(0, i);
4956  proc.value_vectorized_indexed(temp1[i_],
4957  global_vector_ptr,
4958  ind);
4959  }
4960  else
4961  {
4962  if (integrate == false)
4963  for (unsigned int i = 0; i < dofs_per_face; ++i)
4964  temp1[i] = VectorizedArrayType();
4965 
4966  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4967  for (unsigned int i = 0; i < dofs_per_face; ++i)
4968  proc.value(
4969  temp1[reorientate(n_face_orientations == 1 ? 0 : v,
4970  i)][v],
4971  global_vector_ptr
4972  [indices[v] +
4973  index_array_nodal[n_face_orientations == 1 ? 0 : v]
4974  [i] *
4975  strides[v]]);
4976  }
4977  }
4978  }
4979 
4980  // case 4: contiguous indices without interleaving
4981  else if (n_face_orientations > 1 ||
4982  dof_info.index_storage_variants[dof_access_index][cell] ==
4984  contiguous)
4985  {
4986  Number2_ *vector_ptr = global_vector_ptr + index_offset;
4987 
4988  const bool vectorization_possible =
4989  all_faces_are_same && (sm_ptr == nullptr);
4990 
4991  std::array<Number2_ *, n_lanes> vector_ptrs;
4992  std::array<unsigned int, n_lanes> reordered_indices;
4993 
4994  if (vectorization_possible == false)
4995  {
4996  vector_ptrs = {};
4997  if (n_face_orientations == 1)
4998  {
4999  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5000  if (sm_ptr == nullptr)
5001  {
5002  vector_ptrs[v] = vector_ptr + dof_indices[v];
5003  }
5004  else
5005  {
5006  const auto &temp =
5007  dof_info
5008  .dof_indices_contiguous_sm[dof_access_index]
5009  [cell * n_lanes + v];
5010  vector_ptrs[v] = const_cast<Number2_ *>(
5011  sm_ptr->operator[](temp.first).data() +
5012  temp.second + index_offset);
5013  }
5014  }
5015  else if (n_face_orientations == n_lanes)
5016  {
5017  const auto &cells = fe_eval.get_cell_ids();
5018  for (unsigned int v = 0; v < n_lanes; ++v)
5019  if (cells[v] != numbers::invalid_unsigned_int)
5020  {
5021  if (sm_ptr == nullptr)
5022  {
5023  vector_ptrs[v] =
5024  vector_ptr +
5025  dof_info
5026  .dof_indices_contiguous[dof_access_index]
5027  [cells[v]];
5028  }
5029  else
5030  {
5031  const auto &temp =
5032  dof_info
5033  .dof_indices_contiguous_sm[dof_access_index]
5034  [cells[v]];
5035  vector_ptrs[v] = const_cast<Number2_ *>(
5036  sm_ptr->operator[](temp.first).data() +
5037  temp.second + index_offset);
5038  }
5039  }
5040  }
5041  else
5042  {
5043  Assert(false, ExcNotImplemented());
5044  }
5045  }
5046  else if (n_face_orientations == n_lanes)
5047  {
5048  for (unsigned int v = 0; v < n_lanes; ++v)
5049  reordered_indices[v] =
5050  dof_info.dof_indices_contiguous[dof_access_index]
5051  [fe_eval.get_cell_ids()[v]];
5052  dof_indices = reordered_indices.data();
5053  }
5054 
5055  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
5056  {
5057  if (vectorization_possible)
5058  for (unsigned int i = 0; i < dofs_per_face; ++i)
5059  {
5060  const unsigned int ind1 = index_array_hermite[0][2 * i];
5061  const unsigned int ind2 =
5062  index_array_hermite[0][2 * i + 1];
5063  const unsigned int i_ = reorientate(0, i);
5064 
5065  proc.hermite_grad_vectorized_indexed(
5066  temp1[i_],
5067  temp1[i_ + dofs_per_face],
5068  vector_ptr + ind1,
5069  vector_ptr + ind2,
5070  grad_weight,
5071  dof_indices,
5072  dof_indices);
5073  }
5074  else if (n_face_orientations == 1)
5075  for (unsigned int i = 0; i < dofs_per_face; ++i)
5076  {
5077  const unsigned int ind1 = index_array_hermite[0][2 * i];
5078  const unsigned int ind2 =
5079  index_array_hermite[0][2 * i + 1];
5080  const unsigned int i_ = reorientate(0, i);
5081 
5082  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5083  proc.hermite_grad(temp1[i_][v],
5084  temp1[i_ + dofs_per_face][v],
5085  vector_ptrs[v][ind1],
5086  vector_ptrs[v][ind2],
5087  grad_weight[v]);
5088 
5089  if (integrate == false)
5090  for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
5091  {
5092  temp1[i][v] = 0.0;
5093  temp1[i + dofs_per_face][v] = 0.0;
5094  }
5095  }
5096  else
5097  {
5098  if (integrate == false && n_filled_lanes < n_lanes)
5099  for (unsigned int i = 0; i < dofs_per_face; ++i)
5100  temp1[i] = temp1[i + dofs_per_face] = Number();
5101 
5102  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5103  for (unsigned int i = 0; i < dofs_per_face; ++i)
5104  proc.hermite_grad(
5105  temp1[reorientate(v, i)][v],
5106  temp1[reorientate(v, i) + dofs_per_face][v],
5107  vector_ptrs[v][index_array_hermite[v][2 * i]],
5108  vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
5109  grad_weight[v]);
5110  }
5111  }
5112  else
5113  {
5114  if (vectorization_possible)
5115  for (unsigned int i = 0; i < dofs_per_face; ++i)
5116  {
5117  const unsigned int ind = index_array_nodal[0][i];
5118  const unsigned int i_ = reorientate(0, i);
5119 
5120  proc.value_vectorized_indexed(temp1[i_],
5121  vector_ptr + ind,
5122  dof_indices);
5123  }
5124  else if (n_face_orientations == 1)
5125  for (unsigned int i = 0; i < dofs_per_face; ++i)
5126  {
5127  const unsigned int ind = index_array_nodal[0][i];
5128  const unsigned int i_ = reorientate(0, i);
5129 
5130  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5131  proc.value(temp1[i_][v], vector_ptrs[v][ind]);
5132 
5133  if (integrate == false)
5134  for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
5135  temp1[i_][v] = 0.0;
5136  }
5137  else
5138  {
5139  if (integrate == false && n_filled_lanes < n_lanes)
5140  for (unsigned int i = 0; i < dofs_per_face; ++i)
5141  temp1[i] = Number();
5142 
5143  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5144  for (unsigned int i = 0; i < dofs_per_face; ++i)
5145  proc.value(temp1[reorientate(v, i)][v],
5146  vector_ptrs[v][index_array_nodal[v][i]]);
5147  }
5148  }
5149  }
5150  else
5151  {
5152  // We should not end up here, this should be caught by
5153  // FEFaceEvaluationImplGatherEvaluateSelector::supports()
5154  Assert(false, ExcInternalError());
5155  }
5156  temp1 += 3 * dofs_per_face;
5157  }
5158  }
5159 
5160 
5161 
5162  template <int dim, typename Number2, typename VectorizedArrayType>
5164  {
5165  using Number = typename VectorizedArrayType::value_type;
5166 
5167  template <int fe_degree, int n_q_points_1d>
5168  static bool
5169  run(const unsigned int n_components,
5170  const EvaluationFlags::EvaluationFlags evaluation_flag,
5171  const Number2 * src_ptr,
5172  const std::vector<ArrayView<const Number2>> * sm_ptr,
5174  {
5175  Assert(fe_degree > -1, ExcInternalError());
5176  Assert(fe_eval.get_shape_info().element_type <=
5178  ExcInternalError());
5179 
5180  const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
5181 
5182  VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
5183  VectorizedArrayType *scratch_data =
5184  temp + 3 * n_components * dofs_per_face;
5185 
5187 
5188  if (fe_eval.get_dof_access_index() ==
5190  fe_eval.is_interior_face() == false)
5191  fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
5192  p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
5193  else
5194  fe_face_evaluation_process_and_io<1>(
5195  p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
5196 
5197  const unsigned int subface_index = fe_eval.get_subface_index();
5198 
5199  if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
5200  FEFaceEvaluationImpl<true,
5201  dim,
5202  fe_degree,
5203  n_q_points_1d,
5204  VectorizedArrayType>::
5205  evaluate_in_face(n_components,
5206  evaluation_flag,
5207  fe_eval.get_shape_info().data.front(),
5208  temp,
5209  fe_eval.begin_values(),
5210  fe_eval.begin_gradients(),
5211  fe_eval.begin_hessians(),
5212  scratch_data,
5213  subface_index);
5214  else
5215  FEFaceEvaluationImpl<false,
5216  dim,
5217  fe_degree,
5218  n_q_points_1d,
5219  VectorizedArrayType>::
5220  evaluate_in_face(n_components,
5221  evaluation_flag,
5222  fe_eval.get_shape_info().data.front(),
5223  temp,
5224  fe_eval.begin_values(),
5225  fe_eval.begin_gradients(),
5226  fe_eval.begin_hessians(),
5227  scratch_data,
5228  subface_index);
5229 
5230  // re-orientation for cases not possible with above algorithm
5231  if (subface_index < GeometryInfo<dim>::max_children_per_cell)
5232  {
5233  if (fe_eval.get_dof_access_index() ==
5235  fe_eval.is_interior_face() == false)
5236  {
5237  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5238  {
5239  // the loop breaks once an invalid_unsigned_int is hit for
5240  // all cases except the exterior faces in the ECL loop (where
5241  // some faces might be at the boundaries but others not)
5242  if (fe_eval.get_cell_ids()[v] ==
5244  continue;
5245 
5246  if (fe_eval.get_face_orientation(v) != 0)
5248  dim,
5249  n_components,
5250  v,
5251  evaluation_flag,
5253  fe_eval.get_face_orientation(v), 0),
5254  false,
5255  Utilities::pow(n_q_points_1d, dim - 1),
5256  &temp[0][0],
5257  fe_eval.begin_values(),
5258  fe_eval.begin_gradients(),
5259  fe_eval.begin_hessians());
5260  }
5261  }
5262  else if (fe_eval.get_face_orientation() != 0)
5264  dim,
5265  n_components,
5266  evaluation_flag,
5268  fe_eval.get_face_orientation(), 0),
5269  false,
5270  Utilities::pow(n_q_points_1d, dim - 1),
5271  temp,
5272  fe_eval.begin_values(),
5273  fe_eval.begin_gradients(),
5274  fe_eval.begin_hessians());
5275  }
5276 
5277  return false;
5278  }
5279 
5280  static bool
5282  const EvaluationFlags::EvaluationFlags evaluation_flag,
5284  const Number2 * vector_ptr,
5286  {
5287  const unsigned int fe_degree = shape_info.data.front().fe_degree;
5288  if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries ||
5289  (evaluation_flag & EvaluationFlags::gradients &&
5290  (fe_degree < 2 ||
5291  shape_info.data.front().element_type !=
5293  (evaluation_flag & EvaluationFlags::hessians) ||
5294  vector_ptr == nullptr ||
5295  shape_info.data.front().element_type >
5297  storage <
5299  return false;
5300  else
5301  return true;
5302  }
5303 
5304  private:
5305  template <int fe_degree>
5306  struct Processor
5307  {
5308  static const bool do_integrate = false;
5309  static const int dim_ = dim;
5310  static const int fe_degree_ = fe_degree;
5311  using VectorizedArrayType_ = VectorizedArrayType;
5312  using Number_ = Number;
5313  using Number2_ = const Number2;
5314 
5315  template <typename T0, typename T1, typename T2>
5316  void
5318  T0 & temp_2,
5319  const T1 src_ptr_1,
5320  const T1 src_ptr_2,
5321  const T2 &grad_weight)
5322  {
5323  do_vectorized_read(src_ptr_1, temp_1);
5324  do_vectorized_read(src_ptr_2, temp_2);
5325  temp_2 = grad_weight * (temp_1 - temp_2);
5326  }
5327 
5328  template <typename T1, typename T2>
5329  void
5330  value_vectorized(T1 &temp, const T2 src_ptr)
5331  {
5332  do_vectorized_read(src_ptr, temp);
5333  }
5334 
5335  template <typename T0, typename T1, typename T2, typename T3>
5336  void
5338  T0 & temp_2,
5339  const T1 src_ptr_1,
5340  const T1 src_ptr_2,
5341  const T2 &grad_weight,
5342  const T3 &indices_1,
5343  const T3 &indices_2)
5344  {
5345  do_vectorized_gather(src_ptr_1, indices_1, temp_1);
5346  do_vectorized_gather(src_ptr_2, indices_2, temp_2);
5347  temp_2 = grad_weight * (temp_1 - temp_2);
5348  }
5349 
5350  template <typename T0, typename T1, typename T2>
5351  void
5352  value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
5353  {
5354  do_vectorized_gather(src_ptr, indices, temp);
5355  }
5356 
5357  template <typename T0, typename T1, typename T2>
5358  void
5359  hermite_grad(T0 & temp_1,
5360  T0 & temp_2,
5361  const T1 &src_ptr_1,
5362  const T1 &src_ptr_2,
5363  const T2 &grad_weight)
5364  {
5365  // case 3a)
5366  temp_1 = src_ptr_1;
5367  temp_2 = grad_weight * (temp_1 - src_ptr_2);
5368  }
5369 
5370  template <typename T1, typename T2>
5371  void
5372  value(T1 &temp, const T2 &src_ptr)
5373  {
5374  // case 3b)
5375  temp = src_ptr;
5376  }
5377  };
5378  };
5379 
5380 
5381 
5382  template <int dim, typename Number2, typename VectorizedArrayType>
5384  {
5385  using Number = typename VectorizedArrayType::value_type;
5386 
5387  template <int fe_degree, int n_q_points_1d>
5388  static bool
5389  run(const unsigned int n_components,
5390  const EvaluationFlags::EvaluationFlags integration_flag,
5391  Number2 * dst_ptr,
5392  const std::vector<ArrayView<const Number2>> * sm_ptr,
5394  {
5395  Assert(fe_degree > -1, ExcInternalError());
5396  Assert(fe_eval.get_shape_info().element_type <=
5398  ExcInternalError());
5399 
5400  const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
5401 
5402  VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
5403  VectorizedArrayType *scratch_data =
5404  temp + 3 * n_components * dofs_per_face;
5405 
5406  const unsigned int subface_index = fe_eval.get_subface_index();
5407 
5408  // re-orientation for cases not possible with the io function below
5409  if (subface_index < GeometryInfo<dim>::max_children_per_cell)
5410  {
5411  if (fe_eval.get_dof_access_index() ==
5413  fe_eval.is_interior_face() == false)
5414  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5415  {
5416  // the loop breaks once an invalid_unsigned_int is hit for
5417  // all cases except the exterior faces in the ECL loop (where
5418  // some faces might be at the boundaries but others not)
5419  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
5420  continue;
5421 
5422  if (fe_eval.get_face_orientation(v) != 0)
5424  dim,
5425  n_components,
5426  v,
5427  integration_flag,
5429  fe_eval.get_face_orientation(v), 0),
5430  true,
5431  Utilities::pow(n_q_points_1d, dim - 1),
5432  &temp[0][0],
5433  fe_eval.begin_values(),
5434  fe_eval.begin_gradients(),
5435  fe_eval.begin_hessians());
5436  }
5437  else if (fe_eval.get_face_orientation() != 0)
5439  dim,
5440  n_components,
5441  integration_flag,
5443  fe_eval.get_face_orientation(), 0),
5444  true,
5445  Utilities::pow(n_q_points_1d, dim - 1),
5446  temp,
5447  fe_eval.begin_values(),
5448  fe_eval.begin_gradients(),
5449  fe_eval.begin_hessians());
5450  }
5451 
5452  if (fe_degree > -1 && fe_eval.get_subface_index() >=
5453  GeometryInfo<dim - 1>::max_children_per_cell)
5454  FEFaceEvaluationImpl<true,
5455  dim,
5456  fe_degree,
5457  n_q_points_1d,
5458  VectorizedArrayType>::
5459  integrate_in_face(n_components,
5460  integration_flag,
5461  fe_eval.get_shape_info().data.front(),
5462  temp,
5463  fe_eval.begin_values(),
5464  fe_eval.begin_gradients(),
5465  fe_eval.begin_hessians(),
5466  scratch_data,
5467  subface_index);
5468  else
5469  FEFaceEvaluationImpl<false,
5470  dim,
5471  fe_degree,
5472  n_q_points_1d,
5473  VectorizedArrayType>::
5474  integrate_in_face(n_components,
5475  integration_flag,
5476  fe_eval.get_shape_info().data.front(),
5477  temp,
5478  fe_eval.begin_values(),
5479  fe_eval.begin_gradients(),
5480  fe_eval.begin_hessians(),
5481  scratch_data,
5482  subface_index);
5483 
5485 
5486  if (fe_eval.get_dof_access_index() ==
5488  fe_eval.is_interior_face() == false)
5489  fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
5490  p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
5491  else
5492  fe_face_evaluation_process_and_io<1>(
5493  p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
5494 
5495  return false;
5496  }
5497 
5498  private:
5499  template <int fe_degree>
5500  struct Processor
5501  {
5502  static const bool do_integrate = true;
5503  static const int dim_ = dim;
5504  static const int fe_degree_ = fe_degree;
5505  using VectorizedArrayType_ = VectorizedArrayType;
5506  using Number_ = Number;
5507  using Number2_ = Number2;
5508 
5509  template <typename T0, typename T1, typename T2, typename T3, typename T4>
5510  void
5511  hermite_grad_vectorized(const T0 &temp_1,
5512  const T1 &temp_2,
5513  T2 dst_ptr_1,
5514  T3 dst_ptr_2,
5515  const T4 &grad_weight)
5516  {
5517  // case 1a)
5518  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
5519  const VectorizedArrayType grad = grad_weight * temp_2;
5520  do_vectorized_add(val, dst_ptr_1);
5521  do_vectorized_add(grad, dst_ptr_2);
5522  }
5523 
5524  template <typename T0, typename T1>
5525  void
5526  value_vectorized(const T0 &temp, T1 dst_ptr)
5527  {
5528  // case 1b)
5529  do_vectorized_add(temp, dst_ptr);
5530  }
5531 
5532  template <typename T0, typename T1, typename T2, typename T3>
5533  void
5535  const T0 &temp_2,
5536  T1 dst_ptr_1,
5537  T1 dst_ptr_2,
5538  const T2 &grad_weight,
5539  const T3 &indices_1,
5540  const T3 &indices_2)
5541  {
5542  // case 2a)
5543  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
5544  const VectorizedArrayType grad = grad_weight * temp_2;
5545  do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
5546  do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
5547  }
5548 
5549  template <typename T0, typename T1, typename T2>
5550  void
5551  value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
5552  {
5553  // case 2b)
5554  do_vectorized_scatter_add(temp, indices, dst_ptr);
5555  }
5556 
5557  template <typename T0, typename T1, typename T2>
5558  void
5559  hermite_grad(const T0 &temp_1,
5560  const T0 &temp_2,
5561  T1 & dst_ptr_1,
5562  T1 & dst_ptr_2,
5563  const T2 &grad_weight)
5564  {
5565  // case 3a)
5566  const Number val = temp_1 - grad_weight * temp_2;
5567  const Number grad = grad_weight * temp_2;
5568  dst_ptr_1 += val;
5569  dst_ptr_2 += grad;
5570  }
5571 
5572  template <typename T0, typename T1>
5573  void
5574  value(const T0 &temp, T1 &dst_ptr)
5575  {
5576  // case 3b)
5577  dst_ptr += temp;
5578  }
5579  };
5580  };
5581 
5582 
5583 
5588  template <int dim, typename Number>
5590  {
5591  template <int fe_degree, int = 0>
5592  static bool
5593  run(const unsigned int n_components,
5594  const FEEvaluationData<dim, Number, false> &fe_eval,
5595  const Number * in_array,
5596  Number * out_array,
5597  typename std::enable_if<fe_degree != -1>::type * = nullptr)
5598  {
5599  constexpr unsigned int dofs_per_component =
5600  Utilities::pow(fe_degree + 1, dim);
5601 
5602  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
5603  Assert(fe_eval.get_shape_info().element_type <=
5605  ExcNotImplemented());
5606 
5608  dim,
5609  fe_degree + 1,
5610  fe_degree + 1,
5611  Number>
5612  evaluator(
5615  fe_eval.get_shape_info().data.front().inverse_shape_values_eo);
5616 
5617  for (unsigned int d = 0; d < n_components; ++d)
5618  {
5619  const Number *in = in_array + d * dofs_per_component;
5620  Number * out = out_array + d * dofs_per_component;
5621  // Need to select 'apply' method with hessian slot because values
5622  // assume symmetries that do not exist in the inverse shapes
5623  evaluator.template hessians<0, true, false>(in, out);
5624  if (dim > 1)
5625  evaluator.template hessians<1, true, false>(out, out);
5626  if (dim > 2)
5627  evaluator.template hessians<2, true, false>(out, out);
5628  }
5629  for (unsigned int q = 0; q < dofs_per_component; ++q)
5630  {
5631  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
5632  for (unsigned int d = 0; d < n_components; ++d)
5633  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
5634  }
5635  for (unsigned int d = 0; d < n_components; ++d)
5636  {
5637  Number *out = out_array + d * dofs_per_component;
5638  if (dim > 2)
5639  evaluator.template hessians<2, false, false>(out, out);
5640  if (dim > 1)
5641  evaluator.template hessians<1, false, false>(out, out);
5642  evaluator.template hessians<0, false, false>(out, out);
5643  }
5644  return false;
5645  }
5646 
5647  template <int fe_degree, int = 0>
5648  static bool
5649  run(const unsigned int n_components,
5650  const FEEvaluationData<dim, Number, false> &fe_eval,
5651  const Number * in_array,
5652  Number * out_array,
5653  typename std::enable_if<fe_degree == -1>::type * = nullptr)
5654  {
5655  static_assert(fe_degree == -1, "Only usable for degree -1");
5656  const unsigned int dofs_per_component =
5658 
5659  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
5660 
5662  fe_eval.get_shape_info().data.front().inverse_shape_values,
5665  fe_eval.get_shape_info().data.front().fe_degree + 1,
5666  fe_eval.get_shape_info().data.front().fe_degree + 1);
5667 
5668  for (unsigned int d = 0; d < n_components; ++d)
5669  {
5670  const Number *in = in_array + d * dofs_per_component;
5671  Number * out = out_array + d * dofs_per_component;
5672  // Need to select 'apply' method with hessian slot because values
5673  // assume symmetries that do not exist in the inverse shapes
5674  evaluator.template values<0, true, false>(in, out);
5675  if (dim > 1)
5676  evaluator.template values<1, true, false>(out, out);
5677  if (dim > 2)
5678  evaluator.template values<2, true, false>(out, out);
5679  }
5680  for (unsigned int q = 0; q < dofs_per_component; ++q)
5681  {
5682  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
5683  for (unsigned int d = 0; d < n_components; ++d)
5684  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
5685  }
5686  for (unsigned int d = 0; d < n_components; ++d)
5687  {
5688  Number *out = out_array + d * dofs_per_component;
5689  if (dim > 2)
5690  evaluator.template values<2, false, false>(out, out);
5691  if (dim > 1)
5692  evaluator.template values<1, false, false>(out, out);
5693  evaluator.template values<0, false, false>(out, out);
5694  }
5695  return false;
5696  }
5697  };
5698 
5699 
5700 
5705  template <int dim, typename Number>
5707  {
5708  template <int fe_degree, int = 0>
5709  static bool
5710  run(const unsigned int n_desired_components,
5711  const AlignedVector<Number> &inverse_shape,
5712  const AlignedVector<Number> &inverse_coefficients,
5713  const Number * in_array,
5714  Number * out_array,
5715  typename std::enable_if<fe_degree != -1>::type * = nullptr)
5716  {
5717  constexpr unsigned int dofs_per_component =
5718  Utilities::pow(fe_degree + 1, dim);
5719  Assert(inverse_coefficients.size() > 0 &&
5720  inverse_coefficients.size() % dofs_per_component == 0,
5721  ExcMessage(
5722  "Expected diagonal to be a multiple of scalar dof per cells"));
5723  if (inverse_coefficients.size() != dofs_per_component)
5724  AssertDimension(n_desired_components * dofs_per_component,
5725  inverse_coefficients.size());
5726 
5727  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
5728 
5730  dim,
5731  fe_degree + 1,
5732  fe_degree + 1,
5733  Number>
5734  evaluator(AlignedVector<Number>(),
5736  inverse_shape);
5737 
5738  const unsigned int shift_coefficient =
5739  inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
5740  0;
5741  const Number *inv_coefficient = inverse_coefficients.data();
5742  for (unsigned int d = 0; d < n_desired_components; ++d)
5743  {
5744  const Number *in = in_array + d * dofs_per_component;
5745  Number * out = out_array + d * dofs_per_component;
5746  // Need to select 'apply' method with hessian slot because values
5747  // assume symmetries that do not exist in the inverse shapes
5748  evaluator.template hessians<0, true, false>(in, out);
5749  if (dim > 1)
5750  evaluator.template hessians<1, true, false>(out, out);
5751  if (dim > 2)
5752  evaluator.template hessians<2, true, false>(out, out);
5753 
5754  for (unsigned int q = 0; q < dofs_per_component; ++q)
5755  out[q] *= inv_coefficient[q];
5756 
5757  if (dim > 2)
5758  evaluator.template hessians<2, false, false>(out, out);
5759  if (dim > 1)
5760  evaluator.template hessians<1, false, false>(out, out);
5761  evaluator.template hessians<0, false, false>(out, out);
5762 
5763  inv_coefficient += shift_coefficient;
5764  }
5765  return false;
5766  }
5767 
5771  template <int fe_degree, int = 0>
5772  static bool
5773  run(const unsigned int,
5774  const AlignedVector<Number> &,
5775  const AlignedVector<Number> &,
5776  const Number *,
5777  Number *,
5778  typename std::enable_if<fe_degree == -1>::type * = nullptr)
5779  {
5780  static_assert(fe_degree == -1, "Only usable for degree -1");
5781  Assert(false, ExcNotImplemented());
5782  return false;
5783  }
5784  };
5785 
5786 
5787 
5792  template <int dim, typename Number>
5794  {
5795  template <int fe_degree, int n_q_points_1d>
5796  static bool
5797  run(const unsigned int n_desired_components,
5798  const FEEvaluationData<dim, Number, false> &fe_eval,
5799  const Number * in_array,
5800  Number * out_array)
5801  {
5802  static const bool do_inplace =
5803  fe_degree > -1 && (fe_degree + 1 == n_q_points_1d);
5804 
5805  Assert(fe_eval.get_shape_info().element_type !=
5807  ExcNotImplemented());
5808 
5809  const auto &inverse_shape =
5810  do_inplace ?
5811  fe_eval.get_shape_info().data.front().inverse_shape_values_eo :
5812  fe_eval.get_shape_info().data.front().inverse_shape_values;
5813 
5814  const std::size_t dofs_per_component =
5815  do_inplace ? Utilities::pow(fe_degree + 1, dim) :
5817  const std::size_t n_q_points = do_inplace ?
5818  Utilities::pow(fe_degree + 1, dim) :
5819  fe_eval.get_shape_info().n_q_points;
5820 
5822  dim,
5823  fe_degree + 1,
5824  n_q_points_1d,
5825  Number>
5826  evaluator(AlignedVector<Number>(),
5828  inverse_shape,
5829  fe_eval.get_shape_info().data.front().fe_degree + 1,
5830  fe_eval.get_shape_info().data.front().n_q_points_1d);
5831 
5832  for (unsigned int d = 0; d < n_desired_components; ++d)
5833  {
5834  const Number *in = in_array + d * n_q_points;
5835  Number * out = out_array + d * dofs_per_component;
5836 
5837  auto temp_1 = do_inplace ? out : fe_eval.get_scratch_data().begin();
5838  auto temp_2 = do_inplace ?
5839  out :
5840  (temp_1 + std::max(n_q_points, dofs_per_component));
5841 
5842  if (dim == 3)
5843  {
5844  evaluator.template hessians<2, false, false>(in, temp_1);
5845  evaluator.template hessians<1, false, false>(temp_1, temp_2);
5846  evaluator.template hessians<0, false, false>(temp_2, out);
5847  }
5848  if (dim == 2)
5849  {
5850  evaluator.template hessians<1, false, false>(in, temp_1);
5851  evaluator.template hessians<0, false, false>(temp_1, out);
5852  }
5853  if (dim == 1)
5854  evaluator.template hessians<0, false, false>(in, out);
5855  }
5856  return false;
5857  }
5858  };
5859 
5860 } // end of namespace internal
5861 
5862 
5864 
5865 #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 Number *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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
unsigned int cell_index
Definition: grid_tools.cc:1129
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
EvaluationFlags
The EvaluationFlags enum.
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2050
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type 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)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:462
T fixed_power(const T t)
Definition: utilities.h:1123
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:201
static bool run(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array, typename std::enable_if< fe_degree !=-1 >::type *=nullptr)
static bool run(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array, typename std::enable_if< fe_degree==-1 >::type *=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, typename std::enable_if< fe_degree !=-1 >::type *=nullptr)
static bool run(const unsigned int, const AlignedVector< Number > &, const AlignedVector< Number > &, const Number *, Number *, typename std::enable_if< fe_degree==-1 >::type *=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)
typename VectorizedArrayType::value_type Number