Reference documentation for deal.II version GIT 9042b9283b 2023-12-02 14:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
evaluation_kernels_face.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2023 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_face_h
18 #define dealii_matrix_free_evaluation_kernels_face_h
19 
20 #include <deal.II/base/config.h>
21 
23 #include <deal.II/base/ndarray.h>
24 #include <deal.II/base/utilities.h>
26 
32 
33 
35 
36 
37 namespace internal
38 {
39  template <bool symmetric_evaluate,
40  int dim,
41  int fe_degree,
42  int n_q_points_1d,
43  typename Number>
45  {
46  // We enable a transformation to collocation for derivatives if it gives
47  // correct results (first two conditions), if it is the most efficient
48  // choice in terms of operation counts (third condition) and if we were
49  // able to initialize the fields in shape_info.templates.h from the
50  // polynomials (fourth condition).
51  using Number2 =
53 
54  using Eval = EvaluatorTensorProduct<symmetric_evaluate ? evaluate_evenodd :
56  dim - 1,
57  fe_degree + 1,
58  n_q_points_1d,
59  Number,
60  Number2>;
61 
62  static Eval
65  const unsigned int subface_index,
66  const unsigned int direction)
67  {
68  if (symmetric_evaluate)
69  return Eval(data.shape_values_eo,
70  data.shape_gradients_eo,
71  data.shape_hessians_eo,
72  data.fe_degree + 1,
73  data.n_q_points_1d);
74  else if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
75  return Eval(data.shape_values,
76  data.shape_gradients,
77  data.shape_hessians,
78  data.fe_degree + 1,
79  data.n_q_points_1d);
80  else
81  {
82  const unsigned int index =
83  direction == 0 ? subface_index % 2 : subface_index / 2;
84  return Eval(data.values_within_subface[index],
87  data.fe_degree + 1,
88  data.n_q_points_1d);
89  }
90  }
91 
92  static void
94  const unsigned int n_components,
95  const EvaluationFlags::EvaluationFlags evaluation_flag,
97  Number *values_dofs,
98  Number *values_quad,
99  Number *gradients_quad,
100  Number *hessians_quad,
101  Number *scratch_data,
102  const unsigned int subface_index)
103  {
104  Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
105  Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
106 
107  const std::size_t n_dofs = fe_degree > -1 ?
108  Utilities::pow(fe_degree + 1, dim - 1) :
109  Utilities::pow(data.fe_degree + 1, dim - 1);
110  const std::size_t n_q_points =
111  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
112  Utilities::pow(data.n_q_points_1d, dim - 1);
113 
114  // keep a copy of the original pointer for the case of the Hessians
115  Number *values_dofs_ptr = values_dofs;
116 
117  if ((evaluation_flag & EvaluationFlags::values) != 0u &&
118  ((evaluation_flag & EvaluationFlags::gradients) == 0u))
119  for (unsigned int c = 0; c < n_components; ++c)
120  {
121  switch (dim)
122  {
123  case 3:
124  eval0.template values<0, true, false>(values_dofs,
125  values_quad);
126  eval1.template values<1, true, false>(values_quad,
127  values_quad);
128  break;
129  case 2:
130  eval0.template values<0, true, false>(values_dofs,
131  values_quad);
132  break;
133  case 1:
134  values_quad[0] = values_dofs[0];
135  break;
136  default:
137  Assert(false, ExcNotImplemented());
138  }
139  // Note: we always keep storage of values, 1st and 2nd derivatives
140  // in an array
141  values_dofs += 3 * n_dofs;
142  values_quad += n_q_points;
143  }
144  else if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
145  for (unsigned int c = 0; c < n_components; ++c)
146  {
147  switch (dim)
148  {
149  case 3:
150  if (symmetric_evaluate &&
151  use_collocation_evaluation(fe_degree, n_q_points_1d))
152  {
153  eval0.template values<0, true, false>(values_dofs,
154  values_quad);
155  eval0.template values<1, true, false>(values_quad,
156  values_quad);
158  dim - 1,
159  n_q_points_1d,
160  n_q_points_1d,
161  Number,
162  Number2>
163  eval_grad({}, data.shape_gradients_collocation_eo, {});
164  eval_grad.template gradients<0, true, false, 3>(
165  values_quad, gradients_quad);
166  eval_grad.template gradients<1, true, false, 3>(
167  values_quad, gradients_quad + 1);
168  }
169  else
170  {
171  // grad x
172  eval0.template gradients<0, true, false>(values_dofs,
173  scratch_data);
174  eval1.template values<1, true, false, 3>(scratch_data,
175  gradients_quad);
176 
177  // grad y
178  eval0.template values<0, true, false>(values_dofs,
179  scratch_data);
180  eval1.template gradients<1, true, false, 3>(
181  scratch_data, gradients_quad + 1);
182 
183  if ((evaluation_flag & EvaluationFlags::values) != 0u)
184  eval1.template values<1, true, false>(scratch_data,
185  values_quad);
186  }
187  // grad z
188  eval0.template values<0, true, false>(values_dofs + n_dofs,
189  scratch_data);
190  eval1.template values<1, true, false, 3>(scratch_data,
191  gradients_quad + 2);
192 
193  break;
194  case 2:
195  eval0.template values<0, true, false, 2>(values_dofs + n_dofs,
196  gradients_quad + 1);
197  eval0.template gradients<0, true, false, 2>(values_dofs,
198  gradients_quad);
199  if ((evaluation_flag & EvaluationFlags::values) != 0u)
200  eval0.template values<0, true, false>(values_dofs,
201  values_quad);
202  break;
203  case 1:
204  values_quad[0] = values_dofs[0];
205  gradients_quad[0] = values_dofs[1];
206  break;
207  default:
208  AssertThrow(false, ExcNotImplemented());
209  }
210  values_dofs += 3 * n_dofs;
211  values_quad += n_q_points;
212  gradients_quad += dim * n_q_points;
213  }
214 
215  if ((evaluation_flag & EvaluationFlags::hessians) != 0u)
216  {
217  values_dofs = values_dofs_ptr;
218  for (unsigned int c = 0; c < n_components; ++c)
219  {
220  switch (dim)
221  {
222  case 3:
223  // grad xx
224  eval0.template hessians<0, true, false>(values_dofs,
225  scratch_data);
226  eval1.template values<1, true, false>(scratch_data,
227  hessians_quad);
228 
229  // grad yy
230  eval0.template values<0, true, false>(values_dofs,
231  scratch_data);
232  eval1.template hessians<1, true, false>(scratch_data,
233  hessians_quad +
234  n_q_points);
235 
236  // grad zz
237  eval0.template values<0, true, false>(values_dofs +
238  2 * n_dofs,
239  scratch_data);
240  eval1.template values<1, true, false>(scratch_data,
241  hessians_quad +
242  2 * n_q_points);
243 
244  // grad xy
245  eval0.template gradients<0, true, false>(values_dofs,
246  scratch_data);
247  eval1.template gradients<1, true, false>(scratch_data,
248  hessians_quad +
249  3 * n_q_points);
250 
251  // grad xz
252  eval0.template gradients<0, true, false>(values_dofs +
253  n_dofs,
254  scratch_data);
255  eval1.template values<1, true, false>(scratch_data,
256  hessians_quad +
257  4 * n_q_points);
258 
259  // grad yz
260  eval0.template values<0, true, false>(values_dofs + n_dofs,
261  scratch_data);
262  eval1.template gradients<1, true, false>(scratch_data,
263  hessians_quad +
264  5 * n_q_points);
265 
266  break;
267  case 2:
268  // grad xx
269  eval0.template hessians<0, true, false>(values_dofs,
270  hessians_quad);
271  // grad yy
272  eval0.template values<0, true, false>(
273  values_dofs + 2 * n_dofs, hessians_quad + n_q_points);
274  // grad xy
275  eval0.template gradients<0, true, false>(
276  values_dofs + n_dofs, hessians_quad + 2 * n_q_points);
277  break;
278  case 1:
279  hessians_quad[0] = values_dofs[2];
280  break;
281  default:
282  AssertThrow(false, ExcNotImplemented());
283  }
284  values_dofs += 3 * n_dofs;
285  hessians_quad += dim * (dim + 1) / 2 * n_q_points;
286  }
287  }
288  }
289 
290  static void
292  const unsigned int n_components,
293  const EvaluationFlags::EvaluationFlags integration_flag,
295  Number *values_dofs,
296  Number *values_quad,
297  Number *gradients_quad,
298  Number *hessians_quad,
299  Number *scratch_data,
300  const unsigned int subface_index)
301  {
302  Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
303  Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
304 
305  const std::size_t n_dofs =
306  fe_degree > -1 ?
307  Utilities::pow(fe_degree + 1, dim - 1) :
308  (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
309  const std::size_t n_q_points =
310  fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
311  Utilities::pow(data.n_q_points_1d, dim - 1);
312 
313  // keep a copy of the original pointer for the case of the Hessians
314  Number *values_dofs_ptr = values_dofs;
315 
316  if ((integration_flag & EvaluationFlags::values) != 0u &&
317  (integration_flag & EvaluationFlags::gradients) == 0u)
318  for (unsigned int c = 0; c < n_components; ++c)
319  {
320  switch (dim)
321  {
322  case 3:
323  eval1.template values<1, false, false>(values_quad,
324  values_quad);
325  eval0.template values<0, false, false>(values_quad,
326  values_dofs);
327  break;
328  case 2:
329  eval0.template values<0, false, false>(values_quad,
330  values_dofs);
331  break;
332  case 1:
333  values_dofs[0] = values_quad[0];
334  break;
335  default:
336  Assert(false, ExcNotImplemented());
337  }
338  values_dofs += 3 * n_dofs;
339  values_quad += n_q_points;
340  }
341  else if ((integration_flag & EvaluationFlags::gradients) != 0u)
342  for (unsigned int c = 0; c < n_components; ++c)
343  {
344  switch (dim)
345  {
346  case 3:
347  // grad z
348  eval1.template values<1, false, false, 3>(gradients_quad + 2,
349  scratch_data);
350  eval0.template values<0, false, false>(scratch_data,
351  values_dofs + n_dofs);
352  if (symmetric_evaluate &&
353  use_collocation_evaluation(fe_degree, n_q_points_1d))
354  {
356  dim - 1,
357  n_q_points_1d,
358  n_q_points_1d,
359  Number,
360  Number2>
361  eval_grad({}, data.shape_gradients_collocation_eo, {});
362  if ((integration_flag & EvaluationFlags::values) != 0u)
363  eval_grad.template gradients<1, false, true, 3>(
364  gradients_quad + 1, values_quad);
365  else
366  eval_grad.template gradients<1, false, false, 3>(
367  gradients_quad + 1, values_quad);
368  eval_grad.template gradients<0, false, true, 3>(
369  gradients_quad, values_quad);
370  eval0.template values<1, false, false>(values_quad,
371  values_quad);
372  eval0.template values<0, false, false>(values_quad,
373  values_dofs);
374  }
375  else
376  {
377  if ((integration_flag & EvaluationFlags::values) != 0u)
378  {
379  eval1.template values<1, false, false>(values_quad,
380  scratch_data);
381  eval1.template gradients<1, false, true, 3>(
382  gradients_quad + 1, scratch_data);
383  }
384  else
385  eval1.template gradients<1, false, false, 3>(
386  gradients_quad + 1, scratch_data);
387 
388  // grad y
389  eval0.template values<0, false, false>(scratch_data,
390  values_dofs);
391 
392  // grad x
393  eval1.template values<1, false, false, 3>(gradients_quad,
394  scratch_data);
395  eval0.template gradients<0, false, true>(scratch_data,
396  values_dofs);
397  }
398  break;
399  case 2:
400  eval0.template values<0, false, false, 2>(gradients_quad + 1,
401  values_dofs +
402  n_dofs);
403  eval0.template gradients<0, false, false, 2>(gradients_quad,
404  values_dofs);
405  if ((integration_flag & EvaluationFlags::values) != 0u)
406  eval0.template values<0, false, true>(values_quad,
407  values_dofs);
408  break;
409  case 1:
410  values_dofs[0] = values_quad[0];
411  values_dofs[1] = gradients_quad[0];
412  break;
413  default:
414  AssertThrow(false, ExcNotImplemented());
415  }
416  values_dofs += 3 * n_dofs;
417  values_quad += n_q_points;
418  gradients_quad += dim * n_q_points;
419  }
420 
421  if ((integration_flag & EvaluationFlags::hessians) != 0u)
422  {
423  values_dofs = values_dofs_ptr;
424  for (unsigned int c = 0; c < n_components; ++c)
425  {
426  switch (dim)
427  {
428  case 3:
429  // grad xx
430  eval1.template values<1, false, false>(hessians_quad,
431  scratch_data);
432  if ((integration_flag & (EvaluationFlags::values |
434  eval0.template hessians<0, false, true>(scratch_data,
435  values_dofs);
436  else
437  eval0.template hessians<0, false, false>(scratch_data,
438  values_dofs);
439 
440  // grad yy
441  eval1.template hessians<1, false, false>(hessians_quad +
442  n_q_points,
443  scratch_data);
444  eval0.template values<0, false, true>(scratch_data,
445  values_dofs);
446 
447  // grad zz
448  eval1.template values<1, false, false>(hessians_quad +
449  2 * n_q_points,
450  scratch_data);
451  eval0.template values<0, false, false>(scratch_data,
452  values_dofs +
453  2 * n_dofs);
454 
455  // grad xy
456  eval1.template gradients<1, false, false>(hessians_quad +
457  3 * n_q_points,
458  scratch_data);
459  eval0.template gradients<0, false, true>(scratch_data,
460  values_dofs);
461 
462  // grad xz
463  eval1.template values<1, false, false>(hessians_quad +
464  4 * n_q_points,
465  scratch_data);
466  if ((integration_flag & EvaluationFlags::gradients) != 0u)
467  eval0.template gradients<0, false, true>(scratch_data,
468  values_dofs +
469  n_dofs);
470  else
471  eval0.template gradients<0, false, false>(scratch_data,
472  values_dofs +
473  n_dofs);
474 
475  // grad yz
476  eval1.template gradients<1, false, false>(hessians_quad +
477  5 * n_q_points,
478  scratch_data);
479  eval0.template values<0, false, true>(scratch_data,
480  values_dofs + n_dofs);
481 
482  break;
483  case 2:
484  // grad xx
485  if ((integration_flag & (EvaluationFlags::values |
487  eval0.template hessians<0, false, true>(hessians_quad,
488  values_dofs);
489  else
490  eval0.template hessians<0, false, false>(hessians_quad,
491  values_dofs);
492 
493  // grad yy
494  eval0.template values<0, false, false>(
495  hessians_quad + n_q_points, values_dofs + 2 * n_dofs);
496  // grad xy
497  if ((integration_flag & EvaluationFlags::gradients) != 0u)
498  eval0.template gradients<0, false, true>(
499  hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
500  else
501  eval0.template gradients<0, false, false>(
502  hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
503  break;
504  case 1:
505  values_dofs[2] = hessians_quad[0];
506  if ((integration_flag & EvaluationFlags::values) == 0u)
507  values_dofs[0] = 0;
508  if ((integration_flag & EvaluationFlags::gradients) == 0u)
509  values_dofs[1] = 0;
510  break;
511  default:
512  AssertThrow(false, ExcNotImplemented());
513  }
514  values_dofs += 3 * n_dofs;
515  hessians_quad += dim * (dim + 1) / 2 * n_q_points;
516  }
517  }
518  }
519  };
520 
521 
522 
523  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
525  {
526  using Number2 =
528 
533  template <bool do_integrate>
534  static inline void
536  const EvaluationFlags::EvaluationFlags evaluation_flag,
538  &shape_data,
539  Number *values_dofs_in,
540  Number *values,
541  Number *gradients,
542  Number *scratch_data,
543  const unsigned int subface_index,
544  const unsigned int face_direction)
545  {
546  AssertDimension(shape_data.size(), 2);
547 
548  const int degree = fe_degree != -1 ? fe_degree : shape_data[0].fe_degree;
549  const int n_rows_n = degree + 1;
550  const int n_rows_t = degree;
551  const ::ndarray<int, 3, 3> dofs_per_direction{
552  {{{n_rows_n, n_rows_t, n_rows_t}},
553  {{n_rows_t, n_rows_n, n_rows_t}},
554  {{n_rows_t, n_rows_t, n_rows_n}}}};
555 
556  (void)scratch_data;
557  (void)subface_index;
558  // TODO: This is currently not implemented, but the test
559  // matrix_vector_rt_face_03 apparently works without it -> check
560  // if (subface_index < GeometryInfo<dim - 1>::max_children_per_cell)
561  // Assert(false, ExcNotImplemented());
562 
564  dim - 1,
565  (fe_degree > 0 ? fe_degree : 0),
566  n_q_points_1d,
567  Number,
568  Number2>;
569 
570  std::array<int, dim> values_dofs_offsets = {};
571  for (unsigned int comp = 0; comp < dim - 1; ++comp)
572  {
573  if (dim == 2)
574  values_dofs_offsets[comp + 1] =
575  values_dofs_offsets[comp] +
576  3 * dofs_per_direction[comp][(face_direction + 1) % dim];
577  else
578  values_dofs_offsets[comp + 1] =
579  values_dofs_offsets[comp] +
580  3 * dofs_per_direction[comp][(face_direction + 1) % dim] *
581  dofs_per_direction[comp][(face_direction + 2) % dim];
582  }
583 
584  // Jacobians on faces are reordered to enable simple access with the
585  // regular evaluators; to get the RT Piola transform right, we need to
586  // pass through the values_dofs array in a permuted right order
587  std::array<unsigned int, dim> components;
588  for (unsigned int comp = 0; comp < dim; ++comp)
589  components[comp] = (face_direction + comp + 1) % dim;
590 
591  for (const unsigned int comp : components)
592  {
593  Number *values_dofs = values_dofs_in + values_dofs_offsets[comp];
594 
595  std::array<int, 2> n_blocks{
596  {dofs_per_direction[comp][(face_direction + 1) % dim],
597  (dim > 2 ? dofs_per_direction[comp][(face_direction + 2) % dim] :
598  1)}};
599 
600  if constexpr (dim == 3)
601  {
603  dim - 1,
604  n_q_points_1d,
605  n_q_points_1d,
606  Number,
607  Number2>
608  eval_g({},
609  shape_data[0].shape_gradients_collocation_eo.data(),
610  {});
611  if (!do_integrate)
612  {
614  fe_degree,
615  n_q_points_1d,
616  true>
617  eval;
618  // Evaluate in 3d
619  if (n_blocks[0] == n_rows_n)
620  {
621  eval.template normal<0>(shape_data[0],
622  values_dofs,
623  values);
624  eval.template tangential<1, 0>(shape_data[1],
625  values,
626  values);
627 
628  if (evaluation_flag & EvaluationFlags::gradients)
629  {
630  eval.template normal<0>(shape_data[0],
631  values_dofs +
632  n_blocks[0] * n_blocks[1],
633  scratch_data);
634  eval.template tangential<1, 0, dim>(shape_data[1],
635  scratch_data,
636  gradients + 2);
637  }
638  }
639  else if (n_blocks[1] == n_rows_n)
640  {
641  eval.template normal<1>(shape_data[0],
642  values_dofs,
643  values);
644  eval.template tangential<0, 1>(shape_data[1],
645  values,
646  values);
647 
648  if (evaluation_flag & EvaluationFlags::gradients)
649  {
650  eval.template normal<1>(shape_data[0],
651  values_dofs +
652  n_blocks[0] * n_blocks[1],
653  scratch_data);
654  eval.template tangential<0, 1, dim>(shape_data[1],
655  scratch_data,
656  gradients + 2);
657  }
658  }
659  else
660  {
661  Eval eval(shape_data[1].shape_values_eo.data(), {}, {});
662  eval.template values<0, true, false>(values_dofs, values);
663  eval.template values<1, true, false>(values, values);
664  if (evaluation_flag & EvaluationFlags::gradients)
665  {
666  eval.template values<0, true, false>(values_dofs +
667  n_blocks[0] *
668  n_blocks[1],
669  scratch_data);
670  eval.template values<1, true, false, dim>(
671  scratch_data, gradients + 2);
672  }
673  }
674  if (evaluation_flag & EvaluationFlags::gradients)
675  {
676  eval_g.template gradients<0, true, false, dim>(values,
677  gradients);
678  eval_g.template gradients<1, true, false, dim>(values,
679  gradients +
680  1);
681  }
682  }
683  else
684  {
686  fe_degree,
687  n_q_points_1d,
688  false>
689  eval;
690  // Integrate in 3d
691  if (evaluation_flag & EvaluationFlags::gradients)
692  {
693  if (evaluation_flag & EvaluationFlags::values)
694  eval_g.template gradients<0, false, true, dim>(
695  gradients, values);
696  else
697  eval_g.template gradients<0, false, false, dim>(
698  gradients, values);
699  eval_g.template gradients<1, false, true, dim>(gradients +
700  1,
701  values);
702  }
703  if (n_blocks[0] == n_rows_n)
704  {
705  eval.template tangential<1, 0>(shape_data[1],
706  values,
707  values);
708  eval.template normal<0>(shape_data[0],
709  values,
710  values_dofs);
711 
712  if (evaluation_flag & EvaluationFlags::gradients)
713  {
714  eval.template tangential<1, 0, dim>(shape_data[1],
715  gradients + 2,
716  scratch_data);
717  eval.template normal<0>(shape_data[0],
718  scratch_data,
719  values_dofs +
720  n_blocks[0] * n_blocks[1]);
721  }
722  }
723  else if (n_blocks[1] == n_rows_n)
724  {
725  eval.template tangential<0, 1>(shape_data[1],
726  values,
727  values);
728  eval.template normal<1>(shape_data[0],
729  values,
730  values_dofs);
731 
732  if (evaluation_flag & EvaluationFlags::gradients)
733  {
734  eval.template tangential<0, 1, dim>(shape_data[1],
735  gradients + 2,
736  scratch_data);
737  eval.template normal<1>(shape_data[0],
738  scratch_data,
739  values_dofs +
740  n_blocks[0] * n_blocks[1]);
741  }
742  }
743  else
744  {
745  Eval eval_iso(shape_data[1].shape_values_eo.data(),
746  {},
747  {});
748  eval_iso.template values<1, false, false>(values, values);
749  eval_iso.template values<0, false, false>(values,
750  values_dofs);
751  if (evaluation_flag & EvaluationFlags::gradients)
752  {
753  eval_iso.template values<1, false, false, dim>(
754  gradients + 2, scratch_data);
755  eval_iso.template values<0, false, false>(
756  scratch_data,
757  values_dofs + n_blocks[0] * n_blocks[1]);
758  }
759  }
760  }
761  }
762  else
763  {
765  dim - 1,
766  fe_degree + 1,
767  n_q_points_1d,
768  Number,
769  Number2>;
770  if (!do_integrate)
771  {
772  // Evaluate in 2d
773  if (n_blocks[0] == n_rows_n)
774  {
775  EvalN eval(shape_data[0].shape_values_eo,
776  shape_data[0].shape_gradients_eo,
777  {});
778  eval.template values<0, true, false>(values_dofs, values);
779  if (evaluation_flag & EvaluationFlags::gradients)
780  {
781  eval.template gradients<0, true, false, dim>(
782  values_dofs, gradients);
783  eval.template values<0, true, false, dim>(
784  values_dofs + n_rows_n, gradients + 1);
785  }
786  }
787  else
788  {
789  Eval eval(shape_data[1].shape_values_eo,
790  shape_data[1].shape_gradients_eo,
791  {});
792  eval.template values<0, true, false>(values_dofs, values);
793  if (evaluation_flag & EvaluationFlags::gradients)
794  {
795  eval.template gradients<0, true, false, dim>(
796  values_dofs, gradients);
797  eval.template values<0, true, false, dim>(
798  values_dofs + n_rows_t, gradients + 1);
799  }
800  }
801  }
802  else
803  {
804  // Integrate in 2d
805  if (n_blocks[0] == n_rows_n)
806  {
807  EvalN eval(shape_data[0].shape_values_eo,
808  shape_data[0].shape_gradients_eo,
809  {});
810  if (evaluation_flag & EvaluationFlags::values)
811  eval.template values<0, false, false>(values,
812  values_dofs);
813  if (evaluation_flag & EvaluationFlags::gradients)
814  {
815  if (evaluation_flag & EvaluationFlags::values)
816  eval.template gradients<0, false, true, dim>(
817  gradients, values_dofs);
818  else
819  eval.template gradients<0, false, false, dim>(
820  gradients, values_dofs);
821  eval.template values<0, false, false, dim>(
822  gradients + 1, values_dofs + n_rows_n);
823  }
824  }
825  else
826  {
827  Eval eval(shape_data[1].shape_values_eo,
828  shape_data[1].shape_gradients_eo,
829  {});
830  if (evaluation_flag & EvaluationFlags::values)
831  eval.template values<0, false, false>(values,
832  values_dofs);
833  if (evaluation_flag & EvaluationFlags::gradients)
834  {
835  if (evaluation_flag & EvaluationFlags::values)
836  eval.template gradients<0, false, true, dim>(
837  gradients, values_dofs);
838  else
839  eval.template gradients<0, false, false, dim>(
840  gradients, values_dofs);
841  eval.template values<0, false, false, dim>(
842  gradients + 1, values_dofs + n_rows_t);
843  }
844  }
845  }
846  }
847  values += Utilities::pow(n_q_points_1d, dim - 1);
848  gradients += dim * Utilities::pow(n_q_points_1d, dim - 1);
849  }
850  }
851  };
852 
853 
854 
855  template <int dim, int fe_degree, typename Number>
857  {
858  using Number2 =
860 
861  template <bool do_evaluate, bool add_into_output>
862  static void
863  interpolate(const unsigned int n_components,
865  const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info,
866  const Number *input,
867  Number *output,
868  const unsigned int face_no)
869  {
870  Assert(static_cast<unsigned int>(fe_degree) ==
871  shape_info.data.front().fe_degree ||
872  fe_degree == -1,
873  ExcInternalError());
875  interpolate_raviart_thomas<do_evaluate, add_into_output>(
876  n_components, input, output, flags, face_no, shape_info);
877  else
878  interpolate_generic<do_evaluate, add_into_output>(
879  n_components,
880  input,
881  output,
882  flags,
883  face_no,
884  shape_info.data.front().fe_degree + 1,
885  shape_info.data.front().shape_data_on_face,
886  shape_info.dofs_per_component_on_cell,
887  3 * shape_info.dofs_per_component_on_face);
888  }
889 
893  template <bool do_evaluate, bool add_into_output>
894  static void
896  const unsigned int n_components,
898  const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info,
899  const Number *input,
900  Number *output,
901  const unsigned int face_no)
902  {
903  Assert(static_cast<unsigned int>(fe_degree + 1) ==
904  shape_info.data.front().n_q_points_1d ||
905  fe_degree == -1,
906  ExcInternalError());
907 
908  interpolate_generic<do_evaluate, add_into_output>(
909  n_components,
910  input,
911  output,
912  flags,
913  face_no,
914  shape_info.data.front().quadrature.size(),
915  shape_info.data.front().quadrature_data_on_face,
916  shape_info.n_q_points,
917  shape_info.n_q_points_face);
918  }
919 
920  private:
921  template <bool do_evaluate, bool add_into_output, int face_direction = 0>
922  static void
923  interpolate_generic(const unsigned int n_components,
924  const Number *input,
925  Number *output,
927  const unsigned int face_no,
928  const unsigned int n_points_1d,
929  const std::array<AlignedVector<Number2>, 2> &shape_data,
930  const unsigned int dofs_per_component_on_cell,
931  const unsigned int dofs_per_component_on_face)
932  {
933  if (face_direction == face_no / 2)
934  {
935  constexpr int stride_ = Utilities::pow(fe_degree + 1, face_direction);
936 
937  const int n_rows = fe_degree != -1 ? fe_degree + 1 : n_points_1d;
938  const int stride = Utilities::pow(n_rows, face_direction);
939  const std::array<int, 2> n_blocks{
940  {(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}};
941  std::array<int, 2> steps;
942  if constexpr (face_direction == 0)
943  steps = {{n_rows, 0}};
944  else if constexpr (face_direction == 1 && dim == 2)
945  steps = {{1, 0}};
946  else if constexpr (face_direction == 1)
947  // in 3d, the coordinate system is zx, not xz -> switch indices
948  steps = {{n_rows * n_rows, -n_rows * n_rows * n_rows + 1}};
949  else if constexpr (face_direction == 2)
950  steps = {{1, 0}};
951 
952  for (unsigned int c = 0; c < n_components; ++c)
953  {
954  if (flag & EvaluationFlags::hessians)
955  interpolate_to_face<fe_degree + 1,
956  stride_,
957  do_evaluate,
958  add_into_output,
959  2>(shape_data[face_no % 2].begin(),
960  n_blocks,
961  steps,
962  input,
963  output,
964  n_rows,
965  stride);
966  else if (flag & EvaluationFlags::gradients)
967  interpolate_to_face<fe_degree + 1,
968  stride_,
969  do_evaluate,
970  add_into_output,
971  1>(shape_data[face_no % 2].begin(),
972  n_blocks,
973  steps,
974  input,
975  output,
976  n_rows,
977  stride);
978  else
979  interpolate_to_face<fe_degree + 1,
980  stride_,
981  do_evaluate,
982  add_into_output,
983  0>(shape_data[face_no % 2].begin(),
984  n_blocks,
985  steps,
986  input,
987  output,
988  n_rows,
989  stride);
990  if (do_evaluate)
991  {
992  input += dofs_per_component_on_cell;
993  output += dofs_per_component_on_face;
994  }
995  else
996  {
997  output += dofs_per_component_on_cell;
998  input += dofs_per_component_on_face;
999  }
1000  }
1001  }
1002  else if (face_direction < dim)
1003  {
1004  interpolate_generic<do_evaluate,
1005  add_into_output,
1006  std::min(face_direction + 1, dim - 1)>(
1007  n_components,
1008  input,
1009  output,
1010  flag,
1011  face_no,
1012  n_points_1d,
1013  shape_data,
1014  dofs_per_component_on_cell,
1015  dofs_per_component_on_face);
1016  }
1017  }
1018 
1019  template <bool do_evaluate,
1020  bool add_into_output,
1021  int face_direction = 0,
1022  int max_derivative = 0>
1023  static void
1025  const unsigned int n_components,
1026  const Number *input,
1027  Number *output,
1029  const unsigned int face_no,
1030  const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info)
1031  {
1032  if (dim == 1)
1033  {
1034  // This should never happen since the FE_RaviartThomasNodal is not
1035  // defined for dim = 1. It prevents compiler warnings of infinite
1036  // recursion.
1037  Assert(false, ExcInternalError());
1038  return;
1039  }
1040 
1041  bool increase_max_der = false;
1042  if ((flag & EvaluationFlags::hessians && max_derivative < 2) ||
1043  (flag & EvaluationFlags::gradients && max_derivative < 1))
1044  increase_max_der = true;
1045 
1046  if (face_direction == face_no / 2 && !increase_max_der)
1047  {
1048  constexpr int stride1 = Utilities::pow(fe_degree + 1, face_direction);
1049  constexpr int stride0 = Utilities::pow(fe_degree, face_direction);
1050  constexpr int stride2 = fe_degree * (fe_degree + 1);
1051 
1052  const int degree =
1053  fe_degree != -1 ? fe_degree : shape_info.data[0].fe_degree;
1054  const int n_rows_n = degree + 1;
1055  const int n_rows_t = degree;
1056 
1057  std::array<int, 3> strides{{1, 1, 1}};
1058  if (face_direction > 0)
1059  {
1060  strides[0] =
1061  n_rows_n * Utilities::pow(n_rows_t, face_direction - 1);
1062  strides[1] = n_rows_t * (face_direction == 3 ? n_rows_n : 1);
1063  strides[2] = Utilities::pow(n_rows_t, face_direction);
1064  }
1065  const ::ndarray<int, 3, 3> dofs_per_direction{
1066  {{{n_rows_n, n_rows_t, n_rows_t}},
1067  {{n_rows_t, n_rows_n, n_rows_t}},
1068  {{n_rows_t, n_rows_t, n_rows_n}}}};
1069 
1070  std::array<int, 2> steps, n_blocks;
1071 
1072  if constexpr (face_direction == 0)
1073  steps = {{degree + (face_direction == 0), 0}};
1074  else if constexpr (face_direction == 1 && dim == 2)
1075  steps = {{1, 0}};
1076  else if constexpr (face_direction == 1)
1077  // in 3d, the coordinate system is zx, not xz -> switch indices
1078  steps = {
1079  {n_rows_n * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
1080  else if constexpr (face_direction == 2)
1081  steps = {{1, 0}};
1082 
1083  n_blocks[0] = dofs_per_direction[0][(face_direction + 1) % dim];
1084  n_blocks[1] =
1085  dim > 2 ? dofs_per_direction[0][(face_direction + 2) % dim] : 1;
1086 
1088  (fe_degree != -1 ? (fe_degree + (face_direction == 0)) : 0),
1089  ((face_direction < 2) ? stride1 : stride2),
1090  do_evaluate,
1091  add_into_output,
1092  max_derivative>(shape_info.data[face_direction != 0]
1093  .shape_data_on_face[face_no % 2]
1094  .begin(),
1095  n_blocks,
1096  steps,
1097  input,
1098  output,
1099  degree + (face_direction == 0),
1100  strides[0]);
1101 
1102  if (do_evaluate)
1103  {
1104  input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1105  output += 3 * n_blocks[0] * n_blocks[1];
1106  }
1107  else
1108  {
1109  output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1110  input += 3 * n_blocks[0] * n_blocks[1];
1111  }
1112 
1113  // must only change steps only for face direction 0
1114  if constexpr (face_direction == 0)
1115  steps = {{degree, 0}};
1116 
1117  n_blocks[0] = dofs_per_direction[1][(face_direction + 1) % dim];
1118  n_blocks[1] =
1119  dim > 2 ? dofs_per_direction[1][(face_direction + 2) % dim] : 1;
1120 
1122  (fe_degree != -1 ? (fe_degree + (face_direction == 1)) : 0),
1123  ((face_direction < 2) ? stride0 : stride2),
1124  do_evaluate,
1125  add_into_output,
1126  max_derivative>(shape_info.data[face_direction != 1]
1127  .shape_data_on_face[face_no % 2]
1128  .begin(),
1129  n_blocks,
1130  steps,
1131  input,
1132  output,
1133  degree + (face_direction == 1),
1134  strides[1]);
1135 
1136  if constexpr (dim > 2)
1137  {
1138  if (do_evaluate)
1139  {
1140  input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1141  output += 3 * n_blocks[0] * n_blocks[1];
1142  }
1143  else
1144  {
1145  output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1146  input += 3 * n_blocks[0] * n_blocks[1];
1147  }
1148 
1149  if constexpr (face_direction == 0)
1150  steps = {{degree, 0}};
1151  else if constexpr (face_direction == 1)
1152  // in 3d, the coordinate system is zx, not xz -> switch indices
1153  steps = {
1154  {n_rows_t * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
1155  else if constexpr (face_direction == 2)
1156  steps = {{1, 0}};
1157 
1158  n_blocks[0] = dofs_per_direction[2][(face_direction + 1) % dim];
1159  n_blocks[1] = dofs_per_direction[2][(face_direction + 2) % dim];
1160 
1162  (fe_degree != -1 ? (fe_degree + (face_direction == 2)) : 0),
1163  stride0,
1164  do_evaluate,
1165  add_into_output,
1166  max_derivative>(shape_info.data[face_direction != 2]
1167  .shape_data_on_face[face_no % 2]
1168  .begin(),
1169  n_blocks,
1170  steps,
1171  input,
1172  output,
1173  degree + (face_direction == 2),
1174  strides[2]);
1175  }
1176  }
1177  else if (face_direction == face_no / 2)
1178  {
1179  // Only increase max_derivative
1180  interpolate_raviart_thomas<do_evaluate,
1181  add_into_output,
1182  face_direction,
1183  std::min(max_derivative + 1, 2)>(
1184  n_components, input, output, flag, face_no, shape_info);
1185  }
1186  else if (face_direction < dim)
1187  {
1188  if (increase_max_der)
1189  {
1190  interpolate_raviart_thomas<do_evaluate,
1191  add_into_output,
1192  std::min(face_direction + 1, dim - 1),
1193  std::min(max_derivative + 1, 2)>(
1194  n_components, input, output, flag, face_no, shape_info);
1195  }
1196  else
1197  {
1198  interpolate_raviart_thomas<do_evaluate,
1199  add_into_output,
1200  std::min(face_direction + 1, dim - 1),
1201  max_derivative>(
1202  n_components, input, output, flag, face_no, shape_info);
1203  }
1204  }
1205  }
1206  };
1207 
1208 
1209 
1210  // internal helper function for reading data; base version of different types
1211  template <typename VectorizedArrayType, typename Number2>
1212  void
1213  do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
1214  {
1215  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1216  dst[v] = src_ptr[v];
1217  }
1218 
1219 
1220 
1221  // internal helper function for reading data; specialized version where we
1222  // can use a dedicated load function
1223  template <typename Number, std::size_t width>
1224  void
1226  {
1227  dst.load(src_ptr);
1228  }
1229 
1230 
1231 
1232  // internal helper function for reading data; base version of different types
1233  template <typename VectorizedArrayType, typename Number2>
1234  void
1235  do_vectorized_gather(const Number2 *src_ptr,
1236  const unsigned int *indices,
1237  VectorizedArrayType &dst)
1238  {
1239  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1240  dst[v] = src_ptr[indices[v]];
1241  }
1242 
1243 
1244 
1245  // internal helper function for reading data; specialized version where we
1246  // can use a dedicated gather function
1247  template <typename Number, std::size_t width>
1248  void
1249  do_vectorized_gather(const Number *src_ptr,
1250  const unsigned int *indices,
1252  {
1253  dst.gather(src_ptr, indices);
1254  }
1255 
1256 
1257 
1258  // internal helper function for reading data; base version of different types
1259  template <typename VectorizedArrayType, typename Number2>
1260  void
1261  do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
1262  {
1263  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1264  dst_ptr[v] += src[v];
1265  }
1266 
1267 
1268 
1269  // internal helper function for reading data; specialized version where we
1270  // can use a dedicated load function
1271  template <typename Number, std::size_t width>
1272  void
1274  {
1276  tmp.load(dst_ptr);
1277  (tmp + src).store(dst_ptr);
1278  }
1279 
1280 
1281 
1282  // internal helper function for reading data; base version of different types
1283  template <typename VectorizedArrayType, typename Number2>
1284  void
1285  do_vectorized_scatter_add(const VectorizedArrayType src,
1286  const unsigned int *indices,
1287  Number2 *dst_ptr)
1288  {
1289  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1290  dst_ptr[indices[v]] += src[v];
1291  }
1292 
1293 
1294 
1295  // internal helper function for reading data; specialized version where we
1296  // can use a dedicated gather function
1297  template <typename Number, std::size_t width>
1298  void
1300  const unsigned int *indices,
1301  Number *dst_ptr)
1302  {
1303 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
1304  for (unsigned int v = 0; v < width; ++v)
1305  dst_ptr[indices[v]] += src[v];
1306 #else
1308  tmp.gather(dst_ptr, indices);
1309  (tmp + src).scatter(indices, dst_ptr);
1310 #endif
1311  }
1312 
1313 
1314 
1315  template <typename Number>
1316  void
1317  adjust_for_face_orientation(const unsigned int dim,
1318  const unsigned int n_components,
1320  const unsigned int *orientation,
1321  const bool integrate,
1322  const std::size_t n_q_points,
1323  Number *tmp_values,
1324  Number *values_quad,
1325  Number *gradients_quad,
1326  Number *hessians_quad)
1327  {
1328  for (unsigned int c = 0; c < n_components; ++c)
1329  {
1330  if (flag & EvaluationFlags::values)
1331  {
1332  if (integrate)
1333  for (unsigned int q = 0; q < n_q_points; ++q)
1334  tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
1335  else
1336  for (unsigned int q = 0; q < n_q_points; ++q)
1337  tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
1338  for (unsigned int q = 0; q < n_q_points; ++q)
1339  values_quad[c * n_q_points + q] = tmp_values[q];
1340  }
1341  if (flag & EvaluationFlags::gradients)
1342  for (unsigned int d = 0; d < dim; ++d)
1343  {
1344  if (integrate)
1345  for (unsigned int q = 0; q < n_q_points; ++q)
1346  tmp_values[q] =
1347  gradients_quad[(c * n_q_points + orientation[q]) * dim + d];
1348  else
1349  for (unsigned int q = 0; q < n_q_points; ++q)
1350  tmp_values[orientation[q]] =
1351  gradients_quad[(c * n_q_points + q) * dim + d];
1352  for (unsigned int q = 0; q < n_q_points; ++q)
1353  gradients_quad[(c * n_q_points + q) * dim + d] = tmp_values[q];
1354  }
1355  if (flag & EvaluationFlags::hessians)
1356  {
1357  const unsigned int hdim = (dim * (dim + 1)) / 2;
1358  for (unsigned int d = 0; d < hdim; ++d)
1359  {
1360  if (integrate)
1361  for (unsigned int q = 0; q < n_q_points; ++q)
1362  tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
1363  orientation[q]];
1364  else
1365  for (unsigned int q = 0; q < n_q_points; ++q)
1366  tmp_values[orientation[q]] =
1367  hessians_quad[(c * hdim + d) * n_q_points + q];
1368  for (unsigned int q = 0; q < n_q_points; ++q)
1369  hessians_quad[(c * hdim + d) * n_q_points + q] =
1370  tmp_values[q];
1371  }
1372  }
1373  }
1374  }
1375 
1376 
1377 
1378  template <typename Number, typename VectorizedArrayType>
1379  void
1381  const unsigned int dim,
1382  const unsigned int n_components,
1383  const unsigned int v,
1385  const unsigned int *orientation,
1386  const bool integrate,
1387  const std::size_t n_q_points,
1388  Number *tmp_values,
1389  VectorizedArrayType *values_quad,
1390  VectorizedArrayType *gradients_quad = nullptr,
1391  VectorizedArrayType *hessians_quad = nullptr)
1392  {
1393  for (unsigned int c = 0; c < n_components; ++c)
1394  {
1395  if (flag & EvaluationFlags::values)
1396  {
1397  if (integrate)
1398  for (unsigned int q = 0; q < n_q_points; ++q)
1399  tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v];
1400  else
1401  for (unsigned int q = 0; q < n_q_points; ++q)
1402  tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v];
1403  for (unsigned int q = 0; q < n_q_points; ++q)
1404  values_quad[c * n_q_points + q][v] = tmp_values[q];
1405  }
1406  if (flag & EvaluationFlags::gradients)
1407  for (unsigned int d = 0; d < dim; ++d)
1408  {
1409  Assert(gradients_quad != nullptr, ExcInternalError());
1410  if (integrate)
1411  for (unsigned int q = 0; q < n_q_points; ++q)
1412  tmp_values[q] =
1413  gradients_quad[(c * n_q_points + orientation[q]) * dim + d]
1414  [v];
1415  else
1416  for (unsigned int q = 0; q < n_q_points; ++q)
1417  tmp_values[orientation[q]] =
1418  gradients_quad[(c * n_q_points + q) * dim + d][v];
1419  for (unsigned int q = 0; q < n_q_points; ++q)
1420  gradients_quad[(c * n_q_points + q) * dim + d][v] =
1421  tmp_values[q];
1422  }
1423  if (flag & EvaluationFlags::hessians)
1424  {
1425  Assert(hessians_quad != nullptr, ExcInternalError());
1426  const unsigned int hdim = (dim * (dim + 1)) / 2;
1427  for (unsigned int d = 0; d < hdim; ++d)
1428  {
1429  if (integrate)
1430  for (unsigned int q = 0; q < n_q_points; ++q)
1431  tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
1432  orientation[q]][v];
1433  else
1434  for (unsigned int q = 0; q < n_q_points; ++q)
1435  tmp_values[orientation[q]] =
1436  hessians_quad[(c * hdim + d) * n_q_points + q][v];
1437  for (unsigned int q = 0; q < n_q_points; ++q)
1438  hessians_quad[(c * hdim + d) * n_q_points + q][v] =
1439  tmp_values[q];
1440  }
1441  }
1442  }
1443  }
1444 
1445 
1446 
1447  template <int dim, typename Number>
1449  {
1450  template <int fe_degree, int n_q_points_1d>
1451  static bool
1452  run(const unsigned int n_components,
1453  const EvaluationFlags::EvaluationFlags evaluation_flag,
1454  const Number *values_dofs,
1456  {
1457  const auto &shape_info = fe_eval.get_shape_info();
1458  const auto &shape_data = shape_info.data.front();
1459  using Number2 =
1461 
1462  if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
1463  {
1464  Assert((fe_eval.get_dof_access_index() ==
1466  fe_eval.is_interior_face() == false) == false,
1467  ExcNotImplemented());
1468 
1469  const unsigned int face_no = fe_eval.get_face_no();
1470  const unsigned int face_orientation = fe_eval.get_face_orientation();
1471  const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
1472  const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
1473 
1474  using Eval =
1476 
1477  if (evaluation_flag & EvaluationFlags::values)
1478  {
1479  const auto *const shape_values =
1480  &shape_data.shape_values_face(face_no, face_orientation, 0);
1481 
1482  auto *values_quad_ptr = fe_eval.begin_values();
1483  auto *values_dofs_actual_ptr = values_dofs;
1484 
1485  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
1486  for (unsigned int c = 0; c < n_components; ++c)
1487  {
1488  eval.template values<0, true, false>(values_dofs_actual_ptr,
1489  values_quad_ptr);
1490 
1491  values_quad_ptr += n_q_points;
1492  values_dofs_actual_ptr += n_dofs;
1493  }
1494  }
1495 
1496  if (evaluation_flag & EvaluationFlags::gradients)
1497  {
1498  auto *gradients_quad_ptr = fe_eval.begin_gradients();
1499  const auto *values_dofs_actual_ptr = values_dofs;
1500 
1501  std::array<const Number2 *, dim> shape_gradients;
1502  for (unsigned int d = 0; d < dim; ++d)
1503  shape_gradients[d] = &shape_data.shape_gradients_face(
1504  face_no, face_orientation, d, 0);
1505 
1506  for (unsigned int c = 0; c < n_components; ++c)
1507  {
1508  for (unsigned int d = 0; d < dim; ++d)
1509  {
1510  Eval eval(nullptr,
1511  shape_gradients[d],
1512  nullptr,
1513  n_dofs,
1514  n_q_points);
1515 
1516  eval.template gradients<0, true, false, dim>(
1517  values_dofs_actual_ptr, gradients_quad_ptr + d);
1518  }
1519  gradients_quad_ptr += n_q_points * dim;
1520  values_dofs_actual_ptr += n_dofs;
1521  }
1522  }
1523 
1524  Assert(!(evaluation_flag & EvaluationFlags::hessians),
1525  ExcNotImplemented());
1526 
1527  return true;
1528  }
1529 
1530  const unsigned int dofs_per_face =
1531  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
1532  Utilities::pow(shape_data.fe_degree + 1, dim - 1);
1533 
1534  // Note: we always keep storage of values, 1st and 2nd derivatives in an
1535  // array, so reserve space for all three here
1536  Number *temp = fe_eval.get_scratch_data().begin();
1537  Number *scratch_data = temp + 3 * n_components * dofs_per_face;
1538 
1539  bool use_vectorization = true;
1540 
1541  if (fe_eval.get_dof_access_index() ==
1543  fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
1544  for (unsigned int v = 0; v < Number::size(); ++v)
1545  if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
1546  fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
1547  use_vectorization = false;
1548 
1549  if (use_vectorization == false)
1550  {
1551  for (unsigned int v = 0; v < Number::size(); ++v)
1552  {
1553  // the loop breaks once an invalid_unsigned_int is hit for
1554  // all cases except the exterior faces in the ECL loop (where
1555  // some faces might be at the boundaries but others not)
1556  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
1557  {
1558  for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
1559  ++i)
1560  temp[i][v] = 0;
1561  continue;
1562  }
1563 
1565  template interpolate<true, false>(n_components,
1566  evaluation_flag,
1567  shape_info,
1568  values_dofs,
1569  scratch_data,
1570  fe_eval.get_face_no(v));
1571 
1572  for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
1573  ++i)
1574  temp[i][v] = scratch_data[i][v];
1575  }
1576  }
1577  else
1579  template interpolate<true, false>(n_components,
1580  evaluation_flag,
1581  shape_info,
1582  values_dofs,
1583  temp,
1584  fe_eval.get_face_no());
1585 
1586  const unsigned int subface_index = fe_eval.get_subface_index();
1587  constexpr unsigned int n_q_points_1d_actual =
1588  fe_degree > -1 ? n_q_points_1d : 0;
1589 
1590  if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
1591  {
1593  fe_degree,
1594  n_q_points_1d_actual,
1595  Number>::
1596  template evaluate_or_integrate_in_face<false>(
1597  evaluation_flag,
1598  fe_eval.get_shape_info().data,
1599  temp,
1600  fe_eval.begin_values(),
1601  fe_eval.begin_gradients(),
1602  scratch_data,
1603  subface_index,
1604  fe_eval.get_face_no() / 2);
1605  }
1606  else if (fe_degree > -1 &&
1607  subface_index >= GeometryInfo<dim>::max_children_per_cell &&
1608  shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
1609  FEFaceEvaluationImpl<true,
1610  dim,
1611  fe_degree,
1612  n_q_points_1d_actual,
1613  Number>::evaluate_in_face(n_components,
1614  evaluation_flag,
1615  shape_data,
1616  temp,
1617  fe_eval.begin_values(),
1618  fe_eval
1619  .begin_gradients(),
1620  fe_eval.begin_hessians(),
1621  scratch_data,
1622  subface_index);
1623  else
1624  FEFaceEvaluationImpl<false,
1625  dim,
1626  fe_degree,
1627  n_q_points_1d_actual,
1628  Number>::evaluate_in_face(n_components,
1629  evaluation_flag,
1630  shape_data,
1631  temp,
1632  fe_eval.begin_values(),
1633  fe_eval
1634  .begin_gradients(),
1635  fe_eval.begin_hessians(),
1636  scratch_data,
1637  subface_index);
1638 
1639  if (use_vectorization == false)
1640  {
1641  for (unsigned int v = 0; v < Number::size(); ++v)
1642  {
1643  // the loop breaks once an invalid_unsigned_int is hit for
1644  // all cases except the exterior faces in the ECL loop (where
1645  // some faces might be at the boundaries but others not)
1646  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
1647  continue;
1648 
1649  if (fe_eval.get_face_orientation(v) != 0)
1651  dim,
1652  n_components,
1653  v,
1654  evaluation_flag,
1656  fe_eval.get_face_orientation(v), 0),
1657  false,
1658  shape_info.n_q_points_face,
1659  &temp[0][0],
1660  fe_eval.begin_values(),
1661  fe_eval.begin_gradients(),
1662  fe_eval.begin_hessians());
1663  }
1664  }
1665  else if (fe_eval.get_face_orientation() != 0)
1667  dim,
1668  n_components,
1669  evaluation_flag,
1671  fe_eval.get_face_orientation(), 0),
1672  false,
1673  shape_info.n_q_points_face,
1674  temp,
1675  fe_eval.begin_values(),
1676  fe_eval.begin_gradients(),
1677  fe_eval.begin_hessians());
1678 
1679  return false;
1680  }
1681  };
1682 
1683 
1684 
1685  template <int dim, typename Number>
1687  {
1688  template <int fe_degree, int n_q_points_1d>
1689  static bool
1690  run(const unsigned int n_components,
1691  const EvaluationFlags::EvaluationFlags integration_flag,
1692  Number *values_dofs,
1694  {
1695  const auto &shape_info = fe_eval.get_shape_info();
1696  const auto &shape_data = shape_info.data.front();
1697  using Number2 =
1699 
1700  if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
1701  {
1702  Assert((fe_eval.get_dof_access_index() ==
1704  fe_eval.is_interior_face() == false) == false,
1705  ExcNotImplemented());
1706 
1707  const unsigned int face_no = fe_eval.get_face_no();
1708  const unsigned int face_orientation = fe_eval.get_face_orientation();
1709  const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
1710  const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
1711 
1712  using Eval =
1714 
1715  if (integration_flag & EvaluationFlags::values)
1716  {
1717  const auto *const shape_values =
1718  &shape_data.shape_values_face(face_no, face_orientation, 0);
1719 
1720  auto *values_quad_ptr = fe_eval.begin_values();
1721  auto *values_dofs_actual_ptr = values_dofs;
1722 
1723  Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
1724  for (unsigned int c = 0; c < n_components; ++c)
1725  {
1726  eval.template values<0, false, false>(values_quad_ptr,
1727  values_dofs_actual_ptr);
1728 
1729  values_quad_ptr += n_q_points;
1730  values_dofs_actual_ptr += n_dofs;
1731  }
1732  }
1733 
1734  if (integration_flag & EvaluationFlags::gradients)
1735  {
1736  auto *gradients_quad_ptr = fe_eval.begin_gradients();
1737  auto *values_dofs_actual_ptr = values_dofs;
1738 
1739  std::array<const Number2 *, dim> shape_gradients;
1740  for (unsigned int d = 0; d < dim; ++d)
1741  shape_gradients[d] = &shape_data.shape_gradients_face(
1742  face_no, face_orientation, d, 0);
1743 
1744  for (unsigned int c = 0; c < n_components; ++c)
1745  {
1746  for (unsigned int d = 0; d < dim; ++d)
1747  {
1748  Eval eval(nullptr,
1749  shape_gradients[d],
1750  nullptr,
1751  n_dofs,
1752  n_q_points);
1753 
1754  if (!(integration_flag & EvaluationFlags::values) &&
1755  d == 0)
1756  eval.template gradients<0, false, false, dim>(
1757  gradients_quad_ptr + d, values_dofs_actual_ptr);
1758  else
1759  eval.template gradients<0, false, true, dim>(
1760  gradients_quad_ptr + d, values_dofs_actual_ptr);
1761  }
1762  gradients_quad_ptr += n_q_points * dim;
1763  values_dofs_actual_ptr += n_dofs;
1764  }
1765  }
1766 
1767  Assert(!(integration_flag & EvaluationFlags::hessians),
1768  ExcNotImplemented());
1769 
1770  return true;
1771  }
1772 
1773  const unsigned int dofs_per_face =
1774  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
1775  Utilities::pow(shape_data.fe_degree + 1, dim - 1);
1776 
1777  Number *temp = fe_eval.get_scratch_data().begin();
1778  Number *scratch_data = temp + 3 * n_components * dofs_per_face;
1779 
1780  bool use_vectorization = true;
1781 
1782  if (fe_eval.get_dof_access_index() ==
1784  fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
1785  use_vectorization =
1786  fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
1787  std::all_of(fe_eval.get_cell_ids().begin() + 1,
1788  fe_eval.get_cell_ids().end(),
1789  [&](const auto &v) {
1790  return v == fe_eval.get_cell_ids()[0] ||
1791  v == numbers::invalid_unsigned_int;
1792  });
1793 
1794  if (use_vectorization == false)
1795  {
1796  for (unsigned int v = 0; v < Number::size(); ++v)
1797  {
1798  // the loop breaks once an invalid_unsigned_int is hit for
1799  // all cases except the exterior faces in the ECL loop (where
1800  // some faces might be at the boundaries but others not)
1801  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
1802  continue;
1803 
1804  if (fe_eval.get_face_orientation(v) != 0)
1806  dim,
1807  n_components,
1808  v,
1809  integration_flag,
1811  fe_eval.get_face_orientation(v), 0),
1812  true,
1813  shape_info.n_q_points_face,
1814  &temp[0][0],
1815  fe_eval.begin_values(),
1816  fe_eval.begin_gradients(),
1817  fe_eval.begin_hessians());
1818  }
1819  }
1820  else if (fe_eval.get_face_orientation() != 0)
1822  dim,
1823  n_components,
1824  integration_flag,
1826  fe_eval.get_face_orientation(), 0),
1827  true,
1828  shape_info.n_q_points_face,
1829  temp,
1830  fe_eval.begin_values(),
1831  fe_eval.begin_gradients(),
1832  fe_eval.begin_hessians());
1833 
1834  const unsigned int n_q_points_1d_actual =
1835  fe_degree > -1 ? n_q_points_1d : 0;
1836  const unsigned int subface_index = fe_eval.get_subface_index();
1837 
1838  if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
1839  {
1841  fe_degree,
1842  n_q_points_1d_actual,
1843  Number>::
1844  template evaluate_or_integrate_in_face<true>(
1845  integration_flag,
1846  fe_eval.get_shape_info().data,
1847  temp,
1848  fe_eval.begin_values(),
1849  fe_eval.begin_gradients(),
1850  scratch_data,
1851  subface_index,
1852  fe_eval.get_face_no() / 2);
1853  }
1854  else if (fe_degree > -1 &&
1855  fe_eval.get_subface_index() >=
1857  shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
1859  true,
1860  dim,
1861  fe_degree,
1862  n_q_points_1d_actual,
1863  Number>::integrate_in_face(n_components,
1864  integration_flag,
1865  shape_data,
1866  temp,
1867  fe_eval.begin_values(),
1868  fe_eval.begin_gradients(),
1869  fe_eval.begin_hessians(),
1870  scratch_data,
1871  subface_index);
1872  else
1874  false,
1875  dim,
1876  fe_degree,
1877  n_q_points_1d_actual,
1878  Number>::integrate_in_face(n_components,
1879  integration_flag,
1880  shape_data,
1881  temp,
1882  fe_eval.begin_values(),
1883  fe_eval.begin_gradients(),
1884  fe_eval.begin_hessians(),
1885  scratch_data,
1886  subface_index);
1887 
1888  if (use_vectorization == false)
1889  {
1890  for (unsigned int v = 0; v < Number::size(); ++v)
1891  {
1892  // the loop breaks once an invalid_unsigned_int is hit for
1893  // all cases except the exterior faces in the ECL loop (where
1894  // some faces might be at the boundaries but others not)
1895  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
1896  continue;
1897 
1899  template interpolate<false, false>(n_components,
1900  integration_flag,
1901  shape_info,
1902  values_dofs,
1903  scratch_data,
1904  fe_eval.get_face_no(v));
1905 
1906  for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
1907  ++i)
1908  temp[i][v] = scratch_data[i][v];
1909  }
1910  }
1911  else
1913  template interpolate<false, false>(n_components,
1914  integration_flag,
1915  shape_info,
1916  temp,
1917  values_dofs,
1918  fe_eval.get_face_no());
1919  return false;
1920  }
1921  };
1922 
1923 
1924 
1925  template <int n_face_orientations,
1926  typename Processor,
1927  typename EvaluationData,
1928  const bool check_face_orientations = false>
1929  void
1931  Processor &proc,
1932  const unsigned int n_components,
1933  const EvaluationFlags::EvaluationFlags evaluation_flag,
1934  typename Processor::Number2_ *global_vector_ptr,
1935  const std::vector<ArrayView<const typename Processor::Number2_>> *sm_ptr,
1936  const EvaluationData &fe_eval,
1937  typename Processor::VectorizedArrayType_ *temp1)
1938  {
1939  constexpr int dim = Processor::dim_;
1940  constexpr int fe_degree = Processor::fe_degree_;
1941  using VectorizedArrayType = typename Processor::VectorizedArrayType_;
1942  constexpr int n_lanes = VectorizedArrayType::size();
1943 
1944  using Number = typename Processor::Number_;
1945  using Number2_ = typename Processor::Number2_;
1946 
1947  const auto &shape_data = fe_eval.get_shape_info().data.front();
1948  constexpr bool integrate = Processor::do_integrate;
1949  const unsigned int face_no = fe_eval.get_face_no();
1950  const auto &dof_info = fe_eval.get_dof_info();
1951  const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
1952  const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index =
1953  fe_eval.get_dof_access_index();
1954  AssertIndexRange(cell,
1955  dof_info.index_storage_variants[dof_access_index].size());
1956  constexpr unsigned int dofs_per_face =
1957  Utilities::pow(fe_degree + 1, dim - 1);
1958  const unsigned int subface_index = fe_eval.get_subface_index();
1959 
1960  const unsigned int n_filled_lanes =
1961  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
1962 
1963  bool all_faces_are_same = n_filled_lanes == n_lanes;
1964  if (n_face_orientations == n_lanes)
1965  for (unsigned int v = 1; v < n_lanes; ++v)
1966  if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) ||
1967  fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0))
1968  {
1969  all_faces_are_same = false;
1970  break;
1971  }
1972 
1973  // check for re-orientation ...
1974  std::array<const unsigned int *, n_face_orientations> orientation = {};
1975 
1976  if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same &&
1977  fe_eval.is_interior_face() == 0)
1978  for (unsigned int v = 0; v < n_lanes; ++v)
1979  {
1980  // the loop breaks once an invalid_unsigned_int is hit for
1981  // all cases except the exterior faces in the ECL loop (where
1982  // some faces might be at the boundaries but others not)
1983  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
1984  continue;
1985 
1986  if (shape_data.nodal_at_cell_boundaries &&
1987  fe_eval.get_face_orientation(v) != 0)
1988  {
1989  // ... and in case we detect a re-orientation, go to the other
1990  // version of this function that actually allows for this
1991  if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
1992  check_face_orientations == false)
1993  {
1994  fe_face_evaluation_process_and_io<n_face_orientations,
1995  Processor,
1996  EvaluationData,
1997  true>(proc,
1998  n_components,
1999  evaluation_flag,
2000  global_vector_ptr,
2001  sm_ptr,
2002  fe_eval,
2003  temp1);
2004  return;
2005  }
2006  orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
2007  fe_eval.get_face_orientation(v), 0);
2008  }
2009  }
2010  else if (dim == 3 && fe_eval.get_face_orientation() != 0)
2011  {
2012  // go to the other version of this function
2013  if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
2014  check_face_orientations == false)
2015  {
2016  fe_face_evaluation_process_and_io<n_face_orientations,
2017  Processor,
2018  EvaluationData,
2019  true>(proc,
2020  n_components,
2021  evaluation_flag,
2022  global_vector_ptr,
2023  sm_ptr,
2024  fe_eval,
2025  temp1);
2026  return;
2027  }
2028  for (unsigned int v = 0; v < n_face_orientations; ++v)
2029  orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
2030  fe_eval.get_face_orientation(), 0);
2031  }
2032 
2033  // we know that the gradient weights for the Hermite case on the
2034  // right (side==1) are the negative from the value at the left
2035  // (side==0), so we only read out one of them.
2036  VectorizedArrayType grad_weight =
2037  shape_data
2038  .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) :
2039  (1 + face_no % 2))];
2040 
2041  // face_to_cell_index_hermite
2042  std::array<const unsigned int *, n_face_orientations> index_array_hermite =
2043  {};
2044  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2045  {
2046  if (n_face_orientations == 1)
2047  index_array_hermite[0] =
2048  &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0);
2049  else
2050  {
2051  for (unsigned int v = 0; v < n_lanes; ++v)
2052  {
2053  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2054  continue;
2055 
2056  const auto face_no = fe_eval.get_face_no(v);
2057 
2058  grad_weight[v] =
2059  shape_data.shape_data_on_face[0][fe_degree +
2060  (integrate ?
2061  (2 - (face_no % 2)) :
2062  (1 + (face_no % 2)))][0];
2063 
2064  index_array_hermite[v] =
2065  &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no,
2066  0);
2067  }
2068  }
2069  }
2070 
2071  // face_to_cell_index_nodal
2072  std::array<const unsigned int *, n_face_orientations> index_array_nodal =
2073  {};
2074  if (shape_data.nodal_at_cell_boundaries == true)
2075  {
2076  if (n_face_orientations == 1)
2077  index_array_nodal[0] =
2078  &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0);
2079  else
2080  {
2081  for (unsigned int v = 0; v < n_lanes; ++v)
2082  {
2083  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2084  continue;
2085 
2086  const auto face_no = fe_eval.get_face_no(v);
2087 
2088  index_array_nodal[v] =
2089  &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no,
2090  0);
2091  }
2092  }
2093  }
2094 
2095 
2096  const auto reorientate = [&](const unsigned int v, const unsigned int i) {
2097  return (!check_face_orientations || orientation[v] == nullptr) ?
2098  i :
2099  orientation[v][i];
2100  };
2101 
2102  const unsigned int cell_index =
2103  dof_access_index == MatrixFreeFunctions::DoFInfo::dof_access_cell ?
2104  fe_eval.get_cell_ids()[0] :
2105  cell * n_lanes;
2106  const unsigned int *dof_indices =
2107  &dof_info.dof_indices_contiguous[dof_access_index][cell_index];
2108 
2109  for (unsigned int comp = 0; comp < n_components; ++comp)
2110  {
2111  const std::size_t index_offset =
2112  dof_info.component_dof_indices_offset
2113  [fe_eval.get_active_fe_index()]
2114  [fe_eval.get_first_selected_component()] +
2115  comp * Utilities::pow(fe_degree + 1, dim);
2116 
2117  // case 1: contiguous and interleaved indices
2118  if (n_face_orientations == 1 &&
2119  dof_info.index_storage_variants[dof_access_index][cell] ==
2121  interleaved_contiguous)
2122  {
2124  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2125  n_lanes);
2126  Number2_ *vector_ptr =
2127  global_vector_ptr + dof_indices[0] + index_offset * n_lanes;
2128 
2129  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2130  {
2131  for (unsigned int i = 0; i < dofs_per_face; ++i)
2132  {
2133  Assert(n_face_orientations == 1, ExcNotImplemented());
2134 
2135  const unsigned int ind1 = index_array_hermite[0][2 * i];
2136  const unsigned int ind2 = index_array_hermite[0][2 * i + 1];
2137  const unsigned int i_ = reorientate(0, i);
2138  proc.hermite_grad_vectorized(temp1[i_],
2139  temp1[i_ + dofs_per_face],
2140  vector_ptr + ind1 * n_lanes,
2141  vector_ptr + ind2 * n_lanes,
2142  grad_weight);
2143  }
2144  }
2145  else
2146  {
2147  for (unsigned int i = 0; i < dofs_per_face; ++i)
2148  {
2149  Assert(n_face_orientations == 1, ExcNotImplemented());
2150 
2151  const unsigned int i_ = reorientate(0, i);
2152  const unsigned int ind = index_array_nodal[0][i];
2153  proc.value_vectorized(temp1[i_],
2154  vector_ptr + ind * n_lanes);
2155  }
2156  }
2157  }
2158 
2159  // case 2: contiguous and interleaved indices with fixed stride
2160  else if (n_face_orientations == 1 &&
2161  dof_info.index_storage_variants[dof_access_index][cell] ==
2163  interleaved_contiguous_strided)
2164  {
2166  dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2167  n_lanes);
2168  Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes;
2169  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2170  {
2171  for (unsigned int i = 0; i < dofs_per_face; ++i)
2172  {
2173  Assert(n_face_orientations == 1, ExcNotImplemented());
2174 
2175  const unsigned int i_ = reorientate(0, i);
2176  const unsigned int ind1 =
2177  index_array_hermite[0][2 * i] * n_lanes;
2178  const unsigned int ind2 =
2179  index_array_hermite[0][2 * i + 1] * n_lanes;
2180  proc.hermite_grad_vectorized_indexed(
2181  temp1[i_],
2182  temp1[i_ + dofs_per_face],
2183  vector_ptr + ind1,
2184  vector_ptr + ind2,
2185  grad_weight,
2186  dof_indices,
2187  dof_indices);
2188  }
2189  }
2190  else
2191  {
2192  for (unsigned int i = 0; i < dofs_per_face; ++i)
2193  {
2194  Assert(n_face_orientations == 1, ExcNotImplemented());
2195 
2196  const unsigned int i_ = reorientate(0, i);
2197  const unsigned int ind = index_array_nodal[0][i] * n_lanes;
2198  proc.value_vectorized_indexed(temp1[i_],
2199  vector_ptr + ind,
2200  dof_indices);
2201  }
2202  }
2203  }
2204 
2205  // case 3: contiguous and interleaved indices with mixed stride
2206  else if (n_face_orientations == 1 &&
2207  dof_info.index_storage_variants[dof_access_index][cell] ==
2209  interleaved_contiguous_mixed_strides)
2210  {
2211  const unsigned int *strides =
2212  &dof_info.dof_indices_interleave_strides[dof_access_index]
2213  [cell * n_lanes];
2214  unsigned int indices[n_lanes];
2215  for (unsigned int v = 0; v < n_lanes; ++v)
2216  indices[v] = dof_indices[v] + index_offset * strides[v];
2217  const unsigned int n_filled_lanes =
2218  dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
2219 
2220  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2221  {
2222  if (n_filled_lanes == n_lanes)
2223  for (unsigned int i = 0; i < dofs_per_face; ++i)
2224  {
2225  Assert(n_face_orientations == 1, ExcNotImplemented());
2226 
2227  const unsigned int i_ = reorientate(0, i);
2228  unsigned int ind1[n_lanes];
2230  for (unsigned int v = 0; v < n_lanes; ++v)
2231  ind1[v] = indices[v] +
2232  index_array_hermite[0][2 * i] * strides[v];
2233  unsigned int ind2[n_lanes];
2235  for (unsigned int v = 0; v < n_lanes; ++v)
2236  ind2[v] =
2237  indices[v] +
2238  // TODO
2239  index_array_hermite[0][2 * i + 1] * strides[v];
2240  proc.hermite_grad_vectorized_indexed(
2241  temp1[i_],
2242  temp1[i_ + dofs_per_face],
2243  global_vector_ptr,
2244  global_vector_ptr,
2245  grad_weight,
2246  ind1,
2247  ind2);
2248  }
2249  else
2250  {
2251  if (integrate == false)
2252  for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
2253  temp1[i] = VectorizedArrayType();
2254 
2255  for (unsigned int v = 0; v < n_filled_lanes; ++v)
2256  for (unsigned int i = 0; i < dofs_per_face; ++i)
2257  {
2258  const unsigned int i_ =
2259  reorientate(n_face_orientations == 1 ? 0 : v, i);
2260  proc.hermite_grad(
2261  temp1[i_][v],
2262  temp1[i_ + dofs_per_face][v],
2263  global_vector_ptr
2264  [indices[v] +
2265  index_array_hermite
2266  [n_face_orientations == 1 ? 0 : v][2 * i] *
2267  strides[v]],
2268  global_vector_ptr
2269  [indices[v] +
2270  index_array_hermite[n_face_orientations == 1 ?
2271  0 :
2272  v][2 * i + 1] *
2273  strides[v]],
2274  grad_weight[n_face_orientations == 1 ? 0 : v]);
2275  }
2276  }
2277  }
2278  else
2279  {
2280  if (n_filled_lanes == n_lanes)
2281  for (unsigned int i = 0; i < dofs_per_face; ++i)
2282  {
2283  Assert(n_face_orientations == 1, ExcInternalError());
2284  unsigned int ind[n_lanes];
2286  for (unsigned int v = 0; v < n_lanes; ++v)
2287  ind[v] =
2288  indices[v] + index_array_nodal[0][i] * strides[v];
2289  const unsigned int i_ = reorientate(0, i);
2290  proc.value_vectorized_indexed(temp1[i_],
2291  global_vector_ptr,
2292  ind);
2293  }
2294  else
2295  {
2296  if (integrate == false)
2297  for (unsigned int i = 0; i < dofs_per_face; ++i)
2298  temp1[i] = VectorizedArrayType();
2299 
2300  for (unsigned int v = 0; v < n_filled_lanes; ++v)
2301  for (unsigned int i = 0; i < dofs_per_face; ++i)
2302  proc.value(
2303  temp1[reorientate(n_face_orientations == 1 ? 0 : v,
2304  i)][v],
2305  global_vector_ptr
2306  [indices[v] +
2307  index_array_nodal[n_face_orientations == 1 ? 0 : v]
2308  [i] *
2309  strides[v]]);
2310  }
2311  }
2312  }
2313 
2314  // case 4: contiguous indices without interleaving
2315  else if (n_face_orientations > 1 ||
2316  dof_info.index_storage_variants[dof_access_index][cell] ==
2318  contiguous)
2319  {
2320  Number2_ *vector_ptr = global_vector_ptr + index_offset;
2321 
2322  const bool vectorization_possible =
2323  all_faces_are_same && (sm_ptr == nullptr);
2324 
2325  std::array<Number2_ *, n_lanes> vector_ptrs;
2326  std::array<unsigned int, n_lanes> reordered_indices;
2327 
2328  if (vectorization_possible == false)
2329  {
2330  vector_ptrs = {};
2331  if (n_face_orientations == 1)
2332  {
2333  for (unsigned int v = 0; v < n_filled_lanes; ++v)
2334  if (sm_ptr == nullptr)
2335  {
2336  vector_ptrs[v] = vector_ptr + dof_indices[v];
2337  }
2338  else
2339  {
2340  const auto &temp =
2341  dof_info
2342  .dof_indices_contiguous_sm[dof_access_index]
2343  [cell * n_lanes + v];
2344  vector_ptrs[v] = const_cast<Number2_ *>(
2345  sm_ptr->operator[](temp.first).data() +
2346  temp.second + index_offset);
2347  }
2348  }
2349  else if (n_face_orientations == n_lanes)
2350  {
2351  const auto &cells = fe_eval.get_cell_ids();
2352  for (unsigned int v = 0; v < n_lanes; ++v)
2353  if (cells[v] != numbers::invalid_unsigned_int)
2354  {
2355  if (sm_ptr == nullptr)
2356  {
2357  vector_ptrs[v] =
2358  vector_ptr +
2359  dof_info
2360  .dof_indices_contiguous[dof_access_index]
2361  [cells[v]];
2362  }
2363  else
2364  {
2365  const auto &temp =
2366  dof_info
2367  .dof_indices_contiguous_sm[dof_access_index]
2368  [cells[v]];
2369  vector_ptrs[v] = const_cast<Number2_ *>(
2370  sm_ptr->operator[](temp.first).data() +
2371  temp.second + index_offset);
2372  }
2373  }
2374  }
2375  else
2376  {
2377  Assert(false, ExcNotImplemented());
2378  }
2379  }
2380  else if (n_face_orientations == n_lanes)
2381  {
2382  for (unsigned int v = 0; v < n_lanes; ++v)
2383  reordered_indices[v] =
2384  dof_info.dof_indices_contiguous[dof_access_index]
2385  [fe_eval.get_cell_ids()[v]];
2386  dof_indices = reordered_indices.data();
2387  }
2388 
2389  if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2390  {
2391  if (vectorization_possible)
2392  for (unsigned int i = 0; i < dofs_per_face; ++i)
2393  {
2394  const unsigned int ind1 = index_array_hermite[0][2 * i];
2395  const unsigned int ind2 =
2396  index_array_hermite[0][2 * i + 1];
2397  const unsigned int i_ = reorientate(0, i);
2398 
2399  proc.hermite_grad_vectorized_indexed(
2400  temp1[i_],
2401  temp1[i_ + dofs_per_face],
2402  vector_ptr + ind1,
2403  vector_ptr + ind2,
2404  grad_weight,
2405  dof_indices,
2406  dof_indices);
2407  }
2408  else if (n_face_orientations == 1)
2409  for (unsigned int i = 0; i < dofs_per_face; ++i)
2410  {
2411  const unsigned int ind1 = index_array_hermite[0][2 * i];
2412  const unsigned int ind2 =
2413  index_array_hermite[0][2 * i + 1];
2414  const unsigned int i_ = reorientate(0, i);
2415 
2416  for (unsigned int v = 0; v < n_filled_lanes; ++v)
2417  proc.hermite_grad(temp1[i_][v],
2418  temp1[i_ + dofs_per_face][v],
2419  vector_ptrs[v][ind1],
2420  vector_ptrs[v][ind2],
2421  grad_weight[v]);
2422 
2423  if (integrate == false)
2424  for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
2425  {
2426  temp1[i][v] = 0.0;
2427  temp1[i + dofs_per_face][v] = 0.0;
2428  }
2429  }
2430  else
2431  {
2432  if (integrate == false && n_filled_lanes < n_lanes)
2433  for (unsigned int i = 0; i < dofs_per_face; ++i)
2434  temp1[i] = temp1[i + dofs_per_face] = Number();
2435 
2436  for (unsigned int v = 0; v < n_filled_lanes; ++v)
2437  for (unsigned int i = 0; i < dofs_per_face; ++i)
2438  proc.hermite_grad(
2439  temp1[reorientate(v, i)][v],
2440  temp1[reorientate(v, i) + dofs_per_face][v],
2441  vector_ptrs[v][index_array_hermite[v][2 * i]],
2442  vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
2443  grad_weight[v]);
2444  }
2445  }
2446  else
2447  {
2448  if (vectorization_possible)
2449  for (unsigned int i = 0; i < dofs_per_face; ++i)
2450  {
2451  const unsigned int ind = index_array_nodal[0][i];
2452  const unsigned int i_ = reorientate(0, i);
2453 
2454  proc.value_vectorized_indexed(temp1[i_],
2455  vector_ptr + ind,
2456  dof_indices);
2457  }
2458  else if (n_face_orientations == 1)
2459  for (unsigned int i = 0; i < dofs_per_face; ++i)
2460  {
2461  const unsigned int ind = index_array_nodal[0][i];
2462  const unsigned int i_ = reorientate(0, i);
2463 
2464  for (unsigned int v = 0; v < n_filled_lanes; ++v)
2465  proc.value(temp1[i_][v], vector_ptrs[v][ind]);
2466 
2467  if (integrate == false)
2468  for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
2469  temp1[i_][v] = 0.0;
2470  }
2471  else
2472  {
2473  if (integrate == false && n_filled_lanes < n_lanes)
2474  for (unsigned int i = 0; i < dofs_per_face; ++i)
2475  temp1[i] = Number();
2476 
2477  for (unsigned int v = 0; v < n_filled_lanes; ++v)
2478  for (unsigned int i = 0; i < dofs_per_face; ++i)
2479  proc.value(temp1[reorientate(v, i)][v],
2480  vector_ptrs[v][index_array_nodal[v][i]]);
2481  }
2482  }
2483  }
2484  else
2485  {
2486  // We should not end up here, this should be caught by
2487  // FEFaceEvaluationImplGatherEvaluateSelector::supports()
2488  Assert(false, ExcInternalError());
2489  }
2490  temp1 += 3 * dofs_per_face;
2491  }
2492  }
2493 
2494 
2495 
2496  template <int dim, typename Number2, typename VectorizedArrayType>
2498  {
2499  using Number = typename VectorizedArrayType::value_type;
2500 
2501  template <int fe_degree, int n_q_points_1d>
2502  static bool
2503  run(const unsigned int n_components,
2504  const EvaluationFlags::EvaluationFlags evaluation_flag,
2505  const Number2 *src_ptr,
2506  const std::vector<ArrayView<const Number2>> *sm_ptr,
2508  {
2509  Assert(fe_degree > -1, ExcInternalError());
2510  Assert(fe_eval.get_shape_info().element_type <=
2512  ExcInternalError());
2513 
2514  const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
2515 
2516  VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
2517  VectorizedArrayType *scratch_data =
2518  temp + 3 * n_components * dofs_per_face;
2519 
2521 
2522  if (fe_eval.get_dof_access_index() ==
2524  fe_eval.is_interior_face() == false)
2525  fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
2526  p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
2527  else
2528  fe_face_evaluation_process_and_io<1>(
2529  p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
2530 
2531  const unsigned int subface_index = fe_eval.get_subface_index();
2532 
2533  if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
2534  FEFaceEvaluationImpl<true,
2535  dim,
2536  fe_degree,
2537  n_q_points_1d,
2538  VectorizedArrayType>::
2539  evaluate_in_face(n_components,
2540  evaluation_flag,
2541  fe_eval.get_shape_info().data.front(),
2542  temp,
2543  fe_eval.begin_values(),
2544  fe_eval.begin_gradients(),
2545  fe_eval.begin_hessians(),
2546  scratch_data,
2547  subface_index);
2548  else
2549  FEFaceEvaluationImpl<false,
2550  dim,
2551  fe_degree,
2552  n_q_points_1d,
2553  VectorizedArrayType>::
2554  evaluate_in_face(n_components,
2555  evaluation_flag,
2556  fe_eval.get_shape_info().data.front(),
2557  temp,
2558  fe_eval.begin_values(),
2559  fe_eval.begin_gradients(),
2560  fe_eval.begin_hessians(),
2561  scratch_data,
2562  subface_index);
2563 
2564  // re-orientation for cases not possible with above algorithm
2565  if (subface_index < GeometryInfo<dim>::max_children_per_cell)
2566  {
2567  if (fe_eval.get_dof_access_index() ==
2569  fe_eval.is_interior_face() == false)
2570  {
2571  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2572  {
2573  // the loop breaks once an invalid_unsigned_int is hit for
2574  // all cases except the exterior faces in the ECL loop (where
2575  // some faces might be at the boundaries but others not)
2576  if (fe_eval.get_cell_ids()[v] ==
2578  continue;
2579 
2580  if (fe_eval.get_face_orientation(v) != 0)
2582  dim,
2583  n_components,
2584  v,
2585  evaluation_flag,
2587  fe_eval.get_face_orientation(v), 0),
2588  false,
2589  Utilities::pow(n_q_points_1d, dim - 1),
2590  &temp[0][0],
2591  fe_eval.begin_values(),
2592  fe_eval.begin_gradients(),
2593  fe_eval.begin_hessians());
2594  }
2595  }
2596  else if (fe_eval.get_face_orientation() != 0)
2598  dim,
2599  n_components,
2600  evaluation_flag,
2602  fe_eval.get_face_orientation(), 0),
2603  false,
2604  Utilities::pow(n_q_points_1d, dim - 1),
2605  temp,
2606  fe_eval.begin_values(),
2607  fe_eval.begin_gradients(),
2608  fe_eval.begin_hessians());
2609  }
2610 
2611  return false;
2612  }
2613 
2614  template <typename Number3>
2615  static bool
2617  const MatrixFreeFunctions::ShapeInfo<Number3> &shape_info,
2618  const Number2 *vector_ptr,
2620  {
2621  const unsigned int fe_degree = shape_info.data.front().fe_degree;
2622  if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries ||
2623  (evaluation_flag & EvaluationFlags::gradients &&
2624  (fe_degree < 2 ||
2625  shape_info.data.front().element_type !=
2627  (evaluation_flag & EvaluationFlags::hessians) ||
2628  vector_ptr == nullptr ||
2629  shape_info.data.front().element_type >
2631  storage <
2633  return false;
2634  else
2635  return true;
2636  }
2637 
2638  private:
2639  template <int fe_degree>
2640  struct Processor
2641  {
2642  static const bool do_integrate = false;
2643  static const int dim_ = dim;
2644  static const int fe_degree_ = fe_degree;
2645  using VectorizedArrayType_ = VectorizedArrayType;
2646  using Number_ = Number;
2647  using Number2_ = const Number2;
2648 
2649  template <typename T0, typename T1, typename T2>
2650  void
2652  T0 &temp_2,
2653  const T1 src_ptr_1,
2654  const T1 src_ptr_2,
2655  const T2 &grad_weight)
2656  {
2657  do_vectorized_read(src_ptr_1, temp_1);
2658  do_vectorized_read(src_ptr_2, temp_2);
2659  temp_2 = grad_weight * (temp_1 - temp_2);
2660  }
2661 
2662  template <typename T1, typename T2>
2663  void
2664  value_vectorized(T1 &temp, const T2 src_ptr)
2665  {
2666  do_vectorized_read(src_ptr, temp);
2667  }
2668 
2669  template <typename T0, typename T1, typename T2, typename T3>
2670  void
2672  T0 &temp_2,
2673  const T1 src_ptr_1,
2674  const T1 src_ptr_2,
2675  const T2 &grad_weight,
2676  const T3 &indices_1,
2677  const T3 &indices_2)
2678  {
2679  do_vectorized_gather(src_ptr_1, indices_1, temp_1);
2680  do_vectorized_gather(src_ptr_2, indices_2, temp_2);
2681  temp_2 = grad_weight * (temp_1 - temp_2);
2682  }
2683 
2684  template <typename T0, typename T1, typename T2>
2685  void
2686  value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
2687  {
2688  do_vectorized_gather(src_ptr, indices, temp);
2689  }
2690 
2691  template <typename T0, typename T1, typename T2>
2692  void
2693  hermite_grad(T0 &temp_1,
2694  T0 &temp_2,
2695  const T1 &src_ptr_1,
2696  const T1 &src_ptr_2,
2697  const T2 &grad_weight)
2698  {
2699  // case 3a)
2700  temp_1 = src_ptr_1;
2701  temp_2 = grad_weight * (temp_1 - src_ptr_2);
2702  }
2703 
2704  template <typename T1, typename T2>
2705  void
2706  value(T1 &temp, const T2 &src_ptr)
2707  {
2708  // case 3b)
2709  temp = src_ptr;
2710  }
2711  };
2712  };
2713 
2714 
2715 
2716  template <int dim, typename Number2, typename VectorizedArrayType>
2718  {
2719  using Number = typename VectorizedArrayType::value_type;
2720 
2721  template <int fe_degree, int n_q_points_1d>
2722  static bool
2723  run(const unsigned int n_components,
2724  const EvaluationFlags::EvaluationFlags integration_flag,
2725  Number2 *dst_ptr,
2726  const std::vector<ArrayView<const Number2>> *sm_ptr,
2728  {
2729  Assert(fe_degree > -1, ExcInternalError());
2730  Assert(fe_eval.get_shape_info().element_type <=
2732  ExcInternalError());
2733 
2734  const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
2735 
2736  VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
2737  VectorizedArrayType *scratch_data =
2738  temp + 3 * n_components * dofs_per_face;
2739 
2740  const unsigned int subface_index = fe_eval.get_subface_index();
2741 
2742  // re-orientation for cases not possible with the io function below
2743  if (subface_index < GeometryInfo<dim>::max_children_per_cell)
2744  {
2745  if (fe_eval.get_dof_access_index() ==
2747  fe_eval.is_interior_face() == false)
2748  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2749  {
2750  // the loop breaks once an invalid_unsigned_int is hit for
2751  // all cases except the exterior faces in the ECL loop (where
2752  // some faces might be at the boundaries but others not)
2753  if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2754  continue;
2755 
2756  if (fe_eval.get_face_orientation(v) != 0)
2758  dim,
2759  n_components,
2760  v,
2761  integration_flag,
2763  fe_eval.get_face_orientation(v), 0),
2764  true,
2765  Utilities::pow(n_q_points_1d, dim - 1),
2766  &temp[0][0],
2767  fe_eval.begin_values(),
2768  fe_eval.begin_gradients(),
2769  fe_eval.begin_hessians());
2770  }
2771  else if (fe_eval.get_face_orientation() != 0)
2773  dim,
2774  n_components,
2775  integration_flag,
2777  fe_eval.get_face_orientation(), 0),
2778  true,
2779  Utilities::pow(n_q_points_1d, dim - 1),
2780  temp,
2781  fe_eval.begin_values(),
2782  fe_eval.begin_gradients(),
2783  fe_eval.begin_hessians());
2784  }
2785 
2786  if (fe_degree > -1 && fe_eval.get_subface_index() >=
2787  GeometryInfo<dim - 1>::max_children_per_cell)
2788  FEFaceEvaluationImpl<true,
2789  dim,
2790  fe_degree,
2791  n_q_points_1d,
2792  VectorizedArrayType>::
2793  integrate_in_face(n_components,
2794  integration_flag,
2795  fe_eval.get_shape_info().data.front(),
2796  temp,
2797  fe_eval.begin_values(),
2798  fe_eval.begin_gradients(),
2799  fe_eval.begin_hessians(),
2800  scratch_data,
2801  subface_index);
2802  else
2803  FEFaceEvaluationImpl<false,
2804  dim,
2805  fe_degree,
2806  n_q_points_1d,
2807  VectorizedArrayType>::
2808  integrate_in_face(n_components,
2809  integration_flag,
2810  fe_eval.get_shape_info().data.front(),
2811  temp,
2812  fe_eval.begin_values(),
2813  fe_eval.begin_gradients(),
2814  fe_eval.begin_hessians(),
2815  scratch_data,
2816  subface_index);
2817 
2819 
2820  if (fe_eval.get_dof_access_index() ==
2822  fe_eval.is_interior_face() == false)
2823  fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
2824  p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
2825  else
2826  fe_face_evaluation_process_and_io<1>(
2827  p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
2828 
2829  return false;
2830  }
2831 
2832  private:
2833  template <int fe_degree>
2834  struct Processor
2835  {
2836  static const bool do_integrate = true;
2837  static const int dim_ = dim;
2838  static const int fe_degree_ = fe_degree;
2839  using VectorizedArrayType_ = VectorizedArrayType;
2840  using Number_ = Number;
2841  using Number2_ = Number2;
2842 
2843  template <typename T0, typename T1, typename T2, typename T3, typename T4>
2844  void
2845  hermite_grad_vectorized(const T0 &temp_1,
2846  const T1 &temp_2,
2847  T2 dst_ptr_1,
2848  T3 dst_ptr_2,
2849  const T4 &grad_weight)
2850  {
2851  // case 1a)
2852  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
2853  const VectorizedArrayType grad = grad_weight * temp_2;
2854  do_vectorized_add(val, dst_ptr_1);
2855  do_vectorized_add(grad, dst_ptr_2);
2856  }
2857 
2858  template <typename T0, typename T1>
2859  void
2860  value_vectorized(const T0 &temp, T1 dst_ptr)
2861  {
2862  // case 1b)
2863  do_vectorized_add(temp, dst_ptr);
2864  }
2865 
2866  template <typename T0, typename T1, typename T2, typename T3>
2867  void
2869  const T0 &temp_2,
2870  T1 dst_ptr_1,
2871  T1 dst_ptr_2,
2872  const T2 &grad_weight,
2873  const T3 &indices_1,
2874  const T3 &indices_2)
2875  {
2876  // case 2a)
2877  const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
2878  const VectorizedArrayType grad = grad_weight * temp_2;
2879  do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
2880  do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
2881  }
2882 
2883  template <typename T0, typename T1, typename T2>
2884  void
2885  value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
2886  {
2887  // case 2b)
2888  do_vectorized_scatter_add(temp, indices, dst_ptr);
2889  }
2890 
2891  template <typename T0, typename T1, typename T2>
2892  void
2893  hermite_grad(const T0 &temp_1,
2894  const T0 &temp_2,
2895  T1 &dst_ptr_1,
2896  T1 &dst_ptr_2,
2897  const T2 &grad_weight)
2898  {
2899  // case 3a)
2900  const Number val = temp_1 - grad_weight * temp_2;
2901  const Number grad = grad_weight * temp_2;
2902  dst_ptr_1 += val;
2903  dst_ptr_2 += grad;
2904  }
2905 
2906  template <typename T0, typename T1>
2907  void
2908  value(const T0 &temp, T1 &dst_ptr)
2909  {
2910  // case 3b)
2911  dst_ptr += temp;
2912  }
2913  };
2914  };
2915 } // end of namespace internal
2916 
2917 
2919 
2920 #endif
iterator begin() const
Definition: array_view.h:703
std::uint8_t get_face_no(const unsigned int v=0) const
internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex get_dof_access_index() const
const Number * begin_values() const
const Number * begin_hessians() const
const std::array< unsigned int, n_lanes > & get_cell_ids() const
unsigned int get_subface_index() const
const Number * begin_gradients() const
bool is_interior_face() const
ArrayView< Number > get_scratch_data() const
const ShapeInfoType & get_shape_info() const
std::uint8_t get_face_orientation(const unsigned int v=0) const
void gather(const Number *base_ptr, const unsigned int *offsets)
void load(const OtherNumber *ptr)
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:144
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
unsigned int cell_index
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
EvaluationFlags
The EvaluationFlags enum.
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T scatter(const MPI_Comm comm, const std::vector< T > &objects_to_send, const unsigned int root_process=0)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:447
T fixed_power(const T t)
Definition: utilities.h:975
void do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
void adjust_for_face_orientation_per_lane(const unsigned int dim, const unsigned int n_components, const unsigned int v, const EvaluationFlags::EvaluationFlags flag, const unsigned int *orientation, const bool integrate, const std::size_t n_q_points, Number *tmp_values, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad=nullptr, VectorizedArrayType *hessians_quad=nullptr)
void do_vectorized_scatter_add(const VectorizedArrayType src, const unsigned int *indices, Number2 *dst_ptr)
std::enable_if_t< contract_onto_face, void > interpolate_to_face(const Number2 *shape_values, const std::array< int, 2 > &n_blocks, const std::array< int, 2 > &steps, const Number *input, Number *DEAL_II_RESTRICT output, const int n_rows_runtime=0, const int stride_runtime=1)
void do_vectorized_gather(const Number2 *src_ptr, const unsigned int *indices, VectorizedArrayType &dst)
void do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
void fe_face_evaluation_process_and_io(Processor &proc, const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, typename Processor::Number2_ *global_vector_ptr, const std::vector< ArrayView< const typename Processor::Number2_ >> *sm_ptr, const EvaluationData &fe_eval, typename Processor::VectorizedArrayType_ *temp1)
void adjust_for_face_orientation(const unsigned int dim, const unsigned int n_components, const EvaluationFlags::EvaluationFlags flag, const unsigned int *orientation, const bool integrate, const std::size_t n_q_points, Number *tmp_values, Number *values_quad, Number *gradients_quad, Number *hessians_quad)
static const unsigned int invalid_unsigned_int
Definition: types.h:221
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
void hermite_grad(T0 &temp_1, T0 &temp_2, const T1 &src_ptr_1, const T1 &src_ptr_2, const T2 &grad_weight)
void value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
void hermite_grad_vectorized(T0 &temp_1, T0 &temp_2, const T1 src_ptr_1, const T1 src_ptr_2, const T2 &grad_weight)
void hermite_grad_vectorized_indexed(T0 &temp_1, T0 &temp_2, const T1 src_ptr_1, const T1 src_ptr_2, const T2 &grad_weight, const T3 &indices_1, const T3 &indices_2)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number2 *src_ptr, const std::vector< ArrayView< const Number2 >> *sm_ptr, FEEvaluationData< dim, VectorizedArrayType, true > &fe_eval)
static bool supports(const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::ShapeInfo< Number3 > &shape_info, const Number2 *vector_ptr, MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage)
void hermite_grad_vectorized(const T0 &temp_1, const T1 &temp_2, T2 dst_ptr_1, T3 dst_ptr_2, const T4 &grad_weight)
void hermite_grad_vectorized_indexed(const T0 &temp_1, const T0 &temp_2, T1 dst_ptr_1, T1 dst_ptr_2, const T2 &grad_weight, const T3 &indices_1, const T3 &indices_2)
void hermite_grad(const T0 &temp_1, const T0 &temp_2, T1 &dst_ptr_1, T1 &dst_ptr_2, const T2 &grad_weight)
void value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number2 *dst_ptr, const std::vector< ArrayView< const Number2 >> *sm_ptr, FEEvaluationData< dim, VectorizedArrayType, true > &fe_eval)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
typename FEEvaluationData< dim, Number, true >::shape_info_number_type Number2
static void evaluate_or_integrate_in_face(const EvaluationFlags::EvaluationFlags evaluation_flag, const std::vector< MatrixFreeFunctions::UnivariateShapeData< Number2 >> &shape_data, Number *values_dofs_in, Number *values, Number *gradients, Number *scratch_data, const unsigned int subface_index, const unsigned int face_direction)
EvaluatorTensorProduct< symmetric_evaluate ? evaluate_evenodd :evaluate_general, dim - 1, fe_degree+1, n_q_points_1d, Number, Number2 > Eval
static void evaluate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data, const unsigned int subface_index)
static Eval create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const unsigned int subface_index, const unsigned int direction)
static void integrate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data, const unsigned int subface_index)
typename FEEvaluationData< dim, Number, true >::shape_info_number_type Number2
typename FEEvaluationData< dim, Number, true >::shape_info_number_type Number2
static void interpolate_quadrature(const unsigned int n_components, const EvaluationFlags::EvaluationFlags flags, const MatrixFreeFunctions::ShapeInfo< Number2 > &shape_info, const Number *input, Number *output, const unsigned int face_no)
static void interpolate_generic(const unsigned int n_components, const Number *input, Number *output, const EvaluationFlags::EvaluationFlags flag, const unsigned int face_no, const unsigned int n_points_1d, const std::array< AlignedVector< Number2 >, 2 > &shape_data, const unsigned int dofs_per_component_on_cell, const unsigned int dofs_per_component_on_face)
static void interpolate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags flags, const MatrixFreeFunctions::ShapeInfo< Number2 > &shape_info, const Number *input, Number *output, const unsigned int face_no)
static void interpolate_raviart_thomas(const unsigned int n_components, const Number *input, Number *output, const EvaluationFlags::EvaluationFlags flag, const unsigned int face_no, const MatrixFreeFunctions::ShapeInfo< Number2 > &shape_info)
std::vector< UnivariateShapeData< Number > > data
Definition: shape_info.h:486
::Table< 2, unsigned int > face_orientations_quad
Definition: shape_info.h:625
AlignedVector< Number > shape_gradients_collocation_eo
Definition: shape_info.h:262
std::array< AlignedVector< Number >, 2 > hessians_within_subface
Definition: shape_info.h:327
std::array< AlignedVector< Number >, 2 > values_within_subface
Definition: shape_info.h:315
std::array< AlignedVector< Number >, 2 > gradients_within_subface
Definition: shape_info.h:321