Reference documentation for deal.II version GIT acb749f4f5 2024-02-21 16:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
evaluation_kernels_face.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2017 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_matrix_free_evaluation_kernels_face_h
18#define dealii_matrix_free_evaluation_kernels_face_h
19
20#include <deal.II/base/config.h>
21
26
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,
72 data.fe_degree + 1,
73 data.n_q_points_1d);
74 else if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
75 return Eval(data.shape_values,
76 data.shape_gradients,
77 data.shape_hessians,
78 data.fe_degree + 1,
79 data.n_q_points_1d);
80 else
81 {
82 const unsigned int index =
83 direction == 0 ? subface_index % 2 : subface_index / 2;
84 return Eval(data.values_within_subface[index],
87 data.fe_degree + 1,
88 data.n_q_points_1d);
89 }
90 }
91
92 static void
94 const unsigned int n_components,
95 const EvaluationFlags::EvaluationFlags evaluation_flag,
97 Number *values_dofs,
98 Number *values_quad,
99 Number *gradients_quad,
100 Number *hessians_quad,
101 Number *scratch_data,
102 const unsigned int subface_index)
103 {
104 Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
105 Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
106
107 const std::size_t n_dofs = fe_degree > -1 ?
108 Utilities::pow(fe_degree + 1, dim - 1) :
109 Utilities::pow(data.fe_degree + 1, dim - 1);
110 const std::size_t n_q_points =
111 fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
112 Utilities::pow(data.n_q_points_1d, dim - 1);
113
114 // keep a copy of the original pointer for the case of the Hessians
115 Number *values_dofs_ptr = values_dofs;
116
117 if ((evaluation_flag & EvaluationFlags::values) != 0u &&
118 ((evaluation_flag & EvaluationFlags::gradients) == 0u))
119 for (unsigned int c = 0; c < n_components; ++c)
120 {
121 switch (dim)
122 {
123 case 3:
124 eval0.template values<0, true, false>(values_dofs,
125 values_quad);
126 eval1.template values<1, true, false>(values_quad,
127 values_quad);
128 break;
129 case 2:
130 eval0.template values<0, true, false>(values_dofs,
131 values_quad);
132 break;
133 case 1:
134 values_quad[0] = values_dofs[0];
135 break;
136 default:
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 interpolate_generic<do_evaluate, add_into_output>(
879 n_components,
880 input,
881 output,
882 flags,
883 face_no,
884 shape_info.data.front().fe_degree + 1,
885 shape_info.data.front().shape_data_on_face,
887 3 * shape_info.dofs_per_component_on_face);
888 }
889
893 template <bool do_evaluate, bool add_into_output>
894 static void
896 const unsigned int n_components,
899 const Number *input,
900 Number *output,
901 const unsigned int face_no)
902 {
903 Assert(static_cast<unsigned int>(fe_degree + 1) ==
904 shape_info.data.front().n_q_points_1d ||
905 fe_degree == -1,
907
908 interpolate_generic<do_evaluate, add_into_output>(
909 n_components,
910 input,
911 output,
912 flags,
913 face_no,
914 shape_info.data.front().quadrature.size(),
915 shape_info.data.front().quadrature_data_on_face,
916 shape_info.n_q_points,
917 shape_info.n_q_points_face);
918 }
919
920 private:
921 template <bool do_evaluate, bool add_into_output, int face_direction = 0>
922 static void
923 interpolate_generic(const unsigned int n_components,
924 const Number *input,
925 Number *output,
927 const unsigned int face_no,
928 const unsigned int n_points_1d,
929 const std::array<AlignedVector<Number2>, 2> &shape_data,
930 const unsigned int dofs_per_component_on_cell,
931 const unsigned int dofs_per_component_on_face)
932 {
933 if (face_direction == face_no / 2)
934 {
935 constexpr int stride_ = Utilities::pow(fe_degree + 1, face_direction);
936
937 const int n_rows = fe_degree != -1 ? fe_degree + 1 : n_points_1d;
938 const int stride = Utilities::pow(n_rows, face_direction);
939 const std::array<int, 2> n_blocks{
940 {(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}};
941 std::array<int, 2> steps;
942 if constexpr (face_direction == 0)
943 steps = {{n_rows, 0}};
944 else if constexpr (face_direction == 1 && dim == 2)
945 steps = {{1, 0}};
946 else if constexpr (face_direction == 1)
947 // in 3d, the coordinate system is zx, not xz -> switch indices
948 steps = {{n_rows * n_rows, -n_rows * n_rows * n_rows + 1}};
949 else if constexpr (face_direction == 2)
950 steps = {{1, 0}};
951
952 for (unsigned int c = 0; c < n_components; ++c)
953 {
954 if (flag & EvaluationFlags::hessians)
955 interpolate_to_face<fe_degree + 1,
956 stride_,
957 do_evaluate,
958 add_into_output,
959 2>(shape_data[face_no % 2].begin(),
960 n_blocks,
961 steps,
962 input,
963 output,
964 n_rows,
965 stride);
966 else if (flag & EvaluationFlags::gradients)
967 interpolate_to_face<fe_degree + 1,
968 stride_,
969 do_evaluate,
970 add_into_output,
971 1>(shape_data[face_no % 2].begin(),
972 n_blocks,
973 steps,
974 input,
975 output,
976 n_rows,
977 stride);
978 else
979 interpolate_to_face<fe_degree + 1,
980 stride_,
981 do_evaluate,
982 add_into_output,
983 0>(shape_data[face_no % 2].begin(),
984 n_blocks,
985 steps,
986 input,
987 output,
988 n_rows,
989 stride);
990 if (do_evaluate)
991 {
992 input += dofs_per_component_on_cell;
993 output += dofs_per_component_on_face;
994 }
995 else
996 {
997 output += dofs_per_component_on_cell;
998 input += dofs_per_component_on_face;
999 }
1000 }
1001 }
1002 else if (face_direction < dim)
1003 {
1004 interpolate_generic<do_evaluate,
1005 add_into_output,
1006 std::min(face_direction + 1, dim - 1)>(
1007 n_components,
1008 input,
1009 output,
1010 flag,
1011 face_no,
1012 n_points_1d,
1013 shape_data,
1014 dofs_per_component_on_cell,
1015 dofs_per_component_on_face);
1016 }
1017 }
1018
1019 template <bool do_evaluate,
1020 bool add_into_output,
1021 int face_direction = 0,
1022 int max_derivative = 0>
1023 static void
1025 const unsigned int n_components,
1026 const Number *input,
1027 Number *output,
1029 const unsigned int face_no,
1031 {
1032 if (dim == 1)
1033 {
1034 // This should never happen since the FE_RaviartThomasNodal is not
1035 // defined for dim = 1. It prevents compiler warnings of infinite
1036 // recursion.
1038 return;
1039 }
1040
1041 bool increase_max_der = false;
1042 if ((flag & EvaluationFlags::hessians && max_derivative < 2) ||
1043 (flag & EvaluationFlags::gradients && max_derivative < 1))
1044 increase_max_der = true;
1045
1046 if (face_direction == face_no / 2 && !increase_max_der)
1047 {
1048 constexpr int stride1 = Utilities::pow(fe_degree + 1, face_direction);
1049 constexpr int stride0 = Utilities::pow(fe_degree, face_direction);
1050 constexpr int stride2 = fe_degree * (fe_degree + 1);
1051
1052 const int degree =
1053 fe_degree != -1 ? fe_degree : shape_info.data[0].fe_degree;
1054 const int n_rows_n = degree + 1;
1055 const int n_rows_t = degree;
1056
1057 std::array<int, 3> strides{{1, 1, 1}};
1058 if (face_direction > 0)
1059 {
1060 strides[0] =
1061 n_rows_n * Utilities::pow(n_rows_t, face_direction - 1);
1062 strides[1] = n_rows_t * (face_direction == 3 ? n_rows_n : 1);
1063 strides[2] = Utilities::pow(n_rows_t, face_direction);
1064 }
1065 const ::ndarray<int, 3, 3> dofs_per_direction{
1066 {{{n_rows_n, n_rows_t, n_rows_t}},
1067 {{n_rows_t, n_rows_n, n_rows_t}},
1068 {{n_rows_t, n_rows_t, n_rows_n}}}};
1069
1070 std::array<int, 2> steps, n_blocks;
1071
1072 if constexpr (face_direction == 0)
1073 steps = {{degree + (face_direction == 0), 0}};
1074 else if constexpr (face_direction == 1 && dim == 2)
1075 steps = {{1, 0}};
1076 else if constexpr (face_direction == 1)
1077 // in 3d, the coordinate system is zx, not xz -> switch indices
1078 steps = {
1079 {n_rows_n * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
1080 else if constexpr (face_direction == 2)
1081 steps = {{1, 0}};
1082
1083 n_blocks[0] = dofs_per_direction[0][(face_direction + 1) % dim];
1084 n_blocks[1] =
1085 dim > 2 ? dofs_per_direction[0][(face_direction + 2) % dim] : 1;
1086
1088 (fe_degree != -1 ? (fe_degree + (face_direction == 0)) : 0),
1089 ((face_direction < 2) ? stride1 : stride2),
1090 do_evaluate,
1091 add_into_output,
1092 max_derivative>(shape_info.data[face_direction != 0]
1093 .shape_data_on_face[face_no % 2]
1094 .begin(),
1095 n_blocks,
1096 steps,
1097 input,
1098 output,
1099 degree + (face_direction == 0),
1100 strides[0]);
1101
1102 if (do_evaluate)
1103 {
1104 input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1105 output += 3 * n_blocks[0] * n_blocks[1];
1106 }
1107 else
1108 {
1109 output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1110 input += 3 * n_blocks[0] * n_blocks[1];
1111 }
1112
1113 // must only change steps only for face direction 0
1114 if constexpr (face_direction == 0)
1115 steps = {{degree, 0}};
1116
1117 n_blocks[0] = dofs_per_direction[1][(face_direction + 1) % dim];
1118 n_blocks[1] =
1119 dim > 2 ? dofs_per_direction[1][(face_direction + 2) % dim] : 1;
1120
1122 (fe_degree != -1 ? (fe_degree + (face_direction == 1)) : 0),
1123 ((face_direction < 2) ? stride0 : stride2),
1124 do_evaluate,
1125 add_into_output,
1126 max_derivative>(shape_info.data[face_direction != 1]
1127 .shape_data_on_face[face_no % 2]
1128 .begin(),
1129 n_blocks,
1130 steps,
1131 input,
1132 output,
1133 degree + (face_direction == 1),
1134 strides[1]);
1135
1136 if constexpr (dim > 2)
1137 {
1138 if (do_evaluate)
1139 {
1140 input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1141 output += 3 * n_blocks[0] * n_blocks[1];
1142 }
1143 else
1144 {
1145 output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1146 input += 3 * n_blocks[0] * n_blocks[1];
1147 }
1148
1149 if constexpr (face_direction == 0)
1150 steps = {{degree, 0}};
1151 else if constexpr (face_direction == 1)
1152 // in 3d, the coordinate system is zx, not xz -> switch indices
1153 steps = {
1154 {n_rows_t * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
1155 else if constexpr (face_direction == 2)
1156 steps = {{1, 0}};
1157
1158 n_blocks[0] = dofs_per_direction[2][(face_direction + 1) % dim];
1159 n_blocks[1] = dofs_per_direction[2][(face_direction + 2) % dim];
1160
1162 (fe_degree != -1 ? (fe_degree + (face_direction == 2)) : 0),
1163 stride0,
1164 do_evaluate,
1165 add_into_output,
1166 max_derivative>(shape_info.data[face_direction != 2]
1167 .shape_data_on_face[face_no % 2]
1168 .begin(),
1169 n_blocks,
1170 steps,
1171 input,
1172 output,
1173 degree + (face_direction == 2),
1174 strides[2]);
1175 }
1176 }
1177 else if (face_direction == face_no / 2)
1178 {
1179 // Only increase max_derivative
1180 interpolate_raviart_thomas<do_evaluate,
1181 add_into_output,
1182 face_direction,
1183 std::min(max_derivative + 1, 2)>(
1184 n_components, input, output, flag, face_no, shape_info);
1185 }
1186 else if (face_direction < dim)
1187 {
1188 if (increase_max_der)
1189 {
1190 interpolate_raviart_thomas<do_evaluate,
1191 add_into_output,
1192 std::min(face_direction + 1, dim - 1),
1193 std::min(max_derivative + 1, 2)>(
1194 n_components, input, output, flag, face_no, shape_info);
1195 }
1196 else
1197 {
1198 interpolate_raviart_thomas<do_evaluate,
1199 add_into_output,
1200 std::min(face_direction + 1, dim - 1),
1201 max_derivative>(
1202 n_components, input, output, flag, face_no, shape_info);
1203 }
1204 }
1205 }
1206 };
1207
1208
1209
1210 // internal helper function for reading data; base version of different types
1211 template <typename VectorizedArrayType, typename Number2>
1212 void
1213 do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
1214 {
1215 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1216 dst[v] = src_ptr[v];
1217 }
1218
1219
1220
1221 // internal helper function for reading data; specialized version where we
1222 // can use a dedicated load function
1223 template <typename Number, std::size_t width>
1224 void
1226 {
1227 dst.load(src_ptr);
1228 }
1229
1230
1231
1232 // internal helper function for reading data; base version of different types
1233 template <typename VectorizedArrayType, typename Number2>
1234 void
1235 do_vectorized_gather(const Number2 *src_ptr,
1236 const unsigned int *indices,
1237 VectorizedArrayType &dst)
1238 {
1239 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1240 dst[v] = src_ptr[indices[v]];
1241 }
1242
1243
1244
1245 // internal helper function for reading data; specialized version where we
1246 // can use a dedicated gather function
1247 template <typename Number, std::size_t width>
1248 void
1249 do_vectorized_gather(const Number *src_ptr,
1250 const unsigned int *indices,
1252 {
1253 dst.gather(src_ptr, indices);
1254 }
1255
1256
1257
1258 // internal helper function for reading data; base version of different types
1259 template <typename VectorizedArrayType, typename Number2>
1260 void
1261 do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
1262 {
1263 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1264 dst_ptr[v] += src[v];
1265 }
1266
1267
1268
1269 // internal helper function for reading data; specialized version where we
1270 // can use a dedicated load function
1271 template <typename Number, std::size_t width>
1272 void
1274 {
1276 tmp.load(dst_ptr);
1277 (tmp + src).store(dst_ptr);
1278 }
1279
1280
1281
1282 // internal helper function for reading data; base version of different types
1283 template <typename VectorizedArrayType, typename Number2>
1284 void
1285 do_vectorized_scatter_add(const VectorizedArrayType src,
1286 const unsigned int *indices,
1287 Number2 *dst_ptr)
1288 {
1289 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1290 dst_ptr[indices[v]] += src[v];
1291 }
1292
1293
1294
1295 // internal helper function for reading data; specialized version where we
1296 // can use a dedicated gather function
1297 template <typename Number, std::size_t width>
1298 void
1300 const unsigned int *indices,
1301 Number *dst_ptr)
1302 {
1303#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
1304 for (unsigned int v = 0; v < width; ++v)
1305 dst_ptr[indices[v]] += src[v];
1306#else
1308 tmp.gather(dst_ptr, indices);
1309 (tmp + src).scatter(indices, dst_ptr);
1310#endif
1311 }
1312
1313
1314
1315 template <typename Number>
1316 void
1317 adjust_for_face_orientation(const unsigned int dim,
1318 const unsigned int n_components,
1320 const unsigned int *orientation,
1321 const bool integrate,
1322 const std::size_t n_q_points,
1323 Number *tmp_values,
1324 Number *values_quad,
1325 Number *gradients_quad,
1326 Number *hessians_quad)
1327 {
1328 for (unsigned int c = 0; c < n_components; ++c)
1329 {
1330 if (flag & EvaluationFlags::values)
1331 {
1332 if (integrate)
1333 for (unsigned int q = 0; q < n_q_points; ++q)
1334 tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
1335 else
1336 for (unsigned int q = 0; q < n_q_points; ++q)
1337 tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
1338 for (unsigned int q = 0; q < n_q_points; ++q)
1339 values_quad[c * n_q_points + q] = tmp_values[q];
1340 }
1341 if (flag & EvaluationFlags::gradients)
1342 for (unsigned int d = 0; d < dim; ++d)
1343 {
1344 if (integrate)
1345 for (unsigned int q = 0; q < n_q_points; ++q)
1346 tmp_values[q] =
1347 gradients_quad[(c * n_q_points + orientation[q]) * dim + d];
1348 else
1349 for (unsigned int q = 0; q < n_q_points; ++q)
1350 tmp_values[orientation[q]] =
1351 gradients_quad[(c * n_q_points + q) * dim + d];
1352 for (unsigned int q = 0; q < n_q_points; ++q)
1353 gradients_quad[(c * n_q_points + q) * dim + d] = tmp_values[q];
1354 }
1355 if (flag & EvaluationFlags::hessians)
1356 {
1357 const unsigned int hdim = (dim * (dim + 1)) / 2;
1358 for (unsigned int d = 0; d < hdim; ++d)
1359 {
1360 if (integrate)
1361 for (unsigned int q = 0; q < n_q_points; ++q)
1362 tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
1363 orientation[q]];
1364 else
1365 for (unsigned int q = 0; q < n_q_points; ++q)
1366 tmp_values[orientation[q]] =
1367 hessians_quad[(c * hdim + d) * n_q_points + q];
1368 for (unsigned int q = 0; q < n_q_points; ++q)
1369 hessians_quad[(c * hdim + d) * n_q_points + q] =
1370 tmp_values[q];
1371 }
1372 }
1373 }
1374 }
1375
1376
1377
1378 template <typename Number, typename VectorizedArrayType>
1379 void
1381 const unsigned int dim,
1382 const unsigned int n_components,
1383 const unsigned int v,
1385 const unsigned int *orientation,
1386 const bool integrate,
1387 const std::size_t n_q_points,
1388 Number *tmp_values,
1389 VectorizedArrayType *values_quad,
1390 VectorizedArrayType *gradients_quad = nullptr,
1391 VectorizedArrayType *hessians_quad = nullptr)
1392 {
1393 for (unsigned int c = 0; c < n_components; ++c)
1394 {
1395 if (flag & EvaluationFlags::values)
1396 {
1397 if (integrate)
1398 for (unsigned int q = 0; q < n_q_points; ++q)
1399 tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v];
1400 else
1401 for (unsigned int q = 0; q < n_q_points; ++q)
1402 tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v];
1403 for (unsigned int q = 0; q < n_q_points; ++q)
1404 values_quad[c * n_q_points + q][v] = tmp_values[q];
1405 }
1406 if (flag & EvaluationFlags::gradients)
1407 for (unsigned int d = 0; d < dim; ++d)
1408 {
1409 Assert(gradients_quad != nullptr, ExcInternalError());
1410 if (integrate)
1411 for (unsigned int q = 0; q < n_q_points; ++q)
1412 tmp_values[q] =
1413 gradients_quad[(c * n_q_points + orientation[q]) * dim + d]
1414 [v];
1415 else
1416 for (unsigned int q = 0; q < n_q_points; ++q)
1417 tmp_values[orientation[q]] =
1418 gradients_quad[(c * n_q_points + q) * dim + d][v];
1419 for (unsigned int q = 0; q < n_q_points; ++q)
1420 gradients_quad[(c * n_q_points + q) * dim + d][v] =
1421 tmp_values[q];
1422 }
1423 if (flag & EvaluationFlags::hessians)
1424 {
1425 Assert(hessians_quad != nullptr, ExcInternalError());
1426 const unsigned int hdim = (dim * (dim + 1)) / 2;
1427 for (unsigned int d = 0; d < hdim; ++d)
1428 {
1429 if (integrate)
1430 for (unsigned int q = 0; q < n_q_points; ++q)
1431 tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
1432 orientation[q]][v];
1433 else
1434 for (unsigned int q = 0; q < n_q_points; ++q)
1435 tmp_values[orientation[q]] =
1436 hessians_quad[(c * hdim + d) * n_q_points + q][v];
1437 for (unsigned int q = 0; q < n_q_points; ++q)
1438 hessians_quad[(c * hdim + d) * n_q_points + q][v] =
1439 tmp_values[q];
1440 }
1441 }
1442 }
1443 }
1444
1445
1446
1447 template <int dim, typename Number>
1449 {
1450 static bool
1451 evaluate_tensor_none(const unsigned int n_components,
1452 const EvaluationFlags::EvaluationFlags evaluation_flag,
1453 const Number *values_dofs,
1455 {
1456 const auto &shape_info = fe_eval.get_shape_info();
1457 const auto &shape_data = shape_info.data.front();
1458 using Number2 =
1460
1461 Assert((fe_eval.get_dof_access_index() ==
1463 fe_eval.is_interior_face() == false) == false,
1465
1466 const unsigned int face_no = fe_eval.get_face_no();
1467 const unsigned int face_orientation = fe_eval.get_face_orientation();
1468 const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
1469 const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
1470
1471 using Eval =
1473
1474 if (evaluation_flag & EvaluationFlags::values)
1475 {
1476 const auto *const shape_values =
1477 &shape_data.shape_values_face(face_no, face_orientation, 0);
1478
1479 auto *values_quad_ptr = fe_eval.begin_values();
1480 auto *values_dofs_actual_ptr = values_dofs;
1481
1482 Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
1483 for (unsigned int c = 0; c < n_components; ++c)
1484 {
1485 eval.template values<0, true, false>(values_dofs_actual_ptr,
1486 values_quad_ptr);
1487
1488 values_quad_ptr += n_q_points;
1489 values_dofs_actual_ptr += n_dofs;
1490 }
1491 }
1492
1493 if (evaluation_flag & EvaluationFlags::gradients)
1494 {
1495 auto *gradients_quad_ptr = fe_eval.begin_gradients();
1496 const auto *values_dofs_actual_ptr = values_dofs;
1497
1498 std::array<const Number2 *, dim> shape_gradients;
1499 for (unsigned int d = 0; d < dim; ++d)
1500 shape_gradients[d] =
1501 &shape_data.shape_gradients_face(face_no, face_orientation, d, 0);
1502
1503 for (unsigned int c = 0; c < n_components; ++c)
1504 {
1505 for (unsigned int d = 0; d < dim; ++d)
1506 {
1507 Eval eval(
1508 nullptr, shape_gradients[d], nullptr, n_dofs, n_q_points);
1509
1510 eval.template gradients<0, true, false, dim>(
1511 values_dofs_actual_ptr, gradients_quad_ptr + d);
1512 }
1513 gradients_quad_ptr += n_q_points * dim;
1514 values_dofs_actual_ptr += n_dofs;
1515 }
1516 }
1517
1518 Assert(!(evaluation_flag & EvaluationFlags::hessians),
1520
1521 return true;
1522 }
1523
1524 template <int fe_degree>
1525#ifndef DEBUG
1527#endif
1528 static void
1529 project_to_face(const unsigned int n_components,
1530 const EvaluationFlags::EvaluationFlags evaluation_flag,
1531 const Number *values_dofs,
1533 const bool use_vectorization,
1534 Number *temp,
1535 Number *scratch_data)
1536 {
1537 const auto &shape_info = fe_eval.get_shape_info();
1538
1539 if (use_vectorization == false)
1540 {
1541 const auto &shape_data = shape_info.data.front();
1542
1543 const unsigned int dofs_per_comp_face =
1544 fe_degree > -1 ?
1545 Utilities::pow(fe_degree + 1, dim - 1) :
1546 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1547 const unsigned int dofs_per_face = n_components * dofs_per_comp_face;
1548
1549 for (unsigned int v = 0; v < Number::size(); ++v)
1550 {
1551 // the loop breaks once an invalid_unsigned_int is hit for
1552 // all cases except the exterior faces in the ECL loop (where
1553 // some faces might be at the boundaries but others not)
1554 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
1555 {
1556 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
1557 temp[i][v] = 0;
1558 continue;
1559 }
1560
1562 template interpolate<true, false>(n_components,
1563 evaluation_flag,
1564 shape_info,
1565 values_dofs,
1566 scratch_data,
1567 fe_eval.get_face_no(v));
1568
1569 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
1570 temp[i][v] = scratch_data[i][v];
1571 }
1572 }
1573 else
1575 template interpolate<true, false>(n_components,
1576 evaluation_flag,
1577 shape_info,
1578 values_dofs,
1579 temp,
1580 fe_eval.get_face_no());
1581 }
1582
1583
1584 template <int fe_degree, int n_q_points_1d>
1585#ifndef DEBUG
1587#endif
1588 static void
1589 evaluate_in_face(const unsigned int n_components,
1590 const EvaluationFlags::EvaluationFlags evaluation_flag,
1592 Number *temp,
1593 Number *scratch_data)
1594 {
1595 const auto &shape_info = fe_eval.get_shape_info();
1596 const auto &shape_data = shape_info.data.front();
1597
1598 const unsigned int subface_index = fe_eval.get_subface_index();
1599 constexpr unsigned int n_q_points_1d_actual =
1600 fe_degree > -1 ? n_q_points_1d : 0;
1601
1602 if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
1603 {
1605 fe_degree,
1606 n_q_points_1d_actual,
1607 Number>::
1608 template evaluate_or_integrate_in_face<false>(
1609 evaluation_flag,
1610 fe_eval.get_shape_info().data,
1611 temp,
1612 fe_eval.begin_values(),
1613 fe_eval.begin_gradients(),
1614 scratch_data,
1615 subface_index,
1616 fe_eval.get_face_no() / 2);
1617 }
1618 else if (fe_degree > -1 &&
1620 shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
1622 dim,
1623 fe_degree,
1624 n_q_points_1d_actual,
1625 Number>::evaluate_in_face(n_components,
1626 evaluation_flag,
1627 shape_data,
1628 temp,
1629 fe_eval.begin_values(),
1630 fe_eval
1631 .begin_gradients(),
1632 fe_eval.begin_hessians(),
1633 scratch_data,
1634 subface_index);
1635 else
1637 dim,
1638 fe_degree,
1639 n_q_points_1d_actual,
1640 Number>::evaluate_in_face(n_components,
1641 evaluation_flag,
1642 shape_data,
1643 temp,
1644 fe_eval.begin_values(),
1645 fe_eval
1646 .begin_gradients(),
1647 fe_eval.begin_hessians(),
1648 scratch_data,
1649 subface_index);
1650 }
1651
1652#ifndef DEBUG
1654#endif
1655 static void
1657 const unsigned int n_components,
1658 const EvaluationFlags::EvaluationFlags evaluation_flag,
1660 const bool use_vectorization,
1661 Number *temp)
1662 {
1663 const auto &shape_info = fe_eval.get_shape_info();
1664
1665 if (use_vectorization == false)
1666 {
1667 for (unsigned int v = 0; v < Number::size(); ++v)
1668 {
1669 // the loop breaks once an invalid_unsigned_int is hit for
1670 // all cases except the exterior faces in the ECL loop (where
1671 // some faces might be at the boundaries but others not)
1672 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
1673 continue;
1674
1675 if (fe_eval.get_face_orientation(v) != 0)
1677 dim,
1678 n_components,
1679 v,
1680 evaluation_flag,
1681 &shape_info.face_orientations_quad(
1682 fe_eval.get_face_orientation(v), 0),
1683 false,
1684 shape_info.n_q_points_face,
1685 &temp[0][0],
1686 fe_eval.begin_values(),
1687 fe_eval.begin_gradients(),
1688 fe_eval.begin_hessians());
1689 }
1690 }
1691 else if (fe_eval.get_face_orientation() != 0)
1693 dim,
1694 n_components,
1695 evaluation_flag,
1696 &shape_info.face_orientations_quad(fe_eval.get_face_orientation(), 0),
1697 false,
1698 shape_info.n_q_points_face,
1699 temp,
1700 fe_eval.begin_values(),
1701 fe_eval.begin_gradients(),
1702 fe_eval.begin_hessians());
1703 }
1704
1705
1706
1707 template <int fe_degree, int n_q_points_1d>
1708 static bool
1709 evaluate_tensor(const unsigned int n_components,
1710 const EvaluationFlags::EvaluationFlags evaluation_flag,
1711 const Number *values_dofs,
1713 {
1714 const auto &shape_info = fe_eval.get_shape_info();
1715 const auto &shape_data = shape_info.data.front();
1716
1717 const unsigned int dofs_per_comp_face =
1718 fe_degree > -1 ?
1719 Utilities::pow(fe_degree + 1, dim - 1) :
1720 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1721
1722 // Note: we always keep storage of values, 1st and 2nd derivatives in an
1723 // array, so reserve space for all three here
1724 Number *temp = fe_eval.get_scratch_data().begin();
1725 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
1726
1727 bool use_vectorization = true;
1728 if (fe_eval.get_dof_access_index() ==
1730 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
1731 for (unsigned int v = 0; v < Number::size(); ++v)
1732 if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
1733 fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
1734 use_vectorization = false;
1735
1736 project_to_face<fe_degree>(n_components,
1737 evaluation_flag,
1738 values_dofs,
1739 fe_eval,
1740 use_vectorization,
1741 temp,
1742 scratch_data);
1743
1744 evaluate_in_face<fe_degree, n_q_points_1d>(
1745 n_components, evaluation_flag, fe_eval, temp, scratch_data);
1746
1747 if (dim == 3)
1749 n_components, evaluation_flag, fe_eval, use_vectorization, temp);
1750
1751 return false;
1752 }
1753
1754 template <int fe_degree, int n_q_points_1d>
1755 static bool
1756 run(const unsigned int n_components,
1757 const EvaluationFlags::EvaluationFlags evaluation_flag,
1758 const Number *values_dofs,
1760 {
1761 const auto &shape_info = fe_eval.get_shape_info();
1762
1763 if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
1764 return evaluate_tensor_none(n_components,
1765 evaluation_flag,
1766 values_dofs,
1767 fe_eval);
1768 else
1769 return evaluate_tensor<fe_degree, n_q_points_1d>(n_components,
1770 evaluation_flag,
1771 values_dofs,
1772 fe_eval);
1773 }
1774 };
1775
1776
1777
1778 template <int dim, typename Number>
1780 {
1781 template <int fe_degree>
1782 static bool
1783 run(const unsigned int n_components,
1784 const EvaluationFlags::EvaluationFlags evaluation_flag,
1785 const Number *values_dofs,
1787 {
1788 const auto &shape_info = fe_eval.get_shape_info();
1789 const auto &shape_data = shape_info.data.front();
1790
1791 const unsigned int dofs_per_comp_face =
1792 fe_degree > -1 ?
1793 Utilities::pow(fe_degree + 1, dim - 1) :
1794 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1795
1796 // Note: we always keep storage of values, 1st and 2nd derivatives in an
1797 // array, so reserve space for all three here
1798 Number *temp = fe_eval.get_scratch_data().begin();
1799 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
1800
1801 bool use_vectorization = true;
1802 if (fe_eval.get_dof_access_index() ==
1804 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
1805 for (unsigned int v = 0; v < Number::size(); ++v)
1806 if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
1807 fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
1808 use_vectorization = false;
1809
1811 template project_to_face<fe_degree>(n_components,
1812 evaluation_flag,
1813 values_dofs,
1814 fe_eval,
1815 use_vectorization,
1816 temp,
1817 scratch_data);
1818
1819 return false;
1820 }
1821 };
1822
1823
1824
1825 template <int dim, typename Number>
1827 {
1828 template <int fe_degree, int n_q_points_1d>
1829 static bool
1830 run(const unsigned int n_components,
1831 const EvaluationFlags::EvaluationFlags evaluation_flag,
1833 {
1834 const auto &shape_info = fe_eval.get_shape_info();
1835 const auto &shape_data = shape_info.data.front();
1836
1837 const unsigned int dofs_per_comp_face =
1838 fe_degree > -1 ?
1839 Utilities::pow(fe_degree + 1, dim - 1) :
1840 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1841
1842 // Note: we always keep storage of values, 1st and 2nd derivatives in an
1843 // array, so reserve space for all three here
1844 Number *temp = fe_eval.get_scratch_data().begin();
1845 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
1846
1848 template evaluate_in_face<fe_degree, n_q_points_1d>(
1849 n_components, evaluation_flag, fe_eval, temp, scratch_data);
1850
1851 return false;
1852 }
1853 };
1854
1855
1856
1857 template <int dim, typename Number>
1859 {
1860 static bool
1862 const unsigned int n_components,
1863 const EvaluationFlags::EvaluationFlags integration_flag,
1864 Number *values_dofs,
1866 const bool sum_into_values)
1867 {
1868 const auto &shape_info = fe_eval.get_shape_info();
1869 const auto &shape_data = shape_info.data.front();
1870 using Number2 =
1872
1873 Assert((fe_eval.get_dof_access_index() ==
1875 fe_eval.is_interior_face() == false) == false,
1877
1878 const unsigned int face_no = fe_eval.get_face_no();
1879 const unsigned int face_orientation = fe_eval.get_face_orientation();
1880 const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
1881 const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
1882
1883 using Eval =
1885
1886 if (integration_flag & EvaluationFlags::values)
1887 {
1888 const auto *const shape_values =
1889 &shape_data.shape_values_face(face_no, face_orientation, 0);
1890
1891 auto *values_quad_ptr = fe_eval.begin_values();
1892 auto *values_dofs_actual_ptr = values_dofs;
1893
1894 Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
1895 for (unsigned int c = 0; c < n_components; ++c)
1896 {
1897 if (sum_into_values)
1898 eval.template values<0, false, true>(values_quad_ptr,
1899 values_dofs_actual_ptr);
1900 else
1901 eval.template values<0, false, false>(values_quad_ptr,
1902 values_dofs_actual_ptr);
1903 values_quad_ptr += n_q_points;
1904 values_dofs_actual_ptr += n_dofs;
1905 }
1906 }
1907
1908 if (integration_flag & EvaluationFlags::gradients)
1909 {
1910 auto *gradients_quad_ptr = fe_eval.begin_gradients();
1911 auto *values_dofs_actual_ptr = values_dofs;
1912
1913 std::array<const Number2 *, dim> shape_gradients;
1914 for (unsigned int d = 0; d < dim; ++d)
1915 shape_gradients[d] =
1916 &shape_data.shape_gradients_face(face_no, face_orientation, d, 0);
1917
1918 for (unsigned int c = 0; c < n_components; ++c)
1919 {
1920 for (unsigned int d = 0; d < dim; ++d)
1921 {
1922 Eval eval(
1923 nullptr, shape_gradients[d], nullptr, n_dofs, n_q_points);
1924
1925 if (!sum_into_values &&
1926 !(integration_flag & EvaluationFlags::values) && d == 0)
1927 eval.template gradients<0, false, false, dim>(
1928 gradients_quad_ptr + d, values_dofs_actual_ptr);
1929 else
1930 eval.template gradients<0, false, true, dim>(
1931 gradients_quad_ptr + d, values_dofs_actual_ptr);
1932 }
1933 gradients_quad_ptr += n_q_points * dim;
1934 values_dofs_actual_ptr += n_dofs;
1935 }
1936 }
1937
1938 Assert(!(integration_flag & EvaluationFlags::hessians),
1940
1941 return true;
1942 }
1943
1944#ifndef DEBUG
1946#endif
1947 static void
1949 const unsigned int n_components,
1950 const EvaluationFlags::EvaluationFlags integration_flag,
1952 const bool use_vectorization,
1953 Number *temp)
1954 {
1955 const auto &shape_info = fe_eval.get_shape_info();
1956
1957 if (use_vectorization == false)
1958 {
1959 for (unsigned int v = 0; v < Number::size(); ++v)
1960 {
1961 // the loop breaks once an invalid_unsigned_int is hit for
1962 // all cases except the exterior faces in the ECL loop (where
1963 // some faces might be at the boundaries but others not)
1964 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
1965 continue;
1966
1967 if (fe_eval.get_face_orientation(v) != 0)
1969 dim,
1970 n_components,
1971 v,
1972 integration_flag,
1974 fe_eval.get_face_orientation(v), 0),
1975 true,
1976 shape_info.n_q_points_face,
1977 &temp[0][0],
1978 fe_eval.begin_values(),
1979 fe_eval.begin_gradients(),
1980 fe_eval.begin_hessians());
1981 }
1982 }
1983 else if (fe_eval.get_face_orientation() != 0)
1985 dim,
1986 n_components,
1987 integration_flag,
1989 fe_eval.get_face_orientation(), 0),
1990 true,
1991 shape_info.n_q_points_face,
1992 temp,
1993 fe_eval.begin_values(),
1994 fe_eval.begin_gradients(),
1995 fe_eval.begin_hessians());
1996 }
1997
1998 template <int fe_degree, int n_q_points_1d>
1999#ifndef DEBUG
2001#endif
2002 static void
2003 integrate_in_face(const unsigned int n_components,
2004 const EvaluationFlags::EvaluationFlags integration_flag,
2006 Number *temp,
2007 Number *scratch_data)
2008 {
2009 const auto &shape_info = fe_eval.get_shape_info();
2010 const auto &shape_data = shape_info.data.front();
2011
2012 const unsigned int n_q_points_1d_actual =
2013 fe_degree > -1 ? n_q_points_1d : 0;
2014 const unsigned int subface_index = fe_eval.get_subface_index();
2015
2016 if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
2017 {
2019 fe_degree,
2020 n_q_points_1d_actual,
2021 Number>::
2022 template evaluate_or_integrate_in_face<true>(
2023 integration_flag,
2024 fe_eval.get_shape_info().data,
2025 temp,
2026 fe_eval.begin_values(),
2027 fe_eval.begin_gradients(),
2028 scratch_data,
2029 subface_index,
2030 fe_eval.get_face_no() / 2);
2031 }
2032 else if (fe_degree > -1 &&
2033 fe_eval.get_subface_index() >=
2035 shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
2037 true,
2038 dim,
2039 fe_degree,
2040 n_q_points_1d_actual,
2041 Number>::integrate_in_face(n_components,
2042 integration_flag,
2043 shape_data,
2044 temp,
2045 fe_eval.begin_values(),
2046 fe_eval.begin_gradients(),
2047 fe_eval.begin_hessians(),
2048 scratch_data,
2049 subface_index);
2050 else
2052 false,
2053 dim,
2054 fe_degree,
2055 n_q_points_1d_actual,
2056 Number>::integrate_in_face(n_components,
2057 integration_flag,
2058 shape_data,
2059 temp,
2060 fe_eval.begin_values(),
2061 fe_eval.begin_gradients(),
2062 fe_eval.begin_hessians(),
2063 scratch_data,
2064 subface_index);
2065 }
2066
2067 template <int fe_degree>
2068#ifndef DEBUG
2070#endif
2071 static void
2072 collect_from_face(const unsigned int n_components,
2073 const EvaluationFlags::EvaluationFlags integration_flag,
2074 Number *values_dofs,
2076 const bool use_vectorization,
2077 const Number *temp,
2078 Number *scratch_data,
2079 const bool sum_into_values)
2080 {
2081 const auto &shape_info = fe_eval.get_shape_info();
2082 const auto &shape_data = shape_info.data.front();
2083
2084 const unsigned int dofs_per_comp_face =
2085 fe_degree > -1 ?
2086 Utilities::pow(fe_degree + 1, dim - 1) :
2087 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2088 const unsigned int dofs_per_face = n_components * dofs_per_comp_face;
2089
2090 if (use_vectorization == false)
2091 {
2092 for (unsigned int v = 0; v < Number::size(); ++v)
2093 {
2094 // the loop breaks once an invalid_unsigned_int is hit for
2095 // all cases except the exterior faces in the ECL loop (where
2096 // some faces might be at the boundaries but others not)
2097 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2098 continue;
2099
2101 template interpolate<false, false>(n_components,
2102 integration_flag,
2103 shape_info,
2104 temp,
2105 scratch_data,
2106 fe_eval.get_face_no(v));
2107
2108 if (sum_into_values)
2109 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
2110 values_dofs[i][v] += scratch_data[i][v];
2111 else
2112 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
2113 values_dofs[i][v] = scratch_data[i][v];
2114 }
2115 }
2116 else
2117 {
2118 if (sum_into_values)
2120 template interpolate<false, true>(n_components,
2121 integration_flag,
2122 shape_info,
2123 temp,
2124 values_dofs,
2125 fe_eval.get_face_no());
2126 else
2128 template interpolate<false, false>(n_components,
2129 integration_flag,
2130 shape_info,
2131 temp,
2132 values_dofs,
2133 fe_eval.get_face_no());
2134 }
2135 }
2136
2137 template <int fe_degree, int n_q_points_1d>
2138 static bool
2139 integrate_tensor(const unsigned int n_components,
2140 const EvaluationFlags::EvaluationFlags integration_flag,
2141 Number *values_dofs,
2143 const bool sum_into_values)
2144 {
2145 const auto &shape_info = fe_eval.get_shape_info();
2146 const auto &shape_data = shape_info.data.front();
2147
2148 const unsigned int dofs_per_comp_face =
2149 fe_degree > -1 ?
2150 Utilities::pow(fe_degree + 1, dim - 1) :
2151 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2152
2153 Number *temp = fe_eval.get_scratch_data().begin();
2154 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
2155
2156 bool use_vectorization = true;
2157
2158 if (fe_eval.get_dof_access_index() ==
2160 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
2161 use_vectorization =
2163 std::all_of(fe_eval.get_cell_ids().begin() + 1,
2164 fe_eval.get_cell_ids().end(),
2165 [&](const auto &v) {
2166 return v == fe_eval.get_cell_ids()[0] ||
2167 v == numbers::invalid_unsigned_int;
2168 });
2169
2170 if (dim == 3)
2172 n_components, integration_flag, fe_eval, use_vectorization, temp);
2173
2174 integrate_in_face<fe_degree, n_q_points_1d>(
2175 n_components, integration_flag, fe_eval, temp, scratch_data);
2176
2177 collect_from_face<fe_degree>(n_components,
2178 integration_flag,
2179 values_dofs,
2180 fe_eval,
2181 use_vectorization,
2182 temp,
2183 scratch_data,
2184 sum_into_values);
2185
2186 return false;
2187 }
2188
2189 template <int fe_degree, int n_q_points_1d>
2190 static bool
2191 run(const unsigned int n_components,
2192 const EvaluationFlags::EvaluationFlags integration_flag,
2193 Number *values_dofs,
2195 const bool sum_into_values)
2196 {
2197 const auto &shape_info = fe_eval.get_shape_info();
2198
2199 if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
2200 return integrate_tensor_none(n_components,
2201 integration_flag,
2202 values_dofs,
2203 fe_eval,
2204 sum_into_values);
2205 else
2206 return integrate_tensor<fe_degree, n_q_points_1d>(n_components,
2207 integration_flag,
2208 values_dofs,
2209 fe_eval,
2210 sum_into_values);
2211 }
2212 };
2213
2214
2215
2216 template <int dim, typename Number>
2218 {
2219 template <int fe_degree>
2220 static bool
2221 run(const unsigned int n_components,
2222 const EvaluationFlags::EvaluationFlags integration_flag,
2223 Number *values_dofs,
2225 const bool sum_into_values)
2226 {
2227 const auto &shape_info = fe_eval.get_shape_info();
2228 const auto &shape_data = shape_info.data.front();
2229
2230 const unsigned int dofs_per_comp_face =
2231 fe_degree > -1 ?
2232 Utilities::pow(fe_degree + 1, dim - 1) :
2233 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2234
2235 Number *temp = fe_eval.get_scratch_data().begin();
2236 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
2237
2238 bool use_vectorization = true;
2239
2240 if (fe_eval.get_dof_access_index() ==
2242 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
2243 use_vectorization =
2245 std::all_of(fe_eval.get_cell_ids().begin() + 1,
2246 fe_eval.get_cell_ids().end(),
2247 [&](const auto &v) {
2248 return v == fe_eval.get_cell_ids()[0] ||
2249 v == numbers::invalid_unsigned_int;
2250 });
2251
2253 template collect_from_face<fe_degree>(n_components,
2254 integration_flag,
2255 values_dofs,
2256 fe_eval,
2257 use_vectorization,
2258 temp,
2259 scratch_data,
2260 sum_into_values);
2261
2262 return false;
2263 }
2264 };
2265
2266
2267
2268 template <int dim, typename Number>
2270 {
2271 template <int fe_degree, int n_q_points_1d>
2272 static bool
2273 run(const unsigned int n_components,
2274 const EvaluationFlags::EvaluationFlags integration_flag,
2275
2277 {
2278 const auto &shape_info = fe_eval.get_shape_info();
2279 const auto &shape_data = shape_info.data.front();
2280
2281 const unsigned int dofs_per_comp_face =
2282 fe_degree > -1 ?
2283 Utilities::pow(fe_degree + 1, dim - 1) :
2284 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2285
2286 Number *temp = fe_eval.get_scratch_data().begin();
2287 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
2288
2290 template integrate_in_face<fe_degree, n_q_points_1d>(
2291 n_components, integration_flag, fe_eval, temp, scratch_data);
2292
2293 return false;
2294 }
2295 };
2296
2297
2298
2299 template <int n_face_orientations,
2300 typename Processor,
2301 typename EvaluationData,
2302 const bool check_face_orientations = false>
2303 void
2305 Processor &proc,
2306 const unsigned int n_components,
2307 const EvaluationFlags::EvaluationFlags evaluation_flag,
2308 typename Processor::Number2_ *global_vector_ptr,
2309 const std::vector<ArrayView<const typename Processor::Number2_>> *sm_ptr,
2310 const EvaluationData &fe_eval,
2311 typename Processor::VectorizedArrayType_ *temp1)
2312 {
2313 constexpr int dim = Processor::dim_;
2314 constexpr int fe_degree = Processor::fe_degree_;
2315 using VectorizedArrayType = typename Processor::VectorizedArrayType_;
2316 constexpr int n_lanes = VectorizedArrayType::size();
2317
2318 using Number = typename Processor::Number_;
2319 using Number2_ = typename Processor::Number2_;
2320
2321 const auto &shape_data = fe_eval.get_shape_info().data.front();
2322 constexpr bool integrate = Processor::do_integrate;
2323 const unsigned int face_no = fe_eval.get_face_no();
2324 const auto &dof_info = fe_eval.get_dof_info();
2325 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
2326 const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index =
2327 fe_eval.get_dof_access_index();
2328 AssertIndexRange(cell,
2329 dof_info.index_storage_variants[dof_access_index].size());
2330 constexpr unsigned int dofs_per_face =
2331 Utilities::pow(fe_degree + 1, dim - 1);
2332 const unsigned int subface_index = fe_eval.get_subface_index();
2333
2334 const unsigned int n_filled_lanes =
2335 dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
2336
2337 bool all_faces_are_same = n_filled_lanes == n_lanes;
2338 if (n_face_orientations == n_lanes)
2339 for (unsigned int v = 1; v < n_lanes; ++v)
2340 if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) ||
2341 fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0))
2342 {
2343 all_faces_are_same = false;
2344 break;
2345 }
2346
2347 // check for re-orientation ...
2348 std::array<const unsigned int *, n_face_orientations> orientation = {};
2349
2350 if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same &&
2351 fe_eval.is_interior_face() == 0)
2352 for (unsigned int v = 0; v < n_lanes; ++v)
2353 {
2354 // the loop breaks once an invalid_unsigned_int is hit for
2355 // all cases except the exterior faces in the ECL loop (where
2356 // some faces might be at the boundaries but others not)
2357 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2358 continue;
2359
2360 if (shape_data.nodal_at_cell_boundaries &&
2361 fe_eval.get_face_orientation(v) != 0)
2362 {
2363 // ... and in case we detect a re-orientation, go to the other
2364 // version of this function that actually allows for this
2365 if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
2366 check_face_orientations == false)
2367 {
2368 fe_face_evaluation_process_and_io<n_face_orientations,
2369 Processor,
2370 EvaluationData,
2371 true>(proc,
2372 n_components,
2373 evaluation_flag,
2374 global_vector_ptr,
2375 sm_ptr,
2376 fe_eval,
2377 temp1);
2378 return;
2379 }
2380 orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
2381 fe_eval.get_face_orientation(v), 0);
2382 }
2383 }
2384 else if (dim == 3 && fe_eval.get_face_orientation() != 0)
2385 {
2386 // go to the other version of this function
2387 if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
2388 check_face_orientations == false)
2389 {
2390 fe_face_evaluation_process_and_io<n_face_orientations,
2391 Processor,
2392 EvaluationData,
2393 true>(proc,
2394 n_components,
2395 evaluation_flag,
2396 global_vector_ptr,
2397 sm_ptr,
2398 fe_eval,
2399 temp1);
2400 return;
2401 }
2402 for (unsigned int v = 0; v < n_face_orientations; ++v)
2403 orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
2404 fe_eval.get_face_orientation(), 0);
2405 }
2406
2407 // we know that the gradient weights for the Hermite case on the
2408 // right (side==1) are the negative from the value at the left
2409 // (side==0), so we only read out one of them.
2410 VectorizedArrayType grad_weight =
2411 shape_data
2412 .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) :
2413 (1 + face_no % 2))];
2414
2415 // face_to_cell_index_hermite
2416 std::array<const unsigned int *, n_face_orientations> index_array_hermite =
2417 {};
2418 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2419 {
2420 if (n_face_orientations == 1)
2421 index_array_hermite[0] =
2422 &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0);
2423 else
2424 {
2425 for (unsigned int v = 0; v < n_lanes; ++v)
2426 {
2427 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2428 continue;
2429
2430 const auto face_no = fe_eval.get_face_no(v);
2431
2432 grad_weight[v] =
2433 shape_data.shape_data_on_face[0][fe_degree +
2434 (integrate ?
2435 (2 - (face_no % 2)) :
2436 (1 + (face_no % 2)))][0];
2437
2438 index_array_hermite[v] =
2439 &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no,
2440 0);
2441 }
2442 }
2443 }
2444
2445 // face_to_cell_index_nodal
2446 std::array<const unsigned int *, n_face_orientations> index_array_nodal =
2447 {};
2448 if (shape_data.nodal_at_cell_boundaries == true)
2449 {
2450 if (n_face_orientations == 1)
2451 index_array_nodal[0] =
2452 &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0);
2453 else
2454 {
2455 for (unsigned int v = 0; v < n_lanes; ++v)
2456 {
2457 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2458 continue;
2459
2460 const auto face_no = fe_eval.get_face_no(v);
2461
2462 index_array_nodal[v] =
2463 &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no,
2464 0);
2465 }
2466 }
2467 }
2468
2469
2470 const auto reorientate = [&](const unsigned int v, const unsigned int i) {
2471 return (!check_face_orientations || orientation[v] == nullptr) ?
2472 i :
2473 orientation[v][i];
2474 };
2475
2476 const unsigned int cell_index =
2478 fe_eval.get_cell_ids()[0] :
2479 cell * n_lanes;
2480 const unsigned int *dof_indices =
2481 &dof_info.dof_indices_contiguous[dof_access_index][cell_index];
2482
2483 for (unsigned int comp = 0; comp < n_components; ++comp)
2484 {
2485 const std::size_t index_offset =
2486 dof_info.component_dof_indices_offset
2487 [fe_eval.get_active_fe_index()]
2488 [fe_eval.get_first_selected_component()] +
2489 comp * Utilities::pow(fe_degree + 1, dim);
2490
2491 // case 1: contiguous and interleaved indices
2492 if (n_face_orientations == 1 &&
2493 dof_info.index_storage_variants[dof_access_index][cell] ==
2495 interleaved_contiguous)
2496 {
2498 dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2499 n_lanes);
2500 Number2_ *vector_ptr =
2501 global_vector_ptr + dof_indices[0] + index_offset * n_lanes;
2502
2503 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2504 {
2505 for (unsigned int i = 0; i < dofs_per_face; ++i)
2506 {
2507 Assert(n_face_orientations == 1, ExcNotImplemented());
2508
2509 const unsigned int ind1 = index_array_hermite[0][2 * i];
2510 const unsigned int ind2 = index_array_hermite[0][2 * i + 1];
2511 const unsigned int i_ = reorientate(0, i);
2512 proc.hermite_grad_vectorized(temp1[i_],
2513 temp1[i_ + dofs_per_face],
2514 vector_ptr + ind1 * n_lanes,
2515 vector_ptr + ind2 * n_lanes,
2516 grad_weight);
2517 }
2518 }
2519 else
2520 {
2521 for (unsigned int i = 0; i < dofs_per_face; ++i)
2522 {
2523 Assert(n_face_orientations == 1, ExcNotImplemented());
2524
2525 const unsigned int i_ = reorientate(0, i);
2526 const unsigned int ind = index_array_nodal[0][i];
2527 proc.value_vectorized(temp1[i_],
2528 vector_ptr + ind * n_lanes);
2529 }
2530 }
2531 }
2532
2533 // case 2: contiguous and interleaved indices with fixed stride
2534 else if (n_face_orientations == 1 &&
2535 dof_info.index_storage_variants[dof_access_index][cell] ==
2537 interleaved_contiguous_strided)
2538 {
2540 dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2541 n_lanes);
2542 Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes;
2543 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2544 {
2545 for (unsigned int i = 0; i < dofs_per_face; ++i)
2546 {
2547 Assert(n_face_orientations == 1, ExcNotImplemented());
2548
2549 const unsigned int i_ = reorientate(0, i);
2550 const unsigned int ind1 =
2551 index_array_hermite[0][2 * i] * n_lanes;
2552 const unsigned int ind2 =
2553 index_array_hermite[0][2 * i + 1] * n_lanes;
2554 proc.hermite_grad_vectorized_indexed(
2555 temp1[i_],
2556 temp1[i_ + dofs_per_face],
2557 vector_ptr + ind1,
2558 vector_ptr + ind2,
2559 grad_weight,
2560 dof_indices,
2561 dof_indices);
2562 }
2563 }
2564 else
2565 {
2566 for (unsigned int i = 0; i < dofs_per_face; ++i)
2567 {
2568 Assert(n_face_orientations == 1, ExcNotImplemented());
2569
2570 const unsigned int i_ = reorientate(0, i);
2571 const unsigned int ind = index_array_nodal[0][i] * n_lanes;
2572 proc.value_vectorized_indexed(temp1[i_],
2573 vector_ptr + ind,
2574 dof_indices);
2575 }
2576 }
2577 }
2578
2579 // case 3: contiguous and interleaved indices with mixed stride
2580 else if (n_face_orientations == 1 &&
2581 dof_info.index_storage_variants[dof_access_index][cell] ==
2583 interleaved_contiguous_mixed_strides)
2584 {
2585 const unsigned int *strides =
2586 &dof_info.dof_indices_interleave_strides[dof_access_index]
2587 [cell * n_lanes];
2588 unsigned int indices[n_lanes];
2589 for (unsigned int v = 0; v < n_lanes; ++v)
2590 indices[v] = dof_indices[v] + index_offset * strides[v];
2591 const unsigned int n_filled_lanes =
2592 dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
2593
2594 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2595 {
2596 if (n_filled_lanes == n_lanes)
2597 for (unsigned int i = 0; i < dofs_per_face; ++i)
2598 {
2599 Assert(n_face_orientations == 1, ExcNotImplemented());
2600
2601 const unsigned int i_ = reorientate(0, i);
2602 unsigned int ind1[n_lanes];
2604 for (unsigned int v = 0; v < n_lanes; ++v)
2605 ind1[v] = indices[v] +
2606 index_array_hermite[0][2 * i] * strides[v];
2607 unsigned int ind2[n_lanes];
2609 for (unsigned int v = 0; v < n_lanes; ++v)
2610 ind2[v] =
2611 indices[v] +
2612 // TODO
2613 index_array_hermite[0][2 * i + 1] * strides[v];
2614 proc.hermite_grad_vectorized_indexed(
2615 temp1[i_],
2616 temp1[i_ + dofs_per_face],
2617 global_vector_ptr,
2618 global_vector_ptr,
2619 grad_weight,
2620 ind1,
2621 ind2);
2622 }
2623 else
2624 {
2625 if (integrate == false)
2626 for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
2627 temp1[i] = VectorizedArrayType();
2628
2629 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2630 for (unsigned int i = 0; i < dofs_per_face; ++i)
2631 {
2632 const unsigned int i_ =
2633 reorientate(n_face_orientations == 1 ? 0 : v, i);
2634 proc.hermite_grad(
2635 temp1[i_][v],
2636 temp1[i_ + dofs_per_face][v],
2637 global_vector_ptr
2638 [indices[v] +
2639 index_array_hermite
2640 [n_face_orientations == 1 ? 0 : v][2 * i] *
2641 strides[v]],
2642 global_vector_ptr
2643 [indices[v] +
2644 index_array_hermite[n_face_orientations == 1 ?
2645 0 :
2646 v][2 * i + 1] *
2647 strides[v]],
2648 grad_weight[n_face_orientations == 1 ? 0 : v]);
2649 }
2650 }
2651 }
2652 else
2653 {
2654 if (n_filled_lanes == n_lanes)
2655 for (unsigned int i = 0; i < dofs_per_face; ++i)
2656 {
2657 Assert(n_face_orientations == 1, ExcInternalError());
2658 unsigned int ind[n_lanes];
2660 for (unsigned int v = 0; v < n_lanes; ++v)
2661 ind[v] =
2662 indices[v] + index_array_nodal[0][i] * strides[v];
2663 const unsigned int i_ = reorientate(0, i);
2664 proc.value_vectorized_indexed(temp1[i_],
2665 global_vector_ptr,
2666 ind);
2667 }
2668 else
2669 {
2670 if (integrate == false)
2671 for (unsigned int i = 0; i < dofs_per_face; ++i)
2672 temp1[i] = VectorizedArrayType();
2673
2674 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2675 for (unsigned int i = 0; i < dofs_per_face; ++i)
2676 proc.value(
2677 temp1[reorientate(n_face_orientations == 1 ? 0 : v,
2678 i)][v],
2679 global_vector_ptr
2680 [indices[v] +
2681 index_array_nodal[n_face_orientations == 1 ? 0 : v]
2682 [i] *
2683 strides[v]]);
2684 }
2685 }
2686 }
2687
2688 // case 4: contiguous indices without interleaving
2689 else if (n_face_orientations > 1 ||
2690 dof_info.index_storage_variants[dof_access_index][cell] ==
2692 contiguous)
2693 {
2694 Number2_ *vector_ptr = global_vector_ptr + index_offset;
2695
2696 const bool vectorization_possible =
2697 all_faces_are_same && (sm_ptr == nullptr);
2698
2699 std::array<Number2_ *, n_lanes> vector_ptrs;
2700 std::array<unsigned int, n_lanes> reordered_indices;
2701
2702 if (vectorization_possible == false)
2703 {
2704 vector_ptrs = {};
2705 if (n_face_orientations == 1)
2706 {
2707 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2708 if (sm_ptr == nullptr)
2709 {
2710 vector_ptrs[v] = vector_ptr + dof_indices[v];
2711 }
2712 else
2713 {
2714 const auto &temp =
2715 dof_info
2716 .dof_indices_contiguous_sm[dof_access_index]
2717 [cell * n_lanes + v];
2718 vector_ptrs[v] = const_cast<Number2_ *>(
2719 sm_ptr->operator[](temp.first).data() +
2720 temp.second + index_offset);
2721 }
2722 }
2723 else if (n_face_orientations == n_lanes)
2724 {
2725 const auto &cells = fe_eval.get_cell_ids();
2726 for (unsigned int v = 0; v < n_lanes; ++v)
2727 if (cells[v] != numbers::invalid_unsigned_int)
2728 {
2729 if (sm_ptr == nullptr)
2730 {
2731 vector_ptrs[v] =
2732 vector_ptr +
2733 dof_info
2734 .dof_indices_contiguous[dof_access_index]
2735 [cells[v]];
2736 }
2737 else
2738 {
2739 const auto &temp =
2740 dof_info
2741 .dof_indices_contiguous_sm[dof_access_index]
2742 [cells[v]];
2743 vector_ptrs[v] = const_cast<Number2_ *>(
2744 sm_ptr->operator[](temp.first).data() +
2745 temp.second + index_offset);
2746 }
2747 }
2748 }
2749 else
2750 {
2752 }
2753 }
2754 else if (n_face_orientations == n_lanes)
2755 {
2756 for (unsigned int v = 0; v < n_lanes; ++v)
2757 reordered_indices[v] =
2758 dof_info.dof_indices_contiguous[dof_access_index]
2759 [fe_eval.get_cell_ids()[v]];
2760 dof_indices = reordered_indices.data();
2761 }
2762
2763 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2764 {
2765 if (vectorization_possible)
2766 for (unsigned int i = 0; i < dofs_per_face; ++i)
2767 {
2768 const unsigned int ind1 = index_array_hermite[0][2 * i];
2769 const unsigned int ind2 =
2770 index_array_hermite[0][2 * i + 1];
2771 const unsigned int i_ = reorientate(0, i);
2772
2773 proc.hermite_grad_vectorized_indexed(
2774 temp1[i_],
2775 temp1[i_ + dofs_per_face],
2776 vector_ptr + ind1,
2777 vector_ptr + ind2,
2778 grad_weight,
2779 dof_indices,
2780 dof_indices);
2781 }
2782 else if (n_face_orientations == 1)
2783 for (unsigned int i = 0; i < dofs_per_face; ++i)
2784 {
2785 const unsigned int ind1 = index_array_hermite[0][2 * i];
2786 const unsigned int ind2 =
2787 index_array_hermite[0][2 * i + 1];
2788 const unsigned int i_ = reorientate(0, i);
2789
2790 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2791 proc.hermite_grad(temp1[i_][v],
2792 temp1[i_ + dofs_per_face][v],
2793 vector_ptrs[v][ind1],
2794 vector_ptrs[v][ind2],
2795 grad_weight[v]);
2796
2797 if (integrate == false)
2798 for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
2799 {
2800 temp1[i][v] = 0.0;
2801 temp1[i + dofs_per_face][v] = 0.0;
2802 }
2803 }
2804 else
2805 {
2806 if (integrate == false && n_filled_lanes < n_lanes)
2807 for (unsigned int i = 0; i < dofs_per_face; ++i)
2808 temp1[i] = temp1[i + dofs_per_face] = Number();
2809
2810 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2811 for (unsigned int i = 0; i < dofs_per_face; ++i)
2812 proc.hermite_grad(
2813 temp1[reorientate(v, i)][v],
2814 temp1[reorientate(v, i) + dofs_per_face][v],
2815 vector_ptrs[v][index_array_hermite[v][2 * i]],
2816 vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
2817 grad_weight[v]);
2818 }
2819 }
2820 else
2821 {
2822 if (vectorization_possible)
2823 for (unsigned int i = 0; i < dofs_per_face; ++i)
2824 {
2825 const unsigned int ind = index_array_nodal[0][i];
2826 const unsigned int i_ = reorientate(0, i);
2827
2828 proc.value_vectorized_indexed(temp1[i_],
2829 vector_ptr + ind,
2830 dof_indices);
2831 }
2832 else if (n_face_orientations == 1)
2833 for (unsigned int i = 0; i < dofs_per_face; ++i)
2834 {
2835 const unsigned int ind = index_array_nodal[0][i];
2836 const unsigned int i_ = reorientate(0, i);
2837
2838 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2839 proc.value(temp1[i_][v], vector_ptrs[v][ind]);
2840
2841 if (integrate == false)
2842 for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
2843 temp1[i_][v] = 0.0;
2844 }
2845 else
2846 {
2847 if (integrate == false && n_filled_lanes < n_lanes)
2848 for (unsigned int i = 0; i < dofs_per_face; ++i)
2849 temp1[i] = Number();
2850
2851 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2852 for (unsigned int i = 0; i < dofs_per_face; ++i)
2853 proc.value(temp1[reorientate(v, i)][v],
2854 vector_ptrs[v][index_array_nodal[v][i]]);
2855 }
2856 }
2857 }
2858 else
2859 {
2860 // We should not end up here, this should be caught by
2861 // FEFaceEvaluationImplGatherEvaluateSelector::supports()
2863 }
2864 temp1 += 3 * dofs_per_face;
2865 }
2866 }
2867
2868
2869
2870 template <int dim, typename Number2, typename VectorizedArrayType>
2872 {
2873 using Number = typename VectorizedArrayType::value_type;
2874
2875 template <int fe_degree, int n_q_points_1d>
2876 static bool
2877 run(const unsigned int n_components,
2878 const EvaluationFlags::EvaluationFlags evaluation_flag,
2879 const Number2 *src_ptr,
2880 const std::vector<ArrayView<const Number2>> *sm_ptr,
2882 {
2883 Assert(fe_degree > -1, ExcInternalError());
2887
2888 const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
2889
2890 VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
2891 VectorizedArrayType *scratch_data =
2892 temp + 3 * n_components * dofs_per_face;
2893
2895
2896 if (fe_eval.get_dof_access_index() ==
2898 fe_eval.is_interior_face() == false)
2899 fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
2900 p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
2901 else
2902 fe_face_evaluation_process_and_io<1>(
2903 p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
2904
2905 const unsigned int subface_index = fe_eval.get_subface_index();
2906
2907 if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
2909 dim,
2910 fe_degree,
2911 n_q_points_1d,
2912 VectorizedArrayType>::
2913 evaluate_in_face(n_components,
2914 evaluation_flag,
2915 fe_eval.get_shape_info().data.front(),
2916 temp,
2917 fe_eval.begin_values(),
2918 fe_eval.begin_gradients(),
2919 fe_eval.begin_hessians(),
2920 scratch_data,
2921 subface_index);
2922 else
2924 dim,
2925 fe_degree,
2926 n_q_points_1d,
2927 VectorizedArrayType>::
2928 evaluate_in_face(n_components,
2929 evaluation_flag,
2930 fe_eval.get_shape_info().data.front(),
2931 temp,
2932 fe_eval.begin_values(),
2933 fe_eval.begin_gradients(),
2934 fe_eval.begin_hessians(),
2935 scratch_data,
2936 subface_index);
2937
2938 // re-orientation for cases not possible with above algorithm
2939 if (subface_index < GeometryInfo<dim>::max_children_per_cell)
2940 {
2941 if (fe_eval.get_dof_access_index() ==
2943 fe_eval.is_interior_face() == false)
2944 {
2945 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
2946 {
2947 // the loop breaks once an invalid_unsigned_int is hit for
2948 // all cases except the exterior faces in the ECL loop (where
2949 // some faces might be at the boundaries but others not)
2950 if (fe_eval.get_cell_ids()[v] ==
2952 continue;
2953
2954 if (fe_eval.get_face_orientation(v) != 0)
2956 dim,
2957 n_components,
2958 v,
2959 evaluation_flag,
2961 fe_eval.get_face_orientation(v), 0),
2962 false,
2963 Utilities::pow(n_q_points_1d, dim - 1),
2964 &temp[0][0],
2965 fe_eval.begin_values(),
2966 fe_eval.begin_gradients(),
2967 fe_eval.begin_hessians());
2968 }
2969 }
2970 else if (fe_eval.get_face_orientation() != 0)
2972 dim,
2973 n_components,
2974 evaluation_flag,
2976 fe_eval.get_face_orientation(), 0),
2977 false,
2978 Utilities::pow(n_q_points_1d, dim - 1),
2979 temp,
2980 fe_eval.begin_values(),
2981 fe_eval.begin_gradients(),
2982 fe_eval.begin_hessians());
2983 }
2984
2985 return false;
2986 }
2987
2988 template <typename Number3>
2989 static bool
2992 const Number2 *vector_ptr,
2994 {
2995 const unsigned int fe_degree = shape_info.data.front().fe_degree;
2996 if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries ||
2997 (evaluation_flag & EvaluationFlags::gradients &&
2998 (fe_degree < 2 ||
2999 shape_info.data.front().element_type !=
3001 (evaluation_flag & EvaluationFlags::hessians) ||
3002 vector_ptr == nullptr ||
3003 shape_info.data.front().element_type >
3005 storage <
3007 return false;
3008 else
3009 return true;
3010 }
3011
3012 private:
3013 template <int fe_degree>
3015 {
3016 static const bool do_integrate = false;
3017 static const int dim_ = dim;
3018 static const int fe_degree_ = fe_degree;
3019 using VectorizedArrayType_ = VectorizedArrayType;
3021 using Number2_ = const Number2;
3022
3023 template <typename T0, typename T1, typename T2>
3024 void
3026 T0 &temp_2,
3027 const T1 src_ptr_1,
3028 const T1 src_ptr_2,
3029 const T2 &grad_weight)
3030 {
3031 do_vectorized_read(src_ptr_1, temp_1);
3032 do_vectorized_read(src_ptr_2, temp_2);
3033 temp_2 = grad_weight * (temp_1 - temp_2);
3034 }
3035
3036 template <typename T1, typename T2>
3037 void
3038 value_vectorized(T1 &temp, const T2 src_ptr)
3039 {
3040 do_vectorized_read(src_ptr, temp);
3041 }
3042
3043 template <typename T0, typename T1, typename T2, typename T3>
3044 void
3046 T0 &temp_2,
3047 const T1 src_ptr_1,
3048 const T1 src_ptr_2,
3049 const T2 &grad_weight,
3050 const T3 &indices_1,
3051 const T3 &indices_2)
3052 {
3053 do_vectorized_gather(src_ptr_1, indices_1, temp_1);
3054 do_vectorized_gather(src_ptr_2, indices_2, temp_2);
3055 temp_2 = grad_weight * (temp_1 - temp_2);
3056 }
3057
3058 template <typename T0, typename T1, typename T2>
3059 void
3060 value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
3061 {
3062 do_vectorized_gather(src_ptr, indices, temp);
3063 }
3064
3065 template <typename T0, typename T1, typename T2>
3066 void
3067 hermite_grad(T0 &temp_1,
3068 T0 &temp_2,
3069 const T1 &src_ptr_1,
3070 const T1 &src_ptr_2,
3071 const T2 &grad_weight)
3072 {
3073 // case 3a)
3074 temp_1 = src_ptr_1;
3075 temp_2 = grad_weight * (temp_1 - src_ptr_2);
3076 }
3077
3078 template <typename T1, typename T2>
3079 void
3080 value(T1 &temp, const T2 &src_ptr)
3081 {
3082 // case 3b)
3083 temp = src_ptr;
3084 }
3085 };
3086 };
3087
3088
3089
3090 template <int dim, typename Number2, typename VectorizedArrayType>
3092 {
3093 using Number = typename VectorizedArrayType::value_type;
3094
3095 template <int fe_degree, int n_q_points_1d>
3096 static bool
3097 run(const unsigned int n_components,
3098 const EvaluationFlags::EvaluationFlags integration_flag,
3099 Number2 *dst_ptr,
3100 const std::vector<ArrayView<const Number2>> *sm_ptr,
3102 {
3103 Assert(fe_degree > -1, ExcInternalError());
3107
3108 const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
3109
3110 VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
3111 VectorizedArrayType *scratch_data =
3112 temp + 3 * n_components * dofs_per_face;
3113
3114 const unsigned int subface_index = fe_eval.get_subface_index();
3115
3116 // re-orientation for cases not possible with the io function below
3117 if (subface_index < GeometryInfo<dim>::max_children_per_cell)
3118 {
3119 if (fe_eval.get_dof_access_index() ==
3121 fe_eval.is_interior_face() == false)
3122 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3123 {
3124 // the loop breaks once an invalid_unsigned_int is hit for
3125 // all cases except the exterior faces in the ECL loop (where
3126 // some faces might be at the boundaries but others not)
3127 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
3128 continue;
3129
3130 if (fe_eval.get_face_orientation(v) != 0)
3132 dim,
3133 n_components,
3134 v,
3135 integration_flag,
3137 fe_eval.get_face_orientation(v), 0),
3138 true,
3139 Utilities::pow(n_q_points_1d, dim - 1),
3140 &temp[0][0],
3141 fe_eval.begin_values(),
3142 fe_eval.begin_gradients(),
3143 fe_eval.begin_hessians());
3144 }
3145 else if (fe_eval.get_face_orientation() != 0)
3147 dim,
3148 n_components,
3149 integration_flag,
3151 fe_eval.get_face_orientation(), 0),
3152 true,
3153 Utilities::pow(n_q_points_1d, dim - 1),
3154 temp,
3155 fe_eval.begin_values(),
3156 fe_eval.begin_gradients(),
3157 fe_eval.begin_hessians());
3158 }
3159
3160 if (fe_degree > -1 && fe_eval.get_subface_index() >=
3161 GeometryInfo<dim - 1>::max_children_per_cell)
3163 dim,
3164 fe_degree,
3165 n_q_points_1d,
3166 VectorizedArrayType>::
3167 integrate_in_face(n_components,
3168 integration_flag,
3169 fe_eval.get_shape_info().data.front(),
3170 temp,
3171 fe_eval.begin_values(),
3172 fe_eval.begin_gradients(),
3173 fe_eval.begin_hessians(),
3174 scratch_data,
3175 subface_index);
3176 else
3178 dim,
3179 fe_degree,
3180 n_q_points_1d,
3181 VectorizedArrayType>::
3182 integrate_in_face(n_components,
3183 integration_flag,
3184 fe_eval.get_shape_info().data.front(),
3185 temp,
3186 fe_eval.begin_values(),
3187 fe_eval.begin_gradients(),
3188 fe_eval.begin_hessians(),
3189 scratch_data,
3190 subface_index);
3191
3193
3194 if (fe_eval.get_dof_access_index() ==
3196 fe_eval.is_interior_face() == false)
3197 fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
3198 p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
3199 else
3200 fe_face_evaluation_process_and_io<1>(
3201 p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
3202
3203 return false;
3204 }
3205
3206 private:
3207 template <int fe_degree>
3209 {
3210 static const bool do_integrate = true;
3211 static const int dim_ = dim;
3212 static const int fe_degree_ = fe_degree;
3213 using VectorizedArrayType_ = VectorizedArrayType;
3215 using Number2_ = Number2;
3216
3217 template <typename T0, typename T1, typename T2, typename T3, typename T4>
3218 void
3219 hermite_grad_vectorized(const T0 &temp_1,
3220 const T1 &temp_2,
3221 T2 dst_ptr_1,
3222 T3 dst_ptr_2,
3223 const T4 &grad_weight)
3224 {
3225 // case 1a)
3226 const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
3227 const VectorizedArrayType grad = grad_weight * temp_2;
3228 do_vectorized_add(val, dst_ptr_1);
3229 do_vectorized_add(grad, dst_ptr_2);
3230 }
3231
3232 template <typename T0, typename T1>
3233 void
3234 value_vectorized(const T0 &temp, T1 dst_ptr)
3235 {
3236 // case 1b)
3237 do_vectorized_add(temp, dst_ptr);
3238 }
3239
3240 template <typename T0, typename T1, typename T2, typename T3>
3241 void
3243 const T0 &temp_2,
3244 T1 dst_ptr_1,
3245 T1 dst_ptr_2,
3246 const T2 &grad_weight,
3247 const T3 &indices_1,
3248 const T3 &indices_2)
3249 {
3250 // case 2a)
3251 const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
3252 const VectorizedArrayType grad = grad_weight * temp_2;
3253 do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
3254 do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
3255 }
3256
3257 template <typename T0, typename T1, typename T2>
3258 void
3259 value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
3260 {
3261 // case 2b)
3262 do_vectorized_scatter_add(temp, indices, dst_ptr);
3263 }
3264
3265 template <typename T0, typename T1, typename T2>
3266 void
3267 hermite_grad(const T0 &temp_1,
3268 const T0 &temp_2,
3269 T1 &dst_ptr_1,
3270 T1 &dst_ptr_2,
3271 const T2 &grad_weight)
3272 {
3273 // case 3a)
3274 const Number val = temp_1 - grad_weight * temp_2;
3275 const Number grad = grad_weight * temp_2;
3276 dst_ptr_1 += val;
3277 dst_ptr_2 += grad;
3278 }
3279
3280 template <typename T0, typename T1>
3281 void
3282 value(const T0 &temp, T1 &dst_ptr)
3283 {
3284 // case 3b)
3285 dst_ptr += temp;
3286 }
3287 };
3288 };
3289} // end of namespace internal
3290
3291
3293
3294#endif
iterator begin() const
Definition array_view.h:703
std::uint8_t get_face_no(const unsigned int v=0) const
internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex get_dof_access_index() const
const 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:110
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:144
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
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)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
EvaluationFlags
The EvaluationFlags enum.
constexpr T fixed_power(const T t)
Definition utilities.h:940
constexpr T pow(const T base, const int iexp)
Definition utilities.h:964
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)
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)
static const unsigned int invalid_unsigned_int
Definition types.h:221
::VectorizedArray< Number, width > min(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, 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_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 bool integrate_tensor(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:486
::Table< 2, unsigned int > face_orientations_quad
Definition shape_info.h:625
std::array< AlignedVector< Number >, 2 > hessians_within_subface
Definition shape_info.h:327
std::array< AlignedVector< Number >, 2 > values_within_subface
Definition shape_info.h:315
std::array< AlignedVector< Number >, 2 > gradients_within_subface
Definition shape_info.h:321