Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3075-gc235bd4825 2025-04-15 08:10:00+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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
evaluation_kernels_face.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_matrix_free_evaluation_kernels_face_h
17#define dealii_matrix_free_evaluation_kernels_face_h
18
19#include <deal.II/base/config.h>
20
25
32
33
35
36
37namespace 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],
85 data.gradients_within_subface[index],
86 data.hessians_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:
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:
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:
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:
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:
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:
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 // DEAL_II_NOT_IMPLEMENTED();
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,
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,
875 interpolate_raviart_thomas<do_evaluate, add_into_output>(
876 n_components, input, output, flags, face_no, shape_info);
877 else
878 {
879 const unsigned int fe_degree_ = shape_info.data.front().fe_degree;
880
881 interpolate_generic<do_evaluate, add_into_output>(
882 n_components,
883 input,
884 output,
885 flags,
886 face_no,
887 fe_degree_ + 1,
888 shape_info.data.front().shape_data_on_face,
889 Utilities::pow(fe_degree_ + 1, dim),
890 3 * Utilities::pow(fe_degree_ + 1, dim - 1));
891 }
892 }
893
897 template <bool do_evaluate, bool add_into_output>
898 static void
900 const unsigned int n_components,
903 const Number *input,
904 Number *output,
905 const unsigned int face_no)
906 {
907 Assert(static_cast<unsigned int>(fe_degree + 1) ==
908 shape_info.data.front().n_q_points_1d ||
909 fe_degree == -1,
911
912 interpolate_generic<do_evaluate, add_into_output>(
913 n_components,
914 input,
915 output,
916 flags,
917 face_no,
918 shape_info.data.front().quadrature.size(),
919 shape_info.data.front().quadrature_data_on_face,
920 shape_info.n_q_points,
921 shape_info.n_q_points_face);
922 }
923
924 private:
925 template <bool do_evaluate, bool add_into_output, int face_direction = 0>
926 static void
927 interpolate_generic(const unsigned int n_components,
928 const Number *input,
929 Number *output,
931 const unsigned int face_no,
932 const unsigned int n_points_1d,
933 const std::array<AlignedVector<Number2>, 2> &shape_data,
934 const unsigned int dofs_per_component_on_cell,
935 const unsigned int dofs_per_component_on_face)
936 {
937 if (face_direction == face_no / 2)
938 {
939 constexpr int stride_ = Utilities::pow(fe_degree + 1, face_direction);
940
941 const int n_rows = fe_degree != -1 ? fe_degree + 1 : n_points_1d;
942 const int stride = Utilities::pow(n_rows, face_direction);
943 const std::array<int, 2> n_blocks{
944 {(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}};
945 std::array<int, 2> steps;
946 if constexpr (face_direction == 0)
947 steps = {{n_rows, 0}};
948 else if constexpr (face_direction == 1 && dim == 2)
949 steps = {{1, 0}};
950 else if constexpr (face_direction == 1)
951 // in 3d, the coordinate system is zx, not xz -> switch indices
952 steps = {{n_rows * n_rows, -n_rows * n_rows * n_rows + 1}};
953 else if constexpr (face_direction == 2)
954 steps = {{1, 0}};
955
956 for (unsigned int c = 0; c < n_components; ++c)
957 {
958 if (flag & EvaluationFlags::hessians)
959 interpolate_to_face<fe_degree + 1,
960 stride_,
961 do_evaluate,
962 add_into_output,
963 2>(shape_data[face_no % 2].begin(),
964 n_blocks,
965 steps,
966 input,
967 output,
968 n_rows,
969 stride);
970 else if (flag & EvaluationFlags::gradients)
971 interpolate_to_face<fe_degree + 1,
972 stride_,
973 do_evaluate,
974 add_into_output,
975 1>(shape_data[face_no % 2].begin(),
976 n_blocks,
977 steps,
978 input,
979 output,
980 n_rows,
981 stride);
982 else
983 interpolate_to_face<fe_degree + 1,
984 stride_,
985 do_evaluate,
986 add_into_output,
987 0>(shape_data[face_no % 2].begin(),
988 n_blocks,
989 steps,
990 input,
991 output,
992 n_rows,
993 stride);
994 if (do_evaluate)
995 {
996 input += dofs_per_component_on_cell;
997 output += dofs_per_component_on_face;
998 }
999 else
1000 {
1001 output += dofs_per_component_on_cell;
1002 input += dofs_per_component_on_face;
1003 }
1004 }
1005 }
1006 else if (face_direction < dim)
1007 {
1008 interpolate_generic<do_evaluate,
1009 add_into_output,
1010 std::min(face_direction + 1, dim - 1)>(
1011 n_components,
1012 input,
1013 output,
1014 flag,
1015 face_no,
1016 n_points_1d,
1017 shape_data,
1018 dofs_per_component_on_cell,
1019 dofs_per_component_on_face);
1020 }
1021 }
1022
1023 template <bool do_evaluate,
1024 bool add_into_output,
1025 int face_direction = 0,
1026 int max_derivative = 0>
1027 static void
1029 const unsigned int n_components,
1030 const Number *input,
1031 Number *output,
1033 const unsigned int face_no,
1035 {
1036 if (dim == 1)
1037 {
1038 // This should never happen since the FE_RaviartThomasNodal is not
1039 // defined for dim = 1. It prevents compiler warnings of infinite
1040 // recursion.
1042 return;
1043 }
1044
1045 bool increase_max_der = false;
1046 if ((flag & EvaluationFlags::hessians && max_derivative < 2) ||
1047 (flag & EvaluationFlags::gradients && max_derivative < 1))
1048 increase_max_der = true;
1049
1050 if (face_direction == face_no / 2 && !increase_max_der)
1051 {
1052 constexpr int stride1 = Utilities::pow(fe_degree + 1, face_direction);
1053 constexpr int stride0 = Utilities::pow(fe_degree, face_direction);
1054 constexpr int stride2 = fe_degree * (fe_degree + 1);
1055
1056 const int degree =
1057 fe_degree != -1 ? fe_degree : shape_info.data[0].fe_degree;
1058 const int n_rows_n = degree + 1;
1059 const int n_rows_t = degree;
1060
1061 std::array<int, 3> strides{{1, 1, 1}};
1062 if (face_direction > 0)
1063 {
1064 strides[0] =
1065 n_rows_n * Utilities::pow(n_rows_t, face_direction - 1);
1066 strides[1] = n_rows_t * (face_direction == 3 ? n_rows_n : 1);
1067 strides[2] = Utilities::pow(n_rows_t, face_direction);
1068 }
1069 const ::ndarray<int, 3, 3> dofs_per_direction{
1070 {{{n_rows_n, n_rows_t, n_rows_t}},
1071 {{n_rows_t, n_rows_n, n_rows_t}},
1072 {{n_rows_t, n_rows_t, n_rows_n}}}};
1073
1074 std::array<int, 2> steps, n_blocks;
1075
1076 if constexpr (face_direction == 0)
1077 steps = {{degree + (face_direction == 0), 0}};
1078 else if constexpr (face_direction == 1 && dim == 2)
1079 steps = {{1, 0}};
1080 else if constexpr (face_direction == 1)
1081 // in 3d, the coordinate system is zx, not xz -> switch indices
1082 steps = {
1083 {n_rows_n * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
1084 else if constexpr (face_direction == 2)
1085 steps = {{1, 0}};
1086
1087 n_blocks[0] = dofs_per_direction[0][(face_direction + 1) % dim];
1088 n_blocks[1] =
1089 dim > 2 ? dofs_per_direction[0][(face_direction + 2) % dim] : 1;
1090
1092 (fe_degree != -1 ? (fe_degree + (face_direction == 0)) : 0),
1093 ((face_direction < 2) ? stride1 : stride2),
1094 do_evaluate,
1095 add_into_output,
1096 max_derivative>(shape_info.data[face_direction != 0]
1097 .shape_data_on_face[face_no % 2]
1098 .begin(),
1099 n_blocks,
1100 steps,
1101 input,
1102 output,
1103 degree + (face_direction == 0),
1104 strides[0]);
1105
1106 if (do_evaluate)
1107 {
1108 input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1109 output += 3 * n_blocks[0] * n_blocks[1];
1110 }
1111 else
1112 {
1113 output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1114 input += 3 * n_blocks[0] * n_blocks[1];
1115 }
1116
1117 // must only change steps only for face direction 0
1118 if constexpr (face_direction == 0)
1119 steps = {{degree, 0}};
1120
1121 n_blocks[0] = dofs_per_direction[1][(face_direction + 1) % dim];
1122 n_blocks[1] =
1123 dim > 2 ? dofs_per_direction[1][(face_direction + 2) % dim] : 1;
1124
1126 (fe_degree != -1 ? (fe_degree + (face_direction == 1)) : 0),
1127 ((face_direction < 2) ? stride0 : stride2),
1128 do_evaluate,
1129 add_into_output,
1130 max_derivative>(shape_info.data[face_direction != 1]
1131 .shape_data_on_face[face_no % 2]
1132 .begin(),
1133 n_blocks,
1134 steps,
1135 input,
1136 output,
1137 degree + (face_direction == 1),
1138 strides[1]);
1139
1140 if constexpr (dim > 2)
1141 {
1142 if (do_evaluate)
1143 {
1144 input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1145 output += 3 * n_blocks[0] * n_blocks[1];
1146 }
1147 else
1148 {
1149 output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1150 input += 3 * n_blocks[0] * n_blocks[1];
1151 }
1152
1153 if constexpr (face_direction == 0)
1154 steps = {{degree, 0}};
1155 else if constexpr (face_direction == 1)
1156 // in 3d, the coordinate system is zx, not xz -> switch indices
1157 steps = {
1158 {n_rows_t * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
1159 else if constexpr (face_direction == 2)
1160 steps = {{1, 0}};
1161
1162 n_blocks[0] = dofs_per_direction[2][(face_direction + 1) % dim];
1163 n_blocks[1] = dofs_per_direction[2][(face_direction + 2) % dim];
1164
1166 (fe_degree != -1 ? (fe_degree + (face_direction == 2)) : 0),
1167 stride0,
1168 do_evaluate,
1169 add_into_output,
1170 max_derivative>(shape_info.data[face_direction != 2]
1171 .shape_data_on_face[face_no % 2]
1172 .begin(),
1173 n_blocks,
1174 steps,
1175 input,
1176 output,
1177 degree + (face_direction == 2),
1178 strides[2]);
1179 }
1180 }
1181 else if (face_direction == face_no / 2)
1182 {
1183 // Only increase max_derivative
1184 interpolate_raviart_thomas<do_evaluate,
1185 add_into_output,
1186 face_direction,
1187 std::min(max_derivative + 1, 2)>(
1188 n_components, input, output, flag, face_no, shape_info);
1189 }
1190 else if (face_direction < dim)
1191 {
1192 if (increase_max_der)
1193 {
1194 interpolate_raviart_thomas<do_evaluate,
1195 add_into_output,
1196 std::min(face_direction + 1, dim - 1),
1197 std::min(max_derivative + 1, 2)>(
1198 n_components, input, output, flag, face_no, shape_info);
1199 }
1200 else
1201 {
1202 interpolate_raviart_thomas<do_evaluate,
1203 add_into_output,
1204 std::min(face_direction + 1, dim - 1),
1205 max_derivative>(
1206 n_components, input, output, flag, face_no, shape_info);
1207 }
1208 }
1209 }
1210 };
1211
1212
1213
1214 // internal helper function for reading data; base version of different types
1215 template <typename VectorizedArrayType, typename Number2>
1216 void
1217 do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
1218 {
1219 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1220 dst[v] = src_ptr[v];
1221 }
1222
1223
1224
1225 // internal helper function for reading data; specialized version where we
1226 // can use a dedicated load function
1227 template <typename Number, std::size_t width>
1228 void
1230 {
1231 dst.load(src_ptr);
1232 }
1233
1234
1235
1236 // internal helper function for reading data; base version of different types
1237 template <typename VectorizedArrayType, typename Number2>
1238 void
1239 do_vectorized_gather(const Number2 *src_ptr,
1240 const unsigned int *indices,
1241 VectorizedArrayType &dst)
1242 {
1243 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1244 dst[v] = src_ptr[indices[v]];
1245 }
1246
1247
1248
1249 // internal helper function for reading data; specialized version where we
1250 // can use a dedicated gather function
1251 template <typename Number, std::size_t width>
1252 void
1253 do_vectorized_gather(const Number *src_ptr,
1254 const unsigned int *indices,
1256 {
1257 dst.gather(src_ptr, indices);
1258 }
1259
1260
1261
1262 // internal helper function for reading data; base version of different types
1263 template <typename VectorizedArrayType, typename Number2>
1264 void
1265 do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
1266 {
1267 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1268 dst_ptr[v] += src[v];
1269 }
1270
1271
1272
1273 // internal helper function for reading data; specialized version where we
1274 // can use a dedicated load function
1275 template <typename Number, std::size_t width>
1276 void
1278 {
1280 tmp.load(dst_ptr);
1281 (tmp + src).store(dst_ptr);
1282 }
1283
1284
1285
1286 // internal helper function for reading data; base version of different types
1287 template <typename VectorizedArrayType, typename Number2>
1288 void
1289 do_vectorized_scatter_add(const VectorizedArrayType src,
1290 const unsigned int *indices,
1291 Number2 *dst_ptr)
1292 {
1293 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1294 dst_ptr[indices[v]] += src[v];
1295 }
1296
1297
1298
1299 // internal helper function for reading data; specialized version where we
1300 // can use a dedicated gather function
1301 template <typename Number, std::size_t width>
1302 void
1304 const unsigned int *indices,
1305 Number *dst_ptr)
1306 {
1307#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
1308 for (unsigned int v = 0; v < width; ++v)
1309 dst_ptr[indices[v]] += src[v];
1310#else
1312 tmp.gather(dst_ptr, indices);
1313 (tmp + src).scatter(indices, dst_ptr);
1314#endif
1315 }
1316
1317
1318
1319 template <typename Number>
1320 void
1321 adjust_for_face_orientation(const unsigned int dim,
1322 const unsigned int n_components,
1324 const unsigned int *orientation,
1325 const bool integrate,
1326 const std::size_t n_q_points,
1327 Number *tmp_values,
1328 Number *values_quad,
1329 Number *gradients_quad,
1330 Number *hessians_quad)
1331 {
1332 for (unsigned int c = 0; c < n_components; ++c)
1333 {
1334 if (flag & EvaluationFlags::values)
1335 {
1336 if (integrate)
1337 for (unsigned int q = 0; q < n_q_points; ++q)
1338 tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
1339 else
1340 for (unsigned int q = 0; q < n_q_points; ++q)
1341 tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
1342 for (unsigned int q = 0; q < n_q_points; ++q)
1343 values_quad[c * n_q_points + q] = tmp_values[q];
1344 }
1345 if (flag & EvaluationFlags::gradients)
1346 for (unsigned int d = 0; d < dim; ++d)
1347 {
1348 if (integrate)
1349 for (unsigned int q = 0; q < n_q_points; ++q)
1350 tmp_values[q] =
1351 gradients_quad[(c * n_q_points + orientation[q]) * dim + d];
1352 else
1353 for (unsigned int q = 0; q < n_q_points; ++q)
1354 tmp_values[orientation[q]] =
1355 gradients_quad[(c * n_q_points + q) * dim + d];
1356 for (unsigned int q = 0; q < n_q_points; ++q)
1357 gradients_quad[(c * n_q_points + q) * dim + d] = tmp_values[q];
1358 }
1359 if (flag & EvaluationFlags::hessians)
1360 {
1361 const unsigned int hdim = (dim * (dim + 1)) / 2;
1362 for (unsigned int d = 0; d < hdim; ++d)
1363 {
1364 if (integrate)
1365 for (unsigned int q = 0; q < n_q_points; ++q)
1366 tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
1367 orientation[q]];
1368 else
1369 for (unsigned int q = 0; q < n_q_points; ++q)
1370 tmp_values[orientation[q]] =
1371 hessians_quad[(c * hdim + d) * n_q_points + q];
1372 for (unsigned int q = 0; q < n_q_points; ++q)
1373 hessians_quad[(c * hdim + d) * n_q_points + q] =
1374 tmp_values[q];
1375 }
1376 }
1377 }
1378 }
1379
1380
1381
1382 template <typename Number, typename VectorizedArrayType>
1383 void
1385 const unsigned int dim,
1386 const unsigned int n_components,
1387 const unsigned int v,
1389 const unsigned int *orientation,
1390 const bool integrate,
1391 const std::size_t n_q_points,
1392 Number *tmp_values,
1393 VectorizedArrayType *values_quad,
1394 VectorizedArrayType *gradients_quad = nullptr,
1395 VectorizedArrayType *hessians_quad = nullptr)
1396 {
1397 for (unsigned int c = 0; c < n_components; ++c)
1398 {
1399 if (flag & EvaluationFlags::values)
1400 {
1401 if (integrate)
1402 for (unsigned int q = 0; q < n_q_points; ++q)
1403 tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v];
1404 else
1405 for (unsigned int q = 0; q < n_q_points; ++q)
1406 tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v];
1407 for (unsigned int q = 0; q < n_q_points; ++q)
1408 values_quad[c * n_q_points + q][v] = tmp_values[q];
1409 }
1410 if (flag & EvaluationFlags::gradients)
1411 for (unsigned int d = 0; d < dim; ++d)
1412 {
1413 Assert(gradients_quad != nullptr, ExcInternalError());
1414 if (integrate)
1415 for (unsigned int q = 0; q < n_q_points; ++q)
1416 tmp_values[q] =
1417 gradients_quad[(c * n_q_points + orientation[q]) * dim + d]
1418 [v];
1419 else
1420 for (unsigned int q = 0; q < n_q_points; ++q)
1421 tmp_values[orientation[q]] =
1422 gradients_quad[(c * n_q_points + q) * dim + d][v];
1423 for (unsigned int q = 0; q < n_q_points; ++q)
1424 gradients_quad[(c * n_q_points + q) * dim + d][v] =
1425 tmp_values[q];
1426 }
1427 if (flag & EvaluationFlags::hessians)
1428 {
1429 Assert(hessians_quad != nullptr, ExcInternalError());
1430 const unsigned int hdim = (dim * (dim + 1)) / 2;
1431 for (unsigned int d = 0; d < hdim; ++d)
1432 {
1433 if (integrate)
1434 for (unsigned int q = 0; q < n_q_points; ++q)
1435 tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
1436 orientation[q]][v];
1437 else
1438 for (unsigned int q = 0; q < n_q_points; ++q)
1439 tmp_values[orientation[q]] =
1440 hessians_quad[(c * hdim + d) * n_q_points + q][v];
1441 for (unsigned int q = 0; q < n_q_points; ++q)
1442 hessians_quad[(c * hdim + d) * n_q_points + q][v] =
1443 tmp_values[q];
1444 }
1445 }
1446 }
1447 }
1448
1449
1450
1451 template <int dim, typename Number>
1453 {
1454 static bool
1455 evaluate_tensor_none(const unsigned int n_components,
1456 const EvaluationFlags::EvaluationFlags evaluation_flag,
1457 const Number *values_dofs,
1459 {
1460 const auto &shape_info = fe_eval.get_shape_info();
1461 const auto &shape_data = shape_info.data.front();
1462 using Number2 =
1464
1465 Assert((fe_eval.get_dof_access_index() ==
1467 fe_eval.is_interior_face() == false) == false,
1469
1470 const unsigned int face_no = fe_eval.get_face_no();
1471 const unsigned int face_orientation = fe_eval.get_face_orientation();
1472 const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
1473 const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
1474
1475 if (evaluation_flag & EvaluationFlags::values)
1476 {
1477 const auto *const shape_values =
1478 &shape_data.shape_values_face(face_no, face_orientation, 0);
1479
1480 auto *out = fe_eval.begin_values();
1481 auto *in = values_dofs;
1482
1483 for (unsigned int c = 0; c < n_components; c += 3)
1484 {
1485 if (c + 1 == n_components)
1488 /*transpose_matrix*/ true,
1489 /*add*/ false,
1490 /*consider_strides*/ false,
1491 Number,
1492 Number2,
1493 /*n_components*/ 1>(
1494 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1495 else if (c + 2 == n_components)
1498 /*transpose_matrix*/ true,
1499 /*add*/ false,
1500 /*consider_strides*/ false,
1501 Number,
1502 Number2,
1503 /*n_components*/ 2>(
1504 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1505 else
1508 /*transpose_matrix*/ true,
1509 /*add*/ false,
1510 /*consider_strides*/ false,
1511 Number,
1512 Number2,
1513 /*n_components*/ 3>(
1514 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1515
1516 out += 3 * n_q_points;
1517 in += 3 * n_dofs;
1518 }
1519 }
1520
1521 if (evaluation_flag & EvaluationFlags::gradients)
1522 {
1523 auto *out = fe_eval.begin_gradients();
1524 const auto *in = values_dofs;
1525
1526 const auto *const shape_gradients =
1527 &shape_data.shape_gradients_face(face_no, face_orientation, 0);
1528
1529 for (unsigned int c = 0; c < n_components; c += 3)
1530 {
1531 if (c + 1 == n_components)
1534 /*transpose_matrix*/ true,
1535 /*add*/ false,
1536 /*consider_strides*/ false,
1537 Number,
1538 Number2,
1539 /*n_components*/ 1>(
1540 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
1541 else if (c + 2 == n_components)
1544 /*transpose_matrix*/ true,
1545 /*add*/ false,
1546 /*consider_strides*/ false,
1547 Number,
1548 Number2,
1549 /*n_components*/ 2>(
1550 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
1551 else
1554 /*transpose_matrix*/ true,
1555 /*add*/ false,
1556 /*consider_strides*/ false,
1557 Number,
1558 Number2,
1559 /*n_components*/ 3>(
1560 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
1561 out += 3 * n_q_points * dim;
1562 in += 3 * n_dofs;
1563 }
1564 }
1565
1566 Assert(!(evaluation_flag & EvaluationFlags::hessians),
1568
1569 return true;
1570 }
1571
1572 template <int fe_degree>
1573#ifndef DEBUG
1575#endif
1576 static void
1577 project_to_face(const unsigned int n_components,
1578 const EvaluationFlags::EvaluationFlags evaluation_flag,
1579 const Number *values_dofs,
1581 const bool use_vectorization,
1582 Number *temp,
1583 Number *scratch_data)
1584 {
1585 const auto &shape_info = fe_eval.get_shape_info();
1586
1587 if (use_vectorization == false)
1588 {
1589 const auto &shape_data = shape_info.data.front();
1590
1591 const unsigned int dofs_per_comp_face =
1592 fe_degree > -1 ?
1593 Utilities::pow(fe_degree + 1, dim - 1) :
1594 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1595 const unsigned int dofs_per_face = n_components * dofs_per_comp_face;
1596
1597 for (unsigned int v = 0; v < Number::size(); ++v)
1598 {
1599 // the loop breaks once an invalid_unsigned_int is hit for
1600 // all cases except the exterior faces in the ECL loop (where
1601 // some faces might be at the boundaries but others not)
1602 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
1603 {
1604 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
1605 temp[i][v] = 0;
1606 continue;
1607 }
1608
1610 template interpolate<true, false>(n_components,
1611 evaluation_flag,
1612 shape_info,
1613 values_dofs,
1614 scratch_data,
1615 fe_eval.get_face_no(v));
1616
1617 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
1618 temp[i][v] = scratch_data[i][v];
1619 }
1620 }
1621 else
1623 template interpolate<true, false>(n_components,
1624 evaluation_flag,
1625 shape_info,
1626 values_dofs,
1627 temp,
1628 fe_eval.get_face_no());
1629 }
1630
1631
1632 template <int fe_degree, int n_q_points_1d>
1633#ifndef DEBUG
1635#endif
1636 static void
1637 evaluate_in_face(const unsigned int n_components,
1638 const EvaluationFlags::EvaluationFlags evaluation_flag,
1640 Number *temp,
1641 Number *scratch_data)
1642 {
1643 const auto &shape_info = fe_eval.get_shape_info();
1644 const auto &shape_data = shape_info.data.front();
1645
1646 const unsigned int subface_index = fe_eval.get_subface_index();
1647 constexpr unsigned int n_q_points_1d_actual =
1648 fe_degree > -1 ? n_q_points_1d : 0;
1649
1650 if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
1651 {
1653 fe_degree,
1654 n_q_points_1d_actual,
1655 Number>::
1656 template evaluate_or_integrate_in_face<false>(
1657 evaluation_flag,
1658 fe_eval.get_shape_info().data,
1659 temp,
1660 fe_eval.begin_values(),
1661 fe_eval.begin_gradients(),
1662 scratch_data,
1663 subface_index,
1664 fe_eval.get_face_no() / 2);
1665 }
1666 else if (fe_degree > -1 &&
1668 shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
1670 dim,
1671 fe_degree,
1672 n_q_points_1d_actual,
1673 Number>::evaluate_in_face(n_components,
1674 evaluation_flag,
1675 shape_data,
1676 temp,
1677 fe_eval.begin_values(),
1678 fe_eval
1679 .begin_gradients(),
1680 fe_eval.begin_hessians(),
1681 scratch_data,
1682 subface_index);
1683 else
1685 dim,
1686 fe_degree,
1687 n_q_points_1d_actual,
1688 Number>::evaluate_in_face(n_components,
1689 evaluation_flag,
1690 shape_data,
1691 temp,
1692 fe_eval.begin_values(),
1693 fe_eval
1694 .begin_gradients(),
1695 fe_eval.begin_hessians(),
1696 scratch_data,
1697 subface_index);
1698 }
1699
1700#ifndef DEBUG
1702#endif
1703 static void
1705 const unsigned int n_components,
1706 const EvaluationFlags::EvaluationFlags evaluation_flag,
1708 const bool use_vectorization,
1709 Number *temp)
1710 {
1711 const auto &shape_info = fe_eval.get_shape_info();
1712
1713 if (use_vectorization == false)
1714 {
1715 for (unsigned int v = 0; v < Number::size(); ++v)
1716 {
1717 // the loop breaks once an invalid_unsigned_int is hit for
1718 // all cases except the exterior faces in the ECL loop (where
1719 // some faces might be at the boundaries but others not)
1720 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
1721 continue;
1722
1723 if (fe_eval.get_face_orientation(v) != 0)
1725 dim,
1726 n_components,
1727 v,
1728 evaluation_flag,
1729 &shape_info.face_orientations_quad(
1730 fe_eval.get_face_orientation(v), 0),
1731 false,
1732 shape_info.n_q_points_face,
1733 &temp[0][0],
1734 fe_eval.begin_values(),
1735 fe_eval.begin_gradients(),
1736 fe_eval.begin_hessians());
1737 }
1738 }
1739 else if (fe_eval.get_face_orientation() != 0)
1741 dim,
1742 n_components,
1743 evaluation_flag,
1744 &shape_info.face_orientations_quad(fe_eval.get_face_orientation(), 0),
1745 false,
1746 shape_info.n_q_points_face,
1747 temp,
1748 fe_eval.begin_values(),
1749 fe_eval.begin_gradients(),
1750 fe_eval.begin_hessians());
1751 }
1752
1753
1754
1755 template <int fe_degree, int n_q_points_1d>
1756 static bool
1757 evaluate_tensor(const unsigned int n_components,
1758 const EvaluationFlags::EvaluationFlags evaluation_flag,
1759 const Number *values_dofs_actual,
1761 {
1762 const auto &shape_info = fe_eval.get_shape_info();
1763 const auto &shape_data = shape_info.data.front();
1764
1765 const unsigned int dofs_per_comp_face =
1766 fe_degree > -1 ?
1767 Utilities::pow(fe_degree + 1, dim - 1) :
1768 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1769
1770 // Note: we always keep storage of values, 1st and 2nd derivatives in an
1771 // array, so reserve space for all three here
1772 Number *temp1 = fe_eval.get_scratch_data().begin();
1773 Number *temp2 = temp1 + 3 * n_components * dofs_per_comp_face;
1774
1775 const Number *values_dofs =
1776 (shape_data.element_type == MatrixFreeFunctions::truncated_tensor) ?
1777 temp2 +
1779 shape_info.n_q_points)) :
1780 values_dofs_actual;
1781
1782 if (shape_data.element_type == MatrixFreeFunctions::truncated_tensor)
1783 embed_truncated_into_full_tensor_product<dim, fe_degree>(
1784 n_components,
1785 const_cast<Number *>(values_dofs),
1786 values_dofs_actual,
1787 fe_eval);
1788
1789 bool use_vectorization = true;
1790 if (fe_eval.get_dof_access_index() ==
1792 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
1793 for (unsigned int v = 0; v < Number::size(); ++v)
1794 if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
1795 fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
1796 use_vectorization = false;
1797
1798 project_to_face<fe_degree>(n_components,
1799 evaluation_flag,
1800 values_dofs,
1801 fe_eval,
1802 use_vectorization,
1803 temp1,
1804 temp2);
1805
1806 evaluate_in_face<fe_degree, n_q_points_1d>(
1807 n_components, evaluation_flag, fe_eval, temp1, temp2);
1808
1809 if (dim == 3)
1811 n_components, evaluation_flag, fe_eval, use_vectorization, temp1);
1812
1813 return false;
1814 }
1815
1816 template <int fe_degree, int n_q_points_1d>
1817 static bool
1818 run(const unsigned int n_components,
1819 const EvaluationFlags::EvaluationFlags evaluation_flag,
1820 const Number *values_dofs,
1822 {
1823 const auto &shape_info = fe_eval.get_shape_info();
1824
1825 if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
1826 return evaluate_tensor_none(n_components,
1827 evaluation_flag,
1828 values_dofs,
1829 fe_eval);
1830 else
1831 return evaluate_tensor<fe_degree, n_q_points_1d>(n_components,
1832 evaluation_flag,
1833 values_dofs,
1834 fe_eval);
1835 }
1836 };
1837
1838
1839
1840 template <int dim, typename Number>
1842 {
1843 template <int fe_degree>
1844 static bool
1845 run(const unsigned int n_components,
1846 const EvaluationFlags::EvaluationFlags evaluation_flag,
1847 const Number *values_dofs,
1849 {
1850 const auto &shape_info = fe_eval.get_shape_info();
1851 const auto &shape_data = shape_info.data.front();
1852
1853 const unsigned int dofs_per_comp_face =
1854 fe_degree > -1 ?
1855 Utilities::pow(fe_degree + 1, dim - 1) :
1856 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1857
1858 // Note: we always keep storage of values, 1st and 2nd derivatives in an
1859 // array, so reserve space for all three here
1860 Number *temp = fe_eval.get_scratch_data().begin();
1861 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
1862
1863 bool use_vectorization = true;
1864 if (fe_eval.get_dof_access_index() ==
1866 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
1867 for (unsigned int v = 0; v < Number::size(); ++v)
1868 if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
1869 fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
1870 use_vectorization = false;
1871
1873 template project_to_face<fe_degree>(n_components,
1874 evaluation_flag,
1875 values_dofs,
1876 fe_eval,
1877 use_vectorization,
1878 temp,
1879 scratch_data);
1880
1881 return false;
1882 }
1883 };
1884
1885
1886
1887 template <int dim, typename Number>
1889 {
1890 template <int fe_degree, int n_q_points_1d>
1891 static bool
1892 run(const unsigned int n_components,
1893 const EvaluationFlags::EvaluationFlags evaluation_flag,
1895 {
1896 const auto &shape_info = fe_eval.get_shape_info();
1897 const auto &shape_data = shape_info.data.front();
1898
1899 const unsigned int dofs_per_comp_face =
1900 fe_degree > -1 ?
1901 Utilities::pow(fe_degree + 1, dim - 1) :
1902 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1903
1904 // Note: we always keep storage of values, 1st and 2nd derivatives in an
1905 // array, so reserve space for all three here
1906 Number *temp = fe_eval.get_scratch_data().begin();
1907 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
1908
1910 template evaluate_in_face<fe_degree, n_q_points_1d>(
1911 n_components, evaluation_flag, fe_eval, temp, scratch_data);
1912
1913 return false;
1914 }
1915 };
1916
1917
1918
1919 template <int dim, typename Number>
1921 {
1922 static bool
1924 const unsigned int n_components,
1925 const EvaluationFlags::EvaluationFlags integration_flag,
1926 Number *values_dofs,
1928 const bool sum_into_values)
1929 {
1930 const auto &shape_info = fe_eval.get_shape_info();
1931 const auto &shape_data = shape_info.data.front();
1932 using Number2 =
1934
1935 Assert((fe_eval.get_dof_access_index() ==
1937 fe_eval.is_interior_face() == false) == false,
1939
1940 const unsigned int face_no = fe_eval.get_face_no();
1941 const unsigned int face_orientation = fe_eval.get_face_orientation();
1942 const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
1943 const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
1944
1945
1946 if (integration_flag & EvaluationFlags::values)
1947 {
1948 const auto *const shape_values =
1949 &shape_data.shape_values_face(face_no, face_orientation, 0);
1950
1951 auto *in = fe_eval.begin_values();
1952 auto *out = values_dofs;
1953
1954 for (unsigned int c = 0; c < n_components; c += 3)
1955 {
1956 if (sum_into_values)
1957 {
1958 if (c + 1 == n_components)
1961 /*transpose_matrix*/ false,
1962 /*add*/ true,
1963 /*consider_strides*/ false,
1964 Number,
1965 Number2,
1966 /*n_components*/ 1>(
1967 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1968 else if (c + 2 == n_components)
1971 /*transpose_matrix*/ false,
1972 /*add*/ true,
1973 /*consider_strides*/ false,
1974 Number,
1975 Number2,
1976 /*n_components*/ 2>(
1977 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1978 else
1981 /*transpose_matrix*/ false,
1982 /*add*/ true,
1983 /*consider_strides*/ false,
1984 Number,
1985 Number2,
1986 /*n_components*/ 3>(
1987 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1988 }
1989 else
1990 {
1991 if (c + 1 == n_components)
1994 /*transpose_matrix*/ false,
1995 /*add*/ false,
1996 /*consider_strides*/ false,
1997 Number,
1998 Number2,
1999 /*n_components*/ 1>(
2000 shape_values, in, out, n_dofs, n_q_points, 1, 1);
2001 else if (c + 2 == n_components)
2004 /*transpose_matrix*/ false,
2005 /*add*/ false,
2006 /*consider_strides*/ false,
2007 Number,
2008 Number2,
2009 /*n_components*/ 2>(
2010 shape_values, in, out, n_dofs, n_q_points, 1, 1);
2011 else
2014 /*transpose_matrix*/ false,
2015 /*add*/ false,
2016 /*consider_strides*/ false,
2017 Number,
2018 Number2,
2019 /*n_components*/ 3>(
2020 shape_values, in, out, n_dofs, n_q_points, 1, 1);
2021 }
2022 in += 3 * n_q_points;
2023 out += 3 * n_dofs;
2024 }
2025 }
2026
2027 if (integration_flag & EvaluationFlags::gradients)
2028 {
2029 auto *in = fe_eval.begin_gradients();
2030 auto *out = values_dofs;
2031
2032 const auto *const shape_gradients =
2033 &shape_data.shape_gradients_face(face_no, face_orientation, 0);
2034
2035 for (unsigned int c = 0; c < n_components; ++c)
2036 {
2037 if (!sum_into_values &&
2038 !(integration_flag & EvaluationFlags::values))
2039 {
2040 if (c + 1 == n_components)
2043 /*transpose_matrix*/ false,
2044 /*add*/ false,
2045 /*consider_strides*/ false,
2046 Number,
2047 Number2,
2048 /*n_components*/ 1>(
2049 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2050 else if (c + 2 == n_components)
2053 /*transpose_matrix*/ false,
2054 /*add*/ false,
2055 /*consider_strides*/ false,
2056 Number,
2057 Number2,
2058 /*n_components*/ 2>(
2059 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2060 else
2063 /*transpose_matrix*/ false,
2064 /*add*/ false,
2065 /*consider_strides*/ false,
2066 Number,
2067 Number2,
2068 /*n_components*/ 3>(
2069 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2070 }
2071 else
2072 {
2073 if (c + 1 == n_components)
2076 /*transpose_matrix*/ false,
2077 /*add*/ true,
2078 /*consider_strides*/ false,
2079 Number,
2080 Number2,
2081 /*n_components*/ 1>(
2082 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2083 else if (c + 2 == n_components)
2086 /*transpose_matrix*/ false,
2087 /*add*/ true,
2088 /*consider_strides*/ false,
2089 Number,
2090 Number2,
2091 /*n_components*/ 2>(
2092 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2093 else
2096 /*transpose_matrix*/ false,
2097 /*add*/ true,
2098 /*consider_strides*/ false,
2099 Number,
2100 Number2,
2101 /*n_components*/ 3>(
2102 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2103 }
2104 in += 3 * n_q_points * dim;
2105 out += 3 * n_dofs;
2106 }
2107 }
2108
2109 Assert(!(integration_flag & EvaluationFlags::hessians),
2111
2112 return true;
2113 }
2114
2115#ifndef DEBUG
2117#endif
2118 static void
2120 const unsigned int n_components,
2121 const EvaluationFlags::EvaluationFlags integration_flag,
2123 const bool use_vectorization,
2124 Number *temp)
2125 {
2126 const auto &shape_info = fe_eval.get_shape_info();
2127
2128 if (use_vectorization == false)
2129 {
2130 for (unsigned int v = 0; v < Number::size(); ++v)
2131 {
2132 // the loop breaks once an invalid_unsigned_int is hit for
2133 // all cases except the exterior faces in the ECL loop (where
2134 // some faces might be at the boundaries but others not)
2135 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2136 continue;
2137
2138 if (fe_eval.get_face_orientation(v) != 0)
2140 dim,
2141 n_components,
2142 v,
2143 integration_flag,
2145 fe_eval.get_face_orientation(v), 0),
2146 true,
2147 shape_info.n_q_points_face,
2148 &temp[0][0],
2149 fe_eval.begin_values(),
2150 fe_eval.begin_gradients(),
2151 fe_eval.begin_hessians());
2152 }
2153 }
2154 else if (fe_eval.get_face_orientation() != 0)
2156 dim,
2157 n_components,
2158 integration_flag,
2160 fe_eval.get_face_orientation(), 0),
2161 true,
2162 shape_info.n_q_points_face,
2163 temp,
2164 fe_eval.begin_values(),
2165 fe_eval.begin_gradients(),
2166 fe_eval.begin_hessians());
2167 }
2168
2169 template <int fe_degree, int n_q_points_1d>
2170#ifndef DEBUG
2172#endif
2173 static void
2174 integrate_in_face(const unsigned int n_components,
2175 const EvaluationFlags::EvaluationFlags integration_flag,
2177 Number *temp,
2178 Number *scratch_data)
2179 {
2180 const auto &shape_info = fe_eval.get_shape_info();
2181 const auto &shape_data = shape_info.data.front();
2182
2183 const unsigned int n_q_points_1d_actual =
2184 fe_degree > -1 ? n_q_points_1d : 0;
2185 const unsigned int subface_index = fe_eval.get_subface_index();
2186
2187 if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
2188 {
2190 fe_degree,
2191 n_q_points_1d_actual,
2192 Number>::
2193 template evaluate_or_integrate_in_face<true>(
2194 integration_flag,
2195 fe_eval.get_shape_info().data,
2196 temp,
2197 fe_eval.begin_values(),
2198 fe_eval.begin_gradients(),
2199 scratch_data,
2200 subface_index,
2201 fe_eval.get_face_no() / 2);
2202 }
2203 else if (fe_degree > -1 &&
2204 fe_eval.get_subface_index() >=
2206 shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
2208 true,
2209 dim,
2210 fe_degree,
2211 n_q_points_1d_actual,
2212 Number>::integrate_in_face(n_components,
2213 integration_flag,
2214 shape_data,
2215 temp,
2216 fe_eval.begin_values(),
2217 fe_eval.begin_gradients(),
2218 fe_eval.begin_hessians(),
2219 scratch_data,
2220 subface_index);
2221 else
2223 false,
2224 dim,
2225 fe_degree,
2226 n_q_points_1d_actual,
2227 Number>::integrate_in_face(n_components,
2228 integration_flag,
2229 shape_data,
2230 temp,
2231 fe_eval.begin_values(),
2232 fe_eval.begin_gradients(),
2233 fe_eval.begin_hessians(),
2234 scratch_data,
2235 subface_index);
2236 }
2237
2238 template <int fe_degree>
2239#ifndef DEBUG
2241#endif
2242 static void
2243 collect_from_face(const unsigned int n_components,
2244 const EvaluationFlags::EvaluationFlags integration_flag,
2245 Number *values_dofs,
2247 const bool use_vectorization,
2248 const Number *temp,
2249 Number *scratch_data,
2250 const bool sum_into_values)
2251 {
2252 const auto &shape_info = fe_eval.get_shape_info();
2253 const auto &shape_data = shape_info.data.front();
2254
2255 const unsigned int dofs_per_comp_face =
2256 fe_degree > -1 ?
2257 Utilities::pow(fe_degree + 1, dim - 1) :
2258 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2259 const unsigned int dofs_per_face = n_components * dofs_per_comp_face;
2260
2261 if (use_vectorization == false)
2262 {
2263 for (unsigned int v = 0; v < Number::size(); ++v)
2264 {
2265 // the loop breaks once an invalid_unsigned_int is hit for
2266 // all cases except the exterior faces in the ECL loop (where
2267 // some faces might be at the boundaries but others not)
2268 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2269 continue;
2270
2272 template interpolate<false, false>(n_components,
2273 integration_flag,
2274 shape_info,
2275 temp,
2276 scratch_data,
2277 fe_eval.get_face_no(v));
2278
2279 if (sum_into_values)
2280 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
2281 values_dofs[i][v] += scratch_data[i][v];
2282 else
2283 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
2284 values_dofs[i][v] = scratch_data[i][v];
2285 }
2286 }
2287 else
2288 {
2289 if (sum_into_values)
2291 template interpolate<false, true>(n_components,
2292 integration_flag,
2293 shape_info,
2294 temp,
2295 values_dofs,
2296 fe_eval.get_face_no());
2297 else
2299 template interpolate<false, false>(n_components,
2300 integration_flag,
2301 shape_info,
2302 temp,
2303 values_dofs,
2304 fe_eval.get_face_no());
2305 }
2306 }
2307
2308 template <int fe_degree, int n_q_points_1d>
2309 static bool
2310 integrate_tensor(const unsigned int n_components,
2311 const EvaluationFlags::EvaluationFlags integration_flag,
2312 Number *values_dofs_actual,
2314 const bool sum_into_values)
2315 {
2316 const auto &shape_info = fe_eval.get_shape_info();
2317 const auto &shape_data = shape_info.data.front();
2318
2319 const unsigned int dofs_per_comp_face =
2320 fe_degree > -1 ?
2321 Utilities::pow(fe_degree + 1, dim - 1) :
2322 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2323
2324 Number *temp1 = fe_eval.get_scratch_data().begin();
2325 Number *temp2 = temp1 + 3 * n_components * dofs_per_comp_face;
2326
2327 // expand dof_values to tensor product for truncated tensor products
2328 Number *values_dofs =
2329 (shape_data.element_type == MatrixFreeFunctions::truncated_tensor) ?
2330 temp2 + 2 * (std::max<std::size_t>(
2332 fe_eval.get_shape_info().n_q_points)) :
2333 values_dofs_actual;
2334
2335 bool use_vectorization = true;
2336
2337 if (fe_eval.get_dof_access_index() ==
2339 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
2340 use_vectorization =
2342 std::all_of(fe_eval.get_cell_ids().begin() + 1,
2343 fe_eval.get_cell_ids().end(),
2344 [&](const auto &v) {
2345 return v == fe_eval.get_cell_ids()[0] ||
2346 v == numbers::invalid_unsigned_int;
2347 });
2348
2349 if (dim == 3)
2351 n_components, integration_flag, fe_eval, use_vectorization, temp1);
2352
2353 integrate_in_face<fe_degree, n_q_points_1d>(
2354 n_components, integration_flag, fe_eval, temp1, temp2);
2355
2356 collect_from_face<fe_degree>(n_components,
2357 integration_flag,
2358 values_dofs,
2359 fe_eval,
2360 use_vectorization,
2361 temp1,
2362 temp2,
2363 sum_into_values);
2364
2365
2366 if (shape_data.element_type == MatrixFreeFunctions::truncated_tensor)
2367 truncate_tensor_product_to_complete_degrees<dim, fe_degree>(
2368 n_components, values_dofs_actual, values_dofs, fe_eval);
2369
2370 return false;
2371 }
2372
2373 template <int fe_degree, int n_q_points_1d>
2374 static bool
2375 run(const unsigned int n_components,
2376 const EvaluationFlags::EvaluationFlags integration_flag,
2377 Number *values_dofs,
2379 const bool sum_into_values)
2380 {
2381 const auto &shape_info = fe_eval.get_shape_info();
2382
2383 if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
2384 return integrate_tensor_none(n_components,
2385 integration_flag,
2386 values_dofs,
2387 fe_eval,
2388 sum_into_values);
2389 else
2390 return integrate_tensor<fe_degree, n_q_points_1d>(n_components,
2391 integration_flag,
2392 values_dofs,
2393 fe_eval,
2394 sum_into_values);
2395 }
2396 };
2397
2398
2399
2400 template <int dim, typename Number>
2402 {
2403 template <int fe_degree>
2404 static bool
2405 run(const unsigned int n_components,
2406 const EvaluationFlags::EvaluationFlags integration_flag,
2407 Number *values_dofs,
2409 const bool sum_into_values)
2410 {
2411 const auto &shape_info = fe_eval.get_shape_info();
2412 const auto &shape_data = shape_info.data.front();
2413
2414 const unsigned int dofs_per_comp_face =
2415 fe_degree > -1 ?
2416 Utilities::pow(fe_degree + 1, dim - 1) :
2417 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2418
2419 Number *temp = fe_eval.get_scratch_data().begin();
2420 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
2421
2422 bool use_vectorization = true;
2423
2424 if (fe_eval.get_dof_access_index() ==
2426 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
2427 use_vectorization =
2429 std::all_of(fe_eval.get_cell_ids().begin() + 1,
2430 fe_eval.get_cell_ids().end(),
2431 [&](const auto &v) {
2432 return v == fe_eval.get_cell_ids()[0] ||
2433 v == numbers::invalid_unsigned_int;
2434 });
2435
2437 template collect_from_face<fe_degree>(n_components,
2438 integration_flag,
2439 values_dofs,
2440 fe_eval,
2441 use_vectorization,
2442 temp,
2443 scratch_data,
2444 sum_into_values);
2445
2446 return false;
2447 }
2448 };
2449
2450
2451
2452 template <int dim, typename Number>
2454 {
2455 template <int fe_degree, int n_q_points_1d>
2456 static bool
2457 run(const unsigned int n_components,
2458 const EvaluationFlags::EvaluationFlags integration_flag,
2459
2461 {
2462 const auto &shape_info = fe_eval.get_shape_info();
2463 const auto &shape_data = shape_info.data.front();
2464
2465 const unsigned int dofs_per_comp_face =
2466 fe_degree > -1 ?
2467 Utilities::pow(fe_degree + 1, dim - 1) :
2468 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2469
2470 Number *temp = fe_eval.get_scratch_data().begin();
2471 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
2472
2474 template integrate_in_face<fe_degree, n_q_points_1d>(
2475 n_components, integration_flag, fe_eval, temp, scratch_data);
2476
2477 return false;
2478 }
2479 };
2480
2481
2482
2483 template <int n_face_orientations,
2484 typename Processor,
2485 typename EvaluationData,
2486 const bool check_face_orientations = false>
2487 void
2489 Processor &proc,
2490 const unsigned int n_components,
2491 const EvaluationFlags::EvaluationFlags evaluation_flag,
2492 typename Processor::Number2_ *global_vector_ptr,
2493 const std::vector<ArrayView<const typename Processor::Number2_>> *sm_ptr,
2494 const EvaluationData &fe_eval,
2495 typename Processor::VectorizedArrayType_ *temp1)
2496 {
2497 constexpr int dim = Processor::dim_;
2498 constexpr int fe_degree = Processor::fe_degree_;
2499 using VectorizedArrayType = typename Processor::VectorizedArrayType_;
2500 constexpr int n_lanes = VectorizedArrayType::size();
2501
2502 using Number = typename Processor::Number_;
2503 using Number2_ = typename Processor::Number2_;
2504
2505 const auto &shape_data = fe_eval.get_shape_info().data.front();
2506 constexpr bool integrate = Processor::do_integrate;
2507 const unsigned int face_no = fe_eval.get_face_no();
2508 const auto &dof_info = fe_eval.get_dof_info();
2509 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
2510 const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index =
2511 fe_eval.get_dof_access_index();
2512 AssertIndexRange(cell,
2513 dof_info.index_storage_variants[dof_access_index].size());
2514 constexpr unsigned int dofs_per_face =
2515 Utilities::pow(fe_degree + 1, dim - 1);
2516 const unsigned int subface_index = fe_eval.get_subface_index();
2517
2518 const unsigned int n_filled_lanes =
2519 dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
2520
2521 bool all_faces_are_same = n_filled_lanes == n_lanes;
2522 if (n_face_orientations == n_lanes)
2523 for (unsigned int v = 1; v < n_lanes; ++v)
2524 if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) ||
2525 fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0))
2526 {
2527 all_faces_are_same = false;
2528 break;
2529 }
2530
2531 // check for re-orientation ...
2532 std::array<const unsigned int *, n_face_orientations> orientation = {};
2533
2534 if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same &&
2535 fe_eval.is_interior_face() == 0)
2536 for (unsigned int v = 0; v < n_lanes; ++v)
2537 {
2538 // the loop breaks once an invalid_unsigned_int is hit for
2539 // all cases except the exterior faces in the ECL loop (where
2540 // some faces might be at the boundaries but others not)
2541 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2542 continue;
2543
2544 if (shape_data.nodal_at_cell_boundaries &&
2545 fe_eval.get_face_orientation(v) != 0)
2546 {
2547 // ... and in case we detect a re-orientation, go to the other
2548 // version of this function that actually allows for this
2549 if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
2550 check_face_orientations == false)
2551 {
2552 fe_face_evaluation_process_and_io<n_face_orientations,
2553 Processor,
2554 EvaluationData,
2555 true>(proc,
2556 n_components,
2557 evaluation_flag,
2558 global_vector_ptr,
2559 sm_ptr,
2560 fe_eval,
2561 temp1);
2562 return;
2563 }
2564 orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
2565 fe_eval.get_face_orientation(v), 0);
2566 }
2567 }
2568 else if (dim == 3 && fe_eval.get_face_orientation() != 0)
2569 {
2570 // go to the other version of this function
2571 if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
2572 check_face_orientations == false)
2573 {
2574 fe_face_evaluation_process_and_io<n_face_orientations,
2575 Processor,
2576 EvaluationData,
2577 true>(proc,
2578 n_components,
2579 evaluation_flag,
2580 global_vector_ptr,
2581 sm_ptr,
2582 fe_eval,
2583 temp1);
2584 return;
2585 }
2586 for (unsigned int v = 0; v < n_face_orientations; ++v)
2587 orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
2588 fe_eval.get_face_orientation(), 0);
2589 }
2590
2591 // we know that the gradient weights for the Hermite case on the
2592 // right (side==1) are the negative from the value at the left
2593 // (side==0), so we only read out one of them.
2594 VectorizedArrayType grad_weight =
2595 shape_data
2596 .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) :
2597 (1 + face_no % 2))];
2598
2599 // face_to_cell_index_hermite
2600 std::array<const unsigned int *, n_face_orientations> index_array_hermite =
2601 {};
2602 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2603 {
2604 if (n_face_orientations == 1)
2605 index_array_hermite[0] =
2606 &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0);
2607 else
2608 {
2609 for (unsigned int v = 0; v < n_lanes; ++v)
2610 {
2611 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2612 continue;
2613
2614 const auto face_no = fe_eval.get_face_no(v);
2615
2616 grad_weight[v] =
2617 shape_data
2618 .shape_data_on_face[0][fe_degree + (integrate ?
2619 (2 - (face_no % 2)) :
2620 (1 + (face_no % 2)))];
2621
2622 index_array_hermite[v] =
2623 &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no,
2624 0);
2625 }
2626 }
2627 }
2628
2629 // face_to_cell_index_nodal
2630 std::array<const unsigned int *, n_face_orientations> index_array_nodal =
2631 {};
2632 if (shape_data.nodal_at_cell_boundaries == true)
2633 {
2634 if (n_face_orientations == 1)
2635 index_array_nodal[0] =
2636 &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0);
2637 else
2638 {
2639 for (unsigned int v = 0; v < n_lanes; ++v)
2640 {
2641 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2642 continue;
2643
2644 const auto face_no = fe_eval.get_face_no(v);
2645
2646 index_array_nodal[v] =
2647 &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no,
2648 0);
2649 }
2650 }
2651 }
2652
2653
2654 const auto reorientate = [&](const unsigned int v, const unsigned int i) {
2655 return (!check_face_orientations || orientation[v] == nullptr) ?
2656 i :
2657 orientation[v][i];
2658 };
2659
2660 const unsigned int cell_index =
2662 fe_eval.get_cell_ids()[0] :
2663 cell * n_lanes;
2664 const unsigned int *dof_indices =
2665 &dof_info.dof_indices_contiguous[dof_access_index][cell_index];
2666
2667 for (unsigned int comp = 0; comp < n_components; ++comp)
2668 {
2669 const std::size_t index_offset =
2670 dof_info.component_dof_indices_offset
2671 [fe_eval.get_active_fe_index()]
2672 [fe_eval.get_first_selected_component()] +
2673 comp * Utilities::pow(fe_degree + 1, dim);
2674
2675 // case 1: contiguous and interleaved indices
2676 if (n_face_orientations == 1 &&
2677 dof_info.index_storage_variants[dof_access_index][cell] ==
2679 interleaved_contiguous)
2680 {
2682 dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2683 n_lanes);
2684 Number2_ *vector_ptr =
2685 global_vector_ptr + dof_indices[0] + index_offset * n_lanes;
2686
2687 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2688 {
2689 for (unsigned int i = 0; i < dofs_per_face; ++i)
2690 {
2691 Assert(n_face_orientations == 1, ExcNotImplemented());
2692
2693 const unsigned int ind1 = index_array_hermite[0][2 * i];
2694 const unsigned int ind2 = index_array_hermite[0][2 * i + 1];
2695 const unsigned int i_ = reorientate(0, i);
2696 proc.hermite_grad_vectorized(temp1[i_],
2697 temp1[i_ + dofs_per_face],
2698 vector_ptr + ind1 * n_lanes,
2699 vector_ptr + ind2 * n_lanes,
2700 grad_weight);
2701 }
2702 }
2703 else
2704 {
2705 for (unsigned int i = 0; i < dofs_per_face; ++i)
2706 {
2707 Assert(n_face_orientations == 1, ExcNotImplemented());
2708
2709 const unsigned int i_ = reorientate(0, i);
2710 const unsigned int ind = index_array_nodal[0][i];
2711 proc.value_vectorized(temp1[i_],
2712 vector_ptr + ind * n_lanes);
2713 }
2714 }
2715 }
2716
2717 // case 2: contiguous and interleaved indices with fixed stride
2718 else if (n_face_orientations == 1 &&
2719 dof_info.index_storage_variants[dof_access_index][cell] ==
2721 interleaved_contiguous_strided)
2722 {
2724 dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2725 n_lanes);
2726 Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes;
2727 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2728 {
2729 for (unsigned int i = 0; i < dofs_per_face; ++i)
2730 {
2731 Assert(n_face_orientations == 1, ExcNotImplemented());
2732
2733 const unsigned int i_ = reorientate(0, i);
2734 const unsigned int ind1 =
2735 index_array_hermite[0][2 * i] * n_lanes;
2736 const unsigned int ind2 =
2737 index_array_hermite[0][2 * i + 1] * n_lanes;
2738 proc.hermite_grad_vectorized_indexed(
2739 temp1[i_],
2740 temp1[i_ + dofs_per_face],
2741 vector_ptr + ind1,
2742 vector_ptr + ind2,
2743 grad_weight,
2744 dof_indices,
2745 dof_indices);
2746 }
2747 }
2748 else
2749 {
2750 for (unsigned int i = 0; i < dofs_per_face; ++i)
2751 {
2752 Assert(n_face_orientations == 1, ExcNotImplemented());
2753
2754 const unsigned int i_ = reorientate(0, i);
2755 const unsigned int ind = index_array_nodal[0][i] * n_lanes;
2756 proc.value_vectorized_indexed(temp1[i_],
2757 vector_ptr + ind,
2758 dof_indices);
2759 }
2760 }
2761 }
2762
2763 // case 3: contiguous and interleaved indices with mixed stride
2764 else if (n_face_orientations == 1 &&
2765 dof_info.index_storage_variants[dof_access_index][cell] ==
2767 interleaved_contiguous_mixed_strides)
2768 {
2769 const unsigned int *strides =
2770 &dof_info.dof_indices_interleave_strides[dof_access_index]
2771 [cell * n_lanes];
2772 unsigned int indices[n_lanes];
2773 for (unsigned int v = 0; v < n_lanes; ++v)
2774 indices[v] = dof_indices[v] + index_offset * strides[v];
2775 const unsigned int n_filled_lanes =
2776 dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
2777
2778 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2779 {
2780 if (n_filled_lanes == n_lanes)
2781 for (unsigned int i = 0; i < dofs_per_face; ++i)
2782 {
2783 Assert(n_face_orientations == 1, ExcNotImplemented());
2784
2785 const unsigned int i_ = reorientate(0, i);
2786 unsigned int ind1[n_lanes];
2788 for (unsigned int v = 0; v < n_lanes; ++v)
2789 ind1[v] = indices[v] +
2790 index_array_hermite[0][2 * i] * strides[v];
2791 unsigned int ind2[n_lanes];
2793 for (unsigned int v = 0; v < n_lanes; ++v)
2794 ind2[v] =
2795 indices[v] +
2796 // TODO
2797 index_array_hermite[0][2 * i + 1] * strides[v];
2798 proc.hermite_grad_vectorized_indexed(
2799 temp1[i_],
2800 temp1[i_ + dofs_per_face],
2801 global_vector_ptr,
2802 global_vector_ptr,
2803 grad_weight,
2804 ind1,
2805 ind2);
2806 }
2807 else
2808 {
2809 if (integrate == false)
2810 for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
2811 temp1[i] = VectorizedArrayType();
2812
2813 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2814 for (unsigned int i = 0; i < dofs_per_face; ++i)
2815 {
2816 const unsigned int i_ =
2817 reorientate(n_face_orientations == 1 ? 0 : v, i);
2818 proc.hermite_grad(
2819 temp1[i_][v],
2820 temp1[i_ + dofs_per_face][v],
2821 global_vector_ptr
2822 [indices[v] +
2823 index_array_hermite
2824 [n_face_orientations == 1 ? 0 : v][2 * i] *
2825 strides[v]],
2826 global_vector_ptr
2827 [indices[v] +
2828 index_array_hermite[n_face_orientations == 1 ?
2829 0 :
2830 v][2 * i + 1] *
2831 strides[v]],
2832 grad_weight[n_face_orientations == 1 ? 0 : v]);
2833 }
2834 }
2835 }
2836 else
2837 {
2838 if (n_filled_lanes == n_lanes)
2839 for (unsigned int i = 0; i < dofs_per_face; ++i)
2840 {
2841 Assert(n_face_orientations == 1, ExcInternalError());
2842 unsigned int ind[n_lanes];
2844 for (unsigned int v = 0; v < n_lanes; ++v)
2845 ind[v] =
2846 indices[v] + index_array_nodal[0][i] * strides[v];
2847 const unsigned int i_ = reorientate(0, i);
2848 proc.value_vectorized_indexed(temp1[i_],
2849 global_vector_ptr,
2850 ind);
2851 }
2852 else
2853 {
2854 if (integrate == false)
2855 for (unsigned int i = 0; i < dofs_per_face; ++i)
2856 temp1[i] = VectorizedArrayType();
2857
2858 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2859 for (unsigned int i = 0; i < dofs_per_face; ++i)
2860 proc.value(
2861 temp1[reorientate(n_face_orientations == 1 ? 0 : v,
2862 i)][v],
2863 global_vector_ptr
2864 [indices[v] +
2865 index_array_nodal[n_face_orientations == 1 ? 0 : v]
2866 [i] *
2867 strides[v]]);
2868 }
2869 }
2870 }
2871
2872 // case 4: contiguous indices without interleaving
2873 else if (n_face_orientations > 1 ||
2874 dof_info.index_storage_variants[dof_access_index][cell] ==
2876 contiguous)
2877 {
2878 Number2_ *vector_ptr = global_vector_ptr + index_offset;
2879
2880 const bool vectorization_possible =
2881 all_faces_are_same && (sm_ptr == nullptr);
2882
2883 std::array<Number2_ *, n_lanes> vector_ptrs{{nullptr}};
2884 std::array<unsigned int, n_lanes> reordered_indices{
2886
2887 if (vectorization_possible == false)
2888 {
2889 if (n_face_orientations == 1)
2890 {
2891 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2892 if (sm_ptr == nullptr)
2893 {
2894 vector_ptrs[v] = vector_ptr + dof_indices[v];
2895 }
2896 else
2897 {
2898 const auto &temp =
2899 dof_info
2900 .dof_indices_contiguous_sm[dof_access_index]
2901 [cell * n_lanes + v];
2902 vector_ptrs[v] = const_cast<Number2_ *>(
2903 sm_ptr->operator[](temp.first).data() +
2904 temp.second + index_offset);
2905 }
2906 }
2907 else if (n_face_orientations == n_lanes)
2908 {
2909 const auto &cells = fe_eval.get_cell_ids();
2910 for (unsigned int v = 0; v < n_lanes; ++v)
2911 if (cells[v] != numbers::invalid_unsigned_int)
2912 {
2913 if (sm_ptr == nullptr)
2914 {
2915 vector_ptrs[v] =
2916 vector_ptr +
2917 dof_info
2918 .dof_indices_contiguous[dof_access_index]
2919 [cells[v]];
2920 }
2921 else
2922 {
2923 const auto &temp =
2924 dof_info
2925 .dof_indices_contiguous_sm[dof_access_index]
2926 [cells[v]];
2927 vector_ptrs[v] = const_cast<Number2_ *>(
2928 sm_ptr->operator[](temp.first).data() +
2929 temp.second + index_offset);
2930 }
2931 }
2932 }
2933 else
2934 {
2936 }
2937 }
2938 else if (n_face_orientations == n_lanes)
2939 {
2940 for (unsigned int v = 0; v < n_lanes; ++v)
2941 reordered_indices[v] =
2942 dof_info.dof_indices_contiguous[dof_access_index]
2943 [fe_eval.get_cell_ids()[v]];
2944 dof_indices = reordered_indices.data();
2945 }
2946
2947 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2948 {
2949 if (vectorization_possible)
2950 for (unsigned int i = 0; i < dofs_per_face; ++i)
2951 {
2952 const unsigned int ind1 = index_array_hermite[0][2 * i];
2953 const unsigned int ind2 =
2954 index_array_hermite[0][2 * i + 1];
2955 const unsigned int i_ = reorientate(0, i);
2956
2957 proc.hermite_grad_vectorized_indexed(
2958 temp1[i_],
2959 temp1[i_ + dofs_per_face],
2960 vector_ptr + ind1,
2961 vector_ptr + ind2,
2962 grad_weight,
2963 dof_indices,
2964 dof_indices);
2965 }
2966 else if (n_face_orientations == 1)
2967 for (unsigned int i = 0; i < dofs_per_face; ++i)
2968 {
2969 const unsigned int ind1 = index_array_hermite[0][2 * i];
2970 const unsigned int ind2 =
2971 index_array_hermite[0][2 * i + 1];
2972 const unsigned int i_ = reorientate(0, i);
2973
2974 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2975 proc.hermite_grad(temp1[i_][v],
2976 temp1[i_ + dofs_per_face][v],
2977 vector_ptrs[v][ind1],
2978 vector_ptrs[v][ind2],
2979 grad_weight[v]);
2980
2981 if (integrate == false)
2982 for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
2983 {
2984 temp1[i][v] = 0.0;
2985 temp1[i + dofs_per_face][v] = 0.0;
2986 }
2987 }
2988 else
2989 {
2990 if (integrate == false && n_filled_lanes < n_lanes)
2991 for (unsigned int i = 0; i < dofs_per_face; ++i)
2992 temp1[i] = temp1[i + dofs_per_face] = Number();
2993
2994 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2995 for (unsigned int i = 0; i < dofs_per_face; ++i)
2996 proc.hermite_grad(
2997 temp1[reorientate(v, i)][v],
2998 temp1[reorientate(v, i) + dofs_per_face][v],
2999 vector_ptrs[v][index_array_hermite[v][2 * i]],
3000 vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
3001 grad_weight[v]);
3002 }
3003 }
3004 else
3005 {
3006 if (vectorization_possible)
3007 for (unsigned int i = 0; i < dofs_per_face; ++i)
3008 {
3009 const unsigned int ind = index_array_nodal[0][i];
3010 const unsigned int i_ = reorientate(0, i);
3011
3012 proc.value_vectorized_indexed(temp1[i_],
3013 vector_ptr + ind,
3014 dof_indices);
3015 }
3016 else
3017 {
3018 if constexpr (n_face_orientations == 1)
3019 for (unsigned int i = 0; i < dofs_per_face; ++i)
3020 {
3021 const unsigned int ind = index_array_nodal[0][i];
3022 const unsigned int i_ = reorientate(0, i);
3023
3024 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3025 proc.value(temp1[i_][v], vector_ptrs[v][ind]);
3026
3027 if constexpr (integrate == false)
3028 for (unsigned int v = n_filled_lanes; v < n_lanes;
3029 ++v)
3030 temp1[i_][v] = 0.0;
3031 }
3032 else
3033 {
3034 if (integrate == false && n_filled_lanes < n_lanes)
3035 for (unsigned int i = 0; i < dofs_per_face; ++i)
3036 temp1[i] = Number();
3037
3038 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3039 for (unsigned int i = 0; i < dofs_per_face; ++i)
3040 proc.value(temp1[reorientate(v, i)][v],
3041 vector_ptrs[v][index_array_nodal[v][i]]);
3042 }
3043 }
3044 }
3045 }
3046 else
3047 {
3048 // We should not end up here, this should be caught by
3049 // FEFaceEvaluationImplGatherEvaluateSelector::supports()
3051 }
3052 temp1 += 3 * dofs_per_face;
3053 }
3054 }
3055
3056
3057
3058 template <int dim, typename Number2, typename VectorizedArrayType>
3060 {
3061 using Number = typename VectorizedArrayType::value_type;
3062
3063 template <int fe_degree, int n_q_points_1d>
3064 static bool
3065 run(const unsigned int n_components,
3066 const EvaluationFlags::EvaluationFlags evaluation_flag,
3067 const Number2 *src_ptr,
3068 const std::vector<ArrayView<const Number2>> *sm_ptr,
3070 {
3071 Assert(fe_degree > -1, ExcInternalError());
3075
3076 const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
3077
3078 VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
3079 VectorizedArrayType *scratch_data =
3080 temp + 3 * n_components * dofs_per_face;
3081
3083
3084 if (fe_eval.get_dof_access_index() ==
3086 fe_eval.is_interior_face() == false)
3087 fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
3088 p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
3089 else
3090 fe_face_evaluation_process_and_io<1>(
3091 p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
3092
3093 const unsigned int subface_index = fe_eval.get_subface_index();
3094
3095 if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
3097 dim,
3098 fe_degree,
3099 n_q_points_1d,
3100 VectorizedArrayType>::
3101 evaluate_in_face(n_components,
3102 evaluation_flag,
3103 fe_eval.get_shape_info().data.front(),
3104 temp,
3105 fe_eval.begin_values(),
3106 fe_eval.begin_gradients(),
3107 fe_eval.begin_hessians(),
3108 scratch_data,
3109 subface_index);
3110 else
3112 dim,
3113 fe_degree,
3114 n_q_points_1d,
3115 VectorizedArrayType>::
3116 evaluate_in_face(n_components,
3117 evaluation_flag,
3118 fe_eval.get_shape_info().data.front(),
3119 temp,
3120 fe_eval.begin_values(),
3121 fe_eval.begin_gradients(),
3122 fe_eval.begin_hessians(),
3123 scratch_data,
3124 subface_index);
3125
3126 // re-orientation for cases not possible with above algorithm
3127 if (subface_index < GeometryInfo<dim>::max_children_per_cell)
3128 {
3129 if (fe_eval.get_dof_access_index() ==
3131 fe_eval.is_interior_face() == false)
3132 {
3133 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3134 {
3135 // the loop breaks once an invalid_unsigned_int is hit for
3136 // all cases except the exterior faces in the ECL loop (where
3137 // some faces might be at the boundaries but others not)
3138 if (fe_eval.get_cell_ids()[v] ==
3140 continue;
3141
3142 if (fe_eval.get_face_orientation(v) != 0)
3144 dim,
3145 n_components,
3146 v,
3147 evaluation_flag,
3149 fe_eval.get_face_orientation(v), 0),
3150 false,
3151 Utilities::pow(n_q_points_1d, dim - 1),
3152 &temp[0][0],
3153 fe_eval.begin_values(),
3154 fe_eval.begin_gradients(),
3155 fe_eval.begin_hessians());
3156 }
3157 }
3158 else if (fe_eval.get_face_orientation() != 0)
3160 dim,
3161 n_components,
3162 evaluation_flag,
3164 fe_eval.get_face_orientation(), 0),
3165 false,
3166 Utilities::pow(n_q_points_1d, dim - 1),
3167 temp,
3168 fe_eval.begin_values(),
3169 fe_eval.begin_gradients(),
3170 fe_eval.begin_hessians());
3171 }
3172
3173 return false;
3174 }
3175
3176 template <typename Number3>
3177 static bool
3180 const Number2 *vector_ptr,
3182 {
3183 const unsigned int fe_degree = shape_info.data.front().fe_degree;
3184 if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries ||
3185 (evaluation_flag & EvaluationFlags::gradients &&
3186 (fe_degree < 2 ||
3187 shape_info.data.front().element_type !=
3189 (evaluation_flag & EvaluationFlags::hessians) ||
3190 vector_ptr == nullptr ||
3191 shape_info.data.front().element_type >
3193 storage <
3195 return false;
3196 else
3197 return true;
3198 }
3199
3200 private:
3201 template <int fe_degree>
3203 {
3204 static const bool do_integrate = false;
3205 static const int dim_ = dim;
3206 static const int fe_degree_ = fe_degree;
3207 using VectorizedArrayType_ = VectorizedArrayType;
3209 using Number2_ = const Number2;
3210
3211 template <typename T0, typename T1, typename T2>
3212 void
3214 T0 &temp_2,
3215 const T1 src_ptr_1,
3216 const T1 src_ptr_2,
3217 const T2 &grad_weight)
3218 {
3219 do_vectorized_read(src_ptr_1, temp_1);
3220 do_vectorized_read(src_ptr_2, temp_2);
3221 temp_2 = grad_weight * (temp_1 - temp_2);
3222 }
3223
3224 template <typename T1, typename T2>
3225 void
3226 value_vectorized(T1 &temp, const T2 src_ptr)
3227 {
3228 do_vectorized_read(src_ptr, temp);
3229 }
3230
3231 template <typename T0, typename T1, typename T2, typename T3>
3232 void
3234 T0 &temp_2,
3235 const T1 src_ptr_1,
3236 const T1 src_ptr_2,
3237 const T2 &grad_weight,
3238 const T3 &indices_1,
3239 const T3 &indices_2)
3240 {
3241 do_vectorized_gather(src_ptr_1, indices_1, temp_1);
3242 do_vectorized_gather(src_ptr_2, indices_2, temp_2);
3243 temp_2 = grad_weight * (temp_1 - temp_2);
3244 }
3245
3246 template <typename T0, typename T1, typename T2>
3247 void
3248 value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
3249 {
3250 do_vectorized_gather(src_ptr, indices, temp);
3251 }
3252
3253 template <typename T0, typename T1, typename T2>
3254 void
3255 hermite_grad(T0 &temp_1,
3256 T0 &temp_2,
3257 const T1 &src_ptr_1,
3258 const T1 &src_ptr_2,
3259 const T2 &grad_weight)
3260 {
3261 // case 3a)
3262 temp_1 = src_ptr_1;
3263 temp_2 = grad_weight * (temp_1 - src_ptr_2);
3264 }
3265
3266 template <typename T1, typename T2>
3267 void
3268 value(T1 &temp, const T2 &src_ptr)
3269 {
3270 // case 3b)
3271 temp = src_ptr;
3272 }
3273 };
3274 };
3275
3276
3277
3278 template <int dim, typename Number2, typename VectorizedArrayType>
3280 {
3281 using Number = typename VectorizedArrayType::value_type;
3282
3283 template <int fe_degree, int n_q_points_1d>
3284 static bool
3285 run(const unsigned int n_components,
3286 const EvaluationFlags::EvaluationFlags integration_flag,
3287 Number2 *dst_ptr,
3288 const std::vector<ArrayView<const Number2>> *sm_ptr,
3290 {
3291 Assert(fe_degree > -1, ExcInternalError());
3295
3296 const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
3297
3298 VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
3299 VectorizedArrayType *scratch_data =
3300 temp + 3 * n_components * dofs_per_face;
3301
3302 const unsigned int subface_index = fe_eval.get_subface_index();
3303
3304 // re-orientation for cases not possible with the io function below
3305 if (subface_index < GeometryInfo<dim>::max_children_per_cell)
3306 {
3307 if (fe_eval.get_dof_access_index() ==
3309 fe_eval.is_interior_face() == false)
3310 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3311 {
3312 // the loop breaks once an invalid_unsigned_int is hit for
3313 // all cases except the exterior faces in the ECL loop (where
3314 // some faces might be at the boundaries but others not)
3315 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
3316 continue;
3317
3318 if (fe_eval.get_face_orientation(v) != 0)
3320 dim,
3321 n_components,
3322 v,
3323 integration_flag,
3325 fe_eval.get_face_orientation(v), 0),
3326 true,
3327 Utilities::pow(n_q_points_1d, dim - 1),
3328 &temp[0][0],
3329 fe_eval.begin_values(),
3330 fe_eval.begin_gradients(),
3331 fe_eval.begin_hessians());
3332 }
3333 else if (fe_eval.get_face_orientation() != 0)
3335 dim,
3336 n_components,
3337 integration_flag,
3339 fe_eval.get_face_orientation(), 0),
3340 true,
3341 Utilities::pow(n_q_points_1d, dim - 1),
3342 temp,
3343 fe_eval.begin_values(),
3344 fe_eval.begin_gradients(),
3345 fe_eval.begin_hessians());
3346 }
3347
3348 if (fe_degree > -1 && fe_eval.get_subface_index() >=
3349 GeometryInfo<dim - 1>::max_children_per_cell)
3351 dim,
3352 fe_degree,
3353 n_q_points_1d,
3354 VectorizedArrayType>::
3355 integrate_in_face(n_components,
3356 integration_flag,
3357 fe_eval.get_shape_info().data.front(),
3358 temp,
3359 fe_eval.begin_values(),
3360 fe_eval.begin_gradients(),
3361 fe_eval.begin_hessians(),
3362 scratch_data,
3363 subface_index);
3364 else
3366 dim,
3367 fe_degree,
3368 n_q_points_1d,
3369 VectorizedArrayType>::
3370 integrate_in_face(n_components,
3371 integration_flag,
3372 fe_eval.get_shape_info().data.front(),
3373 temp,
3374 fe_eval.begin_values(),
3375 fe_eval.begin_gradients(),
3376 fe_eval.begin_hessians(),
3377 scratch_data,
3378 subface_index);
3379
3381
3382 if (fe_eval.get_dof_access_index() ==
3384 fe_eval.is_interior_face() == false)
3385 fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
3386 p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
3387 else
3388 fe_face_evaluation_process_and_io<1>(
3389 p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
3390
3391 return false;
3392 }
3393
3394 private:
3395 template <int fe_degree>
3397 {
3398 static const bool do_integrate = true;
3399 static const int dim_ = dim;
3400 static const int fe_degree_ = fe_degree;
3401 using VectorizedArrayType_ = VectorizedArrayType;
3403 using Number2_ = Number2;
3404
3405 template <typename T0, typename T1, typename T2, typename T3, typename T4>
3406 void
3407 hermite_grad_vectorized(const T0 &temp_1,
3408 const T1 &temp_2,
3409 T2 dst_ptr_1,
3410 T3 dst_ptr_2,
3411 const T4 &grad_weight)
3412 {
3413 // case 1a)
3414 const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
3415 const VectorizedArrayType grad = grad_weight * temp_2;
3416 do_vectorized_add(val, dst_ptr_1);
3417 do_vectorized_add(grad, dst_ptr_2);
3418 }
3419
3420 template <typename T0, typename T1>
3421 void
3422 value_vectorized(const T0 &temp, T1 dst_ptr)
3423 {
3424 // case 1b)
3425 do_vectorized_add(temp, dst_ptr);
3426 }
3427
3428 template <typename T0, typename T1, typename T2, typename T3>
3429 void
3431 const T0 &temp_2,
3432 T1 dst_ptr_1,
3433 T1 dst_ptr_2,
3434 const T2 &grad_weight,
3435 const T3 &indices_1,
3436 const T3 &indices_2)
3437 {
3438 // case 2a)
3439 const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
3440 const VectorizedArrayType grad = grad_weight * temp_2;
3441 do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
3442 do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
3443 }
3444
3445 template <typename T0, typename T1, typename T2>
3446 void
3447 value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
3448 {
3449 // case 2b)
3450 do_vectorized_scatter_add(temp, indices, dst_ptr);
3451 }
3452
3453 template <typename T0, typename T1, typename T2>
3454 void
3455 hermite_grad(const T0 &temp_1,
3456 const T0 &temp_2,
3457 T1 &dst_ptr_1,
3458 T1 &dst_ptr_2,
3459 const T2 &grad_weight)
3460 {
3461 // case 3a)
3462 const Number val = temp_1 - grad_weight * temp_2;
3463 const Number grad = grad_weight * temp_2;
3464 dst_ptr_1 += val;
3465 dst_ptr_2 += grad;
3466 }
3467
3468 template <typename T0, typename T1>
3469 void
3470 value(const T0 &temp, T1 &dst_ptr)
3471 {
3472 // case 3b)
3473 dst_ptr += temp;
3474 }
3475 };
3476 };
3477} // end of namespace internal
3478
3479
3481
3482#endif
iterator begin() const
Definition array_view.h:707
std::uint8_t get_face_no(const unsigned int v=0) const
internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex get_dof_access_index() const
ScalarNumber shape_info_number_type
const ShapeInfoType & get_shape_info() const
const std::array< unsigned int, n_lanes > & get_cell_ids() const
const Number * begin_gradients() const
unsigned int get_subface_index() const
bool is_interior_face() const
ArrayView< Number > get_scratch_data() const
const Number * begin_values() const
std::uint8_t get_face_orientation(const unsigned int v=0) const
const Number * begin_hessians() const
void gather(const Number *base_ptr, const unsigned int *offsets)
void load(const OtherNumber *ptr)
#define DEAL_II_ALWAYS_INLINE
Definition config.h:161
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:208
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
unsigned int cell_index
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
std::vector< index_type > data
Definition mpi.cc:746
EvaluationFlags
The EvaluationFlags enum.
constexpr T fixed_power(const T t)
Definition utilities.h:943
constexpr T pow(const T base, const int iexp)
Definition utilities.h:967
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 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)
std::enable_if_t<(variant==evaluate_general), void > apply_matrix_vector_product(const Number2 *matrix, const Number *in, Number *out)
void do_vectorized_scatter_add(const VectorizedArrayType src, const unsigned int *indices, Number2 *dst_ptr)
void do_vectorized_gather(const Number2 *src_ptr, const unsigned int *indices, VectorizedArrayType &dst)
void do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
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 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)
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, FEEvaluationData< dim, Number, true > &fe_eval)
static void evaluate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, FEEvaluationData< dim, Number, true > &fe_eval, Number *temp, Number *scratch_data)
static void project_to_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool use_vectorization, Number *temp, Number *scratch_data)
static bool evaluate_tensor_none(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static void adjust_quadrature_for_face_orientation(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, FEEvaluationData< dim, Number, true > &fe_eval, const bool use_vectorization, Number *temp)
static bool evaluate_tensor(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs_actual, 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)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, FEEvaluationData< dim, Number, true > &fe_eval)
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 integrate_tensor(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs_actual, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
static bool integrate_tensor_none(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
static void collect_from_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool use_vectorization, const Number *temp, Number *scratch_data, const bool sum_into_values)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
static void integrate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, FEEvaluationData< dim, Number, true > &fe_eval, Number *temp, Number *scratch_data)
static void adjust_quadrature_for_face_orientation(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, FEEvaluationData< dim, Number, true > &fe_eval, const bool use_vectorization, Number *temp)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
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)
typename FEEvaluationData< dim, Number, true >::shape_info_number_type Number2
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:485
::Table< 2, unsigned int > face_orientations_quad
Definition shape_info.h:624