Reference documentation for deal.II version Git 8113613837 2020-09-22 13:45:23 -0500
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
evaluation_kernels.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_evaluation_kernels_h
18 #define dealii_matrix_free_evaluation_kernels_h
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/utilities.h>
24 
30 
31 
33 
34 
35 // forward declaration
36 template <int, typename, bool, typename>
38 
39 
40 
41 namespace internal
42 {
43  // Select evaluator type from element shape function type
44  template <MatrixFreeFunctions::ElementType element, bool is_long>
46  {};
47 
48  template <bool is_long>
49  struct EvaluatorSelector<MatrixFreeFunctions::tensor_general, is_long>
50  {
51  static const EvaluatorVariant variant = evaluate_general;
52  };
53 
54  template <>
55  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, false>
56  {
57  static const EvaluatorVariant variant = evaluate_symmetric;
58  };
59 
60  template <>
61  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, true>
62  {
63  static const EvaluatorVariant variant = evaluate_evenodd;
64  };
65 
66  template <bool is_long>
67  struct EvaluatorSelector<MatrixFreeFunctions::truncated_tensor, is_long>
68  {
69  static const EvaluatorVariant variant = evaluate_general;
70  };
71 
72  template <>
73  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,
74  false>
75  {
76  static const EvaluatorVariant variant = evaluate_general;
77  };
78 
79  template <>
80  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0, true>
81  {
82  static const EvaluatorVariant variant = evaluate_evenodd;
83  };
84 
85  template <bool is_long>
86  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_collocation,
87  is_long>
88  {
89  static const EvaluatorVariant variant = evaluate_evenodd;
90  };
91 
92 
93 
110  template <MatrixFreeFunctions::ElementType type,
111  int dim,
112  int fe_degree,
113  int n_q_points_1d,
114  typename Number>
116  {
117  static void
118  evaluate(const unsigned int n_components,
119  const EvaluationFlags::EvaluationFlags evaluation_flag,
120  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
121  const Number * values_dofs_actual,
122  Number * values_quad,
123  Number * gradients_quad,
124  Number * hessians_quad,
125  Number * scratch_data);
126 
127  static void
128  integrate(const unsigned int n_components,
129  const EvaluationFlags::EvaluationFlags integration_flag,
130  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
131  Number * values_dofs_actual,
132  Number * values_quad,
133  Number * gradients_quad,
134  Number * scratch_data,
135  const bool add_into_values_array);
136  };
137 
138 
139 
140  template <MatrixFreeFunctions::ElementType type,
141  int dim,
142  int fe_degree,
143  int n_q_points_1d,
144  typename Number>
145  inline void
147  const unsigned int n_components,
148  const EvaluationFlags::EvaluationFlags evaluation_flag,
149  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
150  const Number * values_dofs_actual,
151  Number * values_quad,
152  Number * gradients_quad,
153  Number * hessians_quad,
154  Number * scratch_data)
155  {
156  if (evaluation_flag == EvaluationFlags::nothing)
157  return;
158 
159  const EvaluatorVariant variant =
161  using Eval = EvaluatorTensorProduct<variant,
162  dim,
163  fe_degree + 1,
164  n_q_points_1d,
165  Number>;
166  Eval eval(variant == evaluate_evenodd ?
167  shape_info.data.front().shape_values_eo :
168  shape_info.data.front().shape_values,
169  variant == evaluate_evenodd ?
170  shape_info.data.front().shape_gradients_eo :
171  shape_info.data.front().shape_gradients,
172  variant == evaluate_evenodd ?
173  shape_info.data.front().shape_hessians_eo :
174  shape_info.data.front().shape_hessians,
175  shape_info.data.front().fe_degree + 1,
176  shape_info.data.front().n_q_points_1d);
177 
178  const unsigned int temp_size =
179  Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
180  0 :
181  (Eval::n_rows_of_product > Eval::n_columns_of_product ?
182  Eval::n_rows_of_product :
183  Eval::n_columns_of_product);
184  Number *temp1;
185  Number *temp2;
186  if (temp_size == 0)
187  {
188  temp1 = scratch_data;
189  temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
190  shape_info.data.front().fe_degree + 1),
191  Utilities::fixed_power<dim>(
192  shape_info.data.front().n_q_points_1d));
193  }
194  else
195  {
196  temp1 = scratch_data;
197  temp2 = temp1 + temp_size;
198  }
199 
200  const unsigned int n_q_points =
201  temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
202  const unsigned int dofs_per_comp =
204  Utilities::fixed_power<dim>(shape_info.data.front().fe_degree + 1) :
205  shape_info.dofs_per_component_on_cell;
206  const Number *values_dofs = values_dofs_actual;
208  {
209  Number *values_dofs_tmp =
210  scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
211  shape_info.n_q_points));
212  const int degree =
213  fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
214  for (unsigned int c = 0; c < n_components; ++c)
215  for (int i = 0, count_p = 0, count_q = 0;
216  i < (dim > 2 ? degree + 1 : 1);
217  ++i)
218  {
219  for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
220  {
221  for (int k = 0; k < degree + 1 - j - i;
222  ++k, ++count_p, ++count_q)
223  values_dofs_tmp[c * dofs_per_comp + count_q] =
224  values_dofs_actual
225  [c * shape_info.dofs_per_component_on_cell + count_p];
226  for (int k = degree + 1 - j - i; k < degree + 1;
227  ++k, ++count_q)
228  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
229  }
230  for (int j = degree + 1 - i; j < degree + 1; ++j)
231  for (int k = 0; k < degree + 1; ++k, ++count_q)
232  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
233  }
234  values_dofs = values_dofs_tmp;
235  }
236 
237  switch (dim)
238  {
239  case 1:
240  for (unsigned int c = 0; c < n_components; c++)
241  {
242  if (evaluation_flag & EvaluationFlags::values)
243  eval.template values<0, true, false>(values_dofs, values_quad);
244  if (evaluation_flag & EvaluationFlags::gradients)
245  eval.template gradients<0, true, false>(values_dofs,
246  gradients_quad);
247  if (evaluation_flag & EvaluationFlags::hessians)
248  eval.template hessians<0, true, false>(values_dofs,
249  hessians_quad);
250 
251  // advance the next component in 1D array
252  values_dofs += dofs_per_comp;
253  values_quad += n_q_points;
254  gradients_quad += n_q_points;
255  hessians_quad += n_q_points;
256  }
257  break;
258 
259  case 2:
260  for (unsigned int c = 0; c < n_components; c++)
261  {
262  // grad x
263  if (evaluation_flag & EvaluationFlags::gradients)
264  {
265  eval.template gradients<0, true, false>(values_dofs, temp1);
266  eval.template values<1, true, false>(temp1, gradients_quad);
267  }
268  if (evaluation_flag & EvaluationFlags::hessians)
269  {
270  // grad xy
271  if (!(evaluation_flag & EvaluationFlags::gradients))
272  eval.template gradients<0, true, false>(values_dofs, temp1);
273  eval.template gradients<1, true, false>(temp1,
274  hessians_quad +
275  2 * n_q_points);
276 
277  // grad xx
278  eval.template hessians<0, true, false>(values_dofs, temp1);
279  eval.template values<1, true, false>(temp1, hessians_quad);
280  }
281 
282  // grad y
283  eval.template values<0, true, false>(values_dofs, temp1);
284  if (evaluation_flag & EvaluationFlags::gradients)
285  eval.template gradients<1, true, false>(temp1,
286  gradients_quad +
287  n_q_points);
288 
289  // grad yy
290  if (evaluation_flag & EvaluationFlags::hessians)
291  eval.template hessians<1, true, false>(temp1,
292  hessians_quad +
293  n_q_points);
294 
295  // val: can use values applied in x
296  if (evaluation_flag & EvaluationFlags::values)
297  eval.template values<1, true, false>(temp1, values_quad);
298 
299  // advance to the next component in 1D array
300  values_dofs += dofs_per_comp;
301  values_quad += n_q_points;
302  gradients_quad += 2 * n_q_points;
303  hessians_quad += 3 * n_q_points;
304  }
305  break;
306 
307  case 3:
308  for (unsigned int c = 0; c < n_components; c++)
309  {
310  if (evaluation_flag & EvaluationFlags::gradients)
311  {
312  // grad x
313  eval.template gradients<0, true, false>(values_dofs, temp1);
314  eval.template values<1, true, false>(temp1, temp2);
315  eval.template values<2, true, false>(temp2, gradients_quad);
316  }
317 
318  if (evaluation_flag & EvaluationFlags::hessians)
319  {
320  // grad xz
321  if (!(evaluation_flag & EvaluationFlags::gradients))
322  {
323  eval.template gradients<0, true, false>(values_dofs,
324  temp1);
325  eval.template values<1, true, false>(temp1, temp2);
326  }
327  eval.template gradients<2, true, false>(temp2,
328  hessians_quad +
329  4 * n_q_points);
330 
331  // grad xy
332  eval.template gradients<1, true, false>(temp1, temp2);
333  eval.template values<2, true, false>(temp2,
334  hessians_quad +
335  3 * n_q_points);
336 
337  // grad xx
338  eval.template hessians<0, true, false>(values_dofs, temp1);
339  eval.template values<1, true, false>(temp1, temp2);
340  eval.template values<2, true, false>(temp2, hessians_quad);
341  }
342 
343  // grad y
344  eval.template values<0, true, false>(values_dofs, temp1);
345  if (evaluation_flag & EvaluationFlags::gradients)
346  {
347  eval.template gradients<1, true, false>(temp1, temp2);
348  eval.template values<2, true, false>(temp2,
349  gradients_quad +
350  n_q_points);
351  }
352 
353  if (evaluation_flag & EvaluationFlags::hessians)
354  {
355  // grad yz
356  if (!(evaluation_flag & EvaluationFlags::gradients))
357  eval.template gradients<1, true, false>(temp1, temp2);
358  eval.template gradients<2, true, false>(temp2,
359  hessians_quad +
360  5 * n_q_points);
361 
362  // grad yy
363  eval.template hessians<1, true, false>(temp1, temp2);
364  eval.template values<2, true, false>(temp2,
365  hessians_quad +
366  n_q_points);
367  }
368 
369  // grad z: can use the values applied in x direction stored in
370  // temp1
371  eval.template values<1, true, false>(temp1, temp2);
372  if (evaluation_flag & EvaluationFlags::gradients)
373  eval.template gradients<2, true, false>(temp2,
374  gradients_quad +
375  2 * n_q_points);
376 
377  // grad zz: can use the values applied in x and y direction stored
378  // in temp2
379  if (evaluation_flag & EvaluationFlags::hessians)
380  eval.template hessians<2, true, false>(temp2,
381  hessians_quad +
382  2 * n_q_points);
383 
384  // val: can use the values applied in x & y direction stored in
385  // temp2
386  if (evaluation_flag & EvaluationFlags::values)
387  eval.template values<2, true, false>(temp2, values_quad);
388 
389  // advance to the next component in 1D array
390  values_dofs += dofs_per_comp;
391  values_quad += n_q_points;
392  gradients_quad += 3 * n_q_points;
393  hessians_quad += 6 * n_q_points;
394  }
395  break;
396 
397  default:
398  AssertThrow(false, ExcNotImplemented());
399  }
400 
401  // case additional dof for FE_Q_DG0: add values; gradients and second
402  // derivatives evaluate to zero
404  (evaluation_flag & EvaluationFlags::values))
405  {
406  values_quad -= n_components * n_q_points;
407  values_dofs -= n_components * dofs_per_comp;
408  for (unsigned int c = 0; c < n_components; ++c)
409  for (unsigned int q = 0; q < shape_info.n_q_points; ++q)
410  values_quad[c * shape_info.n_q_points + q] +=
411  values_dofs[(c + 1) * shape_info.dofs_per_component_on_cell - 1];
412  }
413  }
414 
415 
416 
417  template <MatrixFreeFunctions::ElementType type,
418  int dim,
419  int fe_degree,
420  int n_q_points_1d,
421  typename Number>
422  inline void
424  const unsigned int n_components,
425  const EvaluationFlags::EvaluationFlags integration_flag,
426  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
427  Number * values_dofs_actual,
428  Number * values_quad,
429  Number * gradients_quad,
430  Number * scratch_data,
431  const bool add_into_values_array)
432  {
433  const EvaluatorVariant variant =
435  using Eval = EvaluatorTensorProduct<variant,
436  dim,
437  fe_degree + 1,
438  n_q_points_1d,
439  Number>;
440  Eval eval(variant == evaluate_evenodd ?
441  shape_info.data.front().shape_values_eo :
442  shape_info.data.front().shape_values,
443  variant == evaluate_evenodd ?
444  shape_info.data.front().shape_gradients_eo :
445  shape_info.data.front().shape_gradients,
446  variant == evaluate_evenodd ?
447  shape_info.data.front().shape_hessians_eo :
448  shape_info.data.front().shape_hessians,
449  shape_info.data.front().fe_degree + 1,
450  shape_info.data.front().n_q_points_1d);
451 
452  const unsigned int temp_size =
453  Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
454  0 :
455  (Eval::n_rows_of_product > Eval::n_columns_of_product ?
456  Eval::n_rows_of_product :
457  Eval::n_columns_of_product);
458  Number *temp1;
459  Number *temp2;
460  if (temp_size == 0)
461  {
462  temp1 = scratch_data;
463  temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
464  shape_info.data.front().fe_degree + 1),
465  Utilities::fixed_power<dim>(
466  shape_info.data.front().n_q_points_1d));
467  }
468  else
469  {
470  temp1 = scratch_data;
471  temp2 = temp1 + temp_size;
472  }
473 
474  const unsigned int n_q_points =
475  temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
476  const unsigned int dofs_per_comp =
478  Utilities::fixed_power<dim>(shape_info.data.front().fe_degree + 1) :
479  shape_info.dofs_per_component_on_cell;
480  // expand dof_values to tensor product for truncated tensor products
481  Number *values_dofs =
483  scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
484  shape_info.n_q_points)) :
485  values_dofs_actual;
486 
487  switch (dim)
488  {
489  case 1:
490  for (unsigned int c = 0; c < n_components; c++)
491  {
492  if (integration_flag & EvaluationFlags::values)
493  {
494  if (add_into_values_array == false)
495  eval.template values<0, false, false>(values_quad,
496  values_dofs);
497  else
498  eval.template values<0, false, true>(values_quad,
499  values_dofs);
500  }
501  if (integration_flag & EvaluationFlags::gradients)
502  {
503  if (integration_flag & EvaluationFlags::values ||
504  add_into_values_array == true)
505  eval.template gradients<0, false, true>(gradients_quad,
506  values_dofs);
507  else
508  eval.template gradients<0, false, false>(gradients_quad,
509  values_dofs);
510  }
511 
512  // advance to the next component in 1D array
513  values_dofs += dofs_per_comp;
514  values_quad += n_q_points;
515  gradients_quad += n_q_points;
516  }
517  break;
518 
519  case 2:
520  for (unsigned int c = 0; c < n_components; c++)
521  {
522  if ((integration_flag & EvaluationFlags::values) &&
523  !(integration_flag & EvaluationFlags::gradients))
524  {
525  eval.template values<1, false, false>(values_quad, temp1);
526  if (add_into_values_array == false)
527  eval.template values<0, false, false>(temp1, values_dofs);
528  else
529  eval.template values<0, false, true>(temp1, values_dofs);
530  }
531  if (integration_flag & EvaluationFlags::gradients)
532  {
533  eval.template gradients<1, false, false>(gradients_quad +
534  n_q_points,
535  temp1);
536  if (integration_flag & EvaluationFlags::values)
537  eval.template values<1, false, true>(values_quad, temp1);
538  if (add_into_values_array == false)
539  eval.template values<0, false, false>(temp1, values_dofs);
540  else
541  eval.template values<0, false, true>(temp1, values_dofs);
542  eval.template values<1, false, false>(gradients_quad, temp1);
543  eval.template gradients<0, false, true>(temp1, values_dofs);
544  }
545 
546  // advance to the next component in 1D array
547  values_dofs += dofs_per_comp;
548  values_quad += n_q_points;
549  gradients_quad += 2 * n_q_points;
550  }
551  break;
552 
553  case 3:
554  for (unsigned int c = 0; c < n_components; c++)
555  {
556  if ((integration_flag & EvaluationFlags::values) &&
557  !(integration_flag & EvaluationFlags::gradients))
558  {
559  eval.template values<2, false, false>(values_quad, temp1);
560  eval.template values<1, false, false>(temp1, temp2);
561  if (add_into_values_array == false)
562  eval.template values<0, false, false>(temp2, values_dofs);
563  else
564  eval.template values<0, false, true>(temp2, values_dofs);
565  }
566  if (integration_flag & EvaluationFlags::gradients)
567  {
568  eval.template gradients<2, false, false>(gradients_quad +
569  2 * n_q_points,
570  temp1);
571  if (integration_flag & EvaluationFlags::values)
572  eval.template values<2, false, true>(values_quad, temp1);
573  eval.template values<1, false, false>(temp1, temp2);
574  eval.template values<2, false, false>(gradients_quad +
575  n_q_points,
576  temp1);
577  eval.template gradients<1, false, true>(temp1, temp2);
578  if (add_into_values_array == false)
579  eval.template values<0, false, false>(temp2, values_dofs);
580  else
581  eval.template values<0, false, true>(temp2, values_dofs);
582  eval.template values<2, false, false>(gradients_quad, temp1);
583  eval.template values<1, false, false>(temp1, temp2);
584  eval.template gradients<0, false, true>(temp2, values_dofs);
585  }
586 
587  // advance to the next component in 1D array
588  values_dofs += dofs_per_comp;
589  values_quad += n_q_points;
590  gradients_quad += 3 * n_q_points;
591  }
592  break;
593 
594  default:
595  AssertThrow(false, ExcNotImplemented());
596  }
597 
598  // case FE_Q_DG0: add values, gradients and second derivatives are zero
600  {
601  values_dofs -= n_components * dofs_per_comp -
602  shape_info.dofs_per_component_on_cell + 1;
603  values_quad -= n_components * n_q_points;
604  if (integration_flag & EvaluationFlags::values)
605  for (unsigned int c = 0; c < n_components; ++c)
606  {
607  values_dofs[0] = values_quad[0];
608  for (unsigned int q = 1; q < shape_info.n_q_points; ++q)
609  values_dofs[0] += values_quad[q];
610  values_dofs += dofs_per_comp;
611  values_quad += n_q_points;
612  }
613  else
614  {
615  for (unsigned int c = 0; c < n_components; ++c)
616  values_dofs[c * shape_info.dofs_per_component_on_cell] = Number();
617  values_dofs += n_components * shape_info.dofs_per_component_on_cell;
618  }
619  }
620 
622  {
623  values_dofs -= dofs_per_comp * n_components;
624  const int degree =
625  fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
626  for (unsigned int c = 0; c < n_components; ++c)
627  for (int i = 0, count_p = 0, count_q = 0;
628  i < (dim > 2 ? degree + 1 : 1);
629  ++i)
630  {
631  for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
632  {
633  for (int k = 0; k < degree + 1 - j - i;
634  ++k, ++count_p, ++count_q)
635  values_dofs_actual[c *
636  shape_info.dofs_per_component_on_cell +
637  count_p] =
638  values_dofs[c * dofs_per_comp + count_q];
639  count_q += j + i;
640  }
641  count_q += i * (degree + 1);
642  }
643  }
644  }
645 
646 
647 
657  template <EvaluatorVariant variant,
658  EvaluatorQuantity quantity,
659  int dim,
660  int basis_size_1,
661  int basis_size_2,
662  typename Number,
663  typename Number2>
665  {
666  static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
667  "The second dimension must not be smaller than the first");
668 
690 #ifndef DEBUG
692 #endif
693  static void
695  const unsigned int n_components,
696  const AlignedVector<Number2> &transformation_matrix,
697  const Number * values_in,
698  Number * values_out,
699  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
700  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
701  {
702  Assert(
703  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
704  ExcMessage("The second dimension must not be smaller than the first"));
705 
707 
708  // we do recursion until dim==1 or dim==2 and we have
709  // basis_size_1==basis_size_2. The latter optimization increases
710  // optimization possibilities for the compiler but does only work for
711  // aliased pointers if the sizes are equal.
712  constexpr int next_dim =
713  (dim > 2 ||
714  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
715  dim - 1 :
716  dim;
717 
718  EvaluatorTensorProduct<variant,
719  dim,
720  basis_size_1,
721  (basis_size_1 == 0 ? 0 : basis_size_2),
722  Number,
723  Number2>
724  eval_val(transformation_matrix,
727  basis_size_1_variable,
728  basis_size_2_variable);
729  const unsigned int np_1 =
730  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
731  const unsigned int np_2 =
732  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
733  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
734  ExcMessage("Cannot transform with 0-point basis"));
735  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
736  ExcMessage("Cannot transform with 0-point basis"));
737 
738  // run loop backwards to ensure correctness if values_in aliases with
739  // values_out in case with basis_size_1 < basis_size_2
740  values_in = values_in + n_components * Utilities::fixed_power<dim>(np_1);
741  values_out =
742  values_out + n_components * Utilities::fixed_power<dim>(np_2);
743  for (unsigned int c = n_components; c != 0; --c)
744  {
745  values_in -= Utilities::fixed_power<dim>(np_1);
746  values_out -= Utilities::fixed_power<dim>(np_2);
747  if (next_dim < dim)
748  for (unsigned int q = np_1; q != 0; --q)
750  variant,
751  quantity,
752  next_dim,
753  basis_size_1,
754  basis_size_2,
755  Number,
756  Number2>::do_forward(1,
757  transformation_matrix,
758  values_in +
759  (q - 1) *
760  Utilities::fixed_power<next_dim>(np_1),
761  values_out +
762  (q - 1) *
763  Utilities::fixed_power<next_dim>(np_2),
764  basis_size_1_variable,
765  basis_size_2_variable);
766 
767  // the recursion stops if dim==1 or if dim==2 and
768  // basis_size_1==basis_size_2 (the latter is used because the
769  // compiler generates nicer code)
770  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
771  {
772  eval_val.template values<0, true, false>(values_in, values_out);
773  eval_val.template values<1, true, false>(values_out, values_out);
774  }
775  else if (dim == 1)
776  eval_val.template values<dim - 1, true, false>(values_in,
777  values_out);
778  else
779  eval_val.template values<dim - 1, true, false>(values_out,
780  values_out);
781  }
782  }
783 
813 #ifndef DEBUG
815 #endif
816  static void
818  const unsigned int n_components,
819  const AlignedVector<Number2> &transformation_matrix,
820  const bool add_into_result,
821  Number * values_in,
822  Number * values_out,
823  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
824  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
825  {
826  Assert(
827  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
828  ExcMessage("The second dimension must not be smaller than the first"));
829  Assert(add_into_result == false || values_in != values_out,
830  ExcMessage(
831  "Input and output cannot alias with each other when "
832  "adding the result of the basis change to existing data"));
833 
834  Assert(quantity == EvaluatorQuantity::value ||
835  quantity == EvaluatorQuantity::hessian,
836  ExcInternalError());
837 
838  constexpr int next_dim =
839  (dim > 2 ||
840  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
841  dim - 1 :
842  dim;
843  EvaluatorTensorProduct<variant,
844  dim,
845  basis_size_1,
846  (basis_size_1 == 0 ? 0 : basis_size_2),
847  Number,
848  Number2>
849  eval_val(transformation_matrix,
850  transformation_matrix,
851  transformation_matrix,
852  basis_size_1_variable,
853  basis_size_2_variable);
854  const unsigned int np_1 =
855  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
856  const unsigned int np_2 =
857  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
858  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
859  ExcMessage("Cannot transform with 0-point basis"));
860  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
861  ExcMessage("Cannot transform with 0-point basis"));
862 
863  for (unsigned int c = 0; c < n_components; ++c)
864  {
865  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
866  {
867  if (quantity == EvaluatorQuantity::value)
868  eval_val.template values<1, false, false>(values_in, values_in);
869  else
870  eval_val.template hessians<1, false, false>(values_in,
871  values_in);
872 
873  if (add_into_result)
874  {
875  if (quantity == EvaluatorQuantity::value)
876  eval_val.template values<0, false, true>(values_in,
877  values_out);
878  else
879  eval_val.template hessians<0, false, true>(values_in,
880  values_out);
881  }
882  else
883  {
884  if (quantity == EvaluatorQuantity::value)
885  eval_val.template values<0, false, false>(values_in,
886  values_out);
887  else
888  eval_val.template hessians<0, false, false>(values_in,
889  values_out);
890  }
891  }
892  else
893  {
894  if (dim == 1 && add_into_result)
895  {
896  if (quantity == EvaluatorQuantity::value)
897  eval_val.template values<0, false, true>(values_in,
898  values_out);
899  else
900  eval_val.template hessians<0, false, true>(values_in,
901  values_out);
902  }
903  else if (dim == 1)
904  {
905  if (quantity == EvaluatorQuantity::value)
906  eval_val.template values<0, false, false>(values_in,
907  values_out);
908  else
909  eval_val.template hessians<0, false, false>(values_in,
910  values_out);
911  }
912  else
913  {
914  if (quantity == EvaluatorQuantity::value)
915  eval_val.template values<dim - 1, false, false>(values_in,
916  values_in);
917  else
918  eval_val.template hessians<dim - 1, false, false>(
919  values_in, values_in);
920  }
921  }
922  if (next_dim < dim)
923  for (unsigned int q = 0; q < np_1; ++q)
925  quantity,
926  next_dim,
927  basis_size_1,
928  basis_size_2,
929  Number,
930  Number2>::
931  do_backward(1,
932  transformation_matrix,
933  add_into_result,
934  values_in +
935  q * Utilities::fixed_power<next_dim>(np_2),
936  values_out +
937  q * Utilities::fixed_power<next_dim>(np_1),
938  basis_size_1_variable,
939  basis_size_2_variable);
940 
941  values_in += Utilities::fixed_power<dim>(np_2);
942  values_out += Utilities::fixed_power<dim>(np_1);
943  }
944  }
945 
965  static void
966  do_mass(const unsigned int n_components,
967  const AlignedVector<Number2> &transformation_matrix,
968  const AlignedVector<Number> & coefficients,
969  const Number * values_in,
970  Number * scratch_data,
971  Number * values_out)
972  {
973  constexpr int next_dim = dim > 1 ? dim - 1 : dim;
974  Number * my_scratch =
975  basis_size_1 != basis_size_2 ? scratch_data : values_out;
976 
977  const unsigned int size_per_component = Utilities::pow(basis_size_2, dim);
978  Assert(coefficients.size() == size_per_component ||
979  coefficients.size() == n_components * size_per_component,
980  ExcDimensionMismatch(coefficients.size(), size_per_component));
981  const unsigned int stride =
982  coefficients.size() == size_per_component ? 0 : 1;
983 
984  for (unsigned int q = basis_size_1; q != 0; --q)
986  variant,
988  next_dim,
989  basis_size_1,
990  basis_size_2,
991  Number,
992  Number2>::do_forward(n_components,
993  transformation_matrix,
994  values_in +
995  (q - 1) *
996  Utilities::pow(basis_size_1, dim - 1),
997  my_scratch +
998  (q - 1) *
999  Utilities::pow(basis_size_2, dim - 1));
1000  EvaluatorTensorProduct<variant,
1001  dim,
1002  basis_size_1,
1003  basis_size_2,
1004  Number,
1005  Number2>
1006  eval_val(transformation_matrix);
1007  const unsigned int n_inner_blocks =
1008  (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
1009  const unsigned int n_blocks = Utilities::pow(basis_size_2, dim - 1);
1010  for (unsigned int ii = 0; ii < n_blocks; ii += n_inner_blocks)
1011  for (unsigned int c = 0; c < n_components; ++c)
1012  {
1013  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1014  eval_val.template values_one_line<dim - 1, true, false>(
1015  my_scratch + i, my_scratch + i);
1016  for (unsigned int q = 0; q < basis_size_2; ++q)
1017  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1018  my_scratch[i + q * n_blocks + c * size_per_component] *=
1019  coefficients[i + q * n_blocks +
1020  c * stride * size_per_component];
1021  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1022  eval_val.template values_one_line<dim - 1, false, false>(
1023  my_scratch + i, my_scratch + i);
1024  }
1025  for (unsigned int q = 0; q < basis_size_1; ++q)
1027  variant,
1029  next_dim,
1030  basis_size_1,
1031  basis_size_2,
1032  Number,
1033  Number2>::do_backward(n_components,
1034  transformation_matrix,
1035  false,
1036  my_scratch +
1037  q * Utilities::pow(basis_size_2, dim - 1),
1038  values_out +
1039  q * Utilities::pow(basis_size_1, dim - 1));
1040  }
1041  };
1042 
1043 
1044 
1057  template <int dim, int fe_degree, typename Number>
1059  {
1060  static void
1061  evaluate(const unsigned int n_components,
1062  const EvaluationFlags::EvaluationFlags evaluation_flag,
1063  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1064  const Number * values_dofs,
1065  Number * values_quad,
1066  Number * gradients_quad,
1067  Number * hessians_quad,
1068  Number * scratch_data);
1069 
1070  static void
1071  integrate(const unsigned int n_components,
1072  const EvaluationFlags::EvaluationFlags integration_flag,
1073  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1074  Number * values_dofs,
1075  Number * values_quad,
1076  Number * gradients_quad,
1077  Number * scratch_data,
1078  const bool add_into_values_array);
1079  };
1080 
1081 
1082 
1083  template <int dim, int fe_degree, typename Number>
1084  inline void
1086  const unsigned int n_components,
1087  const EvaluationFlags::EvaluationFlags evaluation_flag,
1088  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1089  const Number * values_dofs,
1090  Number * values_quad,
1091  Number * gradients_quad,
1092  Number * hessians_quad,
1093  Number *)
1094  {
1096  shape_info.data.front().shape_gradients_collocation_eo.size(),
1097  (fe_degree + 2) / 2 * (fe_degree + 1));
1098 
1100  dim,
1101  fe_degree + 1,
1102  fe_degree + 1,
1103  Number>
1104  eval(AlignedVector<Number>(),
1105  shape_info.data.front().shape_gradients_collocation_eo,
1106  shape_info.data.front().shape_hessians_collocation_eo);
1107  constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
1108 
1109  for (unsigned int c = 0; c < n_components; c++)
1110  {
1111  if (evaluation_flag & EvaluationFlags::values)
1112  for (unsigned int i = 0; i < n_q_points; ++i)
1113  values_quad[i] = values_dofs[i];
1114  if (evaluation_flag &
1116  {
1117  eval.template gradients<0, true, false>(values_dofs,
1118  gradients_quad);
1119  if (dim > 1)
1120  eval.template gradients<1, true, false>(values_dofs,
1121  gradients_quad +
1122  n_q_points);
1123  if (dim > 2)
1124  eval.template gradients<2, true, false>(values_dofs,
1125  gradients_quad +
1126  2 * n_q_points);
1127  }
1128  if (evaluation_flag & EvaluationFlags::hessians)
1129  {
1130  eval.template hessians<0, true, false>(values_dofs, hessians_quad);
1131  if (dim > 1)
1132  {
1133  eval.template gradients<1, true, false>(gradients_quad,
1134  hessians_quad +
1135  dim * n_q_points);
1136  eval.template hessians<1, true, false>(values_dofs,
1137  hessians_quad +
1138  n_q_points);
1139  }
1140  if (dim > 2)
1141  {
1142  eval.template gradients<2, true, false>(gradients_quad,
1143  hessians_quad +
1144  4 * n_q_points);
1145  eval.template gradients<2, true, false>(
1146  gradients_quad + n_q_points, hessians_quad + 5 * n_q_points);
1147  eval.template hessians<2, true, false>(values_dofs,
1148  hessians_quad +
1149  2 * n_q_points);
1150  }
1151  hessians_quad += (dim * (dim + 1)) / 2 * n_q_points;
1152  }
1153  gradients_quad += dim * n_q_points;
1154  values_quad += n_q_points;
1155  values_dofs += n_q_points;
1156  }
1157  }
1158 
1159 
1160 
1161  template <int dim, int fe_degree, typename Number>
1162  inline void
1164  const unsigned int n_components,
1165  const EvaluationFlags::EvaluationFlags integration_flag,
1166  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1167  Number * values_dofs,
1168  Number * values_quad,
1169  Number * gradients_quad,
1170  Number *,
1171  const bool add_into_values_array)
1172  {
1174  shape_info.data.front().shape_gradients_collocation_eo.size(),
1175  (fe_degree + 2) / 2 * (fe_degree + 1));
1176 
1178  dim,
1179  fe_degree + 1,
1180  fe_degree + 1,
1181  Number>
1182  eval(AlignedVector<Number>(),
1183  shape_info.data.front().shape_gradients_collocation_eo,
1184  shape_info.data.front().shape_hessians_collocation_eo);
1185  constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
1186 
1187  for (unsigned int c = 0; c < n_components; c++)
1188  {
1189  if (integration_flag & EvaluationFlags::values)
1190  {
1191  if (add_into_values_array == false)
1192  for (unsigned int i = 0; i < n_q_points; ++i)
1193  values_dofs[i] = values_quad[i];
1194  else
1195  for (unsigned int i = 0; i < n_q_points; ++i)
1196  values_dofs[i] += values_quad[i];
1197  }
1198  if (integration_flag & EvaluationFlags::gradients)
1199  {
1200  if (integration_flag & EvaluationFlags::values ||
1201  add_into_values_array == true)
1202  eval.template gradients<0, false, true>(gradients_quad,
1203  values_dofs);
1204  else
1205  eval.template gradients<0, false, false>(gradients_quad,
1206  values_dofs);
1207  if (dim > 1)
1208  eval.template gradients<1, false, true>(gradients_quad +
1209  n_q_points,
1210  values_dofs);
1211  if (dim > 2)
1212  eval.template gradients<2, false, true>(gradients_quad +
1213  2 * n_q_points,
1214  values_dofs);
1215  }
1216  gradients_quad += dim * n_q_points;
1217  values_quad += n_q_points;
1218  values_dofs += n_q_points;
1219  }
1220  }
1221 
1222 
1223 
1234  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1236  {
1237  static void
1238  evaluate(const unsigned int n_components,
1239  const EvaluationFlags::EvaluationFlags evaluation_flag,
1240  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1241  const Number * values_dofs,
1242  Number * values_quad,
1243  Number * gradients_quad,
1244  Number * hessians_quad,
1245  Number * scratch_data);
1246 
1247  static void
1248  integrate(const unsigned int n_components,
1249  const EvaluationFlags::EvaluationFlags evaluation_flag,
1250  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1251  Number * values_dofs,
1252  Number * values_quad,
1253  Number * gradients_quad,
1254  Number * scratch_data,
1255  const bool add_into_values_array);
1256  };
1257 
1258 
1259 
1260  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1261  inline void
1263  dim,
1264  fe_degree,
1265  n_q_points_1d,
1266  Number>::evaluate(const unsigned int n_components,
1267  const EvaluationFlags::EvaluationFlags evaluation_flag,
1268  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1269  const Number * values_dofs,
1270  Number * values_quad,
1271  Number *gradients_quad,
1272  Number *hessians_quad,
1273  Number *)
1274  {
1275  Assert(n_q_points_1d > fe_degree,
1276  ExcMessage("You lose information when going to a collocation space "
1277  "of lower degree, so the evaluation results would be "
1278  "wrong. Thus, this class does not permit the desired "
1279  "operation."));
1280  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
1281 
1282  for (unsigned int c = 0; c < n_components; c++)
1283  {
1287  dim,
1288  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1289  n_q_points_1d,
1290  Number,
1291  Number>::do_forward(1,
1292  shape_info.data.front().shape_values_eo,
1293  values_dofs,
1294  values_quad);
1295 
1296  // apply derivatives in the collocation space
1297  if (evaluation_flag &
1300  1,
1301  evaluation_flag &
1303  shape_info,
1304  values_quad,
1305  nullptr,
1306  gradients_quad,
1307  hessians_quad,
1308  nullptr);
1309 
1310  values_dofs += shape_info.dofs_per_component_on_cell;
1311  values_quad += n_q_points;
1312  gradients_quad += dim * n_q_points;
1313  hessians_quad += (dim * (dim + 1)) / 2 * n_q_points;
1314  }
1315  }
1316 
1317 
1318 
1319  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1320  inline void
1322  dim,
1323  fe_degree,
1324  n_q_points_1d,
1325  Number>::integrate(const unsigned int n_components,
1326  const EvaluationFlags::EvaluationFlags integration_flag,
1327  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1328  Number *values_dofs,
1329  Number *values_quad,
1330  Number *gradients_quad,
1331  Number *,
1332  const bool add_into_values_array)
1333  {
1334  Assert(n_q_points_1d > fe_degree,
1335  ExcMessage("You lose information when going to a collocation space "
1336  "of lower degree, so the evaluation results would be "
1337  "wrong. Thus, this class does not permit the desired "
1338  "operation."));
1340  shape_info.data.front().shape_gradients_collocation_eo.size(),
1341  (n_q_points_1d + 1) / 2 * n_q_points_1d);
1342  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
1343 
1344  for (unsigned int c = 0; c < n_components; c++)
1345  {
1346  // apply derivatives in collocation space
1347  if (integration_flag & EvaluationFlags::gradients)
1349  integrate(1,
1350  integration_flag & EvaluationFlags::gradients,
1351  shape_info,
1352  values_quad,
1353  nullptr,
1354  gradients_quad,
1355  nullptr,
1356  /*add_into_values_array=*/integration_flag &
1358 
1359  // transform back to the original space
1363  dim,
1364  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1365  n_q_points_1d,
1366  Number,
1367  Number>::do_backward(1,
1368  shape_info.data.front().shape_values_eo,
1369  add_into_values_array,
1370  values_quad,
1371  values_dofs);
1372  gradients_quad += dim * n_q_points;
1373  values_quad += n_q_points;
1374  values_dofs += shape_info.dofs_per_component_on_cell;
1375  }
1376  }
1377 
1378 
1379 
1395  template <int dim, typename Number>
1397  {
1398  template <int fe_degree, int n_q_points_1d>
1399  static bool
1400  run(const unsigned int n_components,
1401  const EvaluationFlags::EvaluationFlags evaluation_flag,
1403  Number *values_dofs_actual,
1404  Number *values_quad,
1405  Number *gradients_quad,
1406  Number *hessians_quad,
1407  Number *scratch_data)
1408  {
1409  // We enable a transformation to collocation for derivatives if it gives
1410  // correct results (first condition), if it is the most efficient choice
1411  // in terms of operation counts (second condition) and if we were able to
1412  // initialize the fields in shape_info.templates.h from the polynomials
1413  // (third condition).
1414  static constexpr bool use_collocation =
1415  n_q_points_1d > fe_degree && n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
1416  n_q_points_1d < 200;
1417 
1418  if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
1419  shape_info.element_type ==
1421  {
1423  evaluate(n_components,
1424  evaluation_flag,
1425  shape_info,
1426  values_dofs_actual,
1427  values_quad,
1428  gradients_quad,
1429  hessians_quad,
1430  scratch_data);
1431  }
1432  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
1433  // shape_info.h for more details
1434  else if (fe_degree >= 0 && use_collocation &&
1435  shape_info.element_type <=
1437  {
1439  dim,
1440  fe_degree,
1441  n_q_points_1d,
1442  Number>::evaluate(n_components,
1443  evaluation_flag,
1444  shape_info,
1445  values_dofs_actual,
1446  values_quad,
1447  gradients_quad,
1448  hessians_quad,
1449  scratch_data);
1450  }
1451  else if (fe_degree >= 0 &&
1452  shape_info.element_type <=
1454  {
1457  dim,
1458  fe_degree,
1459  n_q_points_1d,
1460  Number>::evaluate(n_components,
1461  evaluation_flag,
1462  shape_info,
1463  values_dofs_actual,
1464  values_quad,
1465  gradients_quad,
1466  hessians_quad,
1467  scratch_data);
1468  }
1469  else if (shape_info.element_type ==
1471  {
1474  dim,
1475  fe_degree,
1476  n_q_points_1d,
1477  Number>::evaluate(n_components,
1478  evaluation_flag,
1479  shape_info,
1480  values_dofs_actual,
1481  values_quad,
1482  gradients_quad,
1483  hessians_quad,
1484  scratch_data);
1485  }
1486  else if (shape_info.element_type ==
1488  {
1491  dim,
1492  fe_degree,
1493  n_q_points_1d,
1494  Number>::evaluate(n_components,
1495  evaluation_flag,
1496  shape_info,
1497  values_dofs_actual,
1498  values_quad,
1499  gradients_quad,
1500  hessians_quad,
1501  scratch_data);
1502  }
1503  else
1504  {
1507  dim,
1508  fe_degree,
1509  n_q_points_1d,
1510  Number>::evaluate(n_components,
1511  evaluation_flag,
1512  shape_info,
1513  values_dofs_actual,
1514  values_quad,
1515  gradients_quad,
1516  hessians_quad,
1517  scratch_data);
1518  }
1519 
1520  return false;
1521  }
1522  };
1523 
1524 
1525 
1541  template <int dim, typename Number>
1543  {
1544  template <int fe_degree, int n_q_points_1d>
1545  static bool
1546  run(const unsigned int n_components,
1547  const EvaluationFlags::EvaluationFlags integration_flag,
1549  Number * values_dofs_actual,
1550  Number * values_quad,
1551  Number * gradients_quad,
1552  Number * scratch_data,
1553  const bool sum_into_values_array)
1554  {
1555  // We enable a transformation to collocation for derivatives if it gives
1556  // correct results (first condition), if it is the most efficient choice
1557  // in terms of operation counts (second condition) and if we were able to
1558  // initialize the fields in shape_info.templates.h from the polynomials
1559  // (third condition).
1560  constexpr bool use_collocation = n_q_points_1d > fe_degree &&
1561  n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
1562  n_q_points_1d < 200;
1563 
1564  if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
1565  shape_info.element_type ==
1567  {
1569  integrate(n_components,
1570  integration_flag,
1571  shape_info,
1572  values_dofs_actual,
1573  values_quad,
1574  gradients_quad,
1575  scratch_data,
1576  sum_into_values_array);
1577  }
1578  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
1579  // shape_info.h for more details
1580  else if (fe_degree >= 0 && use_collocation &&
1581  shape_info.element_type <=
1583  {
1585  dim,
1586  fe_degree,
1587  n_q_points_1d,
1588  Number>::integrate(n_components,
1589  integration_flag,
1590  shape_info,
1591  values_dofs_actual,
1592  values_quad,
1593  gradients_quad,
1594  scratch_data,
1595  sum_into_values_array);
1596  }
1597  else if (fe_degree >= 0 &&
1598  shape_info.element_type <=
1600  {
1603  dim,
1604  fe_degree,
1605  n_q_points_1d,
1606  Number>::integrate(n_components,
1607  integration_flag,
1608  shape_info,
1609  values_dofs_actual,
1610  values_quad,
1611  gradients_quad,
1612  scratch_data,
1613  sum_into_values_array);
1614  }
1615  else if (shape_info.element_type ==
1617  {
1620  dim,
1621  fe_degree,
1622  n_q_points_1d,
1623  Number>::integrate(n_components,
1624  integration_flag,
1625  shape_info,
1626  values_dofs_actual,
1627  values_quad,
1628  gradients_quad,
1629  scratch_data,
1630  sum_into_values_array);
1631  }
1632  else if (shape_info.element_type ==
1634  {
1637  dim,
1638  fe_degree,
1639  n_q_points_1d,
1640  Number>::integrate(n_components,
1641  integration_flag,
1642  shape_info,
1643  values_dofs_actual,
1644  values_quad,
1645  gradients_quad,
1646  scratch_data,
1647  sum_into_values_array);
1648  }
1649  else
1650  {
1653  dim,
1654  fe_degree,
1655  n_q_points_1d,
1656  Number>::integrate(n_components,
1657  integration_flag,
1658  shape_info,
1659  values_dofs_actual,
1660  values_quad,
1661  gradients_quad,
1662  scratch_data,
1663  sum_into_values_array);
1664  }
1665 
1666  return false;
1667  }
1668  };
1669 
1670 
1671 
1672  template <bool symmetric_evaluate,
1673  int dim,
1674  int fe_degree,
1675  int n_q_points_1d,
1676  typename Number>
1678  {
1679  // We enable a transformation to collocation for derivatives if it gives
1680  // correct results (first two conditions), if it is the most efficient
1681  // choice in terms of operation counts (third condition) and if we were
1682  // able to initialize the fields in shape_info.templates.h from the
1683  // polynomials (fourth condition).
1684  static constexpr bool use_collocation =
1685  symmetric_evaluate &&
1686  n_q_points_1d > fe_degree &&n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
1687  n_q_points_1d < 200;
1688 
1689  static void
1690  evaluate_in_face(const unsigned int n_components,
1692  Number * values_dofs,
1693  Number * values_quad,
1694  Number * gradients_quad,
1695  Number * scratch_data,
1696  const bool evaluate_val,
1697  const bool evaluate_grad,
1698  const unsigned int subface_index)
1699  {
1700  const AlignedVector<Number> &val1 =
1701  symmetric_evaluate ?
1702  data.data.front().shape_values_eo :
1703  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1704  data.data.front().shape_values :
1705  data.data.front().values_within_subface[subface_index % 2]);
1706  const AlignedVector<Number> &val2 =
1707  symmetric_evaluate ?
1708  data.data.front().shape_values_eo :
1709  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1710  data.data.front().shape_values :
1711  data.data.front().values_within_subface[subface_index / 2]);
1712 
1713  const AlignedVector<Number> &grad1 =
1714  symmetric_evaluate ?
1715  data.data.front().shape_gradients_eo :
1716  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1717  data.data.front().shape_gradients :
1718  data.data.front().gradients_within_subface[subface_index % 2]);
1719  const AlignedVector<Number> &grad2 =
1720  symmetric_evaluate ?
1721  data.data.front().shape_gradients_eo :
1722  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1723  data.data.front().shape_gradients :
1724  data.data.front().gradients_within_subface[subface_index / 2]);
1725 
1726  using Eval =
1727  internal::EvaluatorTensorProduct<symmetric_evaluate ?
1730  dim - 1,
1731  fe_degree + 1,
1732  n_q_points_1d,
1733  Number>;
1734  Eval eval1(val1,
1735  grad1,
1737  data.data.front().fe_degree + 1,
1738  data.data.front().n_q_points_1d);
1739  Eval eval2(val2,
1740  grad2,
1742  data.data.front().fe_degree + 1,
1743  data.data.front().n_q_points_1d);
1744 
1745  const unsigned int size_deg =
1746  fe_degree > -1 ?
1747  Utilities::pow(fe_degree + 1, dim - 1) :
1748  (dim > 1 ?
1749  Utilities::fixed_power<dim - 1>(data.data.front().fe_degree + 1) :
1750  1);
1751 
1752  const unsigned int n_q_points = fe_degree > -1 ?
1753  Utilities::pow(n_q_points_1d, dim - 1) :
1754  data.n_q_points_face;
1755 
1756  if (evaluate_grad == false)
1757  for (unsigned int c = 0; c < n_components; ++c)
1758  {
1759  switch (dim)
1760  {
1761  case 3:
1762  eval1.template values<0, true, false>(values_dofs,
1763  values_quad);
1764  eval2.template values<1, true, false>(values_quad,
1765  values_quad);
1766  break;
1767  case 2:
1768  eval1.template values<0, true, false>(values_dofs,
1769  values_quad);
1770  break;
1771  case 1:
1772  values_quad[0] = values_dofs[0];
1773  break;
1774  default:
1775  Assert(false, ExcNotImplemented());
1776  }
1777  values_dofs += 2 * size_deg;
1778  values_quad += n_q_points;
1779  }
1780  else
1781  for (unsigned int c = 0; c < n_components; ++c)
1782  {
1783  switch (dim)
1784  {
1785  case 3:
1786  if (use_collocation)
1787  {
1788  eval1.template values<0, true, false>(values_dofs,
1789  values_quad);
1790  eval1.template values<1, true, false>(values_quad,
1791  values_quad);
1794  dim - 1,
1795  n_q_points_1d,
1796  n_q_points_1d,
1797  Number>
1798  eval_grad(
1800  data.data.front().shape_gradients_collocation_eo,
1802  eval_grad.template gradients<0, true, false>(
1803  values_quad, gradients_quad);
1804  eval_grad.template gradients<1, true, false>(
1805  values_quad, gradients_quad + n_q_points);
1806  }
1807  else
1808  {
1809  eval1.template gradients<0, true, false>(values_dofs,
1810  scratch_data);
1811  eval2.template values<1, true, false>(scratch_data,
1812  gradients_quad);
1813 
1814  eval1.template values<0, true, false>(values_dofs,
1815  scratch_data);
1816  eval2.template gradients<1, true, false>(scratch_data,
1817  gradients_quad +
1818  n_q_points);
1819 
1820  if (evaluate_val == true)
1821  eval2.template values<1, true, false>(scratch_data,
1822  values_quad);
1823  }
1824  eval1.template values<0, true, false>(values_dofs + size_deg,
1825  scratch_data);
1826  eval2.template values<1, true, false>(
1827  scratch_data, gradients_quad + (dim - 1) * n_q_points);
1828 
1829  break;
1830  case 2:
1831  eval1.template values<0, true, false>(values_dofs + size_deg,
1832  gradients_quad +
1833  (dim - 1) *
1834  n_q_points);
1835  eval1.template gradients<0, true, false>(values_dofs,
1836  gradients_quad);
1837  if (evaluate_val == true)
1838  eval1.template values<0, true, false>(values_dofs,
1839  values_quad);
1840  break;
1841  case 1:
1842  values_quad[0] = values_dofs[0];
1843  gradients_quad[0] = values_dofs[1];
1844  break;
1845  default:
1846  AssertThrow(false, ExcNotImplemented());
1847  }
1848  values_dofs += 2 * size_deg;
1849  values_quad += n_q_points;
1850  gradients_quad += dim * n_q_points;
1851  }
1852  }
1853 
1854  static void
1855  integrate_in_face(const unsigned int n_components,
1857  Number * values_dofs,
1858  Number * values_quad,
1859  Number * gradients_quad,
1860  Number * scratch_data,
1861  const bool integrate_val,
1862  const bool integrate_grad,
1863  const unsigned int subface_index)
1864  {
1865  const AlignedVector<Number> &val1 =
1866  symmetric_evaluate ?
1867  data.data.front().shape_values_eo :
1868  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1869  data.data.front().shape_values :
1870  data.data.front().values_within_subface[subface_index % 2]);
1871  const AlignedVector<Number> &val2 =
1872  symmetric_evaluate ?
1873  data.data.front().shape_values_eo :
1874  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1875  data.data.front().shape_values :
1876  data.data.front().values_within_subface[subface_index / 2]);
1877 
1878  const AlignedVector<Number> &grad1 =
1879  symmetric_evaluate ?
1880  data.data.front().shape_gradients_eo :
1881  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1882  data.data.front().shape_gradients :
1883  data.data.front().gradients_within_subface[subface_index % 2]);
1884  const AlignedVector<Number> &grad2 =
1885  symmetric_evaluate ?
1886  data.data.front().shape_gradients_eo :
1887  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1888  data.data.front().shape_gradients :
1889  data.data.front().gradients_within_subface[subface_index / 2]);
1890 
1891  using Eval =
1892  internal::EvaluatorTensorProduct<symmetric_evaluate ?
1895  dim - 1,
1896  fe_degree + 1,
1897  n_q_points_1d,
1898  Number>;
1899  Eval eval1(val1,
1900  grad1,
1901  val1,
1902  data.data.front().fe_degree + 1,
1903  data.data.front().n_q_points_1d);
1904  Eval eval2(val2,
1905  grad2,
1906  val1,
1907  data.data.front().fe_degree + 1,
1908  data.data.front().n_q_points_1d);
1909 
1910  const unsigned int size_deg =
1911  fe_degree > -1 ?
1912  Utilities::pow(fe_degree + 1, dim - 1) :
1913  (dim > 1 ?
1914  Utilities::fixed_power<dim - 1>(data.data.front().fe_degree + 1) :
1915  1);
1916 
1917  const unsigned int n_q_points = fe_degree > -1 ?
1918  Utilities::pow(n_q_points_1d, dim - 1) :
1919  data.n_q_points_face;
1920 
1921  if (integrate_grad == false)
1922  for (unsigned int c = 0; c < n_components; ++c)
1923  {
1924  switch (dim)
1925  {
1926  case 3:
1927  eval2.template values<1, false, false>(values_quad,
1928  values_quad);
1929  eval1.template values<0, false, false>(values_quad,
1930  values_dofs);
1931  break;
1932  case 2:
1933  eval1.template values<0, false, false>(values_quad,
1934  values_dofs);
1935  break;
1936  case 1:
1937  values_dofs[0] = values_quad[0];
1938  break;
1939  default:
1940  Assert(false, ExcNotImplemented());
1941  }
1942  values_dofs += 2 * size_deg;
1943  values_quad += n_q_points;
1944  }
1945  else
1946  for (unsigned int c = 0; c < n_components; ++c)
1947  {
1948  switch (dim)
1949  {
1950  case 3:
1951  eval2.template values<1, false, false>(gradients_quad +
1952  2 * n_q_points,
1953  gradients_quad +
1954  2 * n_q_points);
1955  eval1.template values<0, false, false>(
1956  gradients_quad + 2 * n_q_points, values_dofs + size_deg);
1957  if (use_collocation)
1958  {
1961  dim - 1,
1962  n_q_points_1d,
1963  n_q_points_1d,
1964  Number>
1965  eval_grad(
1967  data.data.front().shape_gradients_collocation_eo,
1969  if (integrate_val)
1970  eval_grad.template gradients<1, false, true>(
1971  gradients_quad + n_q_points, values_quad);
1972  else
1973  eval_grad.template gradients<1, false, false>(
1974  gradients_quad + n_q_points, values_quad);
1975  eval_grad.template gradients<0, false, true>(
1976  gradients_quad, values_quad);
1977  eval1.template values<1, false, false>(values_quad,
1978  values_quad);
1979  eval1.template values<0, false, false>(values_quad,
1980  values_dofs);
1981  }
1982  else
1983  {
1984  if (integrate_val)
1985  {
1986  eval2.template values<1, false, false>(values_quad,
1987  scratch_data);
1988  eval2.template gradients<1, false, true>(
1989  gradients_quad + n_q_points, scratch_data);
1990  }
1991  else
1992  eval2.template gradients<1, false, false>(
1993  gradients_quad + n_q_points, scratch_data);
1994 
1995  eval1.template values<0, false, false>(scratch_data,
1996  values_dofs);
1997  eval2.template values<1, false, false>(gradients_quad,
1998  scratch_data);
1999  eval1.template gradients<0, false, true>(scratch_data,
2000  values_dofs);
2001  }
2002  break;
2003  case 2:
2004  eval1.template values<0, false, false>(
2005  gradients_quad + n_q_points, values_dofs + size_deg);
2006  eval1.template gradients<0, false, false>(gradients_quad,
2007  values_dofs);
2008  if (integrate_val == true)
2009  eval1.template values<0, false, true>(values_quad,
2010  values_dofs);
2011  break;
2012  case 1:
2013  values_dofs[0] = values_quad[0];
2014  values_dofs[1] = gradients_quad[0];
2015  break;
2016  default:
2017  AssertThrow(false, ExcNotImplemented());
2018  }
2019  values_dofs += 2 * size_deg;
2020  values_quad += n_q_points;
2021  gradients_quad += dim * n_q_points;
2022  }
2023  }
2024  };
2025 
2026 
2027 
2028  template <int dim, int fe_degree, typename Number, bool lex_faces = false>
2030  {
2031  template <bool do_evaluate, bool add_into_output>
2032  static void
2033  interpolate(const unsigned int n_components,
2035  const Number * input,
2036  Number * output,
2037  const bool do_gradients,
2038  const unsigned int face_no)
2039  {
2040  Assert(static_cast<unsigned int>(fe_degree) ==
2041  data.data.front().fe_degree ||
2042  fe_degree == -1,
2043  ExcInternalError());
2044 
2045  interpolate_generic<do_evaluate, add_into_output>(
2046  n_components,
2047  input,
2048  output,
2049  do_gradients,
2050  face_no,
2051  data.data.front().fe_degree + 1,
2052  data.data.front().shape_data_on_face,
2054  2 * data.dofs_per_component_on_face);
2055  }
2056 
2060  template <bool do_evaluate, bool add_into_output>
2061  static void
2062  interpolate_quadrature(const unsigned int n_components,
2064  const Number * input,
2065  Number * output,
2066  const bool do_gradients,
2067  const unsigned int face_no)
2068  {
2069  Assert(static_cast<unsigned int>(fe_degree + 1) ==
2070  data.data.front().quadrature.size() ||
2071  fe_degree == -1,
2072  ExcInternalError());
2073 
2074  interpolate_generic<do_evaluate, add_into_output>(
2075  n_components,
2076  input,
2077  output,
2078  do_gradients,
2079  face_no,
2080  data.data.front().quadrature.size(),
2081  data.data.front().quadrature_data_on_face,
2082  data.n_q_points,
2083  data.n_q_points_face);
2084  }
2085 
2086  private:
2087  template <bool do_evaluate, bool add_into_output, int face_direction = 0>
2088  static void
2089  interpolate_generic(const unsigned int n_components,
2090  const Number * input,
2091  Number * output,
2092  const bool do_gradients,
2093  const unsigned int face_no,
2094  const unsigned int n_points_1d,
2095  const std::array<AlignedVector<Number>, 2> &shape_data,
2096  const unsigned int dofs_per_component_on_cell,
2097  const unsigned int dofs_per_component_on_face)
2098  {
2099  if (face_direction == face_no / 2)
2100  {
2102  dim,
2103  fe_degree + 1,
2104  0,
2105  Number>
2106  evalf(shape_data[face_no % 2],
2109  n_points_1d,
2110  0);
2111 
2112  const unsigned int in_stride = do_evaluate ?
2113  dofs_per_component_on_cell :
2114  dofs_per_component_on_face;
2115  const unsigned int out_stride = do_evaluate ?
2116  dofs_per_component_on_face :
2117  dofs_per_component_on_cell;
2118 
2119  for (unsigned int c = 0; c < n_components; c++)
2120  {
2121  if (do_gradients)
2122  evalf.template apply_face<face_direction,
2123  do_evaluate,
2124  add_into_output,
2125  1,
2126  lex_faces>(input, output);
2127  else
2128  evalf.template apply_face<face_direction,
2129  do_evaluate,
2130  add_into_output,
2131  0,
2132  lex_faces>(input, output);
2133  input += in_stride;
2134  output += out_stride;
2135  }
2136  }
2137  else if (face_direction < dim)
2138  {
2139  interpolate_generic<do_evaluate,
2140  add_into_output,
2141  std::min(face_direction + 1, dim - 1)>(
2142  n_components,
2143  input,
2144  output,
2145  do_gradients,
2146  face_no,
2147  n_points_1d,
2148  shape_data,
2149  dofs_per_component_on_cell,
2150  dofs_per_component_on_face);
2151  }
2152  }
2153  };
2154 
2155 
2156 
2157  // internal helper function for reading data; base version of different types
2158  template <typename VectorizedArrayType, typename Number2>
2159  void
2160  do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
2161  {
2162  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2163  dst[v] = src_ptr[v];
2164  }
2165 
2166 
2167 
2168  // internal helper function for reading data; specialized version where we
2169  // can use a dedicated load function
2170  template <typename Number, unsigned int width>
2171  void
2173  {
2174  dst.load(src_ptr);
2175  }
2176 
2177 
2178 
2179  // internal helper function for reading data; base version of different types
2180  template <typename VectorizedArrayType, typename Number2>
2181  void
2182  do_vectorized_gather(const Number2 * src_ptr,
2183  const unsigned int * indices,
2184  VectorizedArrayType &dst)
2185  {
2186  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2187  dst[v] = src_ptr[indices[v]];
2188  }
2189 
2190 
2191 
2192  // internal helper function for reading data; specialized version where we
2193  // can use a dedicated gather function
2194  template <typename Number, unsigned int width>
2195  void
2196  do_vectorized_gather(const Number * src_ptr,
2197  const unsigned int * indices,
2199  {
2200  dst.gather(src_ptr, indices);
2201  }
2202 
2203 
2204 
2205  // internal helper function for reading data; base version of different types
2206  template <typename VectorizedArrayType, typename Number2>
2207  void
2208  do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
2209  {
2210  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2211  dst_ptr[v] += src[v];
2212  }
2213 
2214 
2215 
2216  // internal helper function for reading data; specialized version where we
2217  // can use a dedicated load function
2218  template <typename Number, unsigned int width>
2219  void
2221  {
2223  tmp.load(dst_ptr);
2224  (tmp + src).store(dst_ptr);
2225  }
2226 
2227 
2228 
2229  // internal helper function for reading data; base version of different types
2230  template <typename VectorizedArrayType, typename Number2>
2231  void
2232  do_vectorized_scatter_add(const VectorizedArrayType src,
2233  const unsigned int * indices,
2234  Number2 * dst_ptr)
2235  {
2236  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2237  dst_ptr[indices[v]] += src[v];
2238  }
2239 
2240 
2241 
2242  // internal helper function for reading data; specialized version where we
2243  // can use a dedicated gather function
2244  template <typename Number, unsigned int width>
2245  void
2247  const unsigned int * indices,
2248  Number * dst_ptr)
2249  {
2250 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
2251  for (unsigned int v = 0; v < width; ++v)
2252  dst_ptr[indices[v]] += src[v];
2253 #else
2255  tmp.gather(dst_ptr, indices);
2256  (tmp + src).scatter(indices, dst_ptr);
2257 #endif
2258  }
2259 
2260 
2261 
2262  template <typename Number>
2263  void
2264  adjust_for_face_orientation(const unsigned int dim,
2265  const unsigned int n_components,
2266  const unsigned int face_orientation,
2267  const Table<2, unsigned int> &orientation_map,
2268  const bool integrate,
2269  const bool values,
2270  const bool gradients,
2271  const unsigned int n_q_points,
2272  Number * tmp_values,
2273  Number * values_quad,
2274  Number * gradients_quad)
2275  {
2276  Assert(face_orientation, ExcInternalError());
2277  const unsigned int *orientation = &orientation_map[face_orientation][0];
2278  for (unsigned int c = 0; c < n_components; ++c)
2279  {
2280  if (values == true)
2281  {
2282  if (integrate)
2283  for (unsigned int q = 0; q < n_q_points; ++q)
2284  tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
2285  else
2286  for (unsigned int q = 0; q < n_q_points; ++q)
2287  tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
2288  for (unsigned int q = 0; q < n_q_points; ++q)
2289  values_quad[c * n_q_points + q] = tmp_values[q];
2290  }
2291  if (gradients == true)
2292  for (unsigned int d = 0; d < dim; ++d)
2293  {
2294  if (integrate)
2295  for (unsigned int q = 0; q < n_q_points; ++q)
2296  tmp_values[q] =
2297  gradients_quad[(c * dim + d) * n_q_points + orientation[q]];
2298  else
2299  for (unsigned int q = 0; q < n_q_points; ++q)
2300  tmp_values[orientation[q]] =
2301  gradients_quad[(c * dim + d) * n_q_points + q];
2302  for (unsigned int q = 0; q < n_q_points; ++q)
2303  gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q];
2304  }
2305  }
2306  }
2307 
2308 
2309 
2310  template <int dim, typename VectorizedArrayType>
2312  {
2313  template <int fe_degree, int n_q_points_1d>
2314  static bool
2315  run(const unsigned int n_components,
2317  const VectorizedArrayType * values_array,
2318  VectorizedArrayType * values_quad,
2319  VectorizedArrayType * gradients_quad,
2320  VectorizedArrayType * scratch_data,
2321  const bool evaluate_values,
2322  const bool evaluate_gradients,
2323  const unsigned int face_no,
2324  const unsigned int subface_index,
2325  const unsigned int face_orientation,
2326  const Table<2, unsigned int> &orientation_map)
2327  {
2328  constexpr unsigned int static_dofs_per_face =
2329  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
2331  const unsigned int dofs_per_face =
2332  fe_degree > -1 ?
2333  static_dofs_per_face :
2334  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
2335 
2336  VectorizedArrayType *temp1 = scratch_data;
2337 
2339  template interpolate<true, false>(
2340  n_components, data, values_array, temp1, evaluate_gradients, face_no);
2341 
2342  const unsigned int n_q_points_1d_actual =
2343  fe_degree > -1 ? n_q_points_1d : 0;
2344  if (fe_degree > -1 &&
2345  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
2346  data.element_type <= MatrixFreeFunctions::tensor_symmetric)
2348  true,
2349  dim,
2350  fe_degree,
2351  n_q_points_1d_actual,
2352  VectorizedArrayType>::evaluate_in_face(n_components,
2353  data,
2354  temp1,
2355  values_quad,
2356  gradients_quad,
2357  scratch_data + 2 *
2358  n_components *
2359  dofs_per_face,
2360  evaluate_values,
2361  evaluate_gradients,
2362  subface_index);
2363  else
2365  false,
2366  dim,
2367  fe_degree,
2368  n_q_points_1d_actual,
2369  VectorizedArrayType>::evaluate_in_face(n_components,
2370  data,
2371  temp1,
2372  values_quad,
2373  gradients_quad,
2374  scratch_data + 2 *
2375  n_components *
2376  dofs_per_face,
2377  evaluate_values,
2378  evaluate_gradients,
2379  subface_index);
2380 
2381  if (face_orientation)
2383  n_components,
2384  face_orientation,
2385  orientation_map,
2386  false,
2387  evaluate_values,
2388  evaluate_gradients,
2389  data.n_q_points_face,
2390  scratch_data,
2391  values_quad,
2392  gradients_quad);
2393 
2394  return false;
2395  }
2396  };
2397 
2398 
2399 
2400  template <int dim, typename VectorizedArrayType>
2402  {
2403  template <int fe_degree, int n_q_points_1d>
2404  static bool
2405  run(const unsigned int n_components,
2407  VectorizedArrayType * values_array,
2408  VectorizedArrayType * values_quad,
2409  VectorizedArrayType * gradients_quad,
2410  VectorizedArrayType * scratch_data,
2411  const bool integrate_values,
2412  const bool integrate_gradients,
2413  const unsigned int face_no,
2414  const unsigned int subface_index,
2415  const unsigned int face_orientation,
2416  const Table<2, unsigned int> &orientation_map)
2417  {
2418  if (face_orientation)
2420  n_components,
2421  face_orientation,
2422  orientation_map,
2423  true,
2424  integrate_values,
2425  integrate_gradients,
2426  data.n_q_points_face,
2427  scratch_data,
2428  values_quad,
2429  gradients_quad);
2430 
2431  constexpr unsigned int static_dofs_per_face =
2432  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
2434  const unsigned int dofs_per_face =
2435  fe_degree > -1 ?
2436  static_dofs_per_face :
2437  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
2438 
2439  VectorizedArrayType *temp1 = scratch_data;
2440 
2441  const unsigned int n_q_points_1d_actual =
2442  fe_degree > -1 ? n_q_points_1d : 0;
2443  if (fe_degree > -1 &&
2447  true,
2448  dim,
2449  fe_degree,
2450  n_q_points_1d_actual,
2451  VectorizedArrayType>::integrate_in_face(n_components,
2452  data,
2453  temp1,
2454  values_quad,
2455  gradients_quad,
2456  scratch_data +
2457  2 * n_components *
2458  dofs_per_face,
2459  integrate_values,
2460  integrate_gradients,
2461  subface_index);
2462  else
2464  false,
2465  dim,
2466  fe_degree,
2467  n_q_points_1d_actual,
2468  VectorizedArrayType>::integrate_in_face(n_components,
2469  data,
2470  temp1,
2471  values_quad,
2472  gradients_quad,
2473  scratch_data +
2474  2 * n_components *
2475  dofs_per_face,
2476  integrate_values,
2477  integrate_gradients,
2478  subface_index);
2479 
2481  template interpolate<false, false>(n_components,
2482  data,
2483  temp1,
2484  values_array,
2485  integrate_gradients,
2486  face_no);
2487  return false;
2488  }
2489  };
2490 
2491 
2492 
2493  template <int dim,
2494  int fe_degree,
2495  int n_q_points_1d,
2496  typename Number,
2497  typename VectorizedArrayType,
2498  std::size_t n_face_orientations,
2499  typename Number2_,
2500  typename Function1a,
2501  typename Function1b,
2502  typename Function2a,
2503  typename Function2b,
2504  typename Function3a,
2505  typename Function3b,
2506  typename Function5,
2507  typename Function0>
2508  static bool
2510  const unsigned int n_components,
2511  const bool integrate,
2512  Number2_ * global_vector_ptr,
2514  const MatrixFreeFunctions::DoFInfo & dof_info,
2515  VectorizedArrayType * values_quad,
2516  VectorizedArrayType * gradients_quad,
2517  VectorizedArrayType * scratch_data,
2518  const bool do_values,
2519  const bool do_gradients,
2520  const unsigned int active_fe_index,
2521  const unsigned int first_selected_component,
2522  const std::array<unsigned int, n_face_orientations> cells,
2523  const std::array<unsigned int, n_face_orientations> face_nos,
2524  const unsigned int subface_index,
2525  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
2526  const std::array<unsigned int, n_face_orientations> face_orientations,
2527  const Table<2, unsigned int> & orientation_map,
2528  const Function1a & function_1a,
2529  const Function1b & function_1b,
2530  const Function2a & function_2a,
2531  const Function2b & function_2b,
2532  const Function3a & function_3a,
2533  const Function3b & function_3b,
2534  const Function5 & function_5,
2535  const Function0 & function_0)
2536  {
2537  const unsigned int cell = cells[0];
2538 
2539  // In the case of integration, we do not need to reshuffle the
2540  // data at the quadrature points to adjust for the face
2541  // orientation if the shape functions are nodal at the cell
2542  // boundaries (and we only requested the integration of the
2543  // values) or Hermite shape functions are used. These cases are
2544  // handled later when the values are written back into the
2545  // glrobal vector.
2546  if (integrate &&
2547  (face_orientations[0] > 0 &&
2548  (subface_index < GeometryInfo<dim>::max_children_per_cell ||
2549  !(((do_gradients == false &&
2550  data.data.front().nodal_at_cell_boundaries == true) ||
2551  (data.element_type ==
2553  fe_degree > 1)) &&
2554  (dof_info.index_storage_variants[dof_access_index][cell] ==
2556  interleaved_contiguous ||
2557  dof_info.index_storage_variants[dof_access_index][cell] ==
2559  interleaved_contiguous_strided ||
2560  dof_info.index_storage_variants[dof_access_index][cell] ==
2562  interleaved_contiguous_mixed_strides ||
2563  dof_info.index_storage_variants[dof_access_index][cell] ==
2565  contiguous)))))
2566  {
2567  AssertDimension(face_orientations.size(), 1);
2569  n_components,
2570  face_orientations[0],
2571  orientation_map,
2572  true,
2573  do_values,
2574  do_gradients,
2575  data.n_q_points_face,
2576  scratch_data,
2577  values_quad,
2578  gradients_quad);
2579  }
2580 
2581  // we know that the gradient weights for the Hermite case on the
2582  // right (side==1) are the negative from the value at the left
2583  // (side==0), so we only read out one of them.
2584  VectorizedArrayType grad_weight =
2585  (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
2587  data.data.front()
2588  .shape_data_on_face[0][fe_degree + (integrate ?
2589  (2 - (face_nos[0] % 2)) :
2590  (1 + (face_nos[0] % 2)))] :
2591  VectorizedArrayType(0.0 /*dummy*/);
2592 
2593  constexpr unsigned int static_dofs_per_component =
2594  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) :
2596  constexpr unsigned int static_dofs_per_face =
2597  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
2599  const unsigned int dofs_per_face =
2600  fe_degree > -1 ? static_dofs_per_face :
2601  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
2602 
2603  // we allocate small amounts of data on the stack to signal the compiler
2604  // that this temporary data is only needed for the calculations but the
2605  // final results can be discarded and need not be written back to
2606  // memory. For large sizes or when the dofs per face is not a
2607  // compile-time constant, however, we want to go to the heap in the
2608  // `scratch_data` variable to not risk a stack overflow.
2609  constexpr unsigned int stack_array_size_threshold = 100;
2610 
2611  VectorizedArrayType *DEAL_II_RESTRICT temp1;
2612  VectorizedArrayType
2613  temp_data[static_dofs_per_face < stack_array_size_threshold ?
2614  2 * dofs_per_face :
2615  1];
2616  if (static_dofs_per_face < stack_array_size_threshold)
2617  temp1 = &temp_data[0];
2618  else
2619  temp1 = scratch_data;
2620 
2621  const unsigned int dummy = 0;
2622 
2623  // re-orientation
2624  std::array<const unsigned int *, n_face_orientations> orientation;
2625  if (n_face_orientations == 1)
2626  orientation[0] = (data.data.front().nodal_at_cell_boundaries == true) ?
2627  &data.face_orientations[face_orientations[0]][0] :
2628  &dummy;
2629  else
2630  {
2631  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2632  {
2633  // the loop breaks once an invalid_unsigned_int is hit for
2634  // all cases except the exterior faces in the ECL loop (where
2635  // some faces might be at the boundaries but others not)
2636  if (cells[v] == numbers::invalid_unsigned_int)
2637  continue;
2638 
2639  orientation[v] =
2640  (data.data.front().nodal_at_cell_boundaries == true) ?
2641  &data.face_orientations[face_orientations[v]][0] :
2642  &dummy;
2643  }
2644  }
2645 
2646  // face_to_cell_index_hermite
2647  std::array<const unsigned int *, n_face_orientations> index_array_hermite;
2648 
2649  if (n_face_orientations == 1)
2650  index_array_hermite[0] =
2651  (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
2652  data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ?
2653  &data.face_to_cell_index_hermite(face_nos[0], 0) :
2654  &dummy;
2655 
2656  if (n_face_orientations > 1 &&
2657  data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
2658  data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite)
2659  {
2660  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2661  {
2662  if (cells[v] == numbers::invalid_unsigned_int)
2663  continue;
2664 
2665  grad_weight[v] =
2666  data.data.front().shape_data_on_face
2667  [0][fe_degree + (integrate ? (2 - (face_nos[v] % 2)) :
2668  (1 + (face_nos[v] % 2)))][v];
2669 
2670  index_array_hermite[v] =
2671  &data.face_to_cell_index_hermite(face_nos[v], 0);
2672  }
2673  }
2674 
2675  // face_to_cell_index_nodal
2676  std::array<const unsigned int *, n_face_orientations> index_array_nodal;
2677 
2678  if (n_face_orientations == 1)
2679  index_array_nodal[0] =
2680  (data.data.front().nodal_at_cell_boundaries == true) ?
2681  &data.face_to_cell_index_nodal(face_nos[0], 0) :
2682  &dummy;
2683 
2684  if (n_face_orientations > 1 &&
2685  (data.data.front().nodal_at_cell_boundaries == true))
2686  {
2687  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2688  {
2689  if (cells[v] == numbers::invalid_unsigned_int)
2690  continue;
2691 
2692  index_array_nodal[v] =
2693  &data.face_to_cell_index_nodal(face_nos[v], 0);
2694  }
2695  }
2696 
2697  const auto reorientate = [&](const unsigned int v, const unsigned int i) {
2698  return (dim < 3 ||
2699  face_orientations[n_face_orientations == 1 ? 0 : v] == 0 ||
2700  subface_index < GeometryInfo<dim>::max_children_per_cell) ?
2701  i :
2702  orientation[v][i];
2703  };
2704 
2705  // this variable keeps track of whether we are able to directly write
2706  // the results into the result (function returns true) or not, requiring
2707  // an additional call to another function
2708  bool accesses_global_vector = true;
2709 
2710  for (unsigned int comp = 0; comp < n_components; ++comp)
2711  {
2712  if (integrate)
2713  function_0(temp1, comp);
2714  if ((do_gradients == false &&
2715  data.data.front().nodal_at_cell_boundaries == true) ||
2716  (data.element_type ==
2718  fe_degree > 1))
2719  {
2720  // case 1: contiguous and interleaved indices
2721  if (n_face_orientations == 1 &&
2722  dof_info.index_storage_variants[dof_access_index][cell] ==
2724  interleaved_contiguous)
2725  {
2726  AssertDimension(n_face_orientations, 1);
2727 
2729  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2730  VectorizedArrayType::size());
2731  Number2_ *vector_ptr =
2732  global_vector_ptr +
2733  dof_info.dof_indices_contiguous[dof_access_index]
2734  [cell *
2735  VectorizedArrayType::size()] +
2736  (dof_info
2737  .component_dof_indices_offset[active_fe_index]
2738  [first_selected_component] +
2739  comp * static_dofs_per_component) *
2740  VectorizedArrayType::size();
2741 
2742  if (fe_degree > 1 && do_gradients == true)
2743  {
2744  for (unsigned int i = 0; i < dofs_per_face; ++i)
2745  {
2746  if (n_face_orientations == 1)
2747  {
2748  const unsigned int ind1 =
2749  index_array_hermite[0][2 * i];
2750  const unsigned int ind2 =
2751  index_array_hermite[0][2 * i + 1];
2752  AssertIndexRange(ind1,
2753  data.dofs_per_component_on_cell);
2754  AssertIndexRange(ind2,
2755  data.dofs_per_component_on_cell);
2756  const unsigned int i_ = reorientate(0, i);
2757  function_1a(temp1[i_],
2758  temp1[i_ + dofs_per_face],
2759  vector_ptr +
2760  ind1 * VectorizedArrayType::size(),
2761  vector_ptr +
2762  ind2 * VectorizedArrayType::size(),
2763  grad_weight);
2764  }
2765  else
2766  {
2767  Assert(false, ExcNotImplemented());
2768  }
2769  }
2770  }
2771  else
2772  {
2773  for (unsigned int i = 0; i < dofs_per_face; ++i)
2774  {
2775  if (n_face_orientations == 1)
2776  {
2777  const unsigned int i_ = reorientate(0, i);
2778  const unsigned int ind = index_array_nodal[0][i];
2779  function_1b(temp1[i_],
2780  vector_ptr +
2781  ind * VectorizedArrayType::size());
2782  }
2783  else
2784  {
2785  Assert(false, ExcNotImplemented());
2786  }
2787  }
2788  }
2789  }
2790 
2791  // case 2: contiguous and interleaved indices with fixed stride
2792  else if (n_face_orientations == 1 &&
2793  dof_info.index_storage_variants[dof_access_index][cell] ==
2795  interleaved_contiguous_strided)
2796  {
2797  AssertDimension(n_face_orientations, 1);
2798 
2800  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2801  VectorizedArrayType::size());
2802  const unsigned int *indices =
2803  &dof_info.dof_indices_contiguous[dof_access_index]
2804  [cell *
2805  VectorizedArrayType::size()];
2806  Number2_ *vector_ptr =
2807  global_vector_ptr +
2808  (comp * static_dofs_per_component +
2809  dof_info
2810  .component_dof_indices_offset[active_fe_index]
2811  [first_selected_component]) *
2812  VectorizedArrayType::size();
2813  if (fe_degree > 1 && do_gradients == true)
2814  {
2815  for (unsigned int i = 0; i < dofs_per_face; ++i)
2816  {
2817  if (n_face_orientations == 1)
2818  {
2819  const unsigned int i_ = reorientate(0, i);
2820  const unsigned int ind1 =
2821  index_array_hermite[0][2 * i] *
2822  VectorizedArrayType::size();
2823  const unsigned int ind2 =
2824  index_array_hermite[0][2 * i + 1] *
2825  VectorizedArrayType::size();
2826  function_2a(temp1[i_],
2827  temp1[i_ + dofs_per_face],
2828  vector_ptr + ind1,
2829  vector_ptr + ind2,
2830  grad_weight,
2831  indices,
2832  indices);
2833  }
2834  else
2835  {
2836  Assert(false, ExcNotImplemented());
2837  }
2838  }
2839  }
2840  else
2841  {
2842  for (unsigned int i = 0; i < dofs_per_face; ++i)
2843  {
2844  if (n_face_orientations == 1)
2845  {
2846  const unsigned int i_ = reorientate(0, i);
2847  const unsigned int ind =
2848  index_array_nodal[0][i] *
2849  VectorizedArrayType::size();
2850  function_2b(temp1[i_], vector_ptr + ind, indices);
2851  }
2852  else
2853  {
2854  Assert(false, ExcNotImplemented());
2855  }
2856  }
2857  }
2858  }
2859 
2860  // case 3: contiguous and interleaved indices with mixed stride
2861  else if (n_face_orientations == 1 &&
2862  dof_info.index_storage_variants[dof_access_index][cell] ==
2864  interleaved_contiguous_mixed_strides)
2865  {
2866  AssertDimension(n_face_orientations, 1);
2867 
2868  const unsigned int *strides =
2869  &dof_info.dof_indices_interleave_strides
2870  [dof_access_index][cell * VectorizedArrayType::size()];
2871  unsigned int indices[VectorizedArrayType::size()];
2872  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2873  indices[v] =
2874  dof_info.dof_indices_contiguous
2875  [dof_access_index]
2876  [cell * VectorizedArrayType::size() + v] +
2877  (dof_info
2878  .component_dof_indices_offset[active_fe_index]
2879  [first_selected_component] +
2880  comp * static_dofs_per_component) *
2881  strides[v];
2882  const unsigned int n_filled_lanes =
2883  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
2884 
2885  if (fe_degree > 1 && do_gradients == true)
2886  {
2887  if (n_filled_lanes == VectorizedArrayType::size())
2888  for (unsigned int i = 0; i < dofs_per_face; ++i)
2889  {
2890  if (n_face_orientations == 1)
2891  {
2892  const unsigned int i_ = reorientate(0, i);
2893  unsigned int ind1[VectorizedArrayType::size()];
2895  for (unsigned int v = 0;
2896  v < VectorizedArrayType::size();
2897  ++v)
2898  ind1[v] =
2899  indices[v] +
2900  index_array_hermite[0 /*TODO*/][2 * i] *
2901  strides[v];
2902  unsigned int ind2[VectorizedArrayType::size()];
2904  for (unsigned int v = 0;
2905  v < VectorizedArrayType::size();
2906  ++v)
2907  ind2[v] =
2908  indices[v] +
2909  index_array_hermite[0 /*TODO*/][2 * i + 1] *
2910  strides[v];
2911  function_2a(temp1[i_],
2912  temp1[i_ + dofs_per_face],
2913  global_vector_ptr,
2914  global_vector_ptr,
2915  grad_weight,
2916  ind1,
2917  ind2);
2918  }
2919  else
2920  {
2921  Assert(false, ExcNotImplemented());
2922  }
2923  }
2924  else
2925  {
2926  if (integrate == false)
2927  for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
2928  temp1[i] = VectorizedArrayType();
2929 
2930  for (unsigned int v = 0; v < n_filled_lanes; ++v)
2931  for (unsigned int i = 0; i < dofs_per_face; ++i)
2932  {
2933  const unsigned int i_ =
2934  reorientate(n_face_orientations == 1 ? 0 : v,
2935  i);
2936  function_3a(
2937  temp1[i_][v],
2938  temp1[i_ + dofs_per_face][v],
2939  global_vector_ptr
2940  [indices[v] +
2941  index_array_hermite
2942  [n_face_orientations == 1 ? 0 : v]
2943  [2 * i] *
2944  strides[v]],
2945  global_vector_ptr
2946  [indices[v] +
2947  index_array_hermite
2948  [n_face_orientations == 1 ? 0 : v]
2949  [2 * i + 1] *
2950  strides[v]],
2951  grad_weight[n_face_orientations == 1 ? 0 : v]);
2952  }
2953  }
2954  }
2955  else
2956  {
2957  if (n_filled_lanes == VectorizedArrayType::size())
2958  for (unsigned int i = 0; i < dofs_per_face; ++i)
2959  {
2960  if (n_face_orientations == 1)
2961  {
2962  unsigned int ind[VectorizedArrayType::size()];
2964  for (unsigned int v = 0;
2965  v < VectorizedArrayType::size();
2966  ++v)
2967  ind[v] = indices[v] +
2968  index_array_nodal[0][i] * strides[v];
2969  const unsigned int i_ = reorientate(0, i);
2970  function_2b(temp1[i_], global_vector_ptr, ind);
2971  }
2972  else
2973  {
2974  Assert(false, ExcNotImplemented());
2975  }
2976  }
2977  else
2978  {
2979  if (integrate == false)
2980  for (unsigned int i = 0; i < dofs_per_face; ++i)
2981  temp1[i] = VectorizedArrayType();
2982 
2983  for (unsigned int v = 0; v < n_filled_lanes; ++v)
2984  for (unsigned int i = 0; i < dofs_per_face; ++i)
2985  function_3b(
2986  temp1[reorientate(
2987  n_face_orientations == 1 ? 0 : v, i)][v],
2988  global_vector_ptr
2989  [indices[v] +
2990  index_array_nodal
2991  [n_face_orientations == 1 ? 0 : v][i] *
2992  strides[v]]);
2993  }
2994  }
2995  }
2996 
2997  // case 4: contiguous indices without interleaving
2998  else if (n_face_orientations > 1 ||
2999  dof_info.index_storage_variants[dof_access_index][cell] ==
3001  contiguous)
3002  {
3003  const unsigned int *indices =
3004  &dof_info.dof_indices_contiguous[dof_access_index]
3005  [cell *
3006  VectorizedArrayType::size()];
3007  Number2_ *vector_ptr =
3008  global_vector_ptr + comp * static_dofs_per_component +
3009  dof_info
3010  .component_dof_indices_offset[active_fe_index]
3011  [first_selected_component];
3012 
3013  if (do_gradients == true &&
3014  data.element_type ==
3016  {
3017  for (unsigned int i = 0; i < dofs_per_face; ++i)
3018  {
3019  if (n_face_orientations == 1 &&
3020  dof_info
3021  .n_vectorization_lanes_filled[dof_access_index]
3022  [cell] ==
3023  VectorizedArrayType::size())
3024  {
3025  const unsigned int ind1 =
3026  index_array_hermite[0][2 * i];
3027  const unsigned int ind2 =
3028  index_array_hermite[0][2 * i + 1];
3029  const unsigned int i_ = reorientate(0, i);
3030 
3031  function_2a(temp1[i_],
3032  temp1[i_ + dofs_per_face],
3033  vector_ptr + ind1,
3034  vector_ptr + ind2,
3035  grad_weight,
3036  indices,
3037  indices);
3038  }
3039  else if (n_face_orientations == 1)
3040  {
3041  const unsigned int ind1 =
3042  index_array_hermite[0][2 * i];
3043  const unsigned int ind2 =
3044  index_array_hermite[0][2 * i + 1];
3045  const unsigned int i_ = reorientate(0, i);
3046 
3047  const unsigned int n_filled_lanes =
3048  dof_info
3049  .n_vectorization_lanes_filled[dof_access_index]
3050  [cell];
3051 
3052  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3053  function_3a(temp1[i_][v],
3054  temp1[i_ + dofs_per_face][v],
3055  vector_ptr[ind1 + indices[v]],
3056  vector_ptr[ind2 + indices[v]],
3057  grad_weight[v]);
3058 
3059  if (integrate == false)
3060  for (unsigned int v = n_filled_lanes;
3061  v < VectorizedArrayType::size();
3062  ++v)
3063  {
3064  temp1[i_][v] = 0.0;
3065  temp1[i_ + dofs_per_face][v] = 0.0;
3066  }
3067  }
3068  else
3069  {
3070  Assert(false, ExcNotImplemented());
3071 
3072  const unsigned int n_filled_lanes =
3073  dof_info
3074  .n_vectorization_lanes_filled[dof_access_index]
3075  [cell];
3076 
3077  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3078  function_3a(
3079  temp1[reorientate(v, i)][v],
3080  temp1[reorientate(v, i) + dofs_per_face][v],
3081  vector_ptr[index_array_hermite[v][2 * i] +
3082  indices[v]],
3083  vector_ptr[index_array_hermite[v][2 * i + 1] +
3084  indices[v]],
3085  grad_weight[v]);
3086  }
3087  }
3088  }
3089  else
3090  {
3091  for (unsigned int i = 0; i < dofs_per_face; ++i)
3092  {
3093  if (n_face_orientations == 1 &&
3094  dof_info
3095  .n_vectorization_lanes_filled[dof_access_index]
3096  [cell] ==
3097  VectorizedArrayType::size())
3098  {
3099  const unsigned int ind = index_array_nodal[0][i];
3100  const unsigned int i_ = reorientate(0, i);
3101 
3102  function_2b(temp1[i_], vector_ptr + ind, indices);
3103  }
3104  else if (n_face_orientations == 1)
3105  {
3106  const unsigned int ind = index_array_nodal[0][i];
3107  const unsigned int i_ = reorientate(0, i);
3108 
3109  const unsigned int n_filled_lanes =
3110  dof_info
3111  .n_vectorization_lanes_filled[dof_access_index]
3112  [cell];
3113 
3114  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3115  function_3b(temp1[i_][v],
3116  vector_ptr[ind + indices[v]]);
3117 
3118  if (integrate == false)
3119  for (unsigned int v = n_filled_lanes;
3120  v < VectorizedArrayType::size();
3121  ++v)
3122  temp1[i_][v] = 0.0;
3123  }
3124  else
3125  {
3126  for (unsigned int v = 0;
3127  v < VectorizedArrayType::size();
3128  ++v)
3129  if (cells[v] != numbers::invalid_unsigned_int)
3130  function_3b(
3131  temp1[reorientate(v, i)][v],
3132  vector_ptr[index_array_nodal[v][i] +
3133  dof_info.dof_indices_contiguous
3134  [dof_access_index][cells[v]]]);
3135  }
3136  }
3137  }
3138  }
3139  else
3140  {
3141  // case 5: default vector access
3142 
3143  // for the integrate_scatter path (integrate == true), we
3144  // need to only prepare the data in this function for all
3145  // components to later call distribute_local_to_global();
3146  // for the gather_evaluate path (integrate == false), we
3147  // instead want to leave early because we need to get the
3148  // vector data from somewhere else
3149  function_5(temp1, comp);
3150  if (integrate)
3151  accesses_global_vector = false;
3152  else
3153  return false;
3154  }
3155  }
3156  else
3157  {
3158  // case 5: default vector access
3159  function_5(temp1, comp);
3160  if (integrate)
3161  accesses_global_vector = false;
3162  else
3163  return false;
3164  }
3165 
3166  if (!integrate)
3167  function_0(temp1, comp);
3168  }
3169 
3170  if (!integrate &&
3171  (face_orientations[0] > 0 &&
3173  {
3174  AssertDimension(face_orientations.size(), 1);
3176  n_components,
3177  face_orientations[0],
3178  orientation_map,
3179  false,
3180  do_values,
3181  do_gradients,
3182  data.n_q_points_face,
3183  scratch_data,
3184  values_quad,
3185  gradients_quad);
3186  }
3187 
3188  return accesses_global_vector;
3189  }
3190 
3191 
3192  template <int dim,
3193  typename Number,
3194  typename VectorizedArrayType,
3195  typename Number2 = Number>
3197  {
3198  template <int fe_degree, int n_q_points_1d, std::size_t n_face_orientations>
3199  static bool
3200  run(const unsigned int n_components,
3201  const Number2 * src_ptr,
3203  const MatrixFreeFunctions::DoFInfo & dof_info,
3204  VectorizedArrayType * values_quad,
3205  VectorizedArrayType *gradients_quad,
3206  VectorizedArrayType *scratch_data,
3207  const bool evaluate_values,
3208  const bool evaluate_gradients,
3209  const unsigned int active_fe_index,
3210  const unsigned int first_selected_component,
3211  const std::array<unsigned int, n_face_orientations> cells,
3212  const std::array<unsigned int, n_face_orientations> face_nos,
3213  const unsigned int subface_index,
3214  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
3215  const std::array<unsigned int, n_face_orientations> face_orientations,
3216  const Table<2, unsigned int> & orientation_map)
3217  {
3218  if (src_ptr == nullptr)
3219  {
3220  return false;
3221  }
3222 
3224  fe_degree,
3225  n_q_points_1d,
3226  Number>( //
3227  n_components,
3228  false /*=evaluate*/,
3229  src_ptr,
3230  data,
3231  dof_info,
3232  values_quad,
3233  gradients_quad,
3234  scratch_data,
3235  evaluate_values,
3236  evaluate_gradients,
3237  active_fe_index,
3238  first_selected_component,
3239  cells,
3240  face_nos,
3241  subface_index,
3242  dof_access_index,
3243  face_orientations,
3244  orientation_map,
3245  [](VectorizedArrayType & temp_1,
3246  VectorizedArrayType & temp_2,
3247  const Number * src_ptr_1,
3248  const Number * src_ptr_2,
3249  const VectorizedArrayType &grad_weight) {
3250  // case 1a)
3251  do_vectorized_read(src_ptr_1, temp_1);
3252  do_vectorized_read(src_ptr_2, temp_2);
3253  temp_2 = grad_weight * (temp_1 - temp_2);
3254  },
3255  [](VectorizedArrayType &temp, const Number *src_ptr) {
3256  // case 1b)
3257  do_vectorized_read(src_ptr, temp);
3258  },
3259  [](VectorizedArrayType & temp_1,
3260  VectorizedArrayType & temp_2,
3261  const Number * src_ptr_1,
3262  const Number * src_ptr_2,
3263  const VectorizedArrayType &grad_weight,
3264  const unsigned int * indices_1,
3265  const unsigned int * indices_2) {
3266  // case 2a)
3267  do_vectorized_gather(src_ptr_1, indices_1, temp_1);
3268  do_vectorized_gather(src_ptr_2, indices_2, temp_2);
3269  temp_2 = grad_weight * (temp_1 - temp_2);
3270  },
3271  [](VectorizedArrayType &temp,
3272  const Number * src_ptr,
3273  const unsigned int * indices) {
3274  // case 2b)
3275  do_vectorized_gather(src_ptr, indices, temp);
3276  },
3277  [](Number & temp_1,
3278  Number & temp_2,
3279  const Number src_ptr_1,
3280  const Number src_ptr_2,
3281  Number grad_weight) {
3282  // case 3a)
3283  temp_1 = src_ptr_1;
3284  temp_2 = grad_weight * (temp_1 - src_ptr_2);
3285  },
3286  [](Number &temp, const Number src_ptr) {
3287  // case 3b)
3288  temp = src_ptr;
3289  },
3290  [&](const VectorizedArrayType *, const unsigned int) {
3291  // case 5)
3292  },
3293  [&](auto &temp1, const unsigned int comp) {
3294  const unsigned int dofs_per_face =
3295  fe_degree > -1 ?
3296  Utilities::pow(fe_degree + 1, dim - 1) :
3297  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
3298  const unsigned int n_q_points =
3299  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
3300  data.n_q_points_face;
3301  if (fe_degree > -1 &&
3302  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
3303  data.element_type <= MatrixFreeFunctions::tensor_symmetric)
3304  FEFaceEvaluationImpl<true,
3305  dim,
3306  fe_degree,
3307  n_q_points_1d,
3308  VectorizedArrayType>::
3309  evaluate_in_face(/* n_components */ 1,
3310  data,
3311  temp1,
3312  values_quad + comp * n_q_points,
3313  gradients_quad + comp * dim * n_q_points,
3314  scratch_data + 2 * dofs_per_face,
3315  evaluate_values,
3316  evaluate_gradients,
3317  subface_index);
3318  else
3319  FEFaceEvaluationImpl<false,
3320  dim,
3321  fe_degree,
3322  n_q_points_1d,
3323  VectorizedArrayType>::
3324  evaluate_in_face(/* n_components */ 1,
3325  data,
3326  temp1,
3327  values_quad + comp * n_q_points,
3328  gradients_quad + comp * dim * n_q_points,
3329  scratch_data + 2 * dofs_per_face,
3330  evaluate_values,
3331  evaluate_gradients,
3332  subface_index);
3333  });
3334  }
3335  };
3336 
3337  template <int dim,
3338  typename Number,
3339  typename VectorizedArrayType,
3340  typename Number2 = Number>
3342  {
3343  template <int fe_degree, int n_q_points_1d, std::size_t n_face_orientations>
3344  static bool
3345  run(const unsigned int n_components,
3346  Number2 * dst_ptr,
3348  const MatrixFreeFunctions::DoFInfo & dof_info,
3349  VectorizedArrayType * values_array,
3350  VectorizedArrayType * values_quad,
3351  VectorizedArrayType *gradients_quad,
3352  VectorizedArrayType *scratch_data,
3353  const bool integrate_values,
3354  const bool integrate_gradients,
3355  const unsigned int active_fe_index,
3356  const unsigned int first_selected_component,
3357  const std::array<unsigned int, n_face_orientations> cells,
3358  const std::array<unsigned int, n_face_orientations> face_nos,
3359  const unsigned int subface_index,
3360  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
3361  const std::array<unsigned int, n_face_orientations> face_orientations,
3362  const Table<2, unsigned int> & orientation_map)
3363  {
3364  if (dst_ptr == nullptr)
3365  {
3366  AssertDimension(face_nos.size(), 1);
3367  AssertDimension(face_orientations.size(), 1);
3368 
3369  // for block vectors simply integrate
3371  template run<fe_degree, n_q_points_1d>(n_components,
3372  data,
3373  values_array,
3374  values_quad,
3375  gradients_quad,
3376  scratch_data,
3377  integrate_values,
3378  integrate_gradients,
3379  face_nos[0],
3380  subface_index,
3381  face_orientations[0],
3382  orientation_map);
3383 
3384  // default vector access
3385  return false;
3386  }
3387 
3389  fe_degree,
3390  n_q_points_1d,
3391  Number>( //
3392  n_components,
3393  true /*=integrate*/,
3394  dst_ptr,
3395  data,
3396  dof_info,
3397  values_quad,
3398  gradients_quad,
3399  scratch_data,
3400  integrate_values,
3401  integrate_gradients,
3402  active_fe_index,
3403  first_selected_component,
3404  cells,
3405  face_nos,
3406  subface_index,
3407  dof_access_index,
3408  face_orientations,
3409  orientation_map,
3410  [](const VectorizedArrayType &temp_1,
3411  const VectorizedArrayType &temp_2,
3412  Number * dst_ptr_1,
3413  Number * dst_ptr_2,
3414  const VectorizedArrayType &grad_weight) {
3415  // case 1a)
3416  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
3417  const VectorizedArrayType grad = grad_weight * temp_2;
3418  do_vectorized_add(val, dst_ptr_1);
3419  do_vectorized_add(grad, dst_ptr_2);
3420  },
3421  [](const auto &temp, auto dst_ptr) {
3422  // case 1b)
3423  do_vectorized_add(temp, dst_ptr);
3424  },
3425  [](VectorizedArrayType & temp_1,
3426  VectorizedArrayType & temp_2,
3427  Number * dst_ptr_1,
3428  Number * dst_ptr_2,
3429  const VectorizedArrayType &grad_weight,
3430  const unsigned int * indices_1,
3431  const unsigned int * indices_2) {
3432  // case 2a)
3433  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
3434  const VectorizedArrayType grad = grad_weight * temp_2;
3435  do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
3436  do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
3437  },
3438  [](const VectorizedArrayType &temp,
3439  Number * dst_ptr,
3440  const unsigned int * indices) {
3441  // case 2b)
3442  do_vectorized_scatter_add(temp, indices, dst_ptr);
3443  },
3444  [](const Number temp_1,
3445  const Number temp_2,
3446  Number & dst_ptr_1,
3447  Number & dst_ptr_2,
3448  const Number grad_weight) {
3449  // case 3a)
3450  const Number val = temp_1 - grad_weight * temp_2;
3451  const Number grad = grad_weight * temp_2;
3452  dst_ptr_1 += val;
3453  dst_ptr_2 += grad;
3454  },
3455  [](const Number temp, Number &dst_ptr) {
3456  // case 3b)
3457  dst_ptr += temp;
3458  },
3459  [&](const VectorizedArrayType *temp1, const unsigned int comp) {
3460  // case 5: default vector access, must be handled separately, just do
3461  // the face-normal interpolation
3462 
3463  AssertDimension(face_nos.size(), 1);
3464 
3466  template interpolate<false, false>(
3467  /* n_components */ 1,
3468  data,
3469  temp1,
3470  values_array + comp * data.dofs_per_component_on_cell,
3471  integrate_gradients,
3472  face_nos[0]);
3473  },
3474  [&](auto &temp1, const unsigned int comp) {
3475  const unsigned int dofs_per_face =
3476  fe_degree > -1 ?
3477  Utilities::pow(fe_degree + 1, dim - 1) :
3478  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
3479  const unsigned int n_q_points =
3480  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
3481  data.n_q_points_face;
3482  if (fe_degree > -1 &&
3483  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
3484  data.element_type <=
3487  dim,
3488  fe_degree,
3489  n_q_points_1d,
3490  VectorizedArrayType>::
3491  integrate_in_face(/* n_components */ 1,
3492  data,
3493  temp1,
3494  values_quad + comp * n_q_points,
3495  gradients_quad + dim * comp * n_q_points,
3496  scratch_data + 2 * dofs_per_face,
3497  integrate_values,
3498  integrate_gradients,
3499  subface_index);
3500  else
3502  dim,
3503  fe_degree,
3504  n_q_points_1d,
3505  VectorizedArrayType>::
3506  integrate_in_face(/* n_components */ 1,
3507  data,
3508  temp1,
3509  values_quad + comp * n_q_points,
3510  gradients_quad + dim * comp * n_q_points,
3511  scratch_data + 2 * dofs_per_face,
3512  integrate_values,
3513  integrate_gradients,
3514  subface_index);
3515  });
3516  }
3517  };
3518 
3519 
3520 
3525  template <int dim, typename Number>
3527  {
3528  template <int fe_degree, int = 0>
3529  static bool
3530  run(const unsigned int n_components,
3531  const FEEvaluationBaseData<dim,
3532  typename Number::value_type,
3533  false,
3534  Number> &fe_eval,
3535  const Number * in_array,
3536  Number * out_array,
3537  typename std::enable_if<fe_degree != -1>::type * = nullptr)
3538  {
3539  constexpr unsigned int dofs_per_component =
3540  Utilities::pow(fe_degree + 1, dim);
3541 
3542  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
3543  Assert(fe_eval.get_shape_info().element_type <=
3545  ExcNotImplemented());
3546 
3548  dim,
3549  fe_degree + 1,
3550  fe_degree + 1,
3551  Number>
3552  evaluator(
3555  fe_eval.get_shape_info().data.front().inverse_shape_values_eo);
3556 
3557  for (unsigned int d = 0; d < n_components; ++d)
3558  {
3559  const Number *in = in_array + d * dofs_per_component;
3560  Number * out = out_array + d * dofs_per_component;
3561  // Need to select 'apply' method with hessian slot because values
3562  // assume symmetries that do not exist in the inverse shapes
3563  evaluator.template hessians<0, true, false>(in, out);
3564  if (dim > 1)
3565  evaluator.template hessians<1, true, false>(out, out);
3566  if (dim > 2)
3567  evaluator.template hessians<2, true, false>(out, out);
3568  }
3569  for (unsigned int q = 0; q < dofs_per_component; ++q)
3570  {
3571  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
3572  for (unsigned int d = 0; d < n_components; ++d)
3573  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
3574  }
3575  for (unsigned int d = 0; d < n_components; ++d)
3576  {
3577  Number *out = out_array + d * dofs_per_component;
3578  if (dim > 2)
3579  evaluator.template hessians<2, false, false>(out, out);
3580  if (dim > 1)
3581  evaluator.template hessians<1, false, false>(out, out);
3582  evaluator.template hessians<0, false, false>(out, out);
3583  }
3584  return false;
3585  }
3586 
3587  template <int fe_degree, int = 0>
3588  static bool
3589  run(const unsigned int n_components,
3590  const FEEvaluationBaseData<dim,
3591  typename Number::value_type,
3592  false,
3593  Number> &fe_eval,
3594  const Number * in_array,
3595  Number * out_array,
3596  typename std::enable_if<fe_degree == -1>::type * = nullptr)
3597  {
3598  static_assert(fe_degree == -1, "Only usable for degree -1");
3599  const unsigned int dofs_per_component =
3600  fe_eval.get_shape_info().dofs_per_component_on_cell;
3601 
3602  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
3603 
3604  internal::
3605  EvaluatorTensorProduct<internal::evaluate_general, dim, 0, 0, Number>
3606  evaluator(fe_eval.get_shape_info().data.front().inverse_shape_values,
3609  fe_eval.get_shape_info().data.front().fe_degree + 1,
3610  fe_eval.get_shape_info().data.front().fe_degree + 1);
3611 
3612  for (unsigned int d = 0; d < n_components; ++d)
3613  {
3614  const Number *in = in_array + d * dofs_per_component;
3615  Number * out = out_array + d * dofs_per_component;
3616  // Need to select 'apply' method with hessian slot because values
3617  // assume symmetries that do not exist in the inverse shapes
3618  evaluator.template values<0, true, false>(in, out);
3619  if (dim > 1)
3620  evaluator.template values<1, true, false>(out, out);
3621  if (dim > 2)
3622  evaluator.template values<2, true, false>(out, out);
3623  }
3624  for (unsigned int q = 0; q < dofs_per_component; ++q)
3625  {
3626  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
3627  for (unsigned int d = 0; d < n_components; ++d)
3628  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
3629  }
3630  for (unsigned int d = 0; d < n_components; ++d)
3631  {
3632  Number *out = out_array + d * dofs_per_component;
3633  if (dim > 2)
3634  evaluator.template values<2, false, false>(out, out);
3635  if (dim > 1)
3636  evaluator.template values<1, false, false>(out, out);
3637  evaluator.template values<0, false, false>(out, out);
3638  }
3639  return false;
3640  }
3641  };
3642 
3643 
3644 
3649  template <int dim, typename Number>
3651  {
3652  template <int fe_degree, int = 0>
3653  static bool
3654  run(const unsigned int n_desired_components,
3655  const AlignedVector<Number> &inverse_shape,
3656  const AlignedVector<Number> &inverse_coefficients,
3657  const Number * in_array,
3658  Number * out_array,
3659  typename std::enable_if<fe_degree != -1>::type * = nullptr)
3660  {
3661  constexpr unsigned int dofs_per_component =
3662  Utilities::pow(fe_degree + 1, dim);
3663  Assert(inverse_coefficients.size() > 0 &&
3664  inverse_coefficients.size() % dofs_per_component == 0,
3665  ExcMessage(
3666  "Expected diagonal to be a multiple of scalar dof per cells"));
3667  if (inverse_coefficients.size() != dofs_per_component)
3668  AssertDimension(n_desired_components * dofs_per_component,
3669  inverse_coefficients.size());
3670 
3671  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
3672 
3674  dim,
3675  fe_degree + 1,
3676  fe_degree + 1,
3677  Number>
3678  evaluator(AlignedVector<Number>(),
3680  inverse_shape);
3681 
3682  const unsigned int shift_coefficient =
3683  inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
3684  0;
3685  const Number *inv_coefficient = inverse_coefficients.data();
3686  for (unsigned int d = 0; d < n_desired_components; ++d)
3687  {
3688  const Number *in = in_array + d * dofs_per_component;
3689  Number * out = out_array + d * dofs_per_component;
3690  // Need to select 'apply' method with hessian slot because values
3691  // assume symmetries that do not exist in the inverse shapes
3692  evaluator.template hessians<0, true, false>(in, out);
3693  if (dim > 1)
3694  evaluator.template hessians<1, true, false>(out, out);
3695  if (dim > 2)
3696  evaluator.template hessians<2, true, false>(out, out);
3697 
3698  for (unsigned int q = 0; q < dofs_per_component; ++q)
3699  out[q] *= inv_coefficient[q];
3700 
3701  if (dim > 2)
3702  evaluator.template hessians<2, false, false>(out, out);
3703  if (dim > 1)
3704  evaluator.template hessians<1, false, false>(out, out);
3705  evaluator.template hessians<0, false, false>(out, out);
3706 
3707  inv_coefficient += shift_coefficient;
3708  }
3709  return false;
3710  }
3711 
3715  template <int fe_degree, int = 0>
3716  static bool
3717  run(const unsigned int,
3718  const AlignedVector<Number> &,
3719  const AlignedVector<Number> &,
3720  const Number *,
3721  Number *,
3722  typename std::enable_if<fe_degree == -1>::type * = nullptr)
3723  {
3724  static_assert(fe_degree == -1, "Only usable for degree -1");
3725  Assert(false, ExcNotImplemented());
3726  return false;
3727  }
3728  };
3729 
3730 
3731 
3736  template <int dim, typename Number>
3738  {
3739  template <int fe_degree, int = 0>
3740  static bool
3741  run(const unsigned int n_desired_components,
3742  const AlignedVector<Number> &inverse_shape,
3743  const Number * in_array,
3744  Number * out_array,
3745  typename std::enable_if<fe_degree != -1>::type * = nullptr)
3746  {
3747  constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim);
3749  dim,
3750  fe_degree + 1,
3751  fe_degree + 1,
3752  Number>
3753  evaluator(AlignedVector<Number>(),
3755  inverse_shape);
3756 
3757  for (unsigned int d = 0; d < n_desired_components; ++d)
3758  {
3759  const Number *in = in_array + d * dofs_per_cell;
3760  Number * out = out_array + d * dofs_per_cell;
3761 
3762  if (dim == 3)
3763  {
3764  evaluator.template hessians<2, false, false>(in, out);
3765  evaluator.template hessians<1, false, false>(out, out);
3766  evaluator.template hessians<0, false, false>(out, out);
3767  }
3768  if (dim == 2)
3769  {
3770  evaluator.template hessians<1, false, false>(in, out);
3771  evaluator.template hessians<0, false, false>(out, out);
3772  }
3773  if (dim == 1)
3774  evaluator.template hessians<0, false, false>(in, out);
3775  }
3776  return false;
3777  }
3778 
3779  template <int fe_degree, int = 0>
3780  static bool
3781  run(const unsigned int,
3782  const AlignedVector<Number> &,
3783  const Number *,
3784  Number *,
3785  typename std::enable_if<fe_degree == -1>::type * = nullptr)
3786  {
3787  static_assert(fe_degree == -1, "Only usable for degree -1");
3788  Assert(false, ExcNotImplemented());
3789  return false;
3790  }
3791  };
3792 
3793 } // end of namespace internal
3794 
3795 
3797 
3798 #endif
#define DEAL_II_RESTRICT
Definition: config.h:95
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 integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs_actual, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool add_into_values_array)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const internal::MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs_actual, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool sum_into_values_array)
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1568
void adjust_for_face_orientation(const unsigned int dim, const unsigned int n_components, const unsigned int face_orientation, const Table< 2, unsigned int > &orientation_map, const bool integrate, const bool values, const bool gradients, const unsigned int n_q_points, Number *tmp_values, Number *values_quad, Number *gradients_quad)
void do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
static void interpolate_quadrature(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< Number > &data, const Number *input, Number *output, const bool do_gradients, const unsigned int face_no)
static void integrate_in_face(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< Number > &data, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool integrate_val, const bool integrate_grad, const unsigned int subface_index)
void do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
static void interpolate_generic(const unsigned int n_components, const Number *input, Number *output, const bool do_gradients, const unsigned int face_no, const unsigned int n_points_1d, const std::array< AlignedVector< Number >, 2 > &shape_data, const unsigned int dofs_per_component_on_cell, const unsigned int dofs_per_component_on_face)
std::array< std::vector< IndexStorageVariants >, 3 > index_storage_variants
Definition: dof_info.h:436
#define AssertIndexRange(index, range)
Definition: exceptions.h:1636
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data)
pointer data()
static bool run(const unsigned int n_components, Number2 *dst_ptr, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, const MatrixFreeFunctions::DoFInfo &dof_info, VectorizedArrayType *values_array, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool integrate_values, const bool integrate_gradients, const unsigned int active_fe_index, const unsigned int first_selected_component, const std::array< unsigned int, n_face_orientations > cells, const std::array< unsigned int, n_face_orientations > face_nos, const unsigned int subface_index, const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, const std::array< unsigned int, n_face_orientations > face_orientations, const Table< 2, unsigned int > &orientation_map)
static void evaluate_in_face(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< Number > &data, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool evaluate_val, const bool evaluate_grad, const unsigned int subface_index)
static bool run(const unsigned int n_components, const FEEvaluationBaseData< dim, typename Number::value_type, false, Number > &fe_eval, const Number *in_array, Number *out_array, typename std::enable_if< fe_degree !=-1 >::type *=nullptr)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1521
static bool run(const unsigned int, const AlignedVector< Number > &, const Number *, Number *, typename std::enable_if< fe_degree==-1 >::type *=nullptr)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:461
T fixed_power(const T t)
Definition: utilities.h:1045
static bool run(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, const VectorizedArrayType *values_array, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool evaluate_values, const bool evaluate_gradients, const unsigned int face_no, const unsigned int subface_index, const unsigned int face_orientation, const Table< 2, unsigned int > &orientation_map)
static bool run(const unsigned int n_desired_components, const AlignedVector< Number > &inverse_shape, const AlignedVector< Number > &inverse_coefficients, const Number *in_array, Number *out_array, typename std::enable_if< fe_degree !=-1 >::type *=nullptr)
static bool run(const unsigned int n_components, const Number2 *src_ptr, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, const MatrixFreeFunctions::DoFInfo &dof_info, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool evaluate_values, const bool evaluate_gradients, const unsigned int active_fe_index, const unsigned int first_selected_component, const std::array< unsigned int, n_face_orientations > cells, const std::array< unsigned int, n_face_orientations > face_nos, const unsigned int subface_index, const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, const std::array< unsigned int, n_face_orientations > face_orientations, const Table< 2, unsigned int > &orientation_map)
static ::ExceptionBase & ExcMessage(std::string arg1)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const Number *values_dofs_actual, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data)
::Table< 2, unsigned int > face_orientations
Definition: shape_info.h:488
void gather(const Number *base_ptr, const unsigned int *offsets)
#define Assert(cond, exc)
Definition: exceptions.h:1411
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static bool run(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, VectorizedArrayType *values_array, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool integrate_values, const bool integrate_gradients, const unsigned int face_no, const unsigned int subface_index, const unsigned int face_orientation, const Table< 2, unsigned int > &orientation_map)
static bool run(const unsigned int n_components, const FEEvaluationBaseData< dim, typename Number::value_type, false, Number > &fe_eval, const Number *in_array, Number *out_array, typename std::enable_if< fe_degree==-1 >::type *=nullptr)
static void do_mass(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const AlignedVector< Number > &coefficients, const Number *values_in, Number *scratch_data, Number *values_out)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
void load(const Number *ptr)
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:94
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
static void do_forward(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const Number *values_in, Number *values_out, const unsigned int basis_size_1_variable=numbers::invalid_unsigned_int, const unsigned int basis_size_2_variable=numbers::invalid_unsigned_int)
static bool run(const unsigned int n_desired_components, const AlignedVector< Number > &inverse_shape, const Number *in_array, Number *out_array, typename std::enable_if< fe_degree !=-1 >::type *=nullptr)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static void interpolate(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< Number > &data, const Number *input, Number *output, const bool do_gradients, const unsigned int face_no)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
T min(const T &t, const MPI_Comm &mpi_communicator)
static bool run(const unsigned int, const AlignedVector< Number > &, const AlignedVector< Number > &, const Number *, Number *, typename std::enable_if< fe_degree==-1 >::type *=nullptr)
size_type size() const
static bool fe_face_evaluation_process_and_io(const unsigned int n_components, const bool integrate, Number2_ *global_vector_ptr, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, const MatrixFreeFunctions::DoFInfo &dof_info, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool do_values, const bool do_gradients, const unsigned int active_fe_index, const unsigned int first_selected_component, const std::array< unsigned int, n_face_orientations > cells, const std::array< unsigned int, n_face_orientations > face_nos, const unsigned int subface_index, const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, const std::array< unsigned int, n_face_orientations > face_orientations, const Table< 2, unsigned int > &orientation_map, const Function1a &function_1a, const Function1b &function_1b, const Function2a &function_2a, const Function2b &function_2b, const Function3a &function_3a, const Function3b &function_3b, const Function5 &function_5, const Function0 &function_0)
static ::ExceptionBase & ExcNotImplemented()
std::vector< UnivariateShapeData< Number > > data
Definition: shape_info.h:363
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const internal::MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs_actual, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data)
EvaluationFlags
The EvaluationFlags enum.
T max(const T &t, const MPI_Comm &mpi_communicator)
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:134
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool add_into_values_array)
static ::ExceptionBase & ExcInternalError()
void do_vectorized_scatter_add(const VectorizedArrayType src, const unsigned int *indices, Number2 *dst_ptr)
void do_vectorized_gather(const Number2 *src_ptr, const unsigned int *indices, VectorizedArrayType &dst)