Reference documentation for deal.II version 9.2.0
\(\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>
23 #include <deal.II/base/vectorization.h>
24 
25 #include <deal.II/matrix_free/dof_info.h>
26 #include <deal.II/matrix_free/shape_info.h>
27 #include <deal.II/matrix_free/tensor_product_kernels.h>
28 
29 
30 DEAL_II_NAMESPACE_OPEN
31 
32 
33 
34 namespace internal
35 {
36  // Select evaluator type from element shape function type
37  template <MatrixFreeFunctions::ElementType element, bool is_long>
38  struct EvaluatorSelector
39  {};
40 
41  template <bool is_long>
42  struct EvaluatorSelector<MatrixFreeFunctions::tensor_general, is_long>
43  {
44  static const EvaluatorVariant variant = evaluate_general;
45  };
46 
47  template <>
48  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, false>
49  {
50  static const EvaluatorVariant variant = evaluate_symmetric;
51  };
52 
53  template <>
54  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, true>
55  {
56  static const EvaluatorVariant variant = evaluate_evenodd;
57  };
58 
59  template <bool is_long>
60  struct EvaluatorSelector<MatrixFreeFunctions::truncated_tensor, is_long>
61  {
62  static const EvaluatorVariant variant = evaluate_general;
63  };
64 
65  template <>
66  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,
67  false>
68  {
69  static const EvaluatorVariant variant = evaluate_general;
70  };
71 
72  template <>
73  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0, true>
74  {
75  static const EvaluatorVariant variant = evaluate_evenodd;
76  };
77 
78  template <bool is_long>
79  struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_collocation,
80  is_long>
81  {
82  static const EvaluatorVariant variant = evaluate_evenodd;
83  };
84 
85 
86 
105  template <MatrixFreeFunctions::ElementType type,
106  int dim,
107  int fe_degree,
108  int n_q_points_1d,
109  int n_components,
110  typename Number>
111  struct FEEvaluationImpl
112  {
113  static void
114  evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
115  const Number * values_dofs_actual,
116  Number * values_quad,
117  Number * gradients_quad,
118  Number * hessians_quad,
119  Number * scratch_data,
120  const bool evaluate_values,
121  const bool evaluate_gradients,
122  const bool evaluate_hessians);
123 
124  static void
125  integrate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
126  Number * values_dofs_actual,
127  Number * values_quad,
128  Number * gradients_quad,
129  Number * scratch_data,
130  const bool integrate_values,
131  const bool integrate_gradients,
132  const bool add_into_values_array);
133  };
134 
135 
136 
137  template <MatrixFreeFunctions::ElementType type,
138  int dim,
139  int fe_degree,
140  int n_q_points_1d,
141  int n_components,
142  typename Number>
143  inline void
144  FEEvaluationImpl<type, dim, fe_degree, n_q_points_1d, n_components, Number>::
145  evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
146  const Number * values_dofs_actual,
147  Number * values_quad,
148  Number * gradients_quad,
149  Number * hessians_quad,
150  Number * scratch_data,
151  const bool evaluate_values,
152  const bool evaluate_gradients,
153  const bool evaluate_hessians)
154  {
155  if (evaluate_values == false && evaluate_gradients == false &&
156  evaluate_hessians == false)
157  return;
158 
159  const EvaluatorVariant variant =
160  EvaluatorSelector<type, (fe_degree + n_q_points_1d > 4)>::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 =
203  (type == MatrixFreeFunctions::truncated_tensor) ?
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;
207  if (type == MatrixFreeFunctions::truncated_tensor)
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  unsigned int count_p = 0, count_q = 0;
215  for (int i = 0; i < (dim > 2 ? degree + 1 : 1); ++i)
216  {
217  for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
218  {
219  for (int k = 0; k < degree + 1 - j - i;
220  ++k, ++count_p, ++count_q)
221  for (unsigned int c = 0; c < n_components; ++c)
222  values_dofs_tmp[c * dofs_per_comp + count_q] =
223  values_dofs_actual
224  [c * shape_info.dofs_per_component_on_cell + count_p];
225  for (int k = degree + 1 - j - i; k < degree + 1; ++k, ++count_q)
226  for (unsigned int c = 0; c < n_components; ++c)
227  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
228  }
229  for (int j = degree + 1 - i; j < degree + 1; ++j)
230  for (int k = 0; k < degree + 1; ++k, ++count_q)
231  for (unsigned int c = 0; c < n_components; ++c)
232  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
233  }
234  AssertDimension(count_q, dofs_per_comp);
235  values_dofs = values_dofs_tmp;
236  }
237 
238  switch (dim)
239  {
240  case 1:
241  for (unsigned int c = 0; c < n_components; c++)
242  {
243  if (evaluate_values == true)
244  eval.template values<0, true, false>(values_dofs, values_quad);
245  if (evaluate_gradients == true)
246  eval.template gradients<0, true, false>(values_dofs,
247  gradients_quad);
248  if (evaluate_hessians == true)
249  eval.template hessians<0, true, false>(values_dofs,
250  hessians_quad);
251 
252  // advance the next component in 1D array
253  values_dofs += dofs_per_comp;
254  values_quad += n_q_points;
255  gradients_quad += n_q_points;
256  hessians_quad += n_q_points;
257  }
258  break;
259 
260  case 2:
261  for (unsigned int c = 0; c < n_components; c++)
262  {
263  // grad x
264  if (evaluate_gradients == true)
265  {
266  eval.template gradients<0, true, false>(values_dofs, temp1);
267  eval.template values<1, true, false>(temp1, gradients_quad);
268  }
269  if (evaluate_hessians == true)
270  {
271  // grad xy
272  if (evaluate_gradients == false)
273  eval.template gradients<0, true, false>(values_dofs, temp1);
274  eval.template gradients<1, true, false>(temp1,
275  hessians_quad +
276  2 * n_q_points);
277 
278  // grad xx
279  eval.template hessians<0, true, false>(values_dofs, temp1);
280  eval.template values<1, true, false>(temp1, hessians_quad);
281  }
282 
283  // grad y
284  eval.template values<0, true, false>(values_dofs, temp1);
285  if (evaluate_gradients == true)
286  eval.template gradients<1, true, false>(temp1,
287  gradients_quad +
288  n_q_points);
289 
290  // grad yy
291  if (evaluate_hessians == true)
292  eval.template hessians<1, true, false>(temp1,
293  hessians_quad +
294  n_q_points);
295 
296  // val: can use values applied in x
297  if (evaluate_values == true)
298  eval.template values<1, true, false>(temp1, values_quad);
299 
300  // advance to the next component in 1D array
301  values_dofs += dofs_per_comp;
302  values_quad += n_q_points;
303  gradients_quad += 2 * n_q_points;
304  hessians_quad += 3 * n_q_points;
305  }
306  break;
307 
308  case 3:
309  for (unsigned int c = 0; c < n_components; c++)
310  {
311  if (evaluate_gradients == true)
312  {
313  // grad x
314  eval.template gradients<0, true, false>(values_dofs, temp1);
315  eval.template values<1, true, false>(temp1, temp2);
316  eval.template values<2, true, false>(temp2, gradients_quad);
317  }
318 
319  if (evaluate_hessians == true)
320  {
321  // grad xz
322  if (evaluate_gradients == false)
323  {
324  eval.template gradients<0, true, false>(values_dofs,
325  temp1);
326  eval.template values<1, true, false>(temp1, temp2);
327  }
328  eval.template gradients<2, true, false>(temp2,
329  hessians_quad +
330  4 * n_q_points);
331 
332  // grad xy
333  eval.template gradients<1, true, false>(temp1, temp2);
334  eval.template values<2, true, false>(temp2,
335  hessians_quad +
336  3 * n_q_points);
337 
338  // grad xx
339  eval.template hessians<0, true, false>(values_dofs, temp1);
340  eval.template values<1, true, false>(temp1, temp2);
341  eval.template values<2, true, false>(temp2, hessians_quad);
342  }
343 
344  // grad y
345  eval.template values<0, true, false>(values_dofs, temp1);
346  if (evaluate_gradients == true)
347  {
348  eval.template gradients<1, true, false>(temp1, temp2);
349  eval.template values<2, true, false>(temp2,
350  gradients_quad +
351  n_q_points);
352  }
353 
354  if (evaluate_hessians == true)
355  {
356  // grad yz
357  if (evaluate_gradients == false)
358  eval.template gradients<1, true, false>(temp1, temp2);
359  eval.template gradients<2, true, false>(temp2,
360  hessians_quad +
361  5 * n_q_points);
362 
363  // grad yy
364  eval.template hessians<1, true, false>(temp1, temp2);
365  eval.template values<2, true, false>(temp2,
366  hessians_quad +
367  n_q_points);
368  }
369 
370  // grad z: can use the values applied in x direction stored in
371  // temp1
372  eval.template values<1, true, false>(temp1, temp2);
373  if (evaluate_gradients == true)
374  eval.template gradients<2, true, false>(temp2,
375  gradients_quad +
376  2 * n_q_points);
377 
378  // grad zz: can use the values applied in x and y direction stored
379  // in temp2
380  if (evaluate_hessians == true)
381  eval.template hessians<2, true, false>(temp2,
382  hessians_quad +
383  2 * n_q_points);
384 
385  // val: can use the values applied in x & y direction stored in
386  // temp2
387  if (evaluate_values == true)
388  eval.template values<2, true, false>(temp2, values_quad);
389 
390  // advance to the next component in 1D array
391  values_dofs += dofs_per_comp;
392  values_quad += n_q_points;
393  gradients_quad += 3 * n_q_points;
394  hessians_quad += 6 * n_q_points;
395  }
396  break;
397 
398  default:
399  AssertThrow(false, ExcNotImplemented());
400  }
401 
402  // case additional dof for FE_Q_DG0: add values; gradients and second
403  // derivatives evaluate to zero
404  if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 &&
405  evaluate_values)
406  {
407  values_quad -= n_components * n_q_points;
408  values_dofs -= n_components * dofs_per_comp;
409  for (unsigned int c = 0; c < n_components; ++c)
410  for (unsigned int q = 0; q < shape_info.n_q_points; ++q)
411  values_quad[c * shape_info.n_q_points + q] +=
412  values_dofs[(c + 1) * shape_info.dofs_per_component_on_cell - 1];
413  }
414  }
415 
416 
417 
418  template <MatrixFreeFunctions::ElementType type,
419  int dim,
420  int fe_degree,
421  int n_q_points_1d,
422  int n_components,
423  typename Number>
424  inline void
425  FEEvaluationImpl<type, dim, fe_degree, n_q_points_1d, n_components, Number>::
426  integrate(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 integrate_values,
432  const bool integrate_gradients,
433  const bool add_into_values_array)
434  {
435  const EvaluatorVariant variant =
436  EvaluatorSelector<type, (fe_degree + n_q_points_1d > 4)>::variant;
437  using Eval = EvaluatorTensorProduct<variant,
438  dim,
439  fe_degree + 1,
440  n_q_points_1d,
441  Number>;
442  Eval eval(variant == evaluate_evenodd ?
443  shape_info.data.front().shape_values_eo :
444  shape_info.data.front().shape_values,
445  variant == evaluate_evenodd ?
446  shape_info.data.front().shape_gradients_eo :
447  shape_info.data.front().shape_gradients,
448  variant == evaluate_evenodd ?
449  shape_info.data.front().shape_hessians_eo :
450  shape_info.data.front().shape_hessians,
451  shape_info.data.front().fe_degree + 1,
452  shape_info.data.front().n_q_points_1d);
453 
454  const unsigned int temp_size =
455  Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
456  0 :
457  (Eval::n_rows_of_product > Eval::n_columns_of_product ?
458  Eval::n_rows_of_product :
459  Eval::n_columns_of_product);
460  Number *temp1;
461  Number *temp2;
462  if (temp_size == 0)
463  {
464  temp1 = scratch_data;
465  temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
466  shape_info.data.front().fe_degree + 1),
467  Utilities::fixed_power<dim>(
468  shape_info.data.front().n_q_points_1d));
469  }
470  else
471  {
472  temp1 = scratch_data;
473  temp2 = temp1 + temp_size;
474  }
475 
476  const unsigned int n_q_points =
477  temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
478  const unsigned int dofs_per_comp =
479  (type == MatrixFreeFunctions::truncated_tensor) ?
480  Utilities::fixed_power<dim>(shape_info.data.front().fe_degree + 1) :
481  shape_info.dofs_per_component_on_cell;
482  // expand dof_values to tensor product for truncated tensor products
483  Number *values_dofs =
484  (type == MatrixFreeFunctions::truncated_tensor) ?
485  scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
486  shape_info.n_q_points)) :
487  values_dofs_actual;
488 
489  switch (dim)
490  {
491  case 1:
492  for (unsigned int c = 0; c < n_components; c++)
493  {
494  if (integrate_values == true)
495  {
496  if (add_into_values_array == false)
497  eval.template values<0, false, false>(values_quad,
498  values_dofs);
499  else
500  eval.template values<0, false, true>(values_quad,
501  values_dofs);
502  }
503  if (integrate_gradients == true)
504  {
505  if (integrate_values == true || add_into_values_array == true)
506  eval.template gradients<0, false, true>(gradients_quad,
507  values_dofs);
508  else
509  eval.template gradients<0, false, false>(gradients_quad,
510  values_dofs);
511  }
512 
513  // advance to the next component in 1D array
514  values_dofs += dofs_per_comp;
515  values_quad += n_q_points;
516  gradients_quad += n_q_points;
517  }
518  break;
519 
520  case 2:
521  for (unsigned int c = 0; c < n_components; c++)
522  {
523  if (integrate_values == true && integrate_gradients == false)
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 (integrate_gradients == true)
532  {
533  eval.template gradients<1, false, false>(gradients_quad +
534  n_q_points,
535  temp1);
536  if (integrate_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 (integrate_values == true && integrate_gradients == false)
557  {
558  eval.template values<2, false, false>(values_quad, temp1);
559  eval.template values<1, false, false>(temp1, temp2);
560  if (add_into_values_array == false)
561  eval.template values<0, false, false>(temp2, values_dofs);
562  else
563  eval.template values<0, false, true>(temp2, values_dofs);
564  }
565  if (integrate_gradients == true)
566  {
567  eval.template gradients<2, false, false>(gradients_quad +
568  2 * n_q_points,
569  temp1);
570  if (integrate_values)
571  eval.template values<2, false, true>(values_quad, temp1);
572  eval.template values<1, false, false>(temp1, temp2);
573  eval.template values<2, false, false>(gradients_quad +
574  n_q_points,
575  temp1);
576  eval.template gradients<1, false, true>(temp1, temp2);
577  if (add_into_values_array == false)
578  eval.template values<0, false, false>(temp2, values_dofs);
579  else
580  eval.template values<0, false, true>(temp2, values_dofs);
581  eval.template values<2, false, false>(gradients_quad, temp1);
582  eval.template values<1, false, false>(temp1, temp2);
583  eval.template gradients<0, false, true>(temp2, values_dofs);
584  }
585 
586  // advance to the next component in 1D array
587  values_dofs += dofs_per_comp;
588  values_quad += n_q_points;
589  gradients_quad += 3 * n_q_points;
590  }
591  break;
592 
593  default:
594  AssertThrow(false, ExcNotImplemented());
595  }
596 
597  // case FE_Q_DG0: add values, gradients and second derivatives are zero
598  if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
599  {
600  values_dofs -= n_components * dofs_per_comp -
601  shape_info.dofs_per_component_on_cell + 1;
602  values_quad -= n_components * n_q_points;
603  if (integrate_values)
604  for (unsigned int c = 0; c < n_components; ++c)
605  {
606  values_dofs[0] = values_quad[0];
607  for (unsigned int q = 1; q < shape_info.n_q_points; ++q)
608  values_dofs[0] += values_quad[q];
609  values_dofs += dofs_per_comp;
610  values_quad += n_q_points;
611  }
612  else
613  {
614  for (unsigned int c = 0; c < n_components; ++c)
615  values_dofs[c * shape_info.dofs_per_component_on_cell] = Number();
616  values_dofs += n_components * shape_info.dofs_per_component_on_cell;
617  }
618  }
619 
620  if (type == MatrixFreeFunctions::truncated_tensor)
621  {
622  values_dofs -= dofs_per_comp * n_components;
623  unsigned int count_p = 0, count_q = 0;
624  const int degree =
625  fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
626  for (int i = 0; i < (dim > 2 ? degree + 1 : 1); ++i)
627  {
628  for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
629  {
630  for (int k = 0; k < degree + 1 - j - i;
631  ++k, ++count_p, ++count_q)
632  {
633  for (unsigned int c = 0; c < n_components; ++c)
634  values_dofs_actual
635  [c * shape_info.dofs_per_component_on_cell + count_p] =
636  values_dofs[c * dofs_per_comp + count_q];
637  }
638  count_q += j + i;
639  }
640  count_q += i * (degree + 1);
641  }
642  AssertDimension(count_q,
643  Utilities::fixed_power<dim>(
644  shape_info.data.front().fe_degree + 1));
645  }
646  }
647 
648 
649 
661  template <EvaluatorVariant variant,
662  int dim,
663  int basis_size_1,
664  int basis_size_2,
665  int n_components,
666  typename Number,
667  typename Number2>
668  struct FEEvaluationImplBasisChange
669  {
670  static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
671  "The second dimension must not be smaller than the first");
672 
694 #ifndef DEBUG
695  DEAL_II_ALWAYS_INLINE
696 #endif
697  static void
698  do_forward(
699  const AlignedVector<Number2> &transformation_matrix,
700  const Number * values_in,
701  Number * values_out,
702  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
703  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
704  {
705  Assert(
706  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
707  ExcMessage("The second dimension must not be smaller than the first"));
708 
709  // we do recursion until dim==1 or dim==2 and we have
710  // basis_size_1==basis_size_2. The latter optimization increases
711  // optimization possibilities for the compiler but does only work for
712  // aliased pointers if the sizes are equal.
713  constexpr int next_dim =
714  (dim > 2 ||
715  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
716  dim - 1 :
717  dim;
718 
719  EvaluatorTensorProduct<variant,
720  dim,
721  basis_size_1,
722  (basis_size_1 == 0 ? 0 : basis_size_2),
723  Number,
724  Number2>
725  eval_val(transformation_matrix,
726  AlignedVector<Number2>(),
727  AlignedVector<Number2>(),
728  basis_size_1_variable,
729  basis_size_2_variable);
730  const unsigned int np_1 =
731  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
732  const unsigned int np_2 =
733  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
734  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
735  ExcMessage("Cannot transform with 0-point basis"));
736  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
737  ExcMessage("Cannot transform with 0-point basis"));
738 
739  // run loop backwards to ensure correctness if values_in aliases with
740  // values_out in case with basis_size_1 < basis_size_2
741  values_in = values_in + n_components * Utilities::fixed_power<dim>(np_1);
742  values_out =
743  values_out + n_components * Utilities::fixed_power<dim>(np_2);
744  for (unsigned int c = n_components; c != 0; --c)
745  {
746  values_in -= Utilities::fixed_power<dim>(np_1);
747  values_out -= Utilities::fixed_power<dim>(np_2);
748  if (next_dim < dim)
749  for (unsigned int q = np_1; q != 0; --q)
750  FEEvaluationImplBasisChange<
751  variant,
752  next_dim,
753  basis_size_1,
754  basis_size_2,
755  1,
756  Number,
757  Number2>::do_forward(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
814  DEAL_II_ALWAYS_INLINE
815 #endif
816  static void
817  do_backward(
818  const AlignedVector<Number2> &transformation_matrix,
819  const bool add_into_result,
820  Number * values_in,
821  Number * values_out,
822  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
823  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
824  {
825  Assert(
826  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
827  ExcMessage("The second dimension must not be smaller than the first"));
828  Assert(add_into_result == false || values_in != values_out,
829  ExcMessage(
830  "Input and output cannot alias with each other when "
831  "adding the result of the basis change to existing data"));
832 
833  constexpr int next_dim =
834  (dim > 2 ||
835  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
836  dim - 1 :
837  dim;
838  EvaluatorTensorProduct<variant,
839  dim,
840  basis_size_1,
841  (basis_size_1 == 0 ? 0 : basis_size_2),
842  Number,
843  Number2>
844  eval_val(transformation_matrix,
845  AlignedVector<Number2>(),
846  AlignedVector<Number2>(),
847  basis_size_1_variable,
848  basis_size_2_variable);
849  const unsigned int np_1 =
850  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
851  const unsigned int np_2 =
852  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
853  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
854  ExcMessage("Cannot transform with 0-point basis"));
855  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
856  ExcMessage("Cannot transform with 0-point basis"));
857 
858  for (unsigned int c = 0; c < n_components; ++c)
859  {
860  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
861  {
862  eval_val.template values<1, false, false>(values_in, values_in);
863  if (add_into_result)
864  eval_val.template values<0, false, true>(values_in, values_out);
865  else
866  eval_val.template values<0, false, false>(values_in,
867  values_out);
868  }
869  else
870  {
871  if (dim == 1 && add_into_result)
872  eval_val.template values<0, false, true>(values_in, values_out);
873  else if (dim == 1)
874  eval_val.template values<0, false, false>(values_in,
875  values_out);
876  else
877  eval_val.template values<dim - 1, false, false>(values_in,
878  values_in);
879  }
880  if (next_dim < dim)
881  for (unsigned int q = 0; q < np_1; ++q)
882  FEEvaluationImplBasisChange<variant,
883  next_dim,
884  basis_size_1,
885  basis_size_2,
886  1,
887  Number,
888  Number2>::
889  do_backward(transformation_matrix,
890  add_into_result,
891  values_in +
892  q * Utilities::fixed_power<next_dim>(np_2),
893  values_out +
894  q * Utilities::fixed_power<next_dim>(np_1),
895  basis_size_1_variable,
896  basis_size_2_variable);
897 
898  values_in += Utilities::fixed_power<dim>(np_2);
899  values_out += Utilities::fixed_power<dim>(np_1);
900  }
901  }
902 
922  static void
923  do_mass(const AlignedVector<Number2> &transformation_matrix,
924  const AlignedVector<Number> & coefficients,
925  const Number * values_in,
926  Number * scratch_data,
927  Number * values_out)
928  {
929  constexpr int next_dim = dim > 1 ? dim - 1 : dim;
930  Number * my_scratch =
931  basis_size_1 != basis_size_2 ? scratch_data : values_out;
932 
933  const unsigned int size_per_component = Utilities::pow(basis_size_2, dim);
934  Assert(coefficients.size() == size_per_component ||
935  coefficients.size() == n_components * size_per_component,
936  ExcDimensionMismatch(coefficients.size(), size_per_component));
937  const unsigned int stride =
938  coefficients.size() == size_per_component ? 0 : 1;
939 
940  for (unsigned int q = basis_size_1; q != 0; --q)
941  FEEvaluationImplBasisChange<
942  variant,
943  next_dim,
944  basis_size_1,
945  basis_size_2,
946  n_components,
947  Number,
948  Number2>::do_forward(transformation_matrix,
949  values_in +
950  (q - 1) *
951  Utilities::pow(basis_size_1, dim - 1),
952  my_scratch +
953  (q - 1) *
954  Utilities::pow(basis_size_2, dim - 1));
955  EvaluatorTensorProduct<variant,
956  dim,
957  basis_size_1,
958  basis_size_2,
959  Number,
960  Number2>
961  eval_val(transformation_matrix);
962  const unsigned int n_inner_blocks =
963  (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
964  const unsigned int n_blocks = Utilities::pow(basis_size_2, dim - 1);
965  for (unsigned int ii = 0; ii < n_blocks; ii += n_inner_blocks)
966  for (unsigned int c = 0; c < n_components; ++c)
967  {
968  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
969  eval_val.template values_one_line<dim - 1, true, false>(
970  my_scratch + i, my_scratch + i);
971  for (unsigned int q = 0; q < basis_size_2; ++q)
972  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
973  my_scratch[i + q * n_blocks + c * size_per_component] *=
974  coefficients[i + q * n_blocks +
975  c * stride * size_per_component];
976  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
977  eval_val.template values_one_line<dim - 1, false, false>(
978  my_scratch + i, my_scratch + i);
979  }
980  for (unsigned int q = 0; q < basis_size_1; ++q)
981  FEEvaluationImplBasisChange<
982  variant,
983  next_dim,
984  basis_size_1,
985  basis_size_2,
986  n_components,
987  Number,
988  Number2>::do_backward(transformation_matrix,
989  false,
990  my_scratch +
991  q * Utilities::pow(basis_size_2, dim - 1),
992  values_out +
993  q * Utilities::pow(basis_size_1, dim - 1));
994  }
995  };
996 
997 
998 
1013  template <int dim, int fe_degree, int n_components, typename Number>
1014  struct FEEvaluationImplCollocation
1015  {
1016  static void
1017  evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1018  const Number * values_dofs,
1019  Number * values_quad,
1020  Number * gradients_quad,
1021  Number * hessians_quad,
1022  Number * scratch_data,
1023  const bool evaluate_values,
1024  const bool evaluate_gradients,
1025  const bool evaluate_hessians);
1026 
1027  static void
1028  integrate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1029  Number * values_dofs,
1030  Number * values_quad,
1031  Number * gradients_quad,
1032  Number * scratch_data,
1033  const bool integrate_values,
1034  const bool integrate_gradients,
1035  const bool add_into_values_array);
1036  };
1037 
1038 
1039 
1040  template <int dim, int fe_degree, int n_components, typename Number>
1041  inline void
1042  FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>::evaluate(
1043  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1044  const Number * values_dofs,
1045  Number * values_quad,
1046  Number * gradients_quad,
1047  Number * hessians_quad,
1048  Number *,
1049  const bool evaluate_values,
1050  const bool evaluate_gradients,
1051  const bool evaluate_hessians)
1052  {
1053  AssertDimension(
1054  shape_info.data.front().shape_gradients_collocation_eo.size(),
1055  (fe_degree + 2) / 2 * (fe_degree + 1));
1056 
1057  EvaluatorTensorProduct<evaluate_evenodd,
1058  dim,
1059  fe_degree + 1,
1060  fe_degree + 1,
1061  Number>
1062  eval(AlignedVector<Number>(),
1063  shape_info.data.front().shape_gradients_collocation_eo,
1064  shape_info.data.front().shape_hessians_collocation_eo);
1065  constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
1066 
1067  for (unsigned int c = 0; c < n_components; c++)
1068  {
1069  if (evaluate_values == true)
1070  for (unsigned int i = 0; i < n_q_points; ++i)
1071  values_quad[i] = values_dofs[i];
1072  if (evaluate_gradients == true || evaluate_hessians == true)
1073  {
1074  eval.template gradients<0, true, false>(values_dofs,
1075  gradients_quad);
1076  if (dim > 1)
1077  eval.template gradients<1, true, false>(values_dofs,
1078  gradients_quad +
1079  n_q_points);
1080  if (dim > 2)
1081  eval.template gradients<2, true, false>(values_dofs,
1082  gradients_quad +
1083  2 * n_q_points);
1084  }
1085  if (evaluate_hessians == true)
1086  {
1087  eval.template hessians<0, true, false>(values_dofs, hessians_quad);
1088  if (dim > 1)
1089  {
1090  eval.template gradients<1, true, false>(gradients_quad,
1091  hessians_quad +
1092  dim * n_q_points);
1093  eval.template hessians<1, true, false>(values_dofs,
1094  hessians_quad +
1095  n_q_points);
1096  }
1097  if (dim > 2)
1098  {
1099  eval.template gradients<2, true, false>(gradients_quad,
1100  hessians_quad +
1101  4 * n_q_points);
1102  eval.template gradients<2, true, false>(
1103  gradients_quad + n_q_points, hessians_quad + 5 * n_q_points);
1104  eval.template hessians<2, true, false>(values_dofs,
1105  hessians_quad +
1106  2 * n_q_points);
1107  }
1108  hessians_quad += (dim * (dim + 1)) / 2 * n_q_points;
1109  }
1110  gradients_quad += dim * n_q_points;
1111  values_quad += n_q_points;
1112  values_dofs += n_q_points;
1113  }
1114  }
1115 
1116 
1117 
1118  template <int dim, int fe_degree, int n_components, typename Number>
1119  inline void
1120  FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>::integrate(
1121  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1122  Number * values_dofs,
1123  Number * values_quad,
1124  Number * gradients_quad,
1125  Number *,
1126  const bool integrate_values,
1127  const bool integrate_gradients,
1128  const bool add_into_values_array)
1129  {
1130  AssertDimension(
1131  shape_info.data.front().shape_gradients_collocation_eo.size(),
1132  (fe_degree + 2) / 2 * (fe_degree + 1));
1133 
1134  EvaluatorTensorProduct<evaluate_evenodd,
1135  dim,
1136  fe_degree + 1,
1137  fe_degree + 1,
1138  Number>
1139  eval(AlignedVector<Number>(),
1140  shape_info.data.front().shape_gradients_collocation_eo,
1141  shape_info.data.front().shape_hessians_collocation_eo);
1142  constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
1143 
1144  for (unsigned int c = 0; c < n_components; c++)
1145  {
1146  if (integrate_values == true && add_into_values_array == false)
1147  for (unsigned int i = 0; i < n_q_points; ++i)
1148  values_dofs[i] = values_quad[i];
1149  else if (integrate_values == true)
1150  for (unsigned int i = 0; i < n_q_points; ++i)
1151  values_dofs[i] += values_quad[i];
1152  if (integrate_gradients == true)
1153  {
1154  if (integrate_values == true || add_into_values_array == true)
1155  eval.template gradients<0, false, true>(gradients_quad,
1156  values_dofs);
1157  else
1158  eval.template gradients<0, false, false>(gradients_quad,
1159  values_dofs);
1160  if (dim > 1)
1161  eval.template gradients<1, false, true>(gradients_quad +
1162  n_q_points,
1163  values_dofs);
1164  if (dim > 2)
1165  eval.template gradients<2, false, true>(gradients_quad +
1166  2 * n_q_points,
1167  values_dofs);
1168  }
1169  gradients_quad += dim * n_q_points;
1170  values_quad += n_q_points;
1171  values_dofs += n_q_points;
1172  }
1173  }
1174 
1175 
1176 
1189  template <int dim,
1190  int fe_degree,
1191  int n_q_points_1d,
1192  int n_components,
1193  typename Number>
1194  struct FEEvaluationImplTransformToCollocation
1195  {
1196  static void
1197  evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1198  const Number * values_dofs,
1199  Number * values_quad,
1200  Number * gradients_quad,
1201  Number * hessians_quad,
1202  Number * scratch_data,
1203  const bool evaluate_values,
1204  const bool evaluate_gradients,
1205  const bool evaluate_hessians);
1206 
1207  static void
1208  integrate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1209  Number * values_dofs,
1210  Number * values_quad,
1211  Number * gradients_quad,
1212  Number * scratch_data,
1213  const bool integrate_values,
1214  const bool integrate_gradients,
1215  const bool add_into_values_array);
1216  };
1217 
1218 
1219 
1220  template <int dim,
1221  int fe_degree,
1222  int n_q_points_1d,
1223  int n_components,
1224  typename Number>
1225  inline void
1226  FEEvaluationImplTransformToCollocation<
1227  dim,
1228  fe_degree,
1229  n_q_points_1d,
1230  n_components,
1231  Number>::evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1232  const Number * values_dofs,
1233  Number * values_quad,
1234  Number *gradients_quad,
1235  Number *hessians_quad,
1236  Number *,
1237  const bool,
1238  const bool evaluate_gradients,
1239  const bool evaluate_hessians)
1240  {
1241  Assert(n_q_points_1d > fe_degree,
1242  ExcMessage("You lose information when going to a collocation space "
1243  "of lower degree, so the evaluation results would be "
1244  "wrong. Thus, this class does not permit the desired "
1245  "operation."));
1246  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
1247 
1248  for (unsigned int c = 0; c < n_components; c++)
1249  {
1250  FEEvaluationImplBasisChange<
1251  evaluate_evenodd,
1252  dim,
1253  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1254  n_q_points_1d,
1255  1,
1256  Number,
1257  Number>::do_forward(shape_info.data.front().shape_values_eo,
1258  values_dofs,
1259  values_quad);
1260 
1261  // apply derivatives in the collocation space
1262  if (evaluate_gradients == true || evaluate_hessians == true)
1263  FEEvaluationImplCollocation<dim, n_q_points_1d - 1, 1, Number>::
1264  evaluate(shape_info,
1265  values_quad,
1266  nullptr,
1267  gradients_quad,
1268  hessians_quad,
1269  nullptr,
1270  false,
1271  evaluate_gradients,
1272  evaluate_hessians);
1273 
1274  values_dofs += shape_info.dofs_per_component_on_cell;
1275  values_quad += n_q_points;
1276  gradients_quad += dim * n_q_points;
1277  hessians_quad += (dim * (dim + 1)) / 2 * n_q_points;
1278  }
1279  }
1280 
1281 
1282 
1283  template <int dim,
1284  int fe_degree,
1285  int n_q_points_1d,
1286  int n_components,
1287  typename Number>
1288  inline void
1289  FEEvaluationImplTransformToCollocation<
1290  dim,
1291  fe_degree,
1292  n_q_points_1d,
1293  n_components,
1294  Number>::integrate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1295  Number *values_dofs,
1296  Number *values_quad,
1297  Number *gradients_quad,
1298  Number *,
1299  const bool integrate_values,
1300  const bool integrate_gradients,
1301  const bool add_into_values_array)
1302  {
1303  Assert(n_q_points_1d > fe_degree,
1304  ExcMessage("You lose information when going to a collocation space "
1305  "of lower degree, so the evaluation results would be "
1306  "wrong. Thus, this class does not permit the desired "
1307  "operation."));
1308  AssertDimension(
1309  shape_info.data.front().shape_gradients_collocation_eo.size(),
1310  (n_q_points_1d + 1) / 2 * n_q_points_1d);
1311  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
1312 
1313  for (unsigned int c = 0; c < n_components; c++)
1314  {
1315  // apply derivatives in collocation space
1316  if (integrate_gradients == true)
1317  FEEvaluationImplCollocation<dim, n_q_points_1d - 1, 1, Number>::
1318  integrate(shape_info,
1319  values_quad,
1320  nullptr,
1321  gradients_quad,
1322  nullptr,
1323  false,
1324  integrate_gradients,
1325  /*add_into_values_array=*/integrate_values);
1326 
1327  // transform back to the original space
1328  FEEvaluationImplBasisChange<
1329  evaluate_evenodd,
1330  dim,
1331  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1332  n_q_points_1d,
1333  1,
1334  Number,
1335  Number>::do_backward(shape_info.data.front().shape_values_eo,
1336  add_into_values_array,
1337  values_quad,
1338  values_dofs);
1339  gradients_quad += dim * n_q_points;
1340  values_quad += n_q_points;
1341  values_dofs += shape_info.dofs_per_component_on_cell;
1342  }
1343  }
1344 
1345 
1346 
1347  template <bool symmetric_evaluate,
1348  int dim,
1349  int fe_degree,
1350  int n_q_points_1d,
1351  int n_components,
1352  typename Number>
1353  struct FEFaceEvaluationImpl
1354  {
1355  // We enable a transformation to collocation for derivatives if it gives
1356  // correct results (first two conditions), if it is the most efficient
1357  // choice in terms of operation counts (third condition) and if we were
1358  // able to initialize the fields in shape_info.templates.h from the
1359  // polynomials (fourth condition).
1360  static constexpr bool use_collocation =
1361  symmetric_evaluate &&
1362  n_q_points_1d > fe_degree &&n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
1363  n_q_points_1d < 200;
1364 
1365  static void
1366  evaluate_in_face(const MatrixFreeFunctions::ShapeInfo<Number> &data,
1367  Number * values_dofs,
1368  Number * values_quad,
1369  Number * gradients_quad,
1370  Number * scratch_data,
1371  const bool evaluate_val,
1372  const bool evaluate_grad,
1373  const unsigned int subface_index)
1374  {
1375  const AlignedVector<Number> &val1 =
1376  symmetric_evaluate ?
1377  data.data.front().shape_values_eo :
1378  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1379  data.data.front().shape_values :
1380  data.data.front().values_within_subface[subface_index % 2]);
1381  const AlignedVector<Number> &val2 =
1382  symmetric_evaluate ?
1383  data.data.front().shape_values_eo :
1384  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1385  data.data.front().shape_values :
1386  data.data.front().values_within_subface[subface_index / 2]);
1387 
1388  const AlignedVector<Number> &grad1 =
1389  symmetric_evaluate ?
1390  data.data.front().shape_gradients_eo :
1391  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1392  data.data.front().shape_gradients :
1393  data.data.front().gradients_within_subface[subface_index % 2]);
1394  const AlignedVector<Number> &grad2 =
1395  symmetric_evaluate ?
1396  data.data.front().shape_gradients_eo :
1397  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1398  data.data.front().shape_gradients :
1399  data.data.front().gradients_within_subface[subface_index / 2]);
1400 
1401  using Eval =
1402  internal::EvaluatorTensorProduct<symmetric_evaluate ?
1403  internal::evaluate_evenodd :
1404  internal::evaluate_general,
1405  dim - 1,
1406  fe_degree + 1,
1407  n_q_points_1d,
1408  Number>;
1409  Eval eval1(val1,
1410  grad1,
1411  AlignedVector<Number>(),
1412  data.data.front().fe_degree + 1,
1413  data.data.front().n_q_points_1d);
1414  Eval eval2(val2,
1415  grad2,
1416  AlignedVector<Number>(),
1417  data.data.front().fe_degree + 1,
1418  data.data.front().n_q_points_1d);
1419 
1420  const unsigned int size_deg =
1421  fe_degree > -1 ?
1422  Utilities::pow(fe_degree + 1, dim - 1) :
1423  (dim > 1 ?
1424  Utilities::fixed_power<dim - 1>(data.data.front().fe_degree + 1) :
1425  1);
1426 
1427  const unsigned int n_q_points = fe_degree > -1 ?
1428  Utilities::pow(n_q_points_1d, dim - 1) :
1429  data.n_q_points_face;
1430 
1431  if (evaluate_grad == false)
1432  for (unsigned int c = 0; c < n_components; ++c)
1433  {
1434  switch (dim)
1435  {
1436  case 3:
1437  eval1.template values<0, true, false>(values_dofs,
1438  values_quad);
1439  eval2.template values<1, true, false>(values_quad,
1440  values_quad);
1441  break;
1442  case 2:
1443  eval1.template values<0, true, false>(values_dofs,
1444  values_quad);
1445  break;
1446  case 1:
1447  values_quad[0] = values_dofs[0];
1448  break;
1449  default:
1450  Assert(false, ExcNotImplemented());
1451  }
1452  values_dofs += 2 * size_deg;
1453  values_quad += n_q_points;
1454  }
1455  else
1456  for (unsigned int c = 0; c < n_components; ++c)
1457  {
1458  switch (dim)
1459  {
1460  case 3:
1461  if (use_collocation)
1462  {
1463  eval1.template values<0, true, false>(values_dofs,
1464  values_quad);
1465  eval1.template values<1, true, false>(values_quad,
1466  values_quad);
1467  internal::EvaluatorTensorProduct<
1468  internal::evaluate_evenodd,
1469  dim - 1,
1470  n_q_points_1d,
1471  n_q_points_1d,
1472  Number>
1473  eval_grad(
1474  AlignedVector<Number>(),
1475  data.data.front().shape_gradients_collocation_eo,
1476  AlignedVector<Number>());
1477  eval_grad.template gradients<0, true, false>(
1478  values_quad, gradients_quad);
1479  eval_grad.template gradients<1, true, false>(
1480  values_quad, gradients_quad + n_q_points);
1481  }
1482  else
1483  {
1484  eval1.template gradients<0, true, false>(values_dofs,
1485  scratch_data);
1486  eval2.template values<1, true, false>(scratch_data,
1487  gradients_quad);
1488 
1489  eval1.template values<0, true, false>(values_dofs,
1490  scratch_data);
1491  eval2.template gradients<1, true, false>(scratch_data,
1492  gradients_quad +
1493  n_q_points);
1494 
1495  if (evaluate_val == true)
1496  eval2.template values<1, true, false>(scratch_data,
1497  values_quad);
1498  }
1499  eval1.template values<0, true, false>(values_dofs + size_deg,
1500  scratch_data);
1501  eval2.template values<1, true, false>(
1502  scratch_data, gradients_quad + (dim - 1) * n_q_points);
1503 
1504  break;
1505  case 2:
1506  eval1.template values<0, true, false>(values_dofs + size_deg,
1507  gradients_quad +
1508  (dim - 1) *
1509  n_q_points);
1510  eval1.template gradients<0, true, false>(values_dofs,
1511  gradients_quad);
1512  if (evaluate_val == true)
1513  eval1.template values<0, true, false>(values_dofs,
1514  values_quad);
1515  break;
1516  case 1:
1517  values_quad[0] = values_dofs[0];
1518  gradients_quad[0] = values_dofs[1];
1519  break;
1520  default:
1521  AssertThrow(false, ExcNotImplemented());
1522  }
1523  values_dofs += 2 * size_deg;
1524  values_quad += n_q_points;
1525  gradients_quad += dim * n_q_points;
1526  }
1527  }
1528 
1529  static void
1530  integrate_in_face(const MatrixFreeFunctions::ShapeInfo<Number> &data,
1531  Number * values_dofs,
1532  Number * values_quad,
1533  Number * gradients_quad,
1534  Number * scratch_data,
1535  const bool integrate_val,
1536  const bool integrate_grad,
1537  const unsigned int subface_index)
1538  {
1539  const AlignedVector<Number> &val1 =
1540  symmetric_evaluate ?
1541  data.data.front().shape_values_eo :
1542  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1543  data.data.front().shape_values :
1544  data.data.front().values_within_subface[subface_index % 2]);
1545  const AlignedVector<Number> &val2 =
1546  symmetric_evaluate ?
1547  data.data.front().shape_values_eo :
1548  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1549  data.data.front().shape_values :
1550  data.data.front().values_within_subface[subface_index / 2]);
1551 
1552  const AlignedVector<Number> &grad1 =
1553  symmetric_evaluate ?
1554  data.data.front().shape_gradients_eo :
1555  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1556  data.data.front().shape_gradients :
1557  data.data.front().gradients_within_subface[subface_index % 2]);
1558  const AlignedVector<Number> &grad2 =
1559  symmetric_evaluate ?
1560  data.data.front().shape_gradients_eo :
1561  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1562  data.data.front().shape_gradients :
1563  data.data.front().gradients_within_subface[subface_index / 2]);
1564 
1565  using Eval =
1566  internal::EvaluatorTensorProduct<symmetric_evaluate ?
1567  internal::evaluate_evenodd :
1568  internal::evaluate_general,
1569  dim - 1,
1570  fe_degree + 1,
1571  n_q_points_1d,
1572  Number>;
1573  Eval eval1(val1,
1574  grad1,
1575  val1,
1576  data.data.front().fe_degree + 1,
1577  data.data.front().n_q_points_1d);
1578  Eval eval2(val2,
1579  grad2,
1580  val1,
1581  data.data.front().fe_degree + 1,
1582  data.data.front().n_q_points_1d);
1583 
1584  const unsigned int size_deg =
1585  fe_degree > -1 ?
1586  Utilities::pow(fe_degree + 1, dim - 1) :
1587  (dim > 1 ?
1588  Utilities::fixed_power<dim - 1>(data.data.front().fe_degree + 1) :
1589  1);
1590 
1591  const unsigned int n_q_points = fe_degree > -1 ?
1592  Utilities::pow(n_q_points_1d, dim - 1) :
1593  data.n_q_points_face;
1594 
1595  if (integrate_grad == false)
1596  for (unsigned int c = 0; c < n_components; ++c)
1597  {
1598  switch (dim)
1599  {
1600  case 3:
1601  eval2.template values<1, false, false>(values_quad,
1602  values_quad);
1603  eval1.template values<0, false, false>(values_quad,
1604  values_dofs);
1605  break;
1606  case 2:
1607  eval1.template values<0, false, false>(values_quad,
1608  values_dofs);
1609  break;
1610  case 1:
1611  values_dofs[0] = values_quad[0];
1612  break;
1613  default:
1614  Assert(false, ExcNotImplemented());
1615  }
1616  values_dofs += 2 * size_deg;
1617  values_quad += n_q_points;
1618  }
1619  else
1620  for (unsigned int c = 0; c < n_components; ++c)
1621  {
1622  switch (dim)
1623  {
1624  case 3:
1625  eval2.template values<1, false, false>(gradients_quad +
1626  2 * n_q_points,
1627  gradients_quad +
1628  2 * n_q_points);
1629  eval1.template values<0, false, false>(
1630  gradients_quad + 2 * n_q_points, values_dofs + size_deg);
1631  if (use_collocation)
1632  {
1633  internal::EvaluatorTensorProduct<
1634  internal::evaluate_evenodd,
1635  dim - 1,
1636  n_q_points_1d,
1637  n_q_points_1d,
1638  Number>
1639  eval_grad(
1640  AlignedVector<Number>(),
1641  data.data.front().shape_gradients_collocation_eo,
1642  AlignedVector<Number>());
1643  if (integrate_val)
1644  eval_grad.template gradients<1, false, true>(
1645  gradients_quad + n_q_points, values_quad);
1646  else
1647  eval_grad.template gradients<1, false, false>(
1648  gradients_quad + n_q_points, values_quad);
1649  eval_grad.template gradients<0, false, true>(
1650  gradients_quad, values_quad);
1651  eval1.template values<1, false, false>(values_quad,
1652  values_quad);
1653  eval1.template values<0, false, false>(values_quad,
1654  values_dofs);
1655  }
1656  else
1657  {
1658  if (integrate_val)
1659  {
1660  eval2.template values<1, false, false>(values_quad,
1661  scratch_data);
1662  eval2.template gradients<1, false, true>(
1663  gradients_quad + n_q_points, scratch_data);
1664  }
1665  else
1666  eval2.template gradients<1, false, false>(
1667  gradients_quad + n_q_points, scratch_data);
1668 
1669  eval1.template values<0, false, false>(scratch_data,
1670  values_dofs);
1671  eval2.template values<1, false, false>(gradients_quad,
1672  scratch_data);
1673  eval1.template gradients<0, false, true>(scratch_data,
1674  values_dofs);
1675  }
1676  break;
1677  case 2:
1678  eval1.template values<0, false, false>(
1679  gradients_quad + n_q_points, values_dofs + size_deg);
1680  eval1.template gradients<0, false, false>(gradients_quad,
1681  values_dofs);
1682  if (integrate_val == true)
1683  eval1.template values<0, false, true>(values_quad,
1684  values_dofs);
1685  break;
1686  case 1:
1687  values_dofs[0] = values_quad[0];
1688  values_dofs[1] = gradients_quad[0];
1689  break;
1690  default:
1691  AssertThrow(false, ExcNotImplemented());
1692  }
1693  values_dofs += 2 * size_deg;
1694  values_quad += n_q_points;
1695  gradients_quad += dim * n_q_points;
1696  }
1697  }
1698  };
1699 
1700 
1701 
1702  template <int dim, int fe_degree, int n_components, typename Number>
1703  struct FEFaceNormalEvaluationImpl
1704  {
1705  template <bool do_evaluate, bool add_into_output>
1706  static void
1707  interpolate(const MatrixFreeFunctions::ShapeInfo<Number> &data,
1708  const Number * input,
1709  Number * output,
1710  const bool do_gradients,
1711  const unsigned int face_no)
1712  {
1713  internal::EvaluatorTensorProduct<internal::evaluate_general,
1714  dim,
1715  fe_degree + 1,
1716  0,
1717  Number>
1718  evalf(data.data.front().shape_data_on_face[face_no % 2],
1719  AlignedVector<Number>(),
1720  AlignedVector<Number>(),
1721  data.data.front().fe_degree + 1,
1722  0);
1723 
1724  const unsigned int in_stride = do_evaluate ?
1725  data.dofs_per_component_on_cell :
1726  2 * data.dofs_per_component_on_face;
1727  const unsigned int out_stride = do_evaluate ?
1728  2 * data.dofs_per_component_on_face :
1729  data.dofs_per_component_on_cell;
1730  const unsigned int face_direction = face_no / 2;
1731  for (unsigned int c = 0; c < n_components; c++)
1732  {
1733  if (do_gradients)
1734  {
1735  if (face_direction == 0)
1736  evalf.template apply_face<0, do_evaluate, add_into_output, 1>(
1737  input, output);
1738  else if (face_direction == 1)
1739  evalf.template apply_face<1, do_evaluate, add_into_output, 1>(
1740  input, output);
1741  else
1742  evalf.template apply_face<2, do_evaluate, add_into_output, 1>(
1743  input, output);
1744  }
1745  else
1746  {
1747  if (face_direction == 0)
1748  evalf.template apply_face<0, do_evaluate, add_into_output, 0>(
1749  input, output);
1750  else if (face_direction == 1)
1751  evalf.template apply_face<1, do_evaluate, add_into_output, 0>(
1752  input, output);
1753  else
1754  evalf.template apply_face<2, do_evaluate, add_into_output, 0>(
1755  input, output);
1756  }
1757  input += in_stride;
1758  output += out_stride;
1759  }
1760  }
1761  };
1762 
1763 
1764 
1765  // internal helper function for reading data; base version of different types
1766  template <typename VectorizedArrayType, typename Number2>
1767  void
1768  do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
1769  {
1770  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1771  dst[v] = src_ptr[v];
1772  }
1773 
1774 
1775 
1776  // internal helper function for reading data; specialized version where we
1777  // can use a dedicated load function
1778  template <typename Number, unsigned int width>
1779  void
1780  do_vectorized_read(const Number *src_ptr, VectorizedArray<Number, width> &dst)
1781  {
1782  dst.load(src_ptr);
1783  }
1784 
1785 
1786 
1787  // internal helper function for reading data; base version of different types
1788  template <typename VectorizedArrayType, typename Number2>
1789  void
1790  do_vectorized_gather(const Number2 * src_ptr,
1791  const unsigned int * indices,
1792  VectorizedArrayType &dst)
1793  {
1794  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1795  dst[v] = src_ptr[indices[v]];
1796  }
1797 
1798 
1799 
1800  // internal helper function for reading data; specialized version where we
1801  // can use a dedicated gather function
1802  template <typename Number, unsigned int width>
1803  void
1804  do_vectorized_gather(const Number * src_ptr,
1805  const unsigned int * indices,
1806  VectorizedArray<Number, width> &dst)
1807  {
1808  dst.gather(src_ptr, indices);
1809  }
1810 
1811 
1812 
1813  // internal helper function for reading data; base version of different types
1814  template <typename VectorizedArrayType, typename Number2>
1815  void
1816  do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
1817  {
1818  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1819  dst_ptr[v] += src[v];
1820  }
1821 
1822 
1823 
1824  // internal helper function for reading data; specialized version where we
1825  // can use a dedicated load function
1826  template <typename Number, unsigned int width>
1827  void
1828  do_vectorized_add(const VectorizedArray<Number, width> src, Number *dst_ptr)
1829  {
1830  VectorizedArray<Number, width> tmp;
1831  tmp.load(dst_ptr);
1832  (tmp + src).store(dst_ptr);
1833  }
1834 
1835 
1836 
1837  // internal helper function for reading data; base version of different types
1838  template <typename VectorizedArrayType, typename Number2>
1839  void
1840  do_vectorized_scatter_add(const VectorizedArrayType src,
1841  const unsigned int * indices,
1842  Number2 * dst_ptr)
1843  {
1844  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1845  dst_ptr[indices[v]] += src[v];
1846  }
1847 
1848 
1849 
1850  // internal helper function for reading data; specialized version where we
1851  // can use a dedicated gather function
1852  template <typename Number, unsigned int width>
1853  void
1854  do_vectorized_scatter_add(const VectorizedArray<Number, width> src,
1855  const unsigned int * indices,
1856  Number * dst_ptr)
1857  {
1858 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
1859  for (unsigned int v = 0; v < width; ++v)
1860  dst_ptr[indices[v]] += src[v];
1861 #else
1862  VectorizedArray<Number, width> tmp;
1863  tmp.gather(dst_ptr, indices);
1864  (tmp + src).scatter(indices, dst_ptr);
1865 #endif
1866  }
1867 
1868 
1869 
1870  template <int dim,
1871  int fe_degree,
1872  int n_q_points_1d,
1873  int n_components,
1874  typename Number,
1875  typename VectorizedArrayType,
1876  typename Number2 = Number>
1877  struct FEFaceEvaluationSelector
1878  {
1879  static void
1880  evaluate(const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
1881  const VectorizedArrayType * values_array,
1882  VectorizedArrayType * values_quad,
1883  VectorizedArrayType * gradients_quad,
1884  VectorizedArrayType * scratch_data,
1885  const bool evaluate_values,
1886  const bool evaluate_gradients,
1887  const unsigned int face_no,
1888  const unsigned int subface_index,
1889  const unsigned int face_orientation,
1890  const Table<2, unsigned int> &orientation_map)
1891  {
1892  constexpr unsigned int static_dofs_per_face =
1893  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
1894  numbers::invalid_unsigned_int;
1895  const unsigned int dofs_per_face =
1896  fe_degree > -1 ?
1897  static_dofs_per_face :
1898  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
1899 
1900  // we allocate small amounts of data on the stack to signal the compiler
1901  // that this temporary data is only needed for the calculations but the
1902  // final results can be discarded and need not be written back to
1903  // memory. For large sizes or when the dofs per face is not a compile-time
1904  // constant, however, we want to go to the heap in the `scratch_data`
1905  // variable to not risk a stack overflow.
1906  constexpr unsigned int stack_array_size_threshold = 100;
1907 
1908  VectorizedArrayType
1909  temp_data[static_dofs_per_face < stack_array_size_threshold ?
1910  n_components * 2 * static_dofs_per_face :
1911  1];
1912  VectorizedArrayType *temp1;
1913  if (static_dofs_per_face < stack_array_size_threshold)
1914  temp1 = &temp_data[0];
1915  else
1916  temp1 = scratch_data;
1917 
1918  FEFaceNormalEvaluationImpl<dim,
1919  fe_degree,
1920  n_components,
1921  VectorizedArrayType>::
1922  template interpolate<true, false>(
1923  data, values_array, temp1, evaluate_gradients, face_no);
1924 
1925  const unsigned int n_q_points_1d_actual =
1926  fe_degree > -1 ? n_q_points_1d : 0;
1927  if (fe_degree > -1 &&
1928  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
1929  data.element_type <= MatrixFreeFunctions::tensor_symmetric)
1930  FEFaceEvaluationImpl<
1931  true,
1932  dim,
1933  fe_degree,
1934  n_q_points_1d_actual,
1935  n_components,
1936  VectorizedArrayType>::evaluate_in_face(data,
1937  temp1,
1938  values_quad,
1939  gradients_quad,
1940  scratch_data + 2 *
1941  n_components *
1942  dofs_per_face,
1943  evaluate_values,
1944  evaluate_gradients,
1945  subface_index);
1946  else
1947  FEFaceEvaluationImpl<
1948  false,
1949  dim,
1950  fe_degree,
1951  n_q_points_1d_actual,
1952  n_components,
1953  VectorizedArrayType>::evaluate_in_face(data,
1954  temp1,
1955  values_quad,
1956  gradients_quad,
1957  scratch_data + 2 *
1958  n_components *
1959  dofs_per_face,
1960  evaluate_values,
1961  evaluate_gradients,
1962  subface_index);
1963 
1964  if (face_orientation)
1965  adjust_for_face_orientation(face_orientation,
1966  orientation_map,
1967  false,
1968  evaluate_values,
1969  evaluate_gradients,
1970  data.n_q_points_face,
1971  scratch_data,
1972  values_quad,
1973  gradients_quad);
1974  }
1975 
1976  static void
1977  integrate(const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
1978  VectorizedArrayType * values_array,
1979  VectorizedArrayType * values_quad,
1980  VectorizedArrayType * gradients_quad,
1981  VectorizedArrayType * scratch_data,
1982  const bool integrate_values,
1983  const bool integrate_gradients,
1984  const unsigned int face_no,
1985  const unsigned int subface_index,
1986  const unsigned int face_orientation,
1987  const Table<2, unsigned int> &orientation_map)
1988  {
1989  if (face_orientation)
1990  adjust_for_face_orientation(face_orientation,
1991  orientation_map,
1992  true,
1993  integrate_values,
1994  integrate_gradients,
1995  data.n_q_points_face,
1996  scratch_data,
1997  values_quad,
1998  gradients_quad);
1999 
2000  constexpr unsigned int static_dofs_per_face =
2001  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
2002  numbers::invalid_unsigned_int;
2003  const unsigned int dofs_per_face =
2004  fe_degree > -1 ?
2005  static_dofs_per_face :
2006  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
2007 
2008  constexpr unsigned int stack_array_size_threshold = 100;
2009 
2010  VectorizedArrayType
2011  temp_data[static_dofs_per_face < stack_array_size_threshold ?
2012  n_components * 2 * static_dofs_per_face :
2013  1];
2014  VectorizedArrayType *temp1;
2015  if (static_dofs_per_face < stack_array_size_threshold)
2016  temp1 = &temp_data[0];
2017  else
2018  temp1 = scratch_data;
2019 
2020  const unsigned int n_q_points_1d_actual =
2021  fe_degree > -1 ? n_q_points_1d : 0;
2022  if (fe_degree > -1 &&
2023  subface_index >= GeometryInfo<dim - 1>::max_children_per_cell &&
2024  data.element_type <= MatrixFreeFunctions::tensor_symmetric)
2025  FEFaceEvaluationImpl<
2026  true,
2027  dim,
2028  fe_degree,
2029  n_q_points_1d_actual,
2030  n_components,
2031  VectorizedArrayType>::integrate_in_face(data,
2032  temp1,
2033  values_quad,
2034  gradients_quad,
2035  scratch_data +
2036  2 * n_components *
2037  dofs_per_face,
2038  integrate_values,
2039  integrate_gradients,
2040  subface_index);
2041  else
2042  FEFaceEvaluationImpl<
2043  false,
2044  dim,
2045  fe_degree,
2046  n_q_points_1d_actual,
2047  n_components,
2048  VectorizedArrayType>::integrate_in_face(data,
2049  temp1,
2050  values_quad,
2051  gradients_quad,
2052  scratch_data +
2053  2 * n_components *
2054  dofs_per_face,
2055  integrate_values,
2056  integrate_gradients,
2057  subface_index);
2058 
2059  FEFaceNormalEvaluationImpl<dim,
2060  fe_degree,
2061  n_components,
2062  VectorizedArrayType>::
2063  template interpolate<false, false>(
2064  data, temp1, values_array, integrate_gradients, face_no);
2065  }
2066 
2067  static bool
2068  gather_evaluate(
2069  const Number2 * src_ptr,
2070  const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
2071  const MatrixFreeFunctions::DoFInfo & dof_info,
2072  VectorizedArrayType * values_quad,
2073  VectorizedArrayType * gradients_quad,
2074  VectorizedArrayType * scratch_data,
2075  const bool evaluate_values,
2076  const bool evaluate_gradients,
2077  const unsigned int active_fe_index,
2078  const unsigned int first_selected_component,
2079  const unsigned int cell,
2080  const unsigned int face_no,
2081  const unsigned int subface_index,
2082  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
2083  const unsigned int face_orientation,
2084  const Table<2, unsigned int> & orientation_map)
2085  {
2086  const unsigned int side = face_no % 2;
2087 
2088  constexpr unsigned int static_dofs_per_component =
2089  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) :
2090  numbers::invalid_unsigned_int;
2091  constexpr unsigned int static_dofs_per_face =
2092  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
2093  numbers::invalid_unsigned_int;
2094  const unsigned int dofs_per_face =
2095  fe_degree > -1 ?
2096  static_dofs_per_face :
2097  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
2098 
2099  constexpr unsigned int stack_array_size_threshold = 100;
2100 
2101  VectorizedArrayType
2102  temp_data[static_dofs_per_face < stack_array_size_threshold ?
2103  n_components * 2 * dofs_per_face :
2104  1];
2105  VectorizedArrayType *__restrict temp1;
2106  if (static_dofs_per_face < stack_array_size_threshold)
2107  temp1 = &temp_data[0];
2108  else
2109  temp1 = scratch_data;
2110 
2111  // case 1: contiguous and interleaved indices
2112  if (((evaluate_gradients == false &&
2113  data.data.front().nodal_at_cell_boundaries == true) ||
2114  (data.element_type ==
2115  MatrixFreeFunctions::tensor_symmetric_hermite &&
2116  fe_degree > 1)) &&
2117  dof_info.index_storage_variants[dof_access_index][cell] ==
2118  MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2119  interleaved_contiguous)
2120  {
2121  AssertDimension(
2122  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2123  VectorizedArrayType::size());
2124  const unsigned int dof_index =
2125  dof_info
2126  .dof_indices_contiguous[dof_access_index]
2127  [cell * VectorizedArrayType::size()] +
2128  dof_info.component_dof_indices_offset[active_fe_index]
2129  [first_selected_component] *
2130  VectorizedArrayType::size();
2131 
2132  if (fe_degree > 1 && evaluate_gradients == true)
2133  {
2134  // we know that the gradient weights for the Hermite case on the
2135  // right (side==1) are the negative from the value at the left
2136  // (side==0), so we only read out one of them.
2137  const VectorizedArrayType grad_weight =
2138  data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
2139  AssertDimension(data.face_to_cell_index_hermite.size(1),
2140  2 * dofs_per_face);
2141  const unsigned int *index_array =
2142  &data.face_to_cell_index_hermite(face_no, 0);
2143  for (unsigned int i = 0; i < dofs_per_face; ++i)
2144  {
2145  const unsigned int ind1 = index_array[2 * i];
2146  const unsigned int ind2 = index_array[2 * i + 1];
2147  AssertIndexRange(ind1, data.dofs_per_component_on_cell);
2148  AssertIndexRange(ind2, data.dofs_per_component_on_cell);
2149  for (unsigned int comp = 0; comp < n_components; ++comp)
2150  {
2151  do_vectorized_read(src_ptr + dof_index +
2152  (ind1 +
2153  comp * static_dofs_per_component) *
2154  VectorizedArrayType::size(),
2155  temp1[i + 2 * comp * dofs_per_face]);
2156  do_vectorized_read(
2157  src_ptr + dof_index +
2158  (ind2 + comp * static_dofs_per_component) *
2159  VectorizedArrayType::size(),
2160  temp1[dofs_per_face + i + 2 * comp * dofs_per_face]);
2161  temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
2162  grad_weight *
2163  (temp1[i + 2 * comp * dofs_per_face] -
2164  temp1[i + dofs_per_face + 2 * comp * dofs_per_face]);
2165  }
2166  }
2167  }
2168  else
2169  {
2170  AssertDimension(data.face_to_cell_index_nodal.size(1),
2171  dofs_per_face);
2172  const unsigned int *index_array =
2173  &data.face_to_cell_index_nodal(face_no, 0);
2174  for (unsigned int i = 0; i < dofs_per_face; ++i)
2175  {
2176  const unsigned int ind = index_array[i];
2177  for (unsigned int comp = 0; comp < n_components; ++comp)
2178  do_vectorized_read(src_ptr + dof_index +
2179  (ind +
2180  comp * static_dofs_per_component) *
2181  VectorizedArrayType::size(),
2182  temp1[i + 2 * comp * dofs_per_face]);
2183  }
2184  }
2185  }
2186 
2187  // case 2: contiguous and interleaved indices with fixed stride
2188  else if (((evaluate_gradients == false &&
2189  data.data.front().nodal_at_cell_boundaries == true) ||
2190  (data.element_type ==
2191  MatrixFreeFunctions::tensor_symmetric_hermite &&
2192  fe_degree > 1)) &&
2193  dof_info.index_storage_variants[dof_access_index][cell] ==
2194  MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2195  interleaved_contiguous_strided)
2196  {
2197  AssertDimension(
2198  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2199  VectorizedArrayType::size());
2200  const unsigned int *indices =
2201  &dof_info
2202  .dof_indices_contiguous[dof_access_index]
2203  [cell * VectorizedArrayType::size()];
2204  if (fe_degree > 1 && evaluate_gradients == true)
2205  {
2206  // we know that the gradient weights for the Hermite case on the
2207  // right (side==1) are the negative from the value at the left
2208  // (side==0), so we only read out one of them.
2209  const VectorizedArrayType grad_weight =
2210  data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
2211  AssertDimension(data.face_to_cell_index_hermite.size(1),
2212  2 * dofs_per_face);
2213 
2214  const unsigned int *index_array =
2215  &data.face_to_cell_index_hermite(face_no, 0);
2216  for (unsigned int i = 0; i < dofs_per_face; ++i)
2217  {
2218  const unsigned int ind1 =
2219  index_array[2 * i] * VectorizedArrayType::size();
2220  const unsigned int ind2 =
2221  index_array[2 * i + 1] * VectorizedArrayType::size();
2222  for (unsigned int comp = 0; comp < n_components; ++comp)
2223  {
2224  do_vectorized_gather(
2225  src_ptr + ind1 +
2226  comp * static_dofs_per_component *
2227  VectorizedArrayType::size() +
2228  dof_info.component_dof_indices_offset
2229  [active_fe_index][first_selected_component] *
2230  VectorizedArrayType::size(),
2231  indices,
2232  temp1[i + 2 * comp * dofs_per_face]);
2233  VectorizedArrayType grad;
2234  do_vectorized_gather(
2235  src_ptr + ind2 +
2236  comp * static_dofs_per_component *
2237  VectorizedArrayType::size() +
2238  dof_info.component_dof_indices_offset
2239  [active_fe_index][first_selected_component] *
2240  VectorizedArrayType::size(),
2241  indices,
2242  grad);
2243  temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
2244  grad_weight *
2245  (temp1[i + 2 * comp * dofs_per_face] - grad);
2246  }
2247  }
2248  }
2249  else
2250  {
2251  AssertDimension(data.face_to_cell_index_nodal.size(1),
2252  dofs_per_face);
2253  const unsigned int *index_array =
2254  &data.face_to_cell_index_nodal(face_no, 0);
2255  for (unsigned int i = 0; i < dofs_per_face; ++i)
2256  {
2257  const unsigned int ind =
2258  index_array[i] * VectorizedArrayType::size();
2259  for (unsigned int comp = 0; comp < n_components; ++comp)
2260  do_vectorized_gather(
2261  src_ptr + ind +
2262  comp * static_dofs_per_component *
2263  VectorizedArrayType::size() +
2264  dof_info.component_dof_indices_offset
2265  [active_fe_index][first_selected_component] *
2266  VectorizedArrayType::size(),
2267  indices,
2268  temp1[i + 2 * comp * dofs_per_face]);
2269  }
2270  }
2271  }
2272 
2273  // case 3: contiguous and interleaved indices with mixed stride
2274  else if (((evaluate_gradients == false &&
2275  data.data.front().nodal_at_cell_boundaries == true) ||
2276  (data.element_type ==
2277  MatrixFreeFunctions::tensor_symmetric_hermite &&
2278  fe_degree > 1)) &&
2279  dof_info.index_storage_variants[dof_access_index][cell] ==
2280  MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2281  interleaved_contiguous_mixed_strides)
2282  {
2283  const unsigned int *strides =
2284  &dof_info.dof_indices_interleave_strides
2285  [dof_access_index][cell * VectorizedArrayType::size()];
2286  unsigned int indices[VectorizedArrayType::size()];
2287  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2288  indices[v] =
2289  dof_info.dof_indices_contiguous
2290  [dof_access_index][cell * VectorizedArrayType::size() + v] +
2291  dof_info.component_dof_indices_offset[active_fe_index]
2292  [first_selected_component] *
2293  strides[v];
2294  const unsigned int nvec =
2295  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
2296 
2297  if (fe_degree > 1 && evaluate_gradients == true)
2298  {
2299  // we know that the gradient weights for the Hermite case on the
2300  // right (side==1) are the negative from the value at the left
2301  // (side==0), so we only read out one of them.
2302  const VectorizedArrayType grad_weight =
2303  data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
2304  AssertDimension(data.face_to_cell_index_hermite.size(1),
2305  2 * dofs_per_face);
2306 
2307  const unsigned int *index_array =
2308  &data.face_to_cell_index_hermite(face_no, 0);
2309  if (nvec == VectorizedArrayType::size())
2310  for (unsigned int comp = 0; comp < n_components; ++comp)
2311  for (unsigned int i = 0; i < dofs_per_face; ++i)
2312  {
2313  unsigned int ind1[VectorizedArrayType::size()];
2314  DEAL_II_OPENMP_SIMD_PRAGMA
2315  for (unsigned int v = 0; v < VectorizedArrayType::size();
2316  ++v)
2317  ind1[v] =
2318  indices[v] + (comp * static_dofs_per_component +
2319  index_array[2 * i]) *
2320  strides[v];
2321  unsigned int ind2[VectorizedArrayType::size()];
2322  DEAL_II_OPENMP_SIMD_PRAGMA
2323  for (unsigned int v = 0; v < VectorizedArrayType::size();
2324  ++v)
2325  ind2[v] =
2326  indices[v] + (comp * static_dofs_per_component +
2327  index_array[2 * i + 1]) *
2328  strides[v];
2329  do_vectorized_gather(src_ptr,
2330  ind1,
2331  temp1[i + 2 * comp * dofs_per_face]);
2332  VectorizedArrayType grad;
2333  do_vectorized_gather(src_ptr, ind2, grad);
2334  temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
2335  grad_weight *
2336  (temp1[i + 2 * comp * dofs_per_face] - grad);
2337  }
2338  else
2339  {
2340  for (unsigned int i = 0; i < n_components * 2 * dofs_per_face;
2341  ++i)
2342  temp1[i] = VectorizedArrayType();
2343  for (unsigned int v = 0; v < nvec; ++v)
2344  for (unsigned int comp = 0; comp < n_components; ++comp)
2345  for (unsigned int i = 0; i < dofs_per_face; ++i)
2346  {
2347  const unsigned int ind1 =
2348  indices[v] + (comp * static_dofs_per_component +
2349  index_array[2 * i]) *
2350  strides[v];
2351  const unsigned int ind2 =
2352  indices[v] + (comp * static_dofs_per_component +
2353  index_array[2 * i + 1]) *
2354  strides[v];
2355  temp1[i + 2 * comp * dofs_per_face][v] =
2356  src_ptr[ind1];
2357  const Number grad = src_ptr[ind2];
2358  temp1[i + dofs_per_face +
2359  2 * comp * dofs_per_face][v] =
2360  grad_weight[0] *
2361  (temp1[i + 2 * comp * dofs_per_face][v] - grad);
2362  }
2363  }
2364  }
2365  else
2366  {
2367  AssertDimension(data.face_to_cell_index_nodal.size(1),
2368  dofs_per_face);
2369  const unsigned int *index_array =
2370  &data.face_to_cell_index_nodal(face_no, 0);
2371  if (nvec == VectorizedArrayType::size())
2372  for (unsigned int comp = 0; comp < n_components; ++comp)
2373  for (unsigned int i = 0; i < dofs_per_face; ++i)
2374  {
2375  unsigned int ind[VectorizedArrayType::size()];
2376  DEAL_II_OPENMP_SIMD_PRAGMA
2377  for (unsigned int v = 0; v < VectorizedArrayType::size();
2378  ++v)
2379  ind[v] =
2380  indices[v] +
2381  (comp * static_dofs_per_component + index_array[i]) *
2382  strides[v];
2383  do_vectorized_gather(src_ptr,
2384  ind,
2385  temp1[i + 2 * comp * dofs_per_face]);
2386  }
2387  else
2388  {
2389  for (unsigned int i = 0; i < n_components * dofs_per_face;
2390  ++i)
2391  temp1[i] = VectorizedArrayType();
2392  for (unsigned int v = 0; v < nvec; ++v)
2393  for (unsigned int comp = 0; comp < n_components; ++comp)
2394  for (unsigned int i = 0; i < dofs_per_face; ++i)
2395  {
2396  const unsigned int ind1 =
2397  indices[v] + (comp * static_dofs_per_component +
2398  index_array[i]) *
2399  strides[v];
2400  temp1[i + 2 * comp * dofs_per_face][v] =
2401  src_ptr[ind1];
2402  }
2403  }
2404  }
2405  }
2406 
2407  // case 4: contiguous indices without interleaving
2408  else if (((evaluate_gradients == false &&
2409  data.data.front().nodal_at_cell_boundaries == true) ||
2410  (data.element_type ==
2411  MatrixFreeFunctions::tensor_symmetric_hermite &&
2412  fe_degree > 1)) &&
2413  dof_info.index_storage_variants[dof_access_index][cell] ==
2414  MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2415  contiguous &&
2416  dof_info.n_vectorization_lanes_filled[dof_access_index][cell] ==
2417  VectorizedArrayType::size())
2418  {
2419  const unsigned int *indices =
2420  &dof_info
2421  .dof_indices_contiguous[dof_access_index]
2422  [cell * VectorizedArrayType::size()];
2423  if (evaluate_gradients == true &&
2424  data.element_type ==
2425  MatrixFreeFunctions::tensor_symmetric_hermite)
2426  {
2427  // we know that the gradient weights for the Hermite case on the
2428  // right (side==1) are the negative from the value at the left
2429  // (side==0), so we only read out one of them.
2430  const VectorizedArrayType grad_weight =
2431  data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
2432  AssertDimension(data.face_to_cell_index_hermite.size(1),
2433  2 * dofs_per_face);
2434 
2435  const unsigned int *index_array =
2436  &data.face_to_cell_index_hermite(face_no, 0);
2437  for (unsigned int i = 0; i < dofs_per_face; ++i)
2438  {
2439  const unsigned int ind1 = index_array[2 * i];
2440  const unsigned int ind2 = index_array[2 * i + 1];
2441  for (unsigned int comp = 0; comp < n_components; ++comp)
2442  {
2443  do_vectorized_gather(
2444  src_ptr + ind1 + comp * static_dofs_per_component +
2445  dof_info.component_dof_indices_offset
2446  [active_fe_index][first_selected_component],
2447  indices,
2448  temp1[i + 2 * comp * dofs_per_face]);
2449  VectorizedArrayType grad;
2450  do_vectorized_gather(
2451  src_ptr + ind2 + comp * static_dofs_per_component +
2452  dof_info.component_dof_indices_offset
2453  [active_fe_index][first_selected_component],
2454  indices,
2455  grad);
2456  temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
2457  grad_weight *
2458  (temp1[i + 2 * comp * dofs_per_face] - grad);
2459  }
2460  }
2461  }
2462  else
2463  {
2464  AssertDimension(data.face_to_cell_index_nodal.size(1),
2465  dofs_per_face);
2466  const unsigned int *index_array =
2467  &data.face_to_cell_index_nodal(face_no, 0);
2468  for (unsigned int i = 0; i < dofs_per_face; ++i)
2469  for (unsigned int comp = 0; comp < n_components; ++comp)
2470  {
2471  const unsigned int ind = index_array[i];
2472  do_vectorized_gather(
2473  src_ptr + ind + comp * static_dofs_per_component +
2474  dof_info.component_dof_indices_offset
2475  [active_fe_index][first_selected_component],
2476  indices,
2477  temp1[i + comp * 2 * dofs_per_face]);
2478  }
2479  }
2480  }
2481 
2482  // case 5: default vector access
2483  else
2484  {
2485  return false;
2486  }
2487 
2488  if (fe_degree > -1 &&
2489  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
2490  data.element_type <= MatrixFreeFunctions::tensor_symmetric)
2491  FEFaceEvaluationImpl<
2492  true,
2493  dim,
2494  fe_degree,
2495  n_q_points_1d,
2496  n_components,
2497  VectorizedArrayType>::evaluate_in_face(data,
2498  temp1,
2499  values_quad,
2500  gradients_quad,
2501  scratch_data + 2 *
2502  n_components *
2503  dofs_per_face,
2504  evaluate_values,
2505  evaluate_gradients,
2506  subface_index);
2507  else
2508  FEFaceEvaluationImpl<
2509  false,
2510  dim,
2511  fe_degree,
2512  n_q_points_1d,
2513  n_components,
2514  VectorizedArrayType>::evaluate_in_face(data,
2515  temp1,
2516  values_quad,
2517  gradients_quad,
2518  scratch_data + 2 *
2519  n_components *
2520  dofs_per_face,
2521  evaluate_values,
2522  evaluate_gradients,
2523  subface_index);
2524 
2525  if (face_orientation)
2526  adjust_for_face_orientation(face_orientation,
2527  orientation_map,
2528  false,
2529  evaluate_values,
2530  evaluate_gradients,
2531  data.n_q_points_face,
2532  scratch_data,
2533  values_quad,
2534  gradients_quad);
2535 
2536  return true;
2537  }
2538 
2539  static bool
2540  integrate_scatter(
2541  Number2 * dst_ptr,
2542  const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
2543  const MatrixFreeFunctions::DoFInfo & dof_info,
2544  VectorizedArrayType * values_array,
2545  VectorizedArrayType * values_quad,
2546  VectorizedArrayType * gradients_quad,
2547  VectorizedArrayType * scratch_data,
2548  const bool integrate_values,
2549  const bool integrate_gradients,
2550  const unsigned int active_fe_index,
2551  const unsigned int first_selected_component,
2552  const unsigned int cell,
2553  const unsigned int face_no,
2554  const unsigned int subface_index,
2555  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
2556  const unsigned int face_orientation,
2557  const Table<2, unsigned int> & orientation_map)
2558  {
2559  if (face_orientation)
2560  adjust_for_face_orientation(face_orientation,
2561  orientation_map,
2562  true,
2563  integrate_values,
2564  integrate_gradients,
2565  data.n_q_points_face,
2566  scratch_data,
2567  values_quad,
2568  gradients_quad);
2569 
2570  const unsigned int side = face_no % 2;
2571  constexpr unsigned int static_dofs_per_component =
2572  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) :
2573  numbers::invalid_unsigned_int;
2574  const unsigned int dofs_per_face =
2575  fe_degree > -1 ?
2576  Utilities::pow(fe_degree + 1, dim - 1) :
2577  Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
2578 
2579  constexpr unsigned int stack_array_size_threshold = 100;
2580 
2581  VectorizedArrayType temp_data[dofs_per_face < stack_array_size_threshold ?
2582  n_components * 2 * dofs_per_face :
2583  1];
2584  VectorizedArrayType *__restrict temp1;
2585  if (dofs_per_face < stack_array_size_threshold)
2586  temp1 = &temp_data[0];
2587  else
2588  temp1 = scratch_data;
2589 
2590  if (fe_degree > -1 &&
2591  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
2592  data.element_type <= internal::MatrixFreeFunctions::tensor_symmetric)
2593  internal::FEFaceEvaluationImpl<
2594  true,
2595  dim,
2596  fe_degree,
2597  n_q_points_1d,
2598  n_components,
2599  VectorizedArrayType>::integrate_in_face(data,
2600  temp1,
2601  values_quad,
2602  gradients_quad,
2603  scratch_data +
2604  2 * n_components *
2605  dofs_per_face,
2606  integrate_values,
2607  integrate_gradients,
2608  subface_index);
2609  else
2610  internal::FEFaceEvaluationImpl<
2611  false,
2612  dim,
2613  fe_degree,
2614  n_q_points_1d,
2615  n_components,
2616  VectorizedArrayType>::integrate_in_face(data,
2617  temp1,
2618  values_quad,
2619  gradients_quad,
2620  scratch_data +
2621  2 * n_components *
2622  dofs_per_face,
2623  integrate_values,
2624  integrate_gradients,
2625  subface_index);
2626 
2627  // case 1: contiguous and interleaved indices
2628  if (((integrate_gradients == false &&
2629  data.data.front().nodal_at_cell_boundaries == true) ||
2630  (data.element_type ==
2631  internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
2632  fe_degree > 1)) &&
2633  dof_info.index_storage_variants[dof_access_index][cell] ==
2634  internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2635  interleaved_contiguous)
2636  {
2637  AssertDimension(
2638  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2639  VectorizedArrayType::size());
2640  const unsigned int dof_index =
2641  dof_info
2642  .dof_indices_contiguous[dof_access_index]
2643  [cell * VectorizedArrayType::size()] +
2644  dof_info.component_dof_indices_offset[active_fe_index]
2645  [first_selected_component] *
2646  VectorizedArrayType::size();
2647 
2648  if (fe_degree > 1 && integrate_gradients == true)
2649  {
2650  // we know that the gradient weights for the Hermite case on the
2651  // right (side==1) are the negative from the value at the left
2652  // (side==0), so we only read out one of them.
2653  const VectorizedArrayType grad_weight =
2654  data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
2655  AssertDimension(data.face_to_cell_index_hermite.size(1),
2656  2 * dofs_per_face);
2657  const unsigned int *index_array =
2658  &data.face_to_cell_index_hermite(face_no, 0);
2659  for (unsigned int i = 0; i < dofs_per_face; ++i)
2660  {
2661  const unsigned int ind1 = index_array[2 * i];
2662  const unsigned int ind2 = index_array[2 * i + 1];
2663  AssertIndexRange(ind1, data.dofs_per_component_on_cell);
2664  AssertIndexRange(ind2, data.dofs_per_component_on_cell);
2665  for (unsigned int comp = 0; comp < n_components; ++comp)
2666  {
2667  VectorizedArrayType val =
2668  temp1[i + 2 * comp * dofs_per_face] -
2669  grad_weight *
2670  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
2671  VectorizedArrayType grad =
2672  grad_weight *
2673  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
2674  do_vectorized_add(val,
2675  dst_ptr + dof_index +
2676  (ind1 +
2677  comp * static_dofs_per_component) *
2678  VectorizedArrayType::size());
2679  do_vectorized_add(grad,
2680  dst_ptr + dof_index +
2681  (ind2 +
2682  comp * static_dofs_per_component) *
2683  VectorizedArrayType::size());
2684  }
2685  }
2686  }
2687  else
2688  {
2689  AssertDimension(data.face_to_cell_index_nodal.size(1),
2690  dofs_per_face);
2691  const unsigned int *index_array =
2692  &data.face_to_cell_index_nodal(face_no, 0);
2693  for (unsigned int i = 0; i < dofs_per_face; ++i)
2694  {
2695  const unsigned int ind = index_array[i];
2696  for (unsigned int comp = 0; comp < n_components; ++comp)
2697  do_vectorized_add(temp1[i + 2 * comp * dofs_per_face],
2698  dst_ptr + dof_index +
2699  (ind +
2700  comp * static_dofs_per_component) *
2701  VectorizedArrayType::size());
2702  }
2703  }
2704  }
2705 
2706  // case 2: contiguous and interleaved indices with fixed stride
2707  else if (((integrate_gradients == false &&
2708  data.data.front().nodal_at_cell_boundaries == true) ||
2709  (data.element_type ==
2710  internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
2711  fe_degree > 1)) &&
2712  dof_info.index_storage_variants[dof_access_index][cell] ==
2713  internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2714  interleaved_contiguous_strided)
2715  {
2716  AssertDimension(
2717  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2718  VectorizedArrayType::size());
2719  const unsigned int *indices =
2720  &dof_info
2721  .dof_indices_contiguous[dof_access_index]
2722  [cell * VectorizedArrayType::size()];
2723  if (fe_degree > 1 && integrate_gradients == true)
2724  {
2725  // we know that the gradient weights for the Hermite case on the
2726  // right (side==1) are the negative from the value at the left
2727  // (side==0), so we only read out one of them.
2728  const VectorizedArrayType grad_weight =
2729  data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
2730  AssertDimension(data.face_to_cell_index_hermite.size(1),
2731  2 * dofs_per_face);
2732 
2733  const unsigned int *index_array =
2734  &data.face_to_cell_index_hermite(face_no, 0);
2735  for (unsigned int i = 0; i < dofs_per_face; ++i)
2736  {
2737  const unsigned int ind1 =
2738  index_array[2 * i] * VectorizedArrayType::size();
2739  const unsigned int ind2 =
2740  index_array[2 * i + 1] * VectorizedArrayType::size();
2741  for (unsigned int comp = 0; comp < n_components; ++comp)
2742  {
2743  VectorizedArrayType val =
2744  temp1[i + 2 * comp * dofs_per_face] -
2745  grad_weight *
2746  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
2747  VectorizedArrayType grad =
2748  grad_weight *
2749  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
2750  do_vectorized_scatter_add(
2751  val,
2752  indices,
2753  dst_ptr + ind1 +
2754  comp * static_dofs_per_component *
2755  VectorizedArrayType::size() +
2756  dof_info.component_dof_indices_offset
2757  [active_fe_index][first_selected_component] *
2758  VectorizedArrayType::size());
2759  do_vectorized_scatter_add(
2760  grad,
2761  indices,
2762  dst_ptr + ind2 +
2763  comp * static_dofs_per_component *
2764  VectorizedArrayType::size() +
2765  dof_info.component_dof_indices_offset
2766  [active_fe_index][first_selected_component] *
2767  VectorizedArrayType::size());
2768  }
2769  }
2770  }
2771  else
2772  {
2773  AssertDimension(data.face_to_cell_index_nodal.size(1),
2774  dofs_per_face);
2775  const unsigned int *index_array =
2776  &data.face_to_cell_index_nodal(face_no, 0);
2777  for (unsigned int i = 0; i < dofs_per_face; ++i)
2778  {
2779  const unsigned int ind =
2780  index_array[i] * VectorizedArrayType::size();
2781  for (unsigned int comp = 0; comp < n_components; ++comp)
2782  do_vectorized_scatter_add(
2783  temp1[i + 2 * comp * dofs_per_face],
2784  indices,
2785  dst_ptr + ind +
2786  comp * static_dofs_per_component *
2787  VectorizedArrayType::size() +
2788  dof_info.component_dof_indices_offset
2789  [active_fe_index][first_selected_component] *
2790  VectorizedArrayType::size());
2791  }
2792  }
2793  }
2794 
2795  // case 3: contiguous and interleaved indices with mixed stride
2796  else if (((integrate_gradients == false &&
2797  data.data.front().nodal_at_cell_boundaries == true) ||
2798  (data.element_type ==
2799  internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
2800  fe_degree > 1)) &&
2801  dof_info.index_storage_variants[dof_access_index][cell] ==
2802  internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2803  interleaved_contiguous_mixed_strides)
2804  {
2805  const unsigned int *strides =
2806  &dof_info.dof_indices_interleave_strides
2807  [dof_access_index][cell * VectorizedArrayType::size()];
2808  unsigned int indices[VectorizedArrayType::size()];
2809  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2810  indices[v] =
2811  dof_info.dof_indices_contiguous
2812  [dof_access_index][cell * VectorizedArrayType::size() + v] +
2813  dof_info.component_dof_indices_offset[active_fe_index]
2814  [first_selected_component] *
2815  strides[v];
2816  const unsigned int n_filled_lanes =
2817  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
2818 
2819  if (fe_degree > 1 && integrate_gradients == true)
2820  {
2821  // we know that the gradient weights for the Hermite case on the
2822  // right (side==1) are the negative from the value at the left
2823  // (side==0), so we only read out one of them.
2824  const VectorizedArrayType grad_weight =
2825  data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
2826  AssertDimension(data.face_to_cell_index_hermite.size(1),
2827  2 * dofs_per_face);
2828 
2829  const unsigned int *index_array =
2830  &data.face_to_cell_index_hermite(face_no, 0);
2831  if (n_filled_lanes == VectorizedArrayType::size())
2832  for (unsigned int comp = 0; comp < n_components; ++comp)
2833  for (unsigned int i = 0; i < dofs_per_face; ++i)
2834  {
2835  unsigned int ind1[VectorizedArrayType::size()];
2836  DEAL_II_OPENMP_SIMD_PRAGMA
2837  for (unsigned int v = 0; v < VectorizedArrayType::size();
2838  ++v)
2839  ind1[v] =
2840  indices[v] + (comp * static_dofs_per_component +
2841  index_array[2 * i]) *
2842  strides[v];
2843  unsigned int ind2[VectorizedArrayType::size()];
2844  DEAL_II_OPENMP_SIMD_PRAGMA
2845  for (unsigned int v = 0; v < VectorizedArrayType::size();
2846  ++v)
2847  ind2[v] =
2848  indices[v] + (comp * static_dofs_per_component +
2849  index_array[2 * i + 1]) *
2850  strides[v];
2851  VectorizedArrayType val =
2852  temp1[i + 2 * comp * dofs_per_face] -
2853  grad_weight *
2854  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
2855  VectorizedArrayType grad =
2856  grad_weight *
2857  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
2858  do_vectorized_scatter_add(val, ind1, dst_ptr);
2859  do_vectorized_scatter_add(grad, ind2, dst_ptr);
2860  }
2861  else
2862  {
2863  for (unsigned int v = 0; v < n_filled_lanes; ++v)
2864  for (unsigned int comp = 0; comp < n_components; ++comp)
2865  for (unsigned int i = 0; i < dofs_per_face; ++i)
2866  {
2867  const unsigned int ind1 =
2868  indices[v] + (comp * static_dofs_per_component +
2869  index_array[2 * i]) *
2870  strides[v];
2871  const unsigned int ind2 =
2872  indices[v] + (comp * static_dofs_per_component +
2873  index_array[2 * i + 1]) *
2874  strides[v];
2875  Number val =
2876  temp1[i + 2 * comp * dofs_per_face][v] -
2877  grad_weight[0] * temp1[i + dofs_per_face +
2878  2 * comp * dofs_per_face][v];
2879  Number grad =
2880  grad_weight[0] * temp1[i + dofs_per_face +
2881  2 * comp * dofs_per_face][v];
2882  dst_ptr[ind1] += val;
2883  dst_ptr[ind2] += grad;
2884  }
2885  }
2886  }
2887  else
2888  {
2889  AssertDimension(data.face_to_cell_index_nodal.size(1),
2890  dofs_per_face);
2891  const unsigned int *index_array =
2892  &data.face_to_cell_index_nodal(face_no, 0);
2893  if (n_filled_lanes == VectorizedArrayType::size())
2894  for (unsigned int comp = 0; comp < n_components; ++comp)
2895  for (unsigned int i = 0; i < dofs_per_face; ++i)
2896  {
2897  unsigned int ind[VectorizedArrayType::size()];
2898  DEAL_II_OPENMP_SIMD_PRAGMA
2899  for (unsigned int v = 0; v < VectorizedArrayType::size();
2900  ++v)
2901  ind[v] =
2902  indices[v] +
2903  (comp * static_dofs_per_component + index_array[i]) *
2904  strides[v];
2905  do_vectorized_scatter_add(
2906  temp1[i + 2 * comp * dofs_per_face], ind, dst_ptr);
2907  }
2908  else
2909  {
2910  for (unsigned int v = 0; v < n_filled_lanes; ++v)
2911  for (unsigned int comp = 0; comp < n_components; ++comp)
2912  for (unsigned int i = 0; i < dofs_per_face; ++i)
2913  {
2914  const unsigned int ind1 =
2915  indices[v] + (comp * static_dofs_per_component +
2916  index_array[i]) *
2917  strides[v];
2918  dst_ptr[ind1] +=
2919  temp1[i + 2 * comp * dofs_per_face][v];
2920  }
2921  }
2922  }
2923  }
2924 
2925  // case 4: contiguous indices without interleaving
2926  else if (((integrate_gradients == false &&
2927  data.data.front().nodal_at_cell_boundaries == true) ||
2928  (data.element_type ==
2929  internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
2930  fe_degree > 1)) &&
2931  dof_info.index_storage_variants[dof_access_index][cell] ==
2932  internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2933  contiguous &&
2934  dof_info.n_vectorization_lanes_filled[dof_access_index][cell] ==
2935  VectorizedArrayType::size())
2936  {
2937  const unsigned int *indices =
2938  &dof_info
2939  .dof_indices_contiguous[dof_access_index]
2940  [cell * VectorizedArrayType::size()];
2941 
2942  if (integrate_gradients == true &&
2943  data.element_type ==
2944  internal::MatrixFreeFunctions::tensor_symmetric_hermite)
2945  {
2946  // we know that the gradient weights for the Hermite case on the
2947  // right (side==1) are the negative from the value at the left
2948  // (side==0), so we only read out one of them.
2949  const VectorizedArrayType grad_weight =
2950  data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
2951  AssertDimension(data.face_to_cell_index_hermite.size(1),
2952  2 * dofs_per_face);
2953  const unsigned int *index_array =
2954  &data.face_to_cell_index_hermite(face_no, 0);
2955  for (unsigned int i = 0; i < dofs_per_face; ++i)
2956  {
2957  const unsigned int ind1 = index_array[2 * i];
2958  const unsigned int ind2 = index_array[2 * i + 1];
2959  for (unsigned int comp = 0; comp < n_components; ++comp)
2960  {
2961  VectorizedArrayType val =
2962  temp1[i + 2 * comp * dofs_per_face] -
2963  grad_weight *
2964  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
2965  VectorizedArrayType grad =
2966  grad_weight *
2967  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
2968  do_vectorized_scatter_add(
2969  val,
2970  indices,
2971  dst_ptr + comp * static_dofs_per_component + ind1 +
2972  dof_info.component_dof_indices_offset
2973  [active_fe_index][first_selected_component]);
2974  do_vectorized_scatter_add(
2975  grad,
2976  indices,
2977  dst_ptr + comp * static_dofs_per_component + ind2 +
2978  dof_info.component_dof_indices_offset
2979  [active_fe_index][first_selected_component]);
2980  }
2981  }
2982  }
2983  else
2984  {
2985  AssertDimension(data.face_to_cell_index_nodal.size(1),
2986  dofs_per_face);
2987  const unsigned int *index_array =
2988  &data.face_to_cell_index_nodal(face_no, 0);
2989  for (unsigned int i = 0; i < dofs_per_face; ++i)
2990  {
2991  const unsigned int ind = index_array[i];
2992  for (unsigned int comp = 0; comp < n_components; ++comp)
2993  do_vectorized_scatter_add(
2994  temp1[i + 2 * comp * dofs_per_face],
2995  indices,
2996  dst_ptr + comp * static_dofs_per_component + ind +
2997  dof_info.component_dof_indices_offset
2998  [active_fe_index][first_selected_component]);
2999  }
3000  }
3001  }
3002 
3003  // case 5: default vector access, must be handled separately, just do
3004  // the face-normal interpolation
3005  else
3006  {
3007  FEFaceNormalEvaluationImpl<dim,
3008  fe_degree,
3009  n_components,
3010  VectorizedArrayType>::
3011  template interpolate<false, false>(
3012  data, temp1, values_array, integrate_gradients, face_no);
3013  return false;
3014  }
3015 
3016  return true;
3017  }
3018 
3019  static void
3020  adjust_for_face_orientation(const unsigned int face_orientation,
3021  const Table<2, unsigned int> &orientation_map,
3022  const bool integrate,
3023  const bool values,
3024  const bool gradients,
3025  const unsigned int n_q_points,
3026  VectorizedArrayType * tmp_values,
3027  VectorizedArrayType * values_quad,
3028  VectorizedArrayType * gradients_quad)
3029  {
3030  Assert(face_orientation, ExcInternalError());
3031  const unsigned int *orientation = &orientation_map[face_orientation][0];
3032  for (unsigned int c = 0; c < n_components; ++c)
3033  {
3034  if (values == true)
3035  {
3036  if (integrate)
3037  for (unsigned int q = 0; q < n_q_points; ++q)
3038  tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
3039  else
3040  for (unsigned int q = 0; q < n_q_points; ++q)
3041  tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
3042  for (unsigned int q = 0; q < n_q_points; ++q)
3043  values_quad[c * n_q_points + q] = tmp_values[q];
3044  }
3045  if (gradients == true)
3046  for (unsigned int d = 0; d < dim; ++d)
3047  {
3048  if (integrate)
3049  for (unsigned int q = 0; q < n_q_points; ++q)
3050  tmp_values[q] = gradients_quad[(c * dim + d) * n_q_points +
3051  orientation[q]];
3052  else
3053  for (unsigned int q = 0; q < n_q_points; ++q)
3054  tmp_values[orientation[q]] =
3055  gradients_quad[(c * dim + d) * n_q_points + q];
3056  for (unsigned int q = 0; q < n_q_points; ++q)
3057  gradients_quad[(c * dim + d) * n_q_points + q] =
3058  tmp_values[q];
3059  }
3060  }
3061  }
3062  };
3063 
3064 
3065 
3069  template <int dim, int fe_degree, int n_components, typename Number>
3070  struct CellwiseInverseMassMatrixImpl
3071  {
3072  template <typename FEEvaluationType>
3073  static void
3074  apply(const FEEvaluationType &fe_eval,
3075  const Number * in_array,
3076  Number * out_array)
3077  {
3078  constexpr unsigned int dofs_per_component =
3079  Utilities::pow(fe_degree + 1, dim);
3080 
3081  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
3082  Assert(fe_eval.get_shape_info().element_type <=
3083  MatrixFreeFunctions::tensor_symmetric,
3084  ExcNotImplemented());
3085 
3086  internal::EvaluatorTensorProduct<internal::evaluate_evenodd,
3087  dim,
3088  fe_degree + 1,
3089  fe_degree + 1,
3090  Number>
3091  evaluator(
3092  AlignedVector<Number>(),
3093  AlignedVector<Number>(),
3094  fe_eval.get_shape_info().data.front().inverse_shape_values_eo);
3095 
3096  for (unsigned int d = 0; d < n_components; ++d)
3097  {
3098  const Number *in = in_array + d * dofs_per_component;
3099  Number * out = out_array + d * dofs_per_component;
3100  // Need to select 'apply' method with hessian slot because values
3101  // assume symmetries that do not exist in the inverse shapes
3102  evaluator.template hessians<0, true, false>(in, out);
3103  if (dim > 1)
3104  evaluator.template hessians<1, true, false>(out, out);
3105  if (dim > 2)
3106  evaluator.template hessians<2, true, false>(out, out);
3107  }
3108  for (unsigned int q = 0; q < dofs_per_component; ++q)
3109  {
3110  const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
3111  for (unsigned int d = 0; d < n_components; ++d)
3112  out_array[q + d * dofs_per_component] *= inverse_JxW_q;
3113  }
3114  for (unsigned int d = 0; d < n_components; ++d)
3115  {
3116  Number *out = out_array + d * dofs_per_component;
3117  if (dim > 2)
3118  evaluator.template hessians<2, false, false>(out, out);
3119  if (dim > 1)
3120  evaluator.template hessians<1, false, false>(out, out);
3121  evaluator.template hessians<0, false, false>(out, out);
3122  }
3123  }
3124 
3125  static void
3126  apply(const AlignedVector<Number> &inverse_shape,
3127  const AlignedVector<Number> &inverse_coefficients,
3128  const unsigned int n_desired_components,
3129  const Number * in_array,
3130  Number * out_array)
3131  {
3132  constexpr unsigned int dofs_per_component =
3133  Utilities::pow(fe_degree + 1, dim);
3134  Assert(inverse_coefficients.size() > 0 &&
3135  inverse_coefficients.size() % dofs_per_component == 0,
3136  ExcMessage(
3137  "Expected diagonal to be a multiple of scalar dof per cells"));
3138  if (inverse_coefficients.size() != dofs_per_component)
3139  AssertDimension(n_desired_components * dofs_per_component,
3140  inverse_coefficients.size());
3141 
3142  Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
3143 
3144  internal::EvaluatorTensorProduct<internal::evaluate_evenodd,
3145  dim,
3146  fe_degree + 1,
3147  fe_degree + 1,
3148  Number>
3149  evaluator(AlignedVector<Number>(),
3150  AlignedVector<Number>(),
3151  inverse_shape);
3152 
3153  const unsigned int shift_coefficient =
3154  inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
3155  0;
3156  const Number *inv_coefficient = inverse_coefficients.data();
3157  for (unsigned int d = 0; d < n_desired_components; ++d)
3158  {
3159  const Number *in = in_array + d * dofs_per_component;
3160  Number * out = out_array + d * dofs_per_component;
3161  // Need to select 'apply' method with hessian slot because values
3162  // assume symmetries that do not exist in the inverse shapes
3163  evaluator.template hessians<0, true, false>(in, out);
3164  if (dim > 1)
3165  evaluator.template hessians<1, true, false>(out, out);
3166  if (dim > 2)
3167  evaluator.template hessians<2, true, false>(out, out);
3168 
3169  for (unsigned int q = 0; q < dofs_per_component; ++q)
3170  out[q] *= inv_coefficient[q];
3171 
3172  if (dim > 2)
3173  evaluator.template hessians<2, false, false>(out, out);
3174  if (dim > 1)
3175  evaluator.template hessians<1, false, false>(out, out);
3176  evaluator.template hessians<0, false, false>(out, out);
3177 
3178  inv_coefficient += shift_coefficient;
3179  }
3180  }
3181 
3182  static void
3183  transform_from_q_points_to_basis(const AlignedVector<Number> &inverse_shape,
3184  const unsigned int n_desired_components,
3185  const Number * in_array,
3186  Number * out_array)
3187  {
3188  constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim);
3189  internal::EvaluatorTensorProduct<internal::evaluate_evenodd,
3190  dim,
3191  fe_degree + 1,
3192  fe_degree + 1,
3193  Number>
3194  evaluator(AlignedVector<Number>(),
3195  AlignedVector<Number>(),
3196  inverse_shape);
3197 
3198  for (unsigned int d = 0; d < n_desired_components; ++d)
3199  {
3200  const Number *in = in_array + d * dofs_per_cell;
3201  Number * out = out_array + d * dofs_per_cell;
3202 
3203  if (dim == 3)
3204  {
3205  evaluator.template hessians<2, false, false>(in, out);
3206  evaluator.template hessians<1, false, false>(out, out);
3207  evaluator.template hessians<0, false, false>(out, out);
3208  }
3209  if (dim == 2)
3210  {
3211  evaluator.template hessians<1, false, false>(in, out);
3212  evaluator.template hessians<0, false, false>(out, out);
3213  }
3214  if (dim == 1)
3215  evaluator.template hessians<0, false, false>(in, out);
3216  }
3217  }
3218  };
3219 
3220 } // end of namespace internal
3221 
3222 
3223 DEAL_II_NAMESPACE_CLOSE
3224 
3225 #endif
internal::MatrixFreeFunctions::DoFInfo::dof_indices_contiguous
std::vector< unsigned int > dof_indices_contiguous[3]
Definition: dof_info.h:475
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:476
internal::MatrixFreeFunctions::ShapeInfo::data
std::vector< UnivariateShapeData< Number > > data
Definition: shape_info.h:355
internal::EvaluatorVariant
EvaluatorVariant
Definition: tensor_product_kernels.h:36
internal::MatrixFreeFunctions::ShapeInfo::element_type
ElementType element_type
Definition: shape_info.h:297
internal::do_vectorized_read
void do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
Definition: evaluation_kernels.h:1768
tensor_product_kernels.h
internal::MatrixFreeFunctions::DoFInfo::component_dof_indices_offset
std::vector< std::vector< unsigned int > > component_dof_indices_offset
Definition: dof_info.h:590
VectorizedArray::gather
void gather(const Number *base_ptr, const unsigned int *offsets)
Definition: vectorization.h:602
internal::FEFaceEvaluationSelector::gather_evaluate
static bool gather_evaluate(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 unsigned int cell, const unsigned int face_no, const unsigned int subface_index, const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, const unsigned int face_orientation, const Table< 2, unsigned int > &orientation_map)
Definition: evaluation_kernels.h:2068
internal::CellwiseInverseMassMatrixImpl::transform_from_q_points_to_basis
static void transform_from_q_points_to_basis(const AlignedVector< Number > &inverse_shape, const unsigned int n_desired_components, const Number *in_array, Number *out_array)
Definition: evaluation_kernels.h:3183
internal::EvaluatorTensorProduct
Definition: tensor_product_kernels.h:96
vectorization.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
internal::MatrixFreeFunctions::ShapeInfo::n_q_points
unsigned int n_q_points
Definition: shape_info.h:379
internal::MatrixFreeFunctions::ShapeInfo
Definition: shape_info.h:290
internal::FEFaceEvaluationImpl
Definition: evaluation_kernels.h:1353
internal::FEFaceEvaluationImpl::use_collocation
static constexpr bool use_collocation
Definition: evaluation_kernels.h:1360
internal::FEFaceEvaluationImpl::evaluate_in_face
static void evaluate_in_face(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)
Definition: evaluation_kernels.h:1366
internal::FEFaceNormalEvaluationImpl::interpolate
static void interpolate(const MatrixFreeFunctions::ShapeInfo< Number > &data, const Number *input, Number *output, const bool do_gradients, const unsigned int face_no)
Definition: evaluation_kernels.h:1707
utilities.h
internal::MatrixFreeFunctions::DoFInfo::dof_indices_interleave_strides
std::vector< unsigned int > dof_indices_interleave_strides[3]
Definition: dof_info.h:485
internal::MatrixFreeFunctions::ShapeInfo::face_to_cell_index_nodal
::Table< 2, unsigned int > face_to_cell_index_nodal
Definition: shape_info.h:432
AlignedVector::size
size_type size() const
internal::MatrixFreeFunctions::tensor_symmetric
@ tensor_symmetric
Definition: shape_info.h:76
GeometryInfo
Definition: geometry_info.h:1224
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
VectorizedArray
Definition: vectorization.h:395
internal::FEEvaluationImpl::integrate
static void integrate(const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs_actual, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool integrate_values, const bool integrate_gradients, const bool add_into_values_array)
Definition: evaluation_kernels.h:426
internal::MatrixFreeFunctions::truncated_tensor
@ truncated_tensor
Definition: shape_info.h:87
internal::FEEvaluationImplTransformToCollocation::evaluate
static void evaluate(const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians)
Definition: evaluation_kernels.h:1231
shape_info.h
DEAL_II_ALWAYS_INLINE
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:99
DoFTools::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
internal::MatrixFreeFunctions::DoFInfo
Definition: dof_info.h:99
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:140
internal::MatrixFreeFunctions::DoFInfo::n_vectorization_lanes_filled
std::vector< unsigned char > n_vectorization_lanes_filled[3]
Definition: dof_info.h:496
internal::FEFaceEvaluationSelector
Definition: evaluation_kernels.h:1877
Table< 2, unsigned int >
internal::evaluate_symmetric
@ evaluate_symmetric
Definition: tensor_product_kernels.h:48
MatrixFreeOperators::BlockHelper::n_blocks
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
Definition: operators.h:49
internal::MatrixFreeFunctions::ElementType
ElementType
Definition: shape_info.h:52
internal::FEFaceEvaluationSelector::integrate_scatter
static bool integrate_scatter(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 unsigned int cell, const unsigned int face_no, const unsigned int subface_index, const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, const unsigned int face_orientation, const Table< 2, unsigned int > &orientation_map)
Definition: evaluation_kernels.h:2540
internal::MatrixFreeFunctions::ShapeInfo::n_q_points_face
unsigned int n_q_points_face
Definition: shape_info.h:390
internal::FEEvaluationImplBasisChange::do_backward
static void do_backward(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)
Definition: evaluation_kernels.h:817
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
internal::MatrixFreeFunctions::DoFInfo::index_storage_variants
std::vector< IndexStorageVariants > index_storage_variants[3]
Definition: dof_info.h:423
dof_info.h
Utilities::fixed_power
T fixed_power(const T t)
Definition: utilities.h:1072
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
internal::CellwiseInverseMassMatrixImpl
Definition: evaluation_kernels.h:3070
internal::do_vectorized_gather
void do_vectorized_gather(const Number2 *src_ptr, const unsigned int *indices, VectorizedArrayType &dst)
Definition: evaluation_kernels.h:1790
internal::FEFaceEvaluationSelector::evaluate
static void evaluate(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)
Definition: evaluation_kernels.h:1880
AlignedVector
Definition: aligned_vector.h:61
internal::MatrixFreeFunctions::ShapeInfo::dofs_per_component_on_face
unsigned int dofs_per_component_on_face
Definition: shape_info.h:395
internal::evaluate_general
@ evaluate_general
Definition: tensor_product_kernels.h:42
internal::FEEvaluationImplBasisChange::do_forward
static void do_forward(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)
Definition: evaluation_kernels.h:698
internal::FEEvaluationImplTransformToCollocation
Definition: evaluation_kernels.h:1194
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
VectorizedArray::load
void load(const Number *ptr)
Definition: vectorization.h:517
internal::MatrixFreeFunctions::tensor_general
@ tensor_general
Definition: shape_info.h:81
internal::MatrixFreeFunctions::tensor_symmetric_hermite
@ tensor_symmetric_hermite
Definition: shape_info.h:69
internal::FEEvaluationImplCollocation
Definition: evaluation_kernels.h:1014
AlignedVector::data
pointer data()
internal::FEFaceEvaluationSelector::adjust_for_face_orientation
static void adjust_for_face_orientation(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, VectorizedArrayType *tmp_values, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad)
Definition: evaluation_kernels.h:3020
internal::FEFaceEvaluationImpl::integrate_in_face
static void integrate_in_face(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)
Definition: evaluation_kernels.h:1530
TableBase::size
size_type size(const unsigned int i) const
internal::FEEvaluationImpl::evaluate
static void evaluate(const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const Number *values_dofs_actual, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians)
Definition: evaluation_kernels.h:145
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
internal::FEEvaluationImplBasisChange
Definition: evaluation_kernels.h:668
internal::MatrixFreeFunctions::tensor_symmetric_collocation
@ tensor_symmetric_collocation
Definition: shape_info.h:62
internal::FEEvaluationImpl
Definition: evaluation_kernels.h:111
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
internal::do_vectorized_add
void do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
Definition: evaluation_kernels.h:1816
internal::FEEvaluationImplTransformToCollocation::integrate
static void integrate(const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool integrate_values, const bool integrate_gradients, const bool add_into_values_array)
Definition: evaluation_kernels.h:1294
internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants
IndexStorageVariants
Definition: dof_info.h:296
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
internal::FEEvaluationImplCollocation::evaluate
static void evaluate(const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians)
Definition: evaluation_kernels.h:1042
internal::CellwiseInverseMassMatrixImpl::apply
static void apply(const AlignedVector< Number > &inverse_shape, const AlignedVector< Number > &inverse_coefficients, const unsigned int n_desired_components, const Number *in_array, Number *out_array)
Definition: evaluation_kernels.h:3126
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
internal::do_vectorized_scatter_add
void do_vectorized_scatter_add(const VectorizedArrayType src, const unsigned int *indices, Number2 *dst_ptr)
Definition: evaluation_kernels.h:1840
config.h
internal::CellwiseInverseMassMatrixImpl::apply
static void apply(const FEEvaluationType &fe_eval, const Number *in_array, Number *out_array)
Definition: evaluation_kernels.h:3074
internal::MatrixFreeFunctions::ShapeInfo::face_to_cell_index_hermite
::Table< 2, unsigned int > face_to_cell_index_hermite
Definition: shape_info.h:472
internal
Definition: aligned_vector.h:369
internal::FEEvaluationImplCollocation::integrate
static void integrate(const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool integrate_values, const bool integrate_gradients, const bool add_into_values_array)
Definition: evaluation_kernels.h:1120
internal::FEFaceEvaluationSelector::integrate
static void integrate(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)
Definition: evaluation_kernels.h:1977
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
internal::evaluate_evenodd
@ evaluate_evenodd
Definition: tensor_product_kernels.h:54
internal::MatrixFreeFunctions::ShapeInfo::dofs_per_component_on_cell
unsigned int dofs_per_component_on_cell
Definition: shape_info.h:385
internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex
DoFAccessIndex
Definition: dof_info.h:383
internal::EvaluatorSelector
Definition: evaluation_kernels.h:38
internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0
@ tensor_symmetric_plus_dg0
Definition: shape_info.h:94
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
internal::FEEvaluationImplBasisChange::do_mass
static void do_mass(const AlignedVector< Number2 > &transformation_matrix, const AlignedVector< Number > &coefficients, const Number *values_in, Number *scratch_data, Number *values_out)
Definition: evaluation_kernels.h:923
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
internal::FEFaceNormalEvaluationImpl
Definition: evaluation_kernels.h:1703
Generated by   doxygen 1.8.17