deal.II version GIT relicensing-6834-g5b78e6bcdf 2026-10-01 11:20:01+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// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2017 - 2025 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13
14#ifndef dealii_matrix_free_evaluation_kernels_face_h
15#define dealii_matrix_free_evaluation_kernels_face_h
16
17#include <deal.II/base/config.h>
18
19#include <deal.II/base/geometry_info.h>
20#include <deal.II/base/ndarray.h>
21#include <deal.II/base/utilities.h>
22#include <deal.II/base/vectorization.h>
23
24#include <deal.II/matrix_free/dof_info.h>
25#include <deal.II/matrix_free/evaluation_flags.h>
26#include <deal.II/matrix_free/evaluation_kernels_common.h>
27#include <deal.II/matrix_free/fe_evaluation_data.h>
28#include <deal.II/matrix_free/shape_info.h>
29#include <deal.II/matrix_free/tensor_product_kernels.h>
30
31
32DEAL_II_NAMESPACE_OPEN
33
34
35namespace internal
36{
37 template <bool symmetric_evaluate,
38 int dim,
39 int fe_degree,
40 int n_q_points_1d,
41 typename Number>
42 struct FEFaceEvaluationImpl
43 {
44 // We enable a transformation to collocation for derivatives if it gives
45 // correct results (first two conditions), if it is the most efficient
46 // choice in terms of operation counts (third condition) and if we were
47 // able to initialize the fields in shape_info.templates.h from the
48 // polynomials (fourth condition).
49 using Number2 =
50 typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
51
52 using Eval = EvaluatorTensorProduct<symmetric_evaluate ? evaluate_evenodd :
53 evaluate_general,
54 dim - 1,
55 fe_degree + 1,
56 n_q_points_1d,
57 Number,
58 Number2>;
59
60 static Eval
61 create_evaluator_tensor_product(
62 const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
63 const unsigned int subface_index,
64 const unsigned int direction)
65 {
66 if (symmetric_evaluate)
67 return Eval(data.shape_values_eo,
68 data.shape_gradients_eo,
69 data.shape_hessians_eo,
70 data.fe_degree + 1,
71 data.n_q_points_1d);
72 else if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
73 return Eval(data.shape_values,
74 data.shape_gradients,
75 data.shape_hessians,
76 data.fe_degree + 1,
77 data.n_q_points_1d);
78 else
79 {
80 const unsigned int index =
81 direction == 0 ? subface_index % 2 : subface_index / 2;
82 return Eval(data.values_within_subface[index],
83 data.gradients_within_subface[index],
84 data.hessians_within_subface[index],
85 data.fe_degree + 1,
86 data.n_q_points_1d);
87 }
88 }
89
90 static void
91 evaluate_in_face(
92 const unsigned int n_components,
93 const EvaluationFlags::EvaluationFlags evaluation_flag,
94 const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
95 Number *values_dofs,
96 Number *values_quad,
97 Number *gradients_quad,
98 Number *hessians_quad,
99 Number *scratch_data,
100 const unsigned int subface_index)
101 {
102 Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
103 Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
104
105 const std::size_t n_dofs = fe_degree > -1 ?
106 Utilities::pow(fe_degree + 1, dim - 1) :
107 Utilities::pow(data.fe_degree + 1, dim - 1);
108 const std::size_t n_q_points =
109 fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
110 Utilities::pow(data.n_q_points_1d, dim - 1);
111
112 // keep a copy of the original pointer for the case of the Hessians
113 Number *values_dofs_ptr = values_dofs;
114
115 if ((evaluation_flag & EvaluationFlags::values) != 0u &&
116 ((evaluation_flag & EvaluationFlags::gradients) == 0u))
117 for (unsigned int c = 0; c < n_components; ++c)
118 {
119 switch (dim)
120 {
121 case 3:
122 eval0.template values<0, true, false>(values_dofs,
123 values_quad);
124 eval1.template values<1, true, false>(values_quad,
125 values_quad);
126 break;
127 case 2:
128 eval0.template values<0, true, false>(values_dofs,
129 values_quad);
130 break;
131 case 1:
132 values_quad[0] = values_dofs[0];
133 break;
134 default:
135 DEAL_II_NOT_IMPLEMENTED();
136 }
137 // Note: we always keep storage of values, 1st and 2nd derivatives
138 // in an array
139 values_dofs += 3 * n_dofs;
140 values_quad += n_q_points;
141 }
142 else if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
143 for (unsigned int c = 0; c < n_components; ++c)
144 {
145 switch (dim)
146 {
147 case 3:
148 if (symmetric_evaluate &&
149 use_collocation_evaluation(fe_degree, n_q_points_1d))
150 {
151 eval0.template values<0, true, false>(values_dofs,
152 values_quad);
153 eval0.template values<1, true, false>(values_quad,
154 values_quad);
155 EvaluatorTensorProduct<evaluate_evenodd,
156 dim - 1,
157 n_q_points_1d,
158 n_q_points_1d,
159 Number,
160 Number2>
161 eval_grad({}, data.shape_gradients_collocation_eo, {});
162 eval_grad.template gradients<0, true, false, 3>(
163 values_quad, gradients_quad);
164 eval_grad.template gradients<1, true, false, 3>(
165 values_quad, gradients_quad + 1);
166 }
167 else
168 {
169 // grad x
170 eval0.template gradients<0, true, false>(values_dofs,
171 scratch_data);
172 eval1.template values<1, true, false, 3>(scratch_data,
173 gradients_quad);
174
175 // grad y
176 eval0.template values<0, true, false>(values_dofs,
177 scratch_data);
178 eval1.template gradients<1, true, false, 3>(
179 scratch_data, gradients_quad + 1);
180
181 if ((evaluation_flag & EvaluationFlags::values) != 0u)
182 eval1.template values<1, true, false>(scratch_data,
183 values_quad);
184 }
185 // grad z
186 eval0.template values<0, true, false>(values_dofs + n_dofs,
187 scratch_data);
188 eval1.template values<1, true, false, 3>(scratch_data,
189 gradients_quad + 2);
190
191 break;
192 case 2:
193 eval0.template values<0, true, false, 2>(values_dofs + n_dofs,
194 gradients_quad + 1);
195 eval0.template gradients<0, true, false, 2>(values_dofs,
196 gradients_quad);
197 if ((evaluation_flag & EvaluationFlags::values) != 0u)
198 eval0.template values<0, true, false>(values_dofs,
199 values_quad);
200 break;
201 case 1:
202 values_quad[0] = values_dofs[0];
203 gradients_quad[0] = values_dofs[1];
204 break;
205 default:
206 AssertThrow(false, ExcNotImplemented());
207 }
208 values_dofs += 3 * n_dofs;
209 values_quad += n_q_points;
210 gradients_quad += dim * n_q_points;
211 }
212
213 if ((evaluation_flag & EvaluationFlags::hessians) != 0u)
214 {
215 values_dofs = values_dofs_ptr;
216 for (unsigned int c = 0; c < n_components; ++c)
217 {
218 switch (dim)
219 {
220 case 3:
221 // grad xx
222 eval0.template hessians<0, true, false>(values_dofs,
223 scratch_data);
224 eval1.template values<1, true, false>(scratch_data,
225 hessians_quad);
226
227 // grad yy
228 eval0.template values<0, true, false>(values_dofs,
229 scratch_data);
230 eval1.template hessians<1, true, false>(scratch_data,
231 hessians_quad +
232 n_q_points);
233
234 // grad zz
235 eval0.template values<0, true, false>(values_dofs +
236 2 * n_dofs,
237 scratch_data);
238 eval1.template values<1, true, false>(scratch_data,
239 hessians_quad +
240 2 * n_q_points);
241
242 // grad xy
243 eval0.template gradients<0, true, false>(values_dofs,
244 scratch_data);
245 eval1.template gradients<1, true, false>(scratch_data,
246 hessians_quad +
247 3 * n_q_points);
248
249 // grad xz
250 eval0.template gradients<0, true, false>(values_dofs +
251 n_dofs,
252 scratch_data);
253 eval1.template values<1, true, false>(scratch_data,
254 hessians_quad +
255 4 * n_q_points);
256
257 // grad yz
258 eval0.template values<0, true, false>(values_dofs + n_dofs,
259 scratch_data);
260 eval1.template gradients<1, true, false>(scratch_data,
261 hessians_quad +
262 5 * n_q_points);
263
264 break;
265 case 2:
266 // grad xx
267 eval0.template hessians<0, true, false>(values_dofs,
268 hessians_quad);
269 // grad yy
270 eval0.template values<0, true, false>(
271 values_dofs + 2 * n_dofs, hessians_quad + n_q_points);
272 // grad xy
273 eval0.template gradients<0, true, false>(
274 values_dofs + n_dofs, hessians_quad + 2 * n_q_points);
275 break;
276 case 1:
277 hessians_quad[0] = values_dofs[2];
278 break;
279 default:
280 AssertThrow(false, ExcNotImplemented());
281 }
282 values_dofs += 3 * n_dofs;
283 hessians_quad += dim * (dim + 1) / 2 * n_q_points;
284 }
285 }
286 }
287
288 static void
289 integrate_in_face(
290 const unsigned int n_components,
291 const EvaluationFlags::EvaluationFlags integration_flag,
292 const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
293 Number *values_dofs,
294 Number *values_quad,
295 Number *gradients_quad,
296 Number *hessians_quad,
297 Number *scratch_data,
298 const unsigned int subface_index)
299 {
300 Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
301 Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
302
303 const std::size_t n_dofs =
304 fe_degree > -1 ?
305 Utilities::pow(fe_degree + 1, dim - 1) :
306 (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
307 const std::size_t n_q_points =
308 fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
309 Utilities::pow(data.n_q_points_1d, dim - 1);
310
311 // keep a copy of the original pointer for the case of the Hessians
312 Number *values_dofs_ptr = values_dofs;
313
314 if ((integration_flag & EvaluationFlags::values) != 0u &&
315 (integration_flag & EvaluationFlags::gradients) == 0u)
316 for (unsigned int c = 0; c < n_components; ++c)
317 {
318 switch (dim)
319 {
320 case 3:
321 eval1.template values<1, false, false>(values_quad,
322 values_quad);
323 eval0.template values<0, false, false>(values_quad,
324 values_dofs);
325 break;
326 case 2:
327 eval0.template values<0, false, false>(values_quad,
328 values_dofs);
329 break;
330 case 1:
331 values_dofs[0] = values_quad[0];
332 break;
333 default:
334 DEAL_II_NOT_IMPLEMENTED();
335 }
336 values_dofs += 3 * n_dofs;
337 values_quad += n_q_points;
338 }
339 else if ((integration_flag & EvaluationFlags::gradients) != 0u)
340 for (unsigned int c = 0; c < n_components; ++c)
341 {
342 switch (dim)
343 {
344 case 3:
345 // grad z
346 eval1.template values<1, false, false, 3>(gradients_quad + 2,
347 scratch_data);
348 eval0.template values<0, false, false>(scratch_data,
349 values_dofs + n_dofs);
350 if (symmetric_evaluate &&
351 use_collocation_evaluation(fe_degree, n_q_points_1d))
352 {
353 EvaluatorTensorProduct<evaluate_evenodd,
354 dim - 1,
355 n_q_points_1d,
356 n_q_points_1d,
357 Number,
358 Number2>
359 eval_grad({}, data.shape_gradients_collocation_eo, {});
360 if ((integration_flag & EvaluationFlags::values) != 0u)
361 eval_grad.template gradients<1, false, true, 3>(
362 gradients_quad + 1, values_quad);
363 else
364 eval_grad.template gradients<1, false, false, 3>(
365 gradients_quad + 1, values_quad);
366 eval_grad.template gradients<0, false, true, 3>(
367 gradients_quad, values_quad);
368 eval0.template values<1, false, false>(values_quad,
369 values_quad);
370 eval0.template values<0, false, false>(values_quad,
371 values_dofs);
372 }
373 else
374 {
375 if ((integration_flag & EvaluationFlags::values) != 0u)
376 {
377 eval1.template values<1, false, false>(values_quad,
378 scratch_data);
379 eval1.template gradients<1, false, true, 3>(
380 gradients_quad + 1, scratch_data);
381 }
382 else
383 eval1.template gradients<1, false, false, 3>(
384 gradients_quad + 1, scratch_data);
385
386 // grad y
387 eval0.template values<0, false, false>(scratch_data,
388 values_dofs);
389
390 // grad x
391 eval1.template values<1, false, false, 3>(gradients_quad,
392 scratch_data);
393 eval0.template gradients<0, false, true>(scratch_data,
394 values_dofs);
395 }
396 break;
397 case 2:
398 eval0.template values<0, false, false, 2>(gradients_quad + 1,
399 values_dofs +
400 n_dofs);
401 eval0.template gradients<0, false, false, 2>(gradients_quad,
402 values_dofs);
403 if ((integration_flag & EvaluationFlags::values) != 0u)
404 eval0.template values<0, false, true>(values_quad,
405 values_dofs);
406 break;
407 case 1:
408 values_dofs[0] = values_quad[0];
409 values_dofs[1] = gradients_quad[0];
410 break;
411 default:
412 AssertThrow(false, ExcNotImplemented());
413 }
414 values_dofs += 3 * n_dofs;
415 values_quad += n_q_points;
416 gradients_quad += dim * n_q_points;
417 }
418
419 if ((integration_flag & EvaluationFlags::hessians) != 0u)
420 {
421 values_dofs = values_dofs_ptr;
422 for (unsigned int c = 0; c < n_components; ++c)
423 {
424 switch (dim)
425 {
426 case 3:
427 // grad xx
428 eval1.template values<1, false, false>(hessians_quad,
429 scratch_data);
430 if ((integration_flag & (EvaluationFlags::values |
431 EvaluationFlags::gradients)) != 0u)
432 eval0.template hessians<0, false, true>(scratch_data,
433 values_dofs);
434 else
435 eval0.template hessians<0, false, false>(scratch_data,
436 values_dofs);
437
438 // grad yy
439 eval1.template hessians<1, false, false>(hessians_quad +
440 n_q_points,
441 scratch_data);
442 eval0.template values<0, false, true>(scratch_data,
443 values_dofs);
444
445 // grad zz
446 eval1.template values<1, false, false>(hessians_quad +
447 2 * n_q_points,
448 scratch_data);
449 eval0.template values<0, false, false>(scratch_data,
450 values_dofs +
451 2 * n_dofs);
452
453 // grad xy
454 eval1.template gradients<1, false, false>(hessians_quad +
455 3 * n_q_points,
456 scratch_data);
457 eval0.template gradients<0, false, true>(scratch_data,
458 values_dofs);
459
460 // grad xz
461 eval1.template values<1, false, false>(hessians_quad +
462 4 * n_q_points,
463 scratch_data);
464 if ((integration_flag & EvaluationFlags::gradients) != 0u)
465 eval0.template gradients<0, false, true>(scratch_data,
466 values_dofs +
467 n_dofs);
468 else
469 eval0.template gradients<0, false, false>(scratch_data,
470 values_dofs +
471 n_dofs);
472
473 // grad yz
474 eval1.template gradients<1, false, false>(hessians_quad +
475 5 * n_q_points,
476 scratch_data);
477 eval0.template values<0, false, true>(scratch_data,
478 values_dofs + n_dofs);
479
480 break;
481 case 2:
482 // grad xx
483 if ((integration_flag & (EvaluationFlags::values |
484 EvaluationFlags::gradients)) != 0u)
485 eval0.template hessians<0, false, true>(hessians_quad,
486 values_dofs);
487 else
488 eval0.template hessians<0, false, false>(hessians_quad,
489 values_dofs);
490
491 // grad yy
492 eval0.template values<0, false, false>(
493 hessians_quad + n_q_points, values_dofs + 2 * n_dofs);
494 // grad xy
495 if ((integration_flag & EvaluationFlags::gradients) != 0u)
496 eval0.template gradients<0, false, true>(
497 hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
498 else
499 eval0.template gradients<0, false, false>(
500 hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
501 break;
502 case 1:
503 values_dofs[2] = hessians_quad[0];
504 if ((integration_flag & EvaluationFlags::values) == 0u)
505 values_dofs[0] = 0;
506 if ((integration_flag & EvaluationFlags::gradients) == 0u)
507 values_dofs[1] = 0;
508 break;
509 default:
510 AssertThrow(false, ExcNotImplemented());
511 }
512 values_dofs += 3 * n_dofs;
513 hessians_quad += dim * (dim + 1) / 2 * n_q_points;
514 }
515 }
516 }
517 };
518
519
520
521 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
522 struct FEFaceEvaluationImplRaviartThomas
523 {
524 using Number2 =
525 typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
526
531 template <bool do_integrate>
532 static inline void
533 evaluate_or_integrate_in_face(
534 const EvaluationFlags::EvaluationFlags evaluation_flag,
535 const std::vector<MatrixFreeFunctions::UnivariateShapeData<Number2>>
536 &shape_data,
537 Number *values_dofs_in,
538 Number *values,
539 Number *gradients,
540 Number *scratch_data,
541 const unsigned int subface_index,
542 const unsigned int face_direction)
543 {
544 AssertDimension(shape_data.size(), 2);
545
546 const int degree = fe_degree != -1 ? fe_degree : shape_data[0].fe_degree;
547 const int n_rows_n = degree + 1;
548 const int n_rows_t = degree;
549 const ::ndarray<int, 3, 3> dofs_per_direction{
550 {{{n_rows_n, n_rows_t, n_rows_t}},
551 {{n_rows_t, n_rows_n, n_rows_t}},
552 {{n_rows_t, n_rows_t, n_rows_n}}}};
553
554 (void)scratch_data;
555 (void)subface_index;
556 // TODO: This is currently not implemented, but the test
557 // matrix_vector_rt_face_03 apparently works without it -> check
558 // if (subface_index < GeometryInfo<dim - 1>::max_children_per_cell)
559 // DEAL_II_NOT_IMPLEMENTED();
560
561 using Eval = EvaluatorTensorProduct<evaluate_evenodd,
562 dim - 1,
563 (fe_degree > 0 ? fe_degree : 0),
564 n_q_points_1d,
565 Number,
566 Number2>;
567
568 std::array<int, dim> values_dofs_offsets = {};
569 for (unsigned int comp = 0; comp < dim - 1; ++comp)
570 {
571 if (dim == 2)
572 values_dofs_offsets[comp + 1] =
573 values_dofs_offsets[comp] +
574 3 * dofs_per_direction[comp][(face_direction + 1) % dim];
575 else
576 values_dofs_offsets[comp + 1] =
577 values_dofs_offsets[comp] +
578 3 * dofs_per_direction[comp][(face_direction + 1) % dim] *
579 dofs_per_direction[comp][(face_direction + 2) % dim];
580 }
581
582 // Jacobians on faces are reordered to enable simple access with the
583 // regular evaluators; to get the RT Piola transform right, we need to
584 // pass through the values_dofs array in a permuted right order
585 std::array<unsigned int, dim> components;
586 for (unsigned int comp = 0; comp < dim; ++comp)
587 components[comp] = (face_direction + comp + 1) % dim;
588
589 for (const unsigned int comp : components)
590 {
591 Number *values_dofs = values_dofs_in + values_dofs_offsets[comp];
592
593 std::array<int, 2> n_blocks{
594 {dofs_per_direction[comp][(face_direction + 1) % dim],
595 (dim > 2 ? dofs_per_direction[comp][(face_direction + 2) % dim] :
596 1)}};
597
598 if constexpr (dim == 3)
599 {
600 EvaluatorTensorProduct<evaluate_evenodd,
601 dim - 1,
602 n_q_points_1d,
603 n_q_points_1d,
604 Number,
605 Number2>
606 eval_g({},
607 shape_data[0].shape_gradients_collocation_eo.data(),
608 {});
609 if (!do_integrate)
610 {
611 EvaluatorTensorProductAnisotropic<dim - 1,
612 fe_degree,
613 n_q_points_1d,
614 true>
615 eval;
616 // Evaluate in 3d
617 if (n_blocks[0] == n_rows_n)
618 {
619 eval.template normal<0>(shape_data[0],
620 values_dofs,
621 values);
622 eval.template tangential<1, 0>(shape_data[1],
623 values,
624 values);
625
626 if (evaluation_flag & EvaluationFlags::gradients)
627 {
628 eval.template normal<0>(shape_data[0],
629 values_dofs +
630 n_blocks[0] * n_blocks[1],
631 scratch_data);
632 eval.template tangential<1, 0, dim>(shape_data[1],
633 scratch_data,
634 gradients + 2);
635 }
636 }
637 else if (n_blocks[1] == n_rows_n)
638 {
639 eval.template normal<1>(shape_data[0],
640 values_dofs,
641 values);
642 eval.template tangential<0, 1>(shape_data[1],
643 values,
644 values);
645
646 if (evaluation_flag & EvaluationFlags::gradients)
647 {
648 eval.template normal<1>(shape_data[0],
649 values_dofs +
650 n_blocks[0] * n_blocks[1],
651 scratch_data);
652 eval.template tangential<0, 1, dim>(shape_data[1],
653 scratch_data,
654 gradients + 2);
655 }
656 }
657 else
658 {
659 Eval eval(shape_data[1].shape_values_eo.data(), {}, {});
660 eval.template values<0, true, false>(values_dofs, values);
661 eval.template values<1, true, false>(values, values);
662 if (evaluation_flag & EvaluationFlags::gradients)
663 {
664 eval.template values<0, true, false>(values_dofs +
665 n_blocks[0] *
666 n_blocks[1],
667 scratch_data);
668 eval.template values<1, true, false, dim>(
669 scratch_data, gradients + 2);
670 }
671 }
672 if (evaluation_flag & EvaluationFlags::gradients)
673 {
674 eval_g.template gradients<0, true, false, dim>(values,
675 gradients);
676 eval_g.template gradients<1, true, false, dim>(values,
677 gradients +
678 1);
679 }
680 }
681 else
682 {
683 EvaluatorTensorProductAnisotropic<dim - 1,
684 fe_degree,
685 n_q_points_1d,
686 false>
687 eval;
688 // Integrate in 3d
689 if (evaluation_flag & EvaluationFlags::gradients)
690 {
691 if (evaluation_flag & EvaluationFlags::values)
692 eval_g.template gradients<0, false, true, dim>(
693 gradients, values);
694 else
695 eval_g.template gradients<0, false, false, dim>(
696 gradients, values);
697 eval_g.template gradients<1, false, true, dim>(gradients +
698 1,
699 values);
700 }
701 if (n_blocks[0] == n_rows_n)
702 {
703 eval.template tangential<1, 0>(shape_data[1],
704 values,
705 values);
706 eval.template normal<0>(shape_data[0],
707 values,
708 values_dofs);
709
710 if (evaluation_flag & EvaluationFlags::gradients)
711 {
712 eval.template tangential<1, 0, dim>(shape_data[1],
713 gradients + 2,
714 scratch_data);
715 eval.template normal<0>(shape_data[0],
716 scratch_data,
717 values_dofs +
718 n_blocks[0] * n_blocks[1]);
719 }
720 }
721 else if (n_blocks[1] == n_rows_n)
722 {
723 eval.template tangential<0, 1>(shape_data[1],
724 values,
725 values);
726 eval.template normal<1>(shape_data[0],
727 values,
728 values_dofs);
729
730 if (evaluation_flag & EvaluationFlags::gradients)
731 {
732 eval.template tangential<0, 1, dim>(shape_data[1],
733 gradients + 2,
734 scratch_data);
735 eval.template normal<1>(shape_data[0],
736 scratch_data,
737 values_dofs +
738 n_blocks[0] * n_blocks[1]);
739 }
740 }
741 else
742 {
743 Eval eval_iso(shape_data[1].shape_values_eo.data(),
744 {},
745 {});
746 eval_iso.template values<1, false, false>(values, values);
747 eval_iso.template values<0, false, false>(values,
748 values_dofs);
749 if (evaluation_flag & EvaluationFlags::gradients)
750 {
751 eval_iso.template values<1, false, false, dim>(
752 gradients + 2, scratch_data);
753 eval_iso.template values<0, false, false>(
754 scratch_data,
755 values_dofs + n_blocks[0] * n_blocks[1]);
756 }
757 }
758 }
759 }
760 else
761 {
762 using EvalN = EvaluatorTensorProduct<evaluate_evenodd,
763 dim - 1,
764 fe_degree + 1,
765 n_q_points_1d,
766 Number,
767 Number2>;
768 if (!do_integrate)
769 {
770 // Evaluate in 2d
771 if (n_blocks[0] == n_rows_n)
772 {
773 EvalN eval(shape_data[0].shape_values_eo,
774 shape_data[0].shape_gradients_eo,
775 {});
776 eval.template values<0, true, false>(values_dofs, values);
777 if (evaluation_flag & EvaluationFlags::gradients)
778 {
779 eval.template gradients<0, true, false, dim>(
780 values_dofs, gradients);
781 eval.template values<0, true, false, dim>(
782 values_dofs + n_rows_n, gradients + 1);
783 }
784 }
785 else
786 {
787 Eval eval(shape_data[1].shape_values_eo,
788 shape_data[1].shape_gradients_eo,
789 {});
790 eval.template values<0, true, false>(values_dofs, values);
791 if (evaluation_flag & EvaluationFlags::gradients)
792 {
793 eval.template gradients<0, true, false, dim>(
794 values_dofs, gradients);
795 eval.template values<0, true, false, dim>(
796 values_dofs + n_rows_t, gradients + 1);
797 }
798 }
799 }
800 else
801 {
802 // Integrate in 2d
803 if (n_blocks[0] == n_rows_n)
804 {
805 EvalN eval(shape_data[0].shape_values_eo,
806 shape_data[0].shape_gradients_eo,
807 {});
808 if (evaluation_flag & EvaluationFlags::values)
809 eval.template values<0, false, false>(values,
810 values_dofs);
811 if (evaluation_flag & EvaluationFlags::gradients)
812 {
813 if (evaluation_flag & EvaluationFlags::values)
814 eval.template gradients<0, false, true, dim>(
815 gradients, values_dofs);
816 else
817 eval.template gradients<0, false, false, dim>(
818 gradients, values_dofs);
819 eval.template values<0, false, false, dim>(
820 gradients + 1, values_dofs + n_rows_n);
821 }
822 }
823 else
824 {
825 Eval eval(shape_data[1].shape_values_eo,
826 shape_data[1].shape_gradients_eo,
827 {});
828 if (evaluation_flag & EvaluationFlags::values)
829 eval.template values<0, false, false>(values,
830 values_dofs);
831 if (evaluation_flag & EvaluationFlags::gradients)
832 {
833 if (evaluation_flag & EvaluationFlags::values)
834 eval.template gradients<0, false, true, dim>(
835 gradients, values_dofs);
836 else
837 eval.template gradients<0, false, false, dim>(
838 gradients, values_dofs);
839 eval.template values<0, false, false, dim>(
840 gradients + 1, values_dofs + n_rows_t);
841 }
842 }
843 }
844 }
845 values += Utilities::pow(n_q_points_1d, dim - 1);
846 gradients += dim * Utilities::pow(n_q_points_1d, dim - 1);
847 }
848 }
849 };
850
851
852
853 template <int dim, int fe_degree, typename Number>
854 struct FEFaceNormalEvaluationImpl
855 {
856 using Number2 =
857 typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
858
859 template <bool do_evaluate, bool add_into_output>
860 static void
861 interpolate(const unsigned int n_components,
862 const EvaluationFlags::EvaluationFlags flags,
863 const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info,
864 const Number *input,
865 Number *output,
866 const unsigned int face_no)
867 {
868 Assert(static_cast<unsigned int>(fe_degree) ==
869 shape_info.data.front().fe_degree ||
870 fe_degree == -1,
871 ExcInternalError());
872 if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
873 interpolate_raviart_thomas<do_evaluate, add_into_output>(
874 n_components, input, output, flags, face_no, shape_info);
875 else
876 {
877 const unsigned int fe_degree_ = shape_info.data.front().fe_degree;
878
879 interpolate_generic<do_evaluate, add_into_output>(
880 n_components,
881 input,
882 output,
883 flags,
884 face_no,
885 fe_degree_ + 1,
886 shape_info.data.front().shape_data_on_face,
887 Utilities::pow(fe_degree_ + 1, dim),
888 3 * Utilities::pow(fe_degree_ + 1, dim - 1));
889 }
890 }
891
895 template <bool do_evaluate, bool add_into_output>
896 static void
897 interpolate_quadrature(
898 const unsigned int n_components,
899 const EvaluationFlags::EvaluationFlags flags,
900 const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info,
901 const Number *input,
902 Number *output,
903 const unsigned int face_no)
904 {
905 Assert(static_cast<unsigned int>(fe_degree + 1) ==
906 shape_info.data.front().n_q_points_1d ||
907 fe_degree == -1,
908 ExcInternalError());
909
910 interpolate_generic<do_evaluate, add_into_output>(
911 n_components,
912 input,
913 output,
914 flags,
915 face_no,
916 shape_info.data.front().quadrature.size(),
917 shape_info.data.front().quadrature_data_on_face,
918 shape_info.n_q_points,
919 shape_info.n_q_points_face);
920 }
921
922 private:
923 template <bool do_evaluate, bool add_into_output, int face_direction = 0>
924 static void
925 interpolate_generic(const unsigned int n_components,
926 const Number *input,
927 Number *output,
928 const EvaluationFlags::EvaluationFlags flag,
929 const unsigned int face_no,
930 const unsigned int n_points_1d,
931 const std::array<AlignedVector<Number2>, 2> &shape_data,
932 const unsigned int dofs_per_component_on_cell,
933 const unsigned int dofs_per_component_on_face)
934 {
935 if (face_direction == face_no / 2)
936 {
937 constexpr int stride_ = Utilities::pow(fe_degree + 1, face_direction);
938
939 const int n_rows = fe_degree != -1 ? fe_degree + 1 : n_points_1d;
940 const int stride = Utilities::pow(n_rows, face_direction);
941 const std::array<int, 2> n_blocks{
942 {(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}};
943 std::array<int, 2> steps;
944 if constexpr (face_direction == 0)
945 steps = {{n_rows, 0}};
946 else if constexpr (face_direction == 1 && dim == 2)
947 steps = {{1, 0}};
948 else if constexpr (face_direction == 1)
949 // in 3d, the coordinate system is zx, not xz -> switch indices
950 steps = {{n_rows * n_rows, -n_rows * n_rows * n_rows + 1}};
951 else if constexpr (face_direction == 2)
952 steps = {{1, 0}};
953
954 for (unsigned int c = 0; c < n_components; ++c)
955 {
956 if (flag & EvaluationFlags::hessians)
957 interpolate_to_face<fe_degree + 1,
958 stride_,
959 do_evaluate,
960 add_into_output,
961 2>(shape_data[face_no % 2].begin(),
962 n_blocks,
963 steps,
964 input,
965 output,
966 n_rows,
967 stride);
968 else if (flag & EvaluationFlags::gradients)
969 interpolate_to_face<fe_degree + 1,
970 stride_,
971 do_evaluate,
972 add_into_output,
973 1>(shape_data[face_no % 2].begin(),
974 n_blocks,
975 steps,
976 input,
977 output,
978 n_rows,
979 stride);
980 else
981 interpolate_to_face<fe_degree + 1,
982 stride_,
983 do_evaluate,
984 add_into_output,
985 0>(shape_data[face_no % 2].begin(),
986 n_blocks,
987 steps,
988 input,
989 output,
990 n_rows,
991 stride);
992 if (do_evaluate)
993 {
994 input += dofs_per_component_on_cell;
995 output += dofs_per_component_on_face;
996 }
997 else
998 {
999 output += dofs_per_component_on_cell;
1000 input += dofs_per_component_on_face;
1001 }
1002 }
1003 }
1004 else if (face_direction < dim)
1005 {
1006 interpolate_generic<do_evaluate,
1007 add_into_output,
1008 std::min(face_direction + 1, dim - 1)>(
1009 n_components,
1010 input,
1011 output,
1012 flag,
1013 face_no,
1014 n_points_1d,
1015 shape_data,
1016 dofs_per_component_on_cell,
1017 dofs_per_component_on_face);
1018 }
1019 }
1020
1021 template <bool do_evaluate,
1022 bool add_into_output,
1023 int face_direction = 0,
1024 int max_derivative = 0>
1025 static void
1026 interpolate_raviart_thomas(
1027 const unsigned int n_components,
1028 const Number *input,
1029 Number *output,
1030 const EvaluationFlags::EvaluationFlags flag,
1031 const unsigned int face_no,
1032 const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info)
1033 {
1034 if (dim == 1)
1035 {
1036 // This should never happen since the FE_RaviartThomasNodal is not
1037 // defined for dim = 1. It prevents compiler warnings of infinite
1038 // recursion.
1039 DEAL_II_ASSERT_UNREACHABLE();
1040 return;
1041 }
1042
1043 bool increase_max_der = false;
1044 if ((flag & EvaluationFlags::hessians && max_derivative < 2) ||
1045 (flag & EvaluationFlags::gradients && max_derivative < 1))
1046 increase_max_der = true;
1047
1048 if (face_direction == face_no / 2 && !increase_max_der)
1049 {
1050 constexpr int stride1 = Utilities::pow(fe_degree + 1, face_direction);
1051 constexpr int stride0 = Utilities::pow(fe_degree, face_direction);
1052 constexpr int stride2 = fe_degree * (fe_degree + 1);
1053
1054 const int degree =
1055 fe_degree != -1 ? fe_degree : shape_info.data[0].fe_degree;
1056 const int n_rows_n = degree + 1;
1057 const int n_rows_t = degree;
1058
1059 std::array<int, 3> strides{{1, 1, 1}};
1060 if (face_direction > 0)
1061 {
1062 strides[0] =
1063 n_rows_n * Utilities::pow(n_rows_t, face_direction - 1);
1064 strides[1] = n_rows_t * (face_direction == 3 ? n_rows_n : 1);
1065 strides[2] = Utilities::pow(n_rows_t, face_direction);
1066 }
1067 const ::ndarray<int, 3, 3> dofs_per_direction{
1068 {{{n_rows_n, n_rows_t, n_rows_t}},
1069 {{n_rows_t, n_rows_n, n_rows_t}},
1070 {{n_rows_t, n_rows_t, n_rows_n}}}};
1071
1072 std::array<int, 2> steps, n_blocks;
1073
1074 if constexpr (face_direction == 0)
1075 steps = {{degree + (face_direction == 0), 0}};
1076 else if constexpr (face_direction == 1 && dim == 2)
1077 steps = {{1, 0}};
1078 else if constexpr (face_direction == 1)
1079 // in 3d, the coordinate system is zx, not xz -> switch indices
1080 steps = {
1081 {n_rows_n * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
1082 else if constexpr (face_direction == 2)
1083 steps = {{1, 0}};
1084
1085 n_blocks[0] = dofs_per_direction[0][(face_direction + 1) % dim];
1086 n_blocks[1] =
1087 dim > 2 ? dofs_per_direction[0][(face_direction + 2) % dim] : 1;
1088
1089 interpolate_to_face<
1090 (fe_degree != -1 ? (fe_degree + (face_direction == 0)) : 0),
1091 ((face_direction < 2) ? stride1 : stride2),
1092 do_evaluate,
1093 add_into_output,
1094 max_derivative>(shape_info.data[face_direction != 0]
1095 .shape_data_on_face[face_no % 2]
1096 .begin(),
1097 n_blocks,
1098 steps,
1099 input,
1100 output,
1101 degree + (face_direction == 0),
1102 strides[0]);
1103
1104 if (do_evaluate)
1105 {
1106 input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1107 output += 3 * n_blocks[0] * n_blocks[1];
1108 }
1109 else
1110 {
1111 output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1112 input += 3 * n_blocks[0] * n_blocks[1];
1113 }
1114
1115 // must only change steps only for face direction 0
1116 if constexpr (face_direction == 0)
1117 steps = {{degree, 0}};
1118
1119 n_blocks[0] = dofs_per_direction[1][(face_direction + 1) % dim];
1120 n_blocks[1] =
1121 dim > 2 ? dofs_per_direction[1][(face_direction + 2) % dim] : 1;
1122
1123 interpolate_to_face<
1124 (fe_degree != -1 ? (fe_degree + (face_direction == 1)) : 0),
1125 ((face_direction < 2) ? stride0 : stride2),
1126 do_evaluate,
1127 add_into_output,
1128 max_derivative>(shape_info.data[face_direction != 1]
1129 .shape_data_on_face[face_no % 2]
1130 .begin(),
1131 n_blocks,
1132 steps,
1133 input,
1134 output,
1135 degree + (face_direction == 1),
1136 strides[1]);
1137
1138 if constexpr (dim > 2)
1139 {
1140 if (do_evaluate)
1141 {
1142 input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1143 output += 3 * n_blocks[0] * n_blocks[1];
1144 }
1145 else
1146 {
1147 output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
1148 input += 3 * n_blocks[0] * n_blocks[1];
1149 }
1150
1151 if constexpr (face_direction == 0)
1152 steps = {{degree, 0}};
1153 else if constexpr (face_direction == 1)
1154 // in 3d, the coordinate system is zx, not xz -> switch indices
1155 steps = {
1156 {n_rows_t * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
1157 else if constexpr (face_direction == 2)
1158 steps = {{1, 0}};
1159
1160 n_blocks[0] = dofs_per_direction[2][(face_direction + 1) % dim];
1161 n_blocks[1] = dofs_per_direction[2][(face_direction + 2) % dim];
1162
1163 interpolate_to_face<
1164 (fe_degree != -1 ? (fe_degree + (face_direction == 2)) : 0),
1165 stride0,
1166 do_evaluate,
1167 add_into_output,
1168 max_derivative>(shape_info.data[face_direction != 2]
1169 .shape_data_on_face[face_no % 2]
1170 .begin(),
1171 n_blocks,
1172 steps,
1173 input,
1174 output,
1175 degree + (face_direction == 2),
1176 strides[2]);
1177 }
1178 }
1179 else if (face_direction == face_no / 2)
1180 {
1181 // Only increase max_derivative
1182 interpolate_raviart_thomas<do_evaluate,
1183 add_into_output,
1184 face_direction,
1185 std::min(max_derivative + 1, 2)>(
1186 n_components, input, output, flag, face_no, shape_info);
1187 }
1188 else if (face_direction < dim)
1189 {
1190 if (increase_max_der)
1191 {
1192 interpolate_raviart_thomas<do_evaluate,
1193 add_into_output,
1194 std::min(face_direction + 1, dim - 1),
1195 std::min(max_derivative + 1, 2)>(
1196 n_components, input, output, flag, face_no, shape_info);
1197 }
1198 else
1199 {
1200 interpolate_raviart_thomas<do_evaluate,
1201 add_into_output,
1202 std::min(face_direction + 1, dim - 1),
1203 max_derivative>(
1204 n_components, input, output, flag, face_no, shape_info);
1205 }
1206 }
1207 }
1208 };
1209
1210
1211
1212 // internal helper function for reading data; base version of different types
1213 template <typename VectorizedArrayType, typename Number2>
1214 void
1215 do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
1216 {
1217 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1218 dst[v] = src_ptr[v];
1219 }
1220
1221
1222
1223 // internal helper function for reading data; specialized version where we
1224 // can use a dedicated load function
1225 template <typename Number, std::size_t width>
1226 void
1227 do_vectorized_read(const Number *src_ptr, VectorizedArray<Number, width> &dst)
1228 {
1229 dst.load(src_ptr);
1230 }
1231
1232
1233
1234 // internal helper function for reading data; base version of different types
1235 template <typename VectorizedArrayType, typename Number2>
1236 void
1237 do_vectorized_gather(const Number2 *src_ptr,
1238 const unsigned int *indices,
1239 VectorizedArrayType &dst)
1240 {
1241 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1242 dst[v] = src_ptr[indices[v]];
1243 }
1244
1245
1246
1247 // internal helper function for reading data; specialized version where we
1248 // can use a dedicated gather function
1249 template <typename Number, std::size_t width>
1250 void
1251 do_vectorized_gather(const Number *src_ptr,
1252 const unsigned int *indices,
1253 VectorizedArray<Number, width> &dst)
1254 {
1255 dst.gather(src_ptr, indices);
1256 }
1257
1258
1259
1260 // internal helper function for reading data; base version of different types
1261 template <typename VectorizedArrayType, typename Number2>
1262 void
1263 do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
1264 {
1265 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1266 dst_ptr[v] += src[v];
1267 }
1268
1269
1270
1271 // internal helper function for reading data; specialized version where we
1272 // can use a dedicated load function
1273 template <typename Number, std::size_t width>
1274 void
1275 do_vectorized_add(const VectorizedArray<Number, width> src, Number *dst_ptr)
1276 {
1277 VectorizedArray<Number, width> tmp;
1278 tmp.load(dst_ptr);
1279 (tmp + src).store(dst_ptr);
1280 }
1281
1282
1283
1284 // internal helper function for reading data; base version of different types
1285 template <typename VectorizedArrayType, typename Number2>
1286 void
1287 do_vectorized_scatter_add(const VectorizedArrayType src,
1288 const unsigned int *indices,
1289 Number2 *dst_ptr)
1290 {
1291 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
1292 dst_ptr[indices[v]] += src[v];
1293 }
1294
1295
1296
1297 // internal helper function for reading data; specialized version where we
1298 // can use a dedicated gather function
1299 template <typename Number, std::size_t width>
1300 void
1301 do_vectorized_scatter_add(const VectorizedArray<Number, width> src,
1302 const unsigned int *indices,
1303 Number *dst_ptr)
1304 {
1305#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
1306 for (unsigned int v = 0; v < width; ++v)
1307 dst_ptr[indices[v]] += src[v];
1308#else
1309 VectorizedArray<Number, width> tmp;
1310 tmp.gather(dst_ptr, indices);
1311 (tmp + src).scatter(indices, dst_ptr);
1312#endif
1313 }
1314
1315
1316
1317 template <typename Number>
1318 void
1319 adjust_for_face_orientation(const unsigned int dim,
1320 const unsigned int n_components,
1321 const EvaluationFlags::EvaluationFlags flag,
1322 const unsigned int *orientation,
1323 const bool integrate,
1324 const std::size_t n_q_points,
1325 Number *tmp_values,
1326 Number *values_quad,
1327 Number *gradients_quad,
1328 Number *hessians_quad)
1329 {
1330 for (unsigned int c = 0; c < n_components; ++c)
1331 {
1332 if (flag & EvaluationFlags::values)
1333 {
1334 if (integrate)
1335 for (unsigned int q = 0; q < n_q_points; ++q)
1336 tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
1337 else
1338 for (unsigned int q = 0; q < n_q_points; ++q)
1339 tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
1340 for (unsigned int q = 0; q < n_q_points; ++q)
1341 values_quad[c * n_q_points + q] = tmp_values[q];
1342 }
1343 if (flag & EvaluationFlags::gradients)
1344 for (unsigned int d = 0; d < dim; ++d)
1345 {
1346 if (integrate)
1347 for (unsigned int q = 0; q < n_q_points; ++q)
1348 tmp_values[q] =
1349 gradients_quad[(c * n_q_points + orientation[q]) * dim + d];
1350 else
1351 for (unsigned int q = 0; q < n_q_points; ++q)
1352 tmp_values[orientation[q]] =
1353 gradients_quad[(c * n_q_points + q) * dim + d];
1354 for (unsigned int q = 0; q < n_q_points; ++q)
1355 gradients_quad[(c * n_q_points + q) * dim + d] = tmp_values[q];
1356 }
1357 if (flag & EvaluationFlags::hessians)
1358 {
1359 const unsigned int hdim = (dim * (dim + 1)) / 2;
1360 for (unsigned int d = 0; d < hdim; ++d)
1361 {
1362 if (integrate)
1363 for (unsigned int q = 0; q < n_q_points; ++q)
1364 tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
1365 orientation[q]];
1366 else
1367 for (unsigned int q = 0; q < n_q_points; ++q)
1368 tmp_values[orientation[q]] =
1369 hessians_quad[(c * hdim + d) * n_q_points + q];
1370 for (unsigned int q = 0; q < n_q_points; ++q)
1371 hessians_quad[(c * hdim + d) * n_q_points + q] =
1372 tmp_values[q];
1373 }
1374 }
1375 }
1376 }
1377
1378
1379
1380 template <typename Number, typename VectorizedArrayType>
1381 void
1382 adjust_for_face_orientation_per_lane(
1383 const unsigned int dim,
1384 const unsigned int n_components,
1385 const unsigned int v,
1386 const EvaluationFlags::EvaluationFlags flag,
1387 const unsigned int *orientation,
1388 const bool integrate,
1389 const std::size_t n_q_points,
1390 Number *tmp_values,
1391 VectorizedArrayType *values_quad,
1392 VectorizedArrayType *gradients_quad = nullptr,
1393 VectorizedArrayType *hessians_quad = nullptr)
1394 {
1395 for (unsigned int c = 0; c < n_components; ++c)
1396 {
1397 if (flag & EvaluationFlags::values)
1398 {
1399 if (integrate)
1400 for (unsigned int q = 0; q < n_q_points; ++q)
1401 tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v];
1402 else
1403 for (unsigned int q = 0; q < n_q_points; ++q)
1404 tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v];
1405 for (unsigned int q = 0; q < n_q_points; ++q)
1406 values_quad[c * n_q_points + q][v] = tmp_values[q];
1407 }
1408 if (flag & EvaluationFlags::gradients)
1409 for (unsigned int d = 0; d < dim; ++d)
1410 {
1411 Assert(gradients_quad != nullptr, ExcInternalError());
1412 if (integrate)
1413 for (unsigned int q = 0; q < n_q_points; ++q)
1414 tmp_values[q] =
1415 gradients_quad[(c * n_q_points + orientation[q]) * dim + d]
1416 [v];
1417 else
1418 for (unsigned int q = 0; q < n_q_points; ++q)
1419 tmp_values[orientation[q]] =
1420 gradients_quad[(c * n_q_points + q) * dim + d][v];
1421 for (unsigned int q = 0; q < n_q_points; ++q)
1422 gradients_quad[(c * n_q_points + q) * dim + d][v] =
1423 tmp_values[q];
1424 }
1425 if (flag & EvaluationFlags::hessians)
1426 {
1427 Assert(hessians_quad != nullptr, ExcInternalError());
1428 const unsigned int hdim = (dim * (dim + 1)) / 2;
1429 for (unsigned int d = 0; d < hdim; ++d)
1430 {
1431 if (integrate)
1432 for (unsigned int q = 0; q < n_q_points; ++q)
1433 tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
1434 orientation[q]][v];
1435 else
1436 for (unsigned int q = 0; q < n_q_points; ++q)
1437 tmp_values[orientation[q]] =
1438 hessians_quad[(c * hdim + d) * n_q_points + q][v];
1439 for (unsigned int q = 0; q < n_q_points; ++q)
1440 hessians_quad[(c * hdim + d) * n_q_points + q][v] =
1441 tmp_values[q];
1442 }
1443 }
1444 }
1445 }
1446
1447
1448
1449 template <int dim, typename Number>
1450 struct FEFaceEvaluationImplEvaluateSelector
1451 {
1452 static bool
1453 evaluate_tensor_none(const unsigned int n_components,
1454 const EvaluationFlags::EvaluationFlags evaluation_flag,
1455 const Number *values_dofs,
1456 FEEvaluationData<dim, Number, true> &fe_eval)
1457 {
1458 const auto &shape_info = fe_eval.get_shape_info();
1459 const auto &shape_data = shape_info.data.front();
1460 using Number2 =
1461 typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
1462
1463 Assert((fe_eval.get_dof_access_index() ==
1464 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
1465 fe_eval.is_interior_face() == false) == false,
1466 ExcNotImplemented());
1467
1468 const unsigned int face_no = fe_eval.get_face_no();
1469 const unsigned int face_orientation = fe_eval.get_face_orientation();
1470 const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
1471 const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
1472
1473 if (evaluation_flag & EvaluationFlags::values)
1474 {
1475 const auto *const shape_values =
1476 &shape_data.shape_values_face(face_no, face_orientation, 0);
1477
1478 auto *out = fe_eval.begin_values();
1479 auto *in = values_dofs;
1480
1481 for (unsigned int c = 0; c < n_components; c += 3)
1482 {
1483 if (c + 1 == n_components)
1484 apply_matrix_vector_product<evaluate_general,
1485 EvaluatorQuantity::value,
1486 /*transpose_matrix*/ true,
1487 /*add*/ false,
1488 /*consider_strides*/ false,
1489 Number,
1490 Number2,
1491 /*n_components*/ 1>(
1492 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1493 else if (c + 2 == n_components)
1494 apply_matrix_vector_product<evaluate_general,
1495 EvaluatorQuantity::value,
1496 /*transpose_matrix*/ true,
1497 /*add*/ false,
1498 /*consider_strides*/ false,
1499 Number,
1500 Number2,
1501 /*n_components*/ 2>(
1502 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1503 else
1504 apply_matrix_vector_product<evaluate_general,
1505 EvaluatorQuantity::value,
1506 /*transpose_matrix*/ true,
1507 /*add*/ false,
1508 /*consider_strides*/ false,
1509 Number,
1510 Number2,
1511 /*n_components*/ 3>(
1512 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1513
1514 out += 3 * n_q_points;
1515 in += 3 * n_dofs;
1516 }
1517 }
1518
1519 if (evaluation_flag & EvaluationFlags::gradients)
1520 {
1521 auto *out = fe_eval.begin_gradients();
1522 const auto *in = values_dofs;
1523
1524 const auto *const shape_gradients =
1525 &shape_data.shape_gradients_face(face_no, face_orientation, 0);
1526
1527 for (unsigned int c = 0; c < n_components; c += 3)
1528 {
1529 if (c + 1 == n_components)
1530 apply_matrix_vector_product<evaluate_general,
1531 EvaluatorQuantity::value,
1532 /*transpose_matrix*/ true,
1533 /*add*/ false,
1534 /*consider_strides*/ false,
1535 Number,
1536 Number2,
1537 /*n_components*/ 1>(
1538 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
1539 else if (c + 2 == n_components)
1540 apply_matrix_vector_product<evaluate_general,
1541 EvaluatorQuantity::value,
1542 /*transpose_matrix*/ true,
1543 /*add*/ false,
1544 /*consider_strides*/ false,
1545 Number,
1546 Number2,
1547 /*n_components*/ 2>(
1548 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
1549 else
1550 apply_matrix_vector_product<evaluate_general,
1551 EvaluatorQuantity::value,
1552 /*transpose_matrix*/ true,
1553 /*add*/ false,
1554 /*consider_strides*/ false,
1555 Number,
1556 Number2,
1557 /*n_components*/ 3>(
1558 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
1559 out += 3 * n_q_points * dim;
1560 in += 3 * n_dofs;
1561 }
1562 }
1563
1564 Assert(!(evaluation_flag & EvaluationFlags::hessians),
1565 ExcNotImplemented());
1566
1567 return true;
1568 }
1569
1570 template <int fe_degree>
1571#ifndef DEBUG
1572 DEAL_II_ALWAYS_INLINE
1573#endif
1574 static void
1575 project_to_face(const unsigned int n_components,
1576 const EvaluationFlags::EvaluationFlags evaluation_flag,
1577 const Number *values_dofs,
1578 FEEvaluationData<dim, Number, true> &fe_eval,
1579 const bool use_vectorization,
1580 Number *temp,
1581 Number *scratch_data)
1582 {
1583 const auto &shape_info = fe_eval.get_shape_info();
1584
1585 if (use_vectorization == false)
1586 {
1587 const auto &shape_data = shape_info.data.front();
1588
1589 const unsigned int dofs_per_comp_face =
1590 fe_degree > -1 ?
1591 Utilities::pow(fe_degree + 1, dim - 1) :
1592 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1593 const unsigned int dofs_per_face = n_components * dofs_per_comp_face;
1594
1595 for (unsigned int v = 0; v < Number::size(); ++v)
1596 {
1597 // the loop breaks once an invalid_unsigned_int is hit for
1598 // all cases except the exterior faces in the ECL loop (where
1599 // some faces might be at the boundaries but others not)
1600 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
1601 {
1602 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
1603 temp[i][v] = 0;
1604 continue;
1605 }
1606
1607 FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
1608 template interpolate<true, false>(n_components,
1609 evaluation_flag,
1610 shape_info,
1611 values_dofs,
1612 scratch_data,
1613 fe_eval.get_face_no(v));
1614
1615 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
1616 temp[i][v] = scratch_data[i][v];
1617 }
1618 }
1619 else
1620 FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
1621 template interpolate<true, false>(n_components,
1622 evaluation_flag,
1623 shape_info,
1624 values_dofs,
1625 temp,
1626 fe_eval.get_face_no());
1627 }
1628
1629
1630 template <int fe_degree, int n_q_points_1d>
1631#ifndef DEBUG
1632 DEAL_II_ALWAYS_INLINE
1633#endif
1634 static void
1635 evaluate_in_face(const unsigned int n_components,
1636 const EvaluationFlags::EvaluationFlags evaluation_flag,
1637 FEEvaluationData<dim, Number, true> &fe_eval,
1638 Number *temp,
1639 Number *scratch_data)
1640 {
1641 const auto &shape_info = fe_eval.get_shape_info();
1642 const auto &shape_data = shape_info.data.front();
1643
1644 const unsigned int subface_index = fe_eval.get_subface_index();
1645 constexpr unsigned int n_q_points_1d_actual =
1646 fe_degree > -1 ? n_q_points_1d : 0;
1647
1648 if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
1649 {
1650 FEFaceEvaluationImplRaviartThomas<dim,
1651 fe_degree,
1652 n_q_points_1d_actual,
1653 Number>::
1654 template evaluate_or_integrate_in_face<false>(
1655 evaluation_flag,
1656 fe_eval.get_shape_info().data,
1657 temp,
1658 fe_eval.begin_values(),
1659 fe_eval.begin_gradients(),
1660 scratch_data,
1661 subface_index,
1662 fe_eval.get_face_no() / 2);
1663 }
1664 else if (fe_degree > -1 &&
1665 subface_index >= GeometryInfo<dim>::max_children_per_cell &&
1666 shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
1667 FEFaceEvaluationImpl<true,
1668 dim,
1669 fe_degree,
1670 n_q_points_1d_actual,
1671 Number>::evaluate_in_face(n_components,
1672 evaluation_flag,
1673 shape_data,
1674 temp,
1675 fe_eval.begin_values(),
1676 fe_eval
1677 .begin_gradients(),
1678 fe_eval.begin_hessians(),
1679 scratch_data,
1680 subface_index);
1681 else
1682 FEFaceEvaluationImpl<false,
1683 dim,
1684 fe_degree,
1685 n_q_points_1d_actual,
1686 Number>::evaluate_in_face(n_components,
1687 evaluation_flag,
1688 shape_data,
1689 temp,
1690 fe_eval.begin_values(),
1691 fe_eval
1692 .begin_gradients(),
1693 fe_eval.begin_hessians(),
1694 scratch_data,
1695 subface_index);
1696 }
1697
1698#ifndef DEBUG
1699 DEAL_II_ALWAYS_INLINE
1700#endif
1701 static void
1702 adjust_quadrature_for_face_orientation(
1703 const unsigned int n_components,
1704 const EvaluationFlags::EvaluationFlags evaluation_flag,
1705 FEEvaluationData<dim, Number, true> &fe_eval,
1706 const bool use_vectorization,
1707 Number *temp)
1708 {
1709 const auto &shape_info = fe_eval.get_shape_info();
1710
1711 if (use_vectorization == false)
1712 {
1713 for (unsigned int v = 0; v < Number::size(); ++v)
1714 {
1715 // the loop breaks once an invalid_unsigned_int is hit for
1716 // all cases except the exterior faces in the ECL loop (where
1717 // some faces might be at the boundaries but others not)
1718 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
1719 continue;
1720
1721 if (fe_eval.get_face_orientation(v) != 0)
1722 adjust_for_face_orientation_per_lane(
1723 dim,
1724 n_components,
1725 v,
1726 evaluation_flag,
1727 &shape_info.face_orientations_quad(
1728 fe_eval.get_face_orientation(v), 0),
1729 false,
1730 shape_info.n_q_points_face,
1731 &temp[0][0],
1732 fe_eval.begin_values(),
1733 fe_eval.begin_gradients(),
1734 fe_eval.begin_hessians());
1735 }
1736 }
1737 else if (fe_eval.get_face_orientation() != 0)
1738 adjust_for_face_orientation(
1739 dim,
1740 n_components,
1741 evaluation_flag,
1742 &shape_info.face_orientations_quad(fe_eval.get_face_orientation(), 0),
1743 false,
1744 shape_info.n_q_points_face,
1745 temp,
1746 fe_eval.begin_values(),
1747 fe_eval.begin_gradients(),
1748 fe_eval.begin_hessians());
1749 }
1750
1751
1752
1753 template <int fe_degree, int n_q_points_1d>
1754 static bool
1755 evaluate_tensor(const unsigned int n_components,
1756 const EvaluationFlags::EvaluationFlags evaluation_flag,
1757 const Number *values_dofs_actual,
1758 FEEvaluationData<dim, Number, true> &fe_eval)
1759 {
1760 const auto &shape_info = fe_eval.get_shape_info();
1761 const auto &shape_data = shape_info.data.front();
1762
1763 const unsigned int dofs_per_comp_face =
1764 fe_degree > -1 ?
1765 Utilities::pow(fe_degree + 1, dim - 1) :
1766 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1767
1768 // Note: we always keep storage of values, 1st and 2nd derivatives in an
1769 // array, so reserve space for all three here
1770 Number *temp1 = fe_eval.get_scratch_data().begin();
1771 Number *temp2 = temp1 + 3 * n_components * dofs_per_comp_face;
1772
1773 const Number *values_dofs =
1774 (shape_data.element_type == MatrixFreeFunctions::truncated_tensor) ?
1775 temp2 +
1776 2 * (std::max(fe_eval.get_shape_info().dofs_per_component_on_cell,
1777 shape_info.n_q_points)) :
1778 values_dofs_actual;
1779
1780 if (shape_data.element_type == MatrixFreeFunctions::truncated_tensor)
1781 embed_truncated_into_full_tensor_product<dim, fe_degree>(
1782 n_components,
1783 const_cast<Number *>(values_dofs),
1784 values_dofs_actual,
1785 fe_eval);
1786
1787 bool use_vectorization = true;
1788 if (fe_eval.get_dof_access_index() ==
1789 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
1790 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
1791 for (unsigned int v = 0; v < Number::size(); ++v)
1792 if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
1793 fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
1794 use_vectorization = false;
1795
1796 project_to_face<fe_degree>(n_components,
1797 evaluation_flag,
1798 values_dofs,
1799 fe_eval,
1800 use_vectorization,
1801 temp1,
1802 temp2);
1803
1804 evaluate_in_face<fe_degree, n_q_points_1d>(
1805 n_components, evaluation_flag, fe_eval, temp1, temp2);
1806
1807 if (dim == 3)
1808 adjust_quadrature_for_face_orientation(
1809 n_components, evaluation_flag, fe_eval, use_vectorization, temp1);
1810
1811 return false;
1812 }
1813
1814 template <int fe_degree, int n_q_points_1d>
1815 static bool
1816 run(const unsigned int n_components,
1817 const EvaluationFlags::EvaluationFlags evaluation_flag,
1818 const Number *values_dofs,
1819 FEEvaluationData<dim, Number, true> &fe_eval)
1820 {
1821 const auto &shape_info = fe_eval.get_shape_info();
1822
1823 if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
1824 return evaluate_tensor_none(n_components,
1825 evaluation_flag,
1826 values_dofs,
1827 fe_eval);
1828 else
1829 return evaluate_tensor<fe_degree, n_q_points_1d>(n_components,
1830 evaluation_flag,
1831 values_dofs,
1832 fe_eval);
1833 }
1834 };
1835
1836
1837
1838 template <int dim, typename Number>
1839 struct FEFaceEvaluationImplProjectToFaceSelector
1840 {
1841 template <int fe_degree>
1842 static bool
1843 run(const unsigned int n_components,
1844 const EvaluationFlags::EvaluationFlags evaluation_flag,
1845 const Number *values_dofs,
1846 FEEvaluationData<dim, Number, true> &fe_eval)
1847 {
1848 const auto &shape_info = fe_eval.get_shape_info();
1849 const auto &shape_data = shape_info.data.front();
1850
1851 const unsigned int dofs_per_comp_face =
1852 fe_degree > -1 ?
1853 Utilities::pow(fe_degree + 1, dim - 1) :
1854 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1855
1856 // Note: we always keep storage of values, 1st and 2nd derivatives in an
1857 // array, so reserve space for all three here
1858 Number *temp = fe_eval.get_scratch_data().begin();
1859 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
1860
1861 bool use_vectorization = true;
1862 if (fe_eval.get_dof_access_index() ==
1863 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
1864 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
1865 for (unsigned int v = 0; v < Number::size(); ++v)
1866 if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
1867 fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
1868 use_vectorization = false;
1869
1870 FEFaceEvaluationImplEvaluateSelector<dim, Number>::
1871 template project_to_face<fe_degree>(n_components,
1872 evaluation_flag,
1873 values_dofs,
1874 fe_eval,
1875 use_vectorization,
1876 temp,
1877 scratch_data);
1878
1879 return false;
1880 }
1881 };
1882
1883
1884
1885 template <int dim, typename Number>
1886 struct FEFaceEvaluationImplEvaluateInFaceSelector
1887 {
1888 template <int fe_degree, int n_q_points_1d>
1889 static bool
1890 run(const unsigned int n_components,
1891 const EvaluationFlags::EvaluationFlags evaluation_flag,
1892 FEEvaluationData<dim, Number, true> &fe_eval)
1893 {
1894 const auto &shape_info = fe_eval.get_shape_info();
1895 const auto &shape_data = shape_info.data.front();
1896
1897 const unsigned int dofs_per_comp_face =
1898 fe_degree > -1 ?
1899 Utilities::pow(fe_degree + 1, dim - 1) :
1900 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
1901
1902 // Note: we always keep storage of values, 1st and 2nd derivatives in an
1903 // array, so reserve space for all three here
1904 Number *temp = fe_eval.get_scratch_data().begin();
1905 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
1906
1907 FEFaceEvaluationImplEvaluateSelector<dim, Number>::
1908 template evaluate_in_face<fe_degree, n_q_points_1d>(
1909 n_components, evaluation_flag, fe_eval, temp, scratch_data);
1910
1911 return false;
1912 }
1913 };
1914
1915
1916
1917 template <int dim, typename Number>
1918 struct FEFaceEvaluationImplIntegrateSelector
1919 {
1920 static bool
1921 integrate_tensor_none(
1922 const unsigned int n_components,
1923 const EvaluationFlags::EvaluationFlags integration_flag,
1924 Number *values_dofs,
1925 FEEvaluationData<dim, Number, true> &fe_eval,
1926 const bool sum_into_values)
1927 {
1928 const auto &shape_info = fe_eval.get_shape_info();
1929 const auto &shape_data = shape_info.data.front();
1930 using Number2 =
1931 typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
1932
1933 Assert((fe_eval.get_dof_access_index() ==
1934 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
1935 fe_eval.is_interior_face() == false) == false,
1936 ExcNotImplemented());
1937
1938 const unsigned int face_no = fe_eval.get_face_no();
1939 const unsigned int face_orientation = fe_eval.get_face_orientation();
1940 const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
1941 const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
1942
1943
1944 if (integration_flag & EvaluationFlags::values)
1945 {
1946 const auto *const shape_values =
1947 &shape_data.shape_values_face(face_no, face_orientation, 0);
1948
1949 auto *in = fe_eval.begin_values();
1950 auto *out = values_dofs;
1951
1952 for (unsigned int c = 0; c < n_components; c += 3)
1953 {
1954 if (sum_into_values)
1955 {
1956 if (c + 1 == n_components)
1957 apply_matrix_vector_product<evaluate_general,
1958 EvaluatorQuantity::value,
1959 /*transpose_matrix*/ false,
1960 /*add*/ true,
1961 /*consider_strides*/ false,
1962 Number,
1963 Number2,
1964 /*n_components*/ 1>(
1965 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1966 else if (c + 2 == n_components)
1967 apply_matrix_vector_product<evaluate_general,
1968 EvaluatorQuantity::value,
1969 /*transpose_matrix*/ false,
1970 /*add*/ true,
1971 /*consider_strides*/ false,
1972 Number,
1973 Number2,
1974 /*n_components*/ 2>(
1975 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1976 else
1977 apply_matrix_vector_product<evaluate_general,
1978 EvaluatorQuantity::value,
1979 /*transpose_matrix*/ false,
1980 /*add*/ true,
1981 /*consider_strides*/ false,
1982 Number,
1983 Number2,
1984 /*n_components*/ 3>(
1985 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1986 }
1987 else
1988 {
1989 if (c + 1 == n_components)
1990 apply_matrix_vector_product<evaluate_general,
1991 EvaluatorQuantity::value,
1992 /*transpose_matrix*/ false,
1993 /*add*/ false,
1994 /*consider_strides*/ false,
1995 Number,
1996 Number2,
1997 /*n_components*/ 1>(
1998 shape_values, in, out, n_dofs, n_q_points, 1, 1);
1999 else if (c + 2 == n_components)
2000 apply_matrix_vector_product<evaluate_general,
2001 EvaluatorQuantity::value,
2002 /*transpose_matrix*/ false,
2003 /*add*/ false,
2004 /*consider_strides*/ false,
2005 Number,
2006 Number2,
2007 /*n_components*/ 2>(
2008 shape_values, in, out, n_dofs, n_q_points, 1, 1);
2009 else
2010 apply_matrix_vector_product<evaluate_general,
2011 EvaluatorQuantity::value,
2012 /*transpose_matrix*/ false,
2013 /*add*/ false,
2014 /*consider_strides*/ false,
2015 Number,
2016 Number2,
2017 /*n_components*/ 3>(
2018 shape_values, in, out, n_dofs, n_q_points, 1, 1);
2019 }
2020 in += 3 * n_q_points;
2021 out += 3 * n_dofs;
2022 }
2023 }
2024
2025 if (integration_flag & EvaluationFlags::gradients)
2026 {
2027 auto *in = fe_eval.begin_gradients();
2028 auto *out = values_dofs;
2029
2030 const auto *const shape_gradients =
2031 &shape_data.shape_gradients_face(face_no, face_orientation, 0);
2032
2033 for (unsigned int c = 0; c < n_components; ++c)
2034 {
2035 if (!sum_into_values &&
2036 !(integration_flag & EvaluationFlags::values))
2037 {
2038 if (c + 1 == n_components)
2039 apply_matrix_vector_product<evaluate_general,
2040 EvaluatorQuantity::value,
2041 /*transpose_matrix*/ false,
2042 /*add*/ false,
2043 /*consider_strides*/ false,
2044 Number,
2045 Number2,
2046 /*n_components*/ 1>(
2047 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2048 else if (c + 2 == n_components)
2049 apply_matrix_vector_product<evaluate_general,
2050 EvaluatorQuantity::value,
2051 /*transpose_matrix*/ false,
2052 /*add*/ false,
2053 /*consider_strides*/ false,
2054 Number,
2055 Number2,
2056 /*n_components*/ 2>(
2057 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2058 else
2059 apply_matrix_vector_product<evaluate_general,
2060 EvaluatorQuantity::value,
2061 /*transpose_matrix*/ false,
2062 /*add*/ false,
2063 /*consider_strides*/ false,
2064 Number,
2065 Number2,
2066 /*n_components*/ 3>(
2067 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2068 }
2069 else
2070 {
2071 if (c + 1 == n_components)
2072 apply_matrix_vector_product<evaluate_general,
2073 EvaluatorQuantity::value,
2074 /*transpose_matrix*/ false,
2075 /*add*/ true,
2076 /*consider_strides*/ false,
2077 Number,
2078 Number2,
2079 /*n_components*/ 1>(
2080 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2081 else if (c + 2 == n_components)
2082 apply_matrix_vector_product<evaluate_general,
2083 EvaluatorQuantity::value,
2084 /*transpose_matrix*/ false,
2085 /*add*/ true,
2086 /*consider_strides*/ false,
2087 Number,
2088 Number2,
2089 /*n_components*/ 2>(
2090 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2091 else
2092 apply_matrix_vector_product<evaluate_general,
2093 EvaluatorQuantity::value,
2094 /*transpose_matrix*/ false,
2095 /*add*/ true,
2096 /*consider_strides*/ false,
2097 Number,
2098 Number2,
2099 /*n_components*/ 3>(
2100 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
2101 }
2102 in += 3 * n_q_points * dim;
2103 out += 3 * n_dofs;
2104 }
2105 }
2106
2107 Assert(!(integration_flag & EvaluationFlags::hessians),
2108 ExcNotImplemented());
2109
2110 return true;
2111 }
2112
2113#ifndef DEBUG
2114 DEAL_II_ALWAYS_INLINE
2115#endif
2116 static void
2117 adjust_quadrature_for_face_orientation(
2118 const unsigned int n_components,
2119 const EvaluationFlags::EvaluationFlags integration_flag,
2120 FEEvaluationData<dim, Number, true> &fe_eval,
2121 const bool use_vectorization,
2122 Number *temp)
2123 {
2124 const auto &shape_info = fe_eval.get_shape_info();
2125
2126 if (use_vectorization == false)
2127 {
2128 for (unsigned int v = 0; v < Number::size(); ++v)
2129 {
2130 // the loop breaks once an invalid_unsigned_int is hit for
2131 // all cases except the exterior faces in the ECL loop (where
2132 // some faces might be at the boundaries but others not)
2133 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2134 continue;
2135
2136 if (fe_eval.get_face_orientation(v) != 0)
2137 adjust_for_face_orientation_per_lane(
2138 dim,
2139 n_components,
2140 v,
2141 integration_flag,
2142 &fe_eval.get_shape_info().face_orientations_quad(
2143 fe_eval.get_face_orientation(v), 0),
2144 true,
2145 shape_info.n_q_points_face,
2146 &temp[0][0],
2147 fe_eval.begin_values(),
2148 fe_eval.begin_gradients(),
2149 fe_eval.begin_hessians());
2150 }
2151 }
2152 else if (fe_eval.get_face_orientation() != 0)
2153 adjust_for_face_orientation(
2154 dim,
2155 n_components,
2156 integration_flag,
2157 &fe_eval.get_shape_info().face_orientations_quad(
2158 fe_eval.get_face_orientation(), 0),
2159 true,
2160 shape_info.n_q_points_face,
2161 temp,
2162 fe_eval.begin_values(),
2163 fe_eval.begin_gradients(),
2164 fe_eval.begin_hessians());
2165 }
2166
2167 template <int fe_degree, int n_q_points_1d>
2168#ifndef DEBUG
2169 DEAL_II_ALWAYS_INLINE
2170#endif
2171 static void
2172 integrate_in_face(const unsigned int n_components,
2173 const EvaluationFlags::EvaluationFlags integration_flag,
2174 FEEvaluationData<dim, Number, true> &fe_eval,
2175 Number *temp,
2176 Number *scratch_data)
2177 {
2178 const auto &shape_info = fe_eval.get_shape_info();
2179 const auto &shape_data = shape_info.data.front();
2180
2181 const unsigned int n_q_points_1d_actual =
2182 fe_degree > -1 ? n_q_points_1d : 0;
2183 const unsigned int subface_index = fe_eval.get_subface_index();
2184
2185 if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
2186 {
2187 FEFaceEvaluationImplRaviartThomas<dim,
2188 fe_degree,
2189 n_q_points_1d_actual,
2190 Number>::
2191 template evaluate_or_integrate_in_face<true>(
2192 integration_flag,
2193 fe_eval.get_shape_info().data,
2194 temp,
2195 fe_eval.begin_values(),
2196 fe_eval.begin_gradients(),
2197 scratch_data,
2198 subface_index,
2199 fe_eval.get_face_no() / 2);
2200 }
2201 else if (fe_degree > -1 &&
2202 fe_eval.get_subface_index() >=
2203 GeometryInfo<dim - 1>::max_children_per_cell &&
2204 shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
2205 FEFaceEvaluationImpl<
2206 true,
2207 dim,
2208 fe_degree,
2209 n_q_points_1d_actual,
2210 Number>::integrate_in_face(n_components,
2211 integration_flag,
2212 shape_data,
2213 temp,
2214 fe_eval.begin_values(),
2215 fe_eval.begin_gradients(),
2216 fe_eval.begin_hessians(),
2217 scratch_data,
2218 subface_index);
2219 else
2220 FEFaceEvaluationImpl<
2221 false,
2222 dim,
2223 fe_degree,
2224 n_q_points_1d_actual,
2225 Number>::integrate_in_face(n_components,
2226 integration_flag,
2227 shape_data,
2228 temp,
2229 fe_eval.begin_values(),
2230 fe_eval.begin_gradients(),
2231 fe_eval.begin_hessians(),
2232 scratch_data,
2233 subface_index);
2234 }
2235
2236 template <int fe_degree>
2237#ifndef DEBUG
2238 DEAL_II_ALWAYS_INLINE
2239#endif
2240 static void
2241 collect_from_face(const unsigned int n_components,
2242 const EvaluationFlags::EvaluationFlags integration_flag,
2243 Number *values_dofs,
2244 FEEvaluationData<dim, Number, true> &fe_eval,
2245 const bool use_vectorization,
2246 const Number *temp,
2247 Number *scratch_data,
2248 const bool sum_into_values)
2249 {
2250 const auto &shape_info = fe_eval.get_shape_info();
2251 const auto &shape_data = shape_info.data.front();
2252
2253 const unsigned int dofs_per_comp_face =
2254 fe_degree > -1 ?
2255 Utilities::pow(fe_degree + 1, dim - 1) :
2256 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2257 const unsigned int dofs_per_face = n_components * dofs_per_comp_face;
2258
2259 if (use_vectorization == false)
2260 {
2261 for (unsigned int v = 0; v < Number::size(); ++v)
2262 {
2263 // the loop breaks once an invalid_unsigned_int is hit for
2264 // all cases except the exterior faces in the ECL loop (where
2265 // some faces might be at the boundaries but others not)
2266 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2267 continue;
2268
2269 FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
2270 template interpolate<false, false>(n_components,
2271 integration_flag,
2272 shape_info,
2273 temp,
2274 scratch_data,
2275 fe_eval.get_face_no(v));
2276
2277 if (sum_into_values)
2278 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
2279 values_dofs[i][v] += scratch_data[i][v];
2280 else
2281 for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
2282 values_dofs[i][v] = scratch_data[i][v];
2283 }
2284 }
2285 else
2286 {
2287 if (sum_into_values)
2288 FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
2289 template interpolate<false, true>(n_components,
2290 integration_flag,
2291 shape_info,
2292 temp,
2293 values_dofs,
2294 fe_eval.get_face_no());
2295 else
2296 FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
2297 template interpolate<false, false>(n_components,
2298 integration_flag,
2299 shape_info,
2300 temp,
2301 values_dofs,
2302 fe_eval.get_face_no());
2303 }
2304 }
2305
2306 template <int fe_degree, int n_q_points_1d>
2307 static bool
2308 integrate_tensor(const unsigned int n_components,
2309 const EvaluationFlags::EvaluationFlags integration_flag,
2310 Number *values_dofs_actual,
2311 FEEvaluationData<dim, Number, true> &fe_eval,
2312 const bool sum_into_values)
2313 {
2314 const auto &shape_info = fe_eval.get_shape_info();
2315 const auto &shape_data = shape_info.data.front();
2316
2317 const unsigned int dofs_per_comp_face =
2318 fe_degree > -1 ?
2319 Utilities::pow(fe_degree + 1, dim - 1) :
2320 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2321
2322 Number *temp1 = fe_eval.get_scratch_data().begin();
2323 Number *temp2 = temp1 + 3 * n_components * dofs_per_comp_face;
2324
2325 // expand dof_values to tensor product for truncated tensor products
2326 Number *values_dofs =
2327 (shape_data.element_type == MatrixFreeFunctions::truncated_tensor) ?
2328 temp2 + 2 * (std::max<std::size_t>(
2329 fe_eval.get_shape_info().dofs_per_component_on_cell,
2330 fe_eval.get_shape_info().n_q_points)) :
2331 values_dofs_actual;
2332
2333 bool use_vectorization = true;
2334
2335 if (fe_eval.get_dof_access_index() ==
2336 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
2337 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
2338 use_vectorization =
2339 fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
2340 std::all_of(fe_eval.get_cell_ids().begin() + 1,
2341 fe_eval.get_cell_ids().end(),
2342 [&](const auto &v) {
2343 return v == fe_eval.get_cell_ids()[0] ||
2344 v == numbers::invalid_unsigned_int;
2345 });
2346
2347 if (dim == 3)
2348 adjust_quadrature_for_face_orientation(
2349 n_components, integration_flag, fe_eval, use_vectorization, temp1);
2350
2351 integrate_in_face<fe_degree, n_q_points_1d>(
2352 n_components, integration_flag, fe_eval, temp1, temp2);
2353
2354 collect_from_face<fe_degree>(n_components,
2355 integration_flag,
2356 values_dofs,
2357 fe_eval,
2358 use_vectorization,
2359 temp1,
2360 temp2,
2361 sum_into_values);
2362
2363
2364 if (shape_data.element_type == MatrixFreeFunctions::truncated_tensor)
2365 truncate_tensor_product_to_complete_degrees<dim, fe_degree>(
2366 n_components, values_dofs_actual, values_dofs, fe_eval);
2367
2368 return false;
2369 }
2370
2371 template <int fe_degree, int n_q_points_1d>
2372 static bool
2373 run(const unsigned int n_components,
2374 const EvaluationFlags::EvaluationFlags integration_flag,
2375 Number *values_dofs,
2376 FEEvaluationData<dim, Number, true> &fe_eval,
2377 const bool sum_into_values)
2378 {
2379 const auto &shape_info = fe_eval.get_shape_info();
2380
2381 if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
2382 return integrate_tensor_none(n_components,
2383 integration_flag,
2384 values_dofs,
2385 fe_eval,
2386 sum_into_values);
2387 else
2388 return integrate_tensor<fe_degree, n_q_points_1d>(n_components,
2389 integration_flag,
2390 values_dofs,
2391 fe_eval,
2392 sum_into_values);
2393 }
2394 };
2395
2396
2397
2398 template <int dim, typename Number>
2399 struct FEFaceEvaluationImplCollectFromFaceSelector
2400 {
2401 template <int fe_degree>
2402 static bool
2403 run(const unsigned int n_components,
2404 const EvaluationFlags::EvaluationFlags integration_flag,
2405 Number *values_dofs,
2406 FEEvaluationData<dim, Number, true> &fe_eval,
2407 const bool sum_into_values)
2408 {
2409 const auto &shape_info = fe_eval.get_shape_info();
2410 const auto &shape_data = shape_info.data.front();
2411
2412 const unsigned int dofs_per_comp_face =
2413 fe_degree > -1 ?
2414 Utilities::pow(fe_degree + 1, dim - 1) :
2415 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2416
2417 Number *temp = fe_eval.get_scratch_data().begin();
2418 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
2419
2420 bool use_vectorization = true;
2421
2422 if (fe_eval.get_dof_access_index() ==
2423 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
2424 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
2425 use_vectorization =
2426 fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
2427 std::all_of(fe_eval.get_cell_ids().begin() + 1,
2428 fe_eval.get_cell_ids().end(),
2429 [&](const auto &v) {
2430 return v == fe_eval.get_cell_ids()[0] ||
2431 v == numbers::invalid_unsigned_int;
2432 });
2433
2434 FEFaceEvaluationImplIntegrateSelector<dim, Number>::
2435 template collect_from_face<fe_degree>(n_components,
2436 integration_flag,
2437 values_dofs,
2438 fe_eval,
2439 use_vectorization,
2440 temp,
2441 scratch_data,
2442 sum_into_values);
2443
2444 return false;
2445 }
2446 };
2447
2448
2449
2450 template <int dim, typename Number>
2451 struct FEFaceEvaluationImplIntegrateInFaceSelector
2452 {
2453 template <int fe_degree, int n_q_points_1d>
2454 static bool
2455 run(const unsigned int n_components,
2456 const EvaluationFlags::EvaluationFlags integration_flag,
2457
2458 FEEvaluationData<dim, Number, true> &fe_eval)
2459 {
2460 const auto &shape_info = fe_eval.get_shape_info();
2461 const auto &shape_data = shape_info.data.front();
2462
2463 const unsigned int dofs_per_comp_face =
2464 fe_degree > -1 ?
2465 Utilities::pow(fe_degree + 1, dim - 1) :
2466 Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
2467
2468 Number *temp = fe_eval.get_scratch_data().begin();
2469 Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
2470
2471 FEFaceEvaluationImplIntegrateSelector<dim, Number>::
2472 template integrate_in_face<fe_degree, n_q_points_1d>(
2473 n_components, integration_flag, fe_eval, temp, scratch_data);
2474
2475 return false;
2476 }
2477 };
2478
2479
2480
2481 template <int n_face_orientations,
2482 typename Processor,
2483 typename EvaluationData,
2484 const bool check_face_orientations = false>
2485 void
2486 fe_face_evaluation_process_and_io(
2487 Processor &proc,
2488 const unsigned int n_components,
2489 const EvaluationFlags::EvaluationFlags evaluation_flag,
2490 typename Processor::Number2_ *global_vector_ptr,
2491 const std::vector<ArrayView<const typename Processor::Number2_>> *sm_ptr,
2492 const EvaluationData &fe_eval,
2493 typename Processor::VectorizedArrayType_ *temp1)
2494 {
2495 constexpr int dim = Processor::dim_;
2496 constexpr int fe_degree = Processor::fe_degree_;
2497 using VectorizedArrayType = typename Processor::VectorizedArrayType_;
2498 constexpr int n_lanes = VectorizedArrayType::size();
2499
2500 using Number = typename Processor::Number_;
2501 using Number2_ = typename Processor::Number2_;
2502
2503 const auto &shape_data = fe_eval.get_shape_info().data.front();
2504 constexpr bool integrate = Processor::do_integrate;
2505 const unsigned int face_no = fe_eval.get_face_no();
2506 const auto &dof_info = fe_eval.get_dof_info();
2507 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
2508 const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index =
2509 fe_eval.get_dof_access_index();
2510 AssertIndexRange(cell,
2511 dof_info.index_storage_variants[dof_access_index].size());
2512 constexpr unsigned int dofs_per_face =
2513 Utilities::pow(fe_degree + 1, dim - 1);
2514 const unsigned int subface_index = fe_eval.get_subface_index();
2515
2516 const unsigned int n_filled_lanes =
2517 dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
2518
2519 bool all_faces_are_same = n_filled_lanes == n_lanes;
2520 if (n_face_orientations == n_lanes)
2521 for (unsigned int v = 1; v < n_lanes; ++v)
2522 if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) ||
2523 fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0))
2524 {
2525 all_faces_are_same = false;
2526 break;
2527 }
2528
2529 // check for re-orientation ...
2530 std::array<const unsigned int *, n_face_orientations> orientation = {};
2531
2532 if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same &&
2533 fe_eval.is_interior_face() == 0)
2534 for (unsigned int v = 0; v < n_lanes; ++v)
2535 {
2536 // the loop breaks once an invalid_unsigned_int is hit for
2537 // all cases except the exterior faces in the ECL loop (where
2538 // some faces might be at the boundaries but others not)
2539 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2540 continue;
2541
2542 if (shape_data.nodal_at_cell_boundaries &&
2543 fe_eval.get_face_orientation(v) != 0)
2544 {
2545 // ... and in case we detect a re-orientation, go to the other
2546 // version of this function that actually allows for this
2547 if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
2548 check_face_orientations == false)
2549 {
2550 fe_face_evaluation_process_and_io<n_face_orientations,
2551 Processor,
2552 EvaluationData,
2553 true>(proc,
2554 n_components,
2555 evaluation_flag,
2556 global_vector_ptr,
2557 sm_ptr,
2558 fe_eval,
2559 temp1);
2560 return;
2561 }
2562 orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
2563 fe_eval.get_face_orientation(v), 0);
2564 }
2565 }
2566 else if (dim == 3 && fe_eval.get_face_orientation() != 0)
2567 {
2568 // go to the other version of this function
2569 if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
2570 check_face_orientations == false)
2571 {
2572 fe_face_evaluation_process_and_io<n_face_orientations,
2573 Processor,
2574 EvaluationData,
2575 true>(proc,
2576 n_components,
2577 evaluation_flag,
2578 global_vector_ptr,
2579 sm_ptr,
2580 fe_eval,
2581 temp1);
2582 return;
2583 }
2584 for (unsigned int v = 0; v < n_face_orientations; ++v)
2585 orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
2586 fe_eval.get_face_orientation(), 0);
2587 }
2588
2589 // we know that the gradient weights for the Hermite case on the
2590 // right (side==1) are the negative from the value at the left
2591 // (side==0), so we only read out one of them.
2592 VectorizedArrayType grad_weight =
2593 shape_data
2594 .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) :
2595 (1 + face_no % 2))];
2596
2597 // face_to_cell_index_hermite
2598 std::array<const unsigned int *, n_face_orientations> index_array_hermite =
2599 {};
2600 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2601 {
2602 if (n_face_orientations == 1)
2603 index_array_hermite[0] =
2604 &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0);
2605 else
2606 {
2607 for (unsigned int v = 0; v < n_lanes; ++v)
2608 {
2609 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2610 continue;
2611
2612 const auto face_no = fe_eval.get_face_no(v);
2613
2614 grad_weight[v] =
2615 shape_data
2616 .shape_data_on_face[0][fe_degree + (integrate ?
2617 (2 - (face_no % 2)) :
2618 (1 + (face_no % 2)))];
2619
2620 index_array_hermite[v] =
2621 &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no,
2622 0);
2623 }
2624 }
2625 }
2626
2627 // face_to_cell_index_nodal
2628 std::array<const unsigned int *, n_face_orientations> index_array_nodal =
2629 {};
2630 if (shape_data.nodal_at_cell_boundaries == true)
2631 {
2632 if (n_face_orientations == 1)
2633 index_array_nodal[0] =
2634 &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0);
2635 else
2636 {
2637 for (unsigned int v = 0; v < n_lanes; ++v)
2638 {
2639 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
2640 continue;
2641
2642 const auto face_no = fe_eval.get_face_no(v);
2643
2644 index_array_nodal[v] =
2645 &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no,
2646 0);
2647 }
2648 }
2649 }
2650
2651
2652 const auto reorientate = [&](const unsigned int v, const unsigned int i) {
2653 return (!check_face_orientations || orientation[v] == nullptr) ?
2654 i :
2655 orientation[v][i];
2656 };
2657
2658 const unsigned int cell_index =
2659 dof_access_index == MatrixFreeFunctions::DoFInfo::dof_access_cell ?
2660 fe_eval.get_cell_ids()[0] :
2661 cell * n_lanes;
2662 const unsigned int *dof_indices =
2663 &dof_info.dof_indices_contiguous[dof_access_index][cell_index];
2664
2665 for (unsigned int comp = 0; comp < n_components; ++comp)
2666 {
2667 const std::size_t index_offset =
2668 dof_info.component_dof_indices_offset
2669 [fe_eval.get_active_fe_index()]
2670 [fe_eval.get_first_selected_component()] +
2671 comp * Utilities::pow(fe_degree + 1, dim);
2672
2673 // case 1: contiguous and interleaved indices
2674 if (n_face_orientations == 1 &&
2675 dof_info.index_storage_variants[dof_access_index][cell] ==
2676 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2677 interleaved_contiguous)
2678 {
2679 AssertDimension(
2680 dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2681 n_lanes);
2682 Number2_ *vector_ptr =
2683 global_vector_ptr + dof_indices[0] + index_offset * n_lanes;
2684
2685 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2686 {
2687 for (unsigned int i = 0; i < dofs_per_face; ++i)
2688 {
2689 Assert(n_face_orientations == 1, ExcNotImplemented());
2690
2691 const unsigned int ind1 = index_array_hermite[0][2 * i];
2692 const unsigned int ind2 = index_array_hermite[0][2 * i + 1];
2693 const unsigned int i_ = reorientate(0, i);
2694 proc.hermite_grad_vectorized(temp1[i_],
2695 temp1[i_ + dofs_per_face],
2696 vector_ptr + ind1 * n_lanes,
2697 vector_ptr + ind2 * n_lanes,
2698 grad_weight);
2699 }
2700 }
2701 else
2702 {
2703 for (unsigned int i = 0; i < dofs_per_face; ++i)
2704 {
2705 Assert(n_face_orientations == 1, ExcNotImplemented());
2706
2707 const unsigned int i_ = reorientate(0, i);
2708 const unsigned int ind = index_array_nodal[0][i];
2709 proc.value_vectorized(temp1[i_],
2710 vector_ptr + ind * n_lanes);
2711 }
2712 }
2713 }
2714
2715 // case 2: contiguous and interleaved indices with fixed stride
2716 else if (n_face_orientations == 1 &&
2717 dof_info.index_storage_variants[dof_access_index][cell] ==
2718 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2719 interleaved_contiguous_strided)
2720 {
2721 AssertDimension(
2722 dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
2723 n_lanes);
2724 Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes;
2725 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2726 {
2727 for (unsigned int i = 0; i < dofs_per_face; ++i)
2728 {
2729 Assert(n_face_orientations == 1, ExcNotImplemented());
2730
2731 const unsigned int i_ = reorientate(0, i);
2732 const unsigned int ind1 =
2733 index_array_hermite[0][2 * i] * n_lanes;
2734 const unsigned int ind2 =
2735 index_array_hermite[0][2 * i + 1] * n_lanes;
2736 proc.hermite_grad_vectorized_indexed(
2737 temp1[i_],
2738 temp1[i_ + dofs_per_face],
2739 vector_ptr + ind1,
2740 vector_ptr + ind2,
2741 grad_weight,
2742 dof_indices,
2743 dof_indices);
2744 }
2745 }
2746 else
2747 {
2748 for (unsigned int i = 0; i < dofs_per_face; ++i)
2749 {
2750 Assert(n_face_orientations == 1, ExcNotImplemented());
2751
2752 const unsigned int i_ = reorientate(0, i);
2753 const unsigned int ind = index_array_nodal[0][i] * n_lanes;
2754 proc.value_vectorized_indexed(temp1[i_],
2755 vector_ptr + ind,
2756 dof_indices);
2757 }
2758 }
2759 }
2760
2761 // case 3: contiguous and interleaved indices with mixed stride
2762 else if (n_face_orientations == 1 &&
2763 dof_info.index_storage_variants[dof_access_index][cell] ==
2764 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2765 interleaved_contiguous_mixed_strides)
2766 {
2767 const unsigned int *strides =
2768 &dof_info.dof_indices_interleave_strides[dof_access_index]
2769 [cell * n_lanes];
2770 unsigned int indices[n_lanes];
2771 for (unsigned int v = 0; v < n_lanes; ++v)
2772 indices[v] = dof_indices[v] + index_offset * strides[v];
2773 const unsigned int n_filled_lanes =
2774 dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
2775
2776 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2777 {
2778 if (n_filled_lanes == n_lanes)
2779 for (unsigned int i = 0; i < dofs_per_face; ++i)
2780 {
2781 Assert(n_face_orientations == 1, ExcNotImplemented());
2782
2783 const unsigned int i_ = reorientate(0, i);
2784 unsigned int ind1[n_lanes];
2785 DEAL_II_OPENMP_SIMD_PRAGMA
2786 for (unsigned int v = 0; v < n_lanes; ++v)
2787 ind1[v] = indices[v] +
2788 index_array_hermite[0][2 * i] * strides[v];
2789 unsigned int ind2[n_lanes];
2790 DEAL_II_OPENMP_SIMD_PRAGMA
2791 for (unsigned int v = 0; v < n_lanes; ++v)
2792 ind2[v] =
2793 indices[v] +
2794 // TODO
2795 index_array_hermite[0][2 * i + 1] * strides[v];
2796 proc.hermite_grad_vectorized_indexed(
2797 temp1[i_],
2798 temp1[i_ + dofs_per_face],
2799 global_vector_ptr,
2800 global_vector_ptr,
2801 grad_weight,
2802 ind1,
2803 ind2);
2804 }
2805 else
2806 {
2807 if (integrate == false)
2808 for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
2809 temp1[i] = VectorizedArrayType();
2810
2811 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2812 for (unsigned int i = 0; i < dofs_per_face; ++i)
2813 {
2814 const unsigned int i_ =
2815 reorientate(n_face_orientations == 1 ? 0 : v, i);
2816 proc.hermite_grad(
2817 temp1[i_][v],
2818 temp1[i_ + dofs_per_face][v],
2819 global_vector_ptr
2820 [indices[v] +
2821 index_array_hermite
2822 [n_face_orientations == 1 ? 0 : v][2 * i] *
2823 strides[v]],
2824 global_vector_ptr
2825 [indices[v] +
2826 index_array_hermite[n_face_orientations == 1 ?
2827 0 :
2828 v][2 * i + 1] *
2829 strides[v]],
2830 grad_weight[n_face_orientations == 1 ? 0 : v]);
2831 }
2832 }
2833 }
2834 else
2835 {
2836 if (n_filled_lanes == n_lanes)
2837 for (unsigned int i = 0; i < dofs_per_face; ++i)
2838 {
2839 Assert(n_face_orientations == 1, ExcInternalError());
2840 unsigned int ind[n_lanes];
2841 DEAL_II_OPENMP_SIMD_PRAGMA
2842 for (unsigned int v = 0; v < n_lanes; ++v)
2843 ind[v] =
2844 indices[v] + index_array_nodal[0][i] * strides[v];
2845 const unsigned int i_ = reorientate(0, i);
2846 proc.value_vectorized_indexed(temp1[i_],
2847 global_vector_ptr,
2848 ind);
2849 }
2850 else
2851 {
2852 if (integrate == false)
2853 for (unsigned int i = 0; i < dofs_per_face; ++i)
2854 temp1[i] = VectorizedArrayType();
2855
2856 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2857 for (unsigned int i = 0; i < dofs_per_face; ++i)
2858 proc.value(
2859 temp1[reorientate(n_face_orientations == 1 ? 0 : v,
2860 i)][v],
2861 global_vector_ptr
2862 [indices[v] +
2863 index_array_nodal[n_face_orientations == 1 ? 0 : v]
2864 [i] *
2865 strides[v]]);
2866 }
2867 }
2868 }
2869
2870 // case 4: contiguous indices without interleaving
2871 else if (n_face_orientations > 1 ||
2872 dof_info.index_storage_variants[dof_access_index][cell] ==
2873 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
2874 contiguous)
2875 {
2876 Number2_ *vector_ptr = global_vector_ptr + index_offset;
2877
2878 const bool vectorization_possible =
2879 all_faces_are_same && (sm_ptr == nullptr);
2880
2881 std::array<Number2_ *, n_lanes> vector_ptrs{{nullptr}};
2882 std::array<unsigned int, n_lanes> reordered_indices{
2883 {numbers::invalid_unsigned_int}};
2884
2885 if (vectorization_possible == false)
2886 {
2887 if (n_face_orientations == 1)
2888 {
2889 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2890 if (sm_ptr == nullptr)
2891 {
2892 vector_ptrs[v] = vector_ptr + dof_indices[v];
2893 }
2894 else
2895 {
2896 const auto &temp =
2897 dof_info
2898 .dof_indices_contiguous_sm[dof_access_index]
2899 [cell * n_lanes + v];
2900 vector_ptrs[v] = const_cast<Number2_ *>(
2901 sm_ptr->operator[](temp.first).data() +
2902 temp.second + index_offset);
2903 }
2904 }
2905 else if (n_face_orientations == n_lanes)
2906 {
2907 const auto &cells = fe_eval.get_cell_ids();
2908 for (unsigned int v = 0; v < n_lanes; ++v)
2909 if (cells[v] != numbers::invalid_unsigned_int)
2910 {
2911 if (sm_ptr == nullptr)
2912 {
2913 vector_ptrs[v] =
2914 vector_ptr +
2915 dof_info
2916 .dof_indices_contiguous[dof_access_index]
2917 [cells[v]];
2918 }
2919 else
2920 {
2921 const auto &temp =
2922 dof_info
2923 .dof_indices_contiguous_sm[dof_access_index]
2924 [cells[v]];
2925 vector_ptrs[v] = const_cast<Number2_ *>(
2926 sm_ptr->operator[](temp.first).data() +
2927 temp.second + index_offset);
2928 }
2929 }
2930 }
2931 else
2932 {
2933 DEAL_II_NOT_IMPLEMENTED();
2934 }
2935 }
2936 else if (n_face_orientations == n_lanes)
2937 {
2938 for (unsigned int v = 0; v < n_lanes; ++v)
2939 reordered_indices[v] =
2940 dof_info.dof_indices_contiguous[dof_access_index]
2941 [fe_eval.get_cell_ids()[v]];
2942 dof_indices = reordered_indices.data();
2943 }
2944
2945 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
2946 {
2947 if (vectorization_possible)
2948 for (unsigned int i = 0; i < dofs_per_face; ++i)
2949 {
2950 const unsigned int ind1 = index_array_hermite[0][2 * i];
2951 const unsigned int ind2 =
2952 index_array_hermite[0][2 * i + 1];
2953 const unsigned int i_ = reorientate(0, i);
2954
2955 proc.hermite_grad_vectorized_indexed(
2956 temp1[i_],
2957 temp1[i_ + dofs_per_face],
2958 vector_ptr + ind1,
2959 vector_ptr + ind2,
2960 grad_weight,
2961 dof_indices,
2962 dof_indices);
2963 }
2964 else if (n_face_orientations == 1)
2965 for (unsigned int i = 0; i < dofs_per_face; ++i)
2966 {
2967 const unsigned int ind1 = index_array_hermite[0][2 * i];
2968 const unsigned int ind2 =
2969 index_array_hermite[0][2 * i + 1];
2970 const unsigned int i_ = reorientate(0, i);
2971
2972 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2973 proc.hermite_grad(temp1[i_][v],
2974 temp1[i_ + dofs_per_face][v],
2975 vector_ptrs[v][ind1],
2976 vector_ptrs[v][ind2],
2977 grad_weight[v]);
2978
2979 if (integrate == false)
2980 for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
2981 {
2982 temp1[i][v] = 0.0;
2983 temp1[i + dofs_per_face][v] = 0.0;
2984 }
2985 }
2986 else
2987 {
2988 if (integrate == false && n_filled_lanes < n_lanes)
2989 for (unsigned int i = 0; i < dofs_per_face; ++i)
2990 temp1[i] = temp1[i + dofs_per_face] = Number();
2991
2992 for (unsigned int v = 0; v < n_filled_lanes; ++v)
2993 for (unsigned int i = 0; i < dofs_per_face; ++i)
2994 proc.hermite_grad(
2995 temp1[reorientate(v, i)][v],
2996 temp1[reorientate(v, i) + dofs_per_face][v],
2997 vector_ptrs[v][index_array_hermite[v][2 * i]],
2998 vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
2999 grad_weight[v]);
3000 }
3001 }
3002 else
3003 {
3004 if (vectorization_possible)
3005 for (unsigned int i = 0; i < dofs_per_face; ++i)
3006 {
3007 const unsigned int ind = index_array_nodal[0][i];
3008 const unsigned int i_ = reorientate(0, i);
3009
3010 proc.value_vectorized_indexed(temp1[i_],
3011 vector_ptr + ind,
3012 dof_indices);
3013 }
3014 else
3015 {
3016 if constexpr (n_face_orientations == 1)
3017 for (unsigned int i = 0; i < dofs_per_face; ++i)
3018 {
3019 const unsigned int ind = index_array_nodal[0][i];
3020 const unsigned int i_ = reorientate(0, i);
3021
3022 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3023 proc.value(temp1[i_][v], vector_ptrs[v][ind]);
3024
3025 if constexpr (integrate == false)
3026 for (unsigned int v = n_filled_lanes; v < n_lanes;
3027 ++v)
3028 temp1[i_][v] = 0.0;
3029 }
3030 else
3031 {
3032 if (integrate == false && n_filled_lanes < n_lanes)
3033 for (unsigned int i = 0; i < dofs_per_face; ++i)
3034 temp1[i] = Number();
3035
3036 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3037 for (unsigned int i = 0; i < dofs_per_face; ++i)
3038 proc.value(temp1[reorientate(v, i)][v],
3039 vector_ptrs[v][index_array_nodal[v][i]]);
3040 }
3041 }
3042 }
3043 }
3044 else
3045 {
3046 // We should not end up here, this should be caught by
3047 // FEFaceEvaluationImplGatherEvaluateSelector::supports()
3048 DEAL_II_ASSERT_UNREACHABLE();
3049 }
3050 temp1 += 3 * dofs_per_face;
3051 }
3052 }
3053
3054
3055
3056 template <int dim, typename Number2, typename VectorizedArrayType>
3057 struct FEFaceEvaluationImplGatherEvaluateSelector
3058 {
3059 using Number = typename VectorizedArrayType::value_type;
3060
3061 template <int fe_degree, int n_q_points_1d>
3062 static bool
3063 run(const unsigned int n_components,
3064 const EvaluationFlags::EvaluationFlags evaluation_flag,
3065 const Number2 *src_ptr,
3066 const std::vector<ArrayView<const Number2>> *sm_ptr,
3067 FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
3068 {
3069 Assert(fe_degree > -1, ExcInternalError());
3070 Assert(fe_eval.get_shape_info().element_type <=
3071 MatrixFreeFunctions::tensor_symmetric_no_collocation,
3072 ExcInternalError());
3073
3074 const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
3075
3076 VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
3077 VectorizedArrayType *scratch_data =
3078 temp + 3 * n_components * dofs_per_face;
3079
3080 Processor<fe_degree> p;
3081
3082 if (fe_eval.get_dof_access_index() ==
3083 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
3084 fe_eval.is_interior_face() == false)
3085 fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
3086 p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
3087 else
3088 fe_face_evaluation_process_and_io<1>(
3089 p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
3090
3091 const unsigned int subface_index = fe_eval.get_subface_index();
3092
3093 if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
3094 FEFaceEvaluationImpl<true,
3095 dim,
3096 fe_degree,
3097 n_q_points_1d,
3098 VectorizedArrayType>::
3099 evaluate_in_face(n_components,
3100 evaluation_flag,
3101 fe_eval.get_shape_info().data.front(),
3102 temp,
3103 fe_eval.begin_values(),
3104 fe_eval.begin_gradients(),
3105 fe_eval.begin_hessians(),
3106 scratch_data,
3107 subface_index);
3108 else
3109 FEFaceEvaluationImpl<false,
3110 dim,
3111 fe_degree,
3112 n_q_points_1d,
3113 VectorizedArrayType>::
3114 evaluate_in_face(n_components,
3115 evaluation_flag,
3116 fe_eval.get_shape_info().data.front(),
3117 temp,
3118 fe_eval.begin_values(),
3119 fe_eval.begin_gradients(),
3120 fe_eval.begin_hessians(),
3121 scratch_data,
3122 subface_index);
3123
3124 // re-orientation for cases not possible with above algorithm
3125 if (subface_index < GeometryInfo<dim>::max_children_per_cell)
3126 {
3127 if (fe_eval.get_dof_access_index() ==
3128 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
3129 fe_eval.is_interior_face() == false)
3130 {
3131 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3132 {
3133 // the loop breaks once an invalid_unsigned_int is hit for
3134 // all cases except the exterior faces in the ECL loop (where
3135 // some faces might be at the boundaries but others not)
3136 if (fe_eval.get_cell_ids()[v] ==
3137 numbers::invalid_unsigned_int)
3138 continue;
3139
3140 if (fe_eval.get_face_orientation(v) != 0)
3141 adjust_for_face_orientation_per_lane(
3142 dim,
3143 n_components,
3144 v,
3145 evaluation_flag,
3146 &fe_eval.get_shape_info().face_orientations_quad(
3147 fe_eval.get_face_orientation(v), 0),
3148 false,
3149 Utilities::pow(n_q_points_1d, dim - 1),
3150 &temp[0][0],
3151 fe_eval.begin_values(),
3152 fe_eval.begin_gradients(),
3153 fe_eval.begin_hessians());
3154 }
3155 }
3156 else if (fe_eval.get_face_orientation() != 0)
3157 adjust_for_face_orientation(
3158 dim,
3159 n_components,
3160 evaluation_flag,
3161 &fe_eval.get_shape_info().face_orientations_quad(
3162 fe_eval.get_face_orientation(), 0),
3163 false,
3164 Utilities::pow(n_q_points_1d, dim - 1),
3165 temp,
3166 fe_eval.begin_values(),
3167 fe_eval.begin_gradients(),
3168 fe_eval.begin_hessians());
3169 }
3170
3171 return false;
3172 }
3173
3174 template <typename Number3>
3175 static bool
3176 supports(const EvaluationFlags::EvaluationFlags evaluation_flag,
3177 const MatrixFreeFunctions::ShapeInfo<Number3> &shape_info,
3178 const Number2 *vector_ptr,
3179 MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage)
3180 {
3181 const unsigned int fe_degree = shape_info.data.front().fe_degree;
3182 if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries ||
3183 (evaluation_flag & EvaluationFlags::gradients &&
3184 (fe_degree < 2 ||
3185 shape_info.data.front().element_type !=
3186 MatrixFreeFunctions::tensor_symmetric_hermite)) ||
3187 (evaluation_flag & EvaluationFlags::hessians) ||
3188 vector_ptr == nullptr ||
3189 shape_info.data.front().element_type >
3190 MatrixFreeFunctions::tensor_symmetric_no_collocation ||
3191 storage <
3192 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous)
3193 return false;
3194 else
3195 return true;
3196 }
3197
3198 private:
3199 template <int fe_degree>
3200 struct Processor
3201 {
3202 static const bool do_integrate = false;
3203 static const int dim_ = dim;
3204 static const int fe_degree_ = fe_degree;
3205 using VectorizedArrayType_ = VectorizedArrayType;
3206 using Number_ = Number;
3207 using Number2_ = const Number2;
3208
3209 template <typename T0, typename T1, typename T2>
3210 void
3211 hermite_grad_vectorized(T0 &temp_1,
3212 T0 &temp_2,
3213 const T1 src_ptr_1,
3214 const T1 src_ptr_2,
3215 const T2 &grad_weight)
3216 {
3217 do_vectorized_read(src_ptr_1, temp_1);
3218 do_vectorized_read(src_ptr_2, temp_2);
3219 temp_2 = grad_weight * (temp_1 - temp_2);
3220 }
3221
3222 template <typename T1, typename T2>
3223 void
3224 value_vectorized(T1 &temp, const T2 src_ptr)
3225 {
3226 do_vectorized_read(src_ptr, temp);
3227 }
3228
3229 template <typename T0, typename T1, typename T2, typename T3>
3230 void
3231 hermite_grad_vectorized_indexed(T0 &temp_1,
3232 T0 &temp_2,
3233 const T1 src_ptr_1,
3234 const T1 src_ptr_2,
3235 const T2 &grad_weight,
3236 const T3 &indices_1,
3237 const T3 &indices_2)
3238 {
3239 do_vectorized_gather(src_ptr_1, indices_1, temp_1);
3240 do_vectorized_gather(src_ptr_2, indices_2, temp_2);
3241 temp_2 = grad_weight * (temp_1 - temp_2);
3242 }
3243
3244 template <typename T0, typename T1, typename T2>
3245 void
3246 value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
3247 {
3248 do_vectorized_gather(src_ptr, indices, temp);
3249 }
3250
3251 template <typename T0, typename T1, typename T2>
3252 void
3253 hermite_grad(T0 &temp_1,
3254 T0 &temp_2,
3255 const T1 &src_ptr_1,
3256 const T1 &src_ptr_2,
3257 const T2 &grad_weight)
3258 {
3259 // case 3a)
3260 temp_1 = src_ptr_1;
3261 temp_2 = grad_weight * (temp_1 - src_ptr_2);
3262 }
3263
3264 template <typename T1, typename T2>
3265 void
3266 value(T1 &temp, const T2 &src_ptr)
3267 {
3268 // case 3b)
3269 temp = src_ptr;
3270 }
3271 };
3272 };
3273
3274
3275
3276 template <int dim, typename Number2, typename VectorizedArrayType>
3277 struct FEFaceEvaluationImplIntegrateScatterSelector
3278 {
3279 using Number = typename VectorizedArrayType::value_type;
3280
3281 template <int fe_degree, int n_q_points_1d>
3282 static bool
3283 run(const unsigned int n_components,
3284 const EvaluationFlags::EvaluationFlags integration_flag,
3285 Number2 *dst_ptr,
3286 const std::vector<ArrayView<const Number2>> *sm_ptr,
3287 FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
3288 {
3289 Assert(fe_degree > -1, ExcInternalError());
3290 Assert(fe_eval.get_shape_info().element_type <=
3291 MatrixFreeFunctions::tensor_symmetric_no_collocation,
3292 ExcInternalError());
3293
3294 const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
3295
3296 VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
3297 VectorizedArrayType *scratch_data =
3298 temp + 3 * n_components * dofs_per_face;
3299
3300 const unsigned int subface_index = fe_eval.get_subface_index();
3301
3302 // re-orientation for cases not possible with the io function below
3303 if (subface_index < GeometryInfo<dim>::max_children_per_cell)
3304 {
3305 if (fe_eval.get_dof_access_index() ==
3306 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
3307 fe_eval.is_interior_face() == false)
3308 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3309 {
3310 // the loop breaks once an invalid_unsigned_int is hit for
3311 // all cases except the exterior faces in the ECL loop (where
3312 // some faces might be at the boundaries but others not)
3313 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
3314 continue;
3315
3316 if (fe_eval.get_face_orientation(v) != 0)
3317 adjust_for_face_orientation_per_lane(
3318 dim,
3319 n_components,
3320 v,
3321 integration_flag,
3322 &fe_eval.get_shape_info().face_orientations_quad(
3323 fe_eval.get_face_orientation(v), 0),
3324 true,
3325 Utilities::pow(n_q_points_1d, dim - 1),
3326 &temp[0][0],
3327 fe_eval.begin_values(),
3328 fe_eval.begin_gradients(),
3329 fe_eval.begin_hessians());
3330 }
3331 else if (fe_eval.get_face_orientation() != 0)
3332 adjust_for_face_orientation(
3333 dim,
3334 n_components,
3335 integration_flag,
3336 &fe_eval.get_shape_info().face_orientations_quad(
3337 fe_eval.get_face_orientation(), 0),
3338 true,
3339 Utilities::pow(n_q_points_1d, dim - 1),
3340 temp,
3341 fe_eval.begin_values(),
3342 fe_eval.begin_gradients(),
3343 fe_eval.begin_hessians());
3344 }
3345
3346 if (fe_degree > -1 && fe_eval.get_subface_index() >=
3347 GeometryInfo<dim - 1>::max_children_per_cell)
3348 FEFaceEvaluationImpl<true,
3349 dim,
3350 fe_degree,
3351 n_q_points_1d,
3352 VectorizedArrayType>::
3353 integrate_in_face(n_components,
3354 integration_flag,
3355 fe_eval.get_shape_info().data.front(),
3356 temp,
3357 fe_eval.begin_values(),
3358 fe_eval.begin_gradients(),
3359 fe_eval.begin_hessians(),
3360 scratch_data,
3361 subface_index);
3362 else
3363 FEFaceEvaluationImpl<false,
3364 dim,
3365 fe_degree,
3366 n_q_points_1d,
3367 VectorizedArrayType>::
3368 integrate_in_face(n_components,
3369 integration_flag,
3370 fe_eval.get_shape_info().data.front(),
3371 temp,
3372 fe_eval.begin_values(),
3373 fe_eval.begin_gradients(),
3374 fe_eval.begin_hessians(),
3375 scratch_data,
3376 subface_index);
3377
3378 Processor<fe_degree> p;
3379
3380 if (fe_eval.get_dof_access_index() ==
3381 MatrixFreeFunctions::DoFInfo::dof_access_cell &&
3382 fe_eval.is_interior_face() == false)
3383 fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
3384 p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
3385 else
3386 fe_face_evaluation_process_and_io<1>(
3387 p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
3388
3389 return false;
3390 }
3391
3392 private:
3393 template <int fe_degree>
3394 struct Processor
3395 {
3396 static const bool do_integrate = true;
3397 static const int dim_ = dim;
3398 static const int fe_degree_ = fe_degree;
3399 using VectorizedArrayType_ = VectorizedArrayType;
3400 using Number_ = Number;
3401 using Number2_ = Number2;
3402
3403 template <typename T0, typename T1, typename T2, typename T3, typename T4>
3404 void
3405 hermite_grad_vectorized(const T0 &temp_1,
3406 const T1 &temp_2,
3407 T2 dst_ptr_1,
3408 T3 dst_ptr_2,
3409 const T4 &grad_weight)
3410 {
3411 // case 1a)
3412 const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
3413 const VectorizedArrayType grad = grad_weight * temp_2;
3414 do_vectorized_add(val, dst_ptr_1);
3415 do_vectorized_add(grad, dst_ptr_2);
3416 }
3417
3418 template <typename T0, typename T1>
3419 void
3420 value_vectorized(const T0 &temp, T1 dst_ptr)
3421 {
3422 // case 1b)
3423 do_vectorized_add(temp, dst_ptr);
3424 }
3425
3426 template <typename T0, typename T1, typename T2, typename T3>
3427 void
3428 hermite_grad_vectorized_indexed(const T0 &temp_1,
3429 const T0 &temp_2,
3430 T1 dst_ptr_1,
3431 T1 dst_ptr_2,
3432 const T2 &grad_weight,
3433 const T3 &indices_1,
3434 const T3 &indices_2)
3435 {
3436 // case 2a)
3437 const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
3438 const VectorizedArrayType grad = grad_weight * temp_2;
3439 do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
3440 do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
3441 }
3442
3443 template <typename T0, typename T1, typename T2>
3444 void
3445 value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
3446 {
3447 // case 2b)
3448 do_vectorized_scatter_add(temp, indices, dst_ptr);
3449 }
3450
3451 template <typename T0, typename T1, typename T2>
3452 void
3453 hermite_grad(const T0 &temp_1,
3454 const T0 &temp_2,
3455 T1 &dst_ptr_1,
3456 T1 &dst_ptr_2,
3457 const T2 &grad_weight)
3458 {
3459 // case 3a)
3460 const Number val = temp_1 - grad_weight * temp_2;
3461 const Number grad = grad_weight * temp_2;
3462 dst_ptr_1 += val;
3463 dst_ptr_2 += grad;
3464 }
3465
3466 template <typename T0, typename T1>
3467 void
3468 value(const T0 &temp, T1 &dst_ptr)
3469 {
3470 // case 3b)
3471 dst_ptr += temp;
3472 }
3473 };
3474 };
3475} // end of namespace internal
3476
3477
3478DEAL_II_NAMESPACE_CLOSE
3479
3480#endif
i_
*  unsigned int i_
Definition ALE_Finite_Strain_Plasticity.h:10366
ArrayView::begin
iterator begin() const
Definition array_view.h:755
FEEvaluationData::get_face_no
std::uint8_t get_face_no(const unsigned int v=0) const
FEEvaluationData::get_dof_access_index
internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex get_dof_access_index() const
FEEvaluationData::shape_info_number_type
ScalarNumber shape_info_number_type
Definition fe_evaluation_data.h:127
FEEvaluationData::get_shape_info
const ShapeInfoType & get_shape_info() const
FEEvaluationData::get_cell_ids
const std::array< unsigned int, n_lanes > & get_cell_ids() const
Definition fe_evaluation_data.h:496
FEEvaluationData::begin_gradients
const Number * begin_gradients() const
FEEvaluationData::get_subface_index
unsigned int get_subface_index() const
FEEvaluationData::is_interior_face
bool is_interior_face() const
FEEvaluationData::get_scratch_data
ArrayView< Number > get_scratch_data() const
FEEvaluationData::begin_values
const Number * begin_values() const
FEEvaluationData::get_face_orientation
std::uint8_t get_face_orientation(const unsigned int v=0) const
FEEvaluationData::begin_hessians
const Number * begin_hessians() const
Function::value
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
VectorizedArray::gather
void gather(const Number *base_ptr, const unsigned int *offsets)
Definition vectorization.h:667
VectorizedArray::load
void load(const OtherNumber *ptr)
Definition vectorization.h:579
DEAL_II_ALWAYS_INLINE
#define DEAL_II_ALWAYS_INLINE
Definition config.h:166
DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:214
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
evaluation_flags.h
evaluation_kernels_common.h
DEAL_II_ASSERT_UNREACHABLE
#define DEAL_II_ASSERT_UNREACHABLE()
Definition exception_macros.h:781
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exception_macros.h:705
DoFCouplingType::contiguous
@ contiguous
fe_evaluation_data.h
geometry_info.h
cell_index
unsigned int cell_index
Definition grid_tools_topology.cc:777
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exception_macros.h:829
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exception_macros.h:901
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
AssertThrow
#define AssertThrow(cond, exc)
Definition exception_macros.h:647
internal::MatrixFreeFunctions::tensor_symmetric_no_collocation
@ tensor_symmetric_no_collocation
Definition shape_info.h:85
internal::MatrixFreeFunctions::tensor_raviart_thomas
@ tensor_raviart_thomas
Definition shape_info.h:110
internal::MatrixFreeFunctions::truncated_tensor
@ truncated_tensor
Definition shape_info.h:96
internal::MatrixFreeFunctions::tensor_none
@ tensor_none
Definition shape_info.h:122
internal::MatrixFreeFunctions::tensor_symmetric
@ tensor_symmetric
Definition shape_info.h:79
internal::MatrixFreeFunctions::tensor_symmetric_hermite
@ tensor_symmetric_hermite
Definition shape_info.h:72
dof_info.h
data
std::vector< index_type > data
Definition mpi.cc:734
EvaluationFlags::EvaluationFlags
EvaluationFlags
The EvaluationFlags enum.
Definition evaluation_flags.h:42
PlasticityLab::Number
*  *  *  RotationFunction< dim, Number >::RotationFunction Number(dim)
Definition ALE_Finite_Strain_Plasticity.h:8733
Utilities::fixed_power
constexpr T fixed_power(const T t)
Definition utilities.h:942
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
internal::do_vectorized_add
void do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
Definition evaluation_kernels_face.h:1263
internal::use_collocation_evaluation
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
Definition tensor_product_kernels.h:2451
internal::adjust_for_face_orientation_per_lane
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)
Definition evaluation_kernels_face.h:1382
internal::fe_face_evaluation_process_and_io
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)
Definition evaluation_kernels_face.h:2486
internal::apply_matrix_vector_product
std::enable_if_t<(variant==evaluate_general), void > apply_matrix_vector_product(const Number2 *matrix, const Number *in, Number *out)
Definition tensor_product_kernels.h:116
internal::do_vectorized_scatter_add
void do_vectorized_scatter_add(const VectorizedArrayType src, const unsigned int *indices, Number2 *dst_ptr)
Definition evaluation_kernels_face.h:1287
internal::do_vectorized_gather
void do_vectorized_gather(const Number2 *src_ptr, const unsigned int *indices, VectorizedArrayType &dst)
Definition evaluation_kernels_face.h:1237
internal::do_vectorized_read
void do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
Definition evaluation_kernels_face.h:1215
internal::interpolate_to_face
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)
Definition tensor_product_kernels.h:2383
internal::adjust_for_face_orientation
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)
Definition evaluation_kernels_face.h:1319
numbers::invalid_unsigned_int
constexpr unsigned int invalid_unsigned_int
Definition types.h:228
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6596
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6579
shape_info.h
internal::FEFaceEvaluationImplCollectFromFaceSelector::run
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)
Definition evaluation_kernels_face.h:2403
internal::FEFaceEvaluationImplEvaluateInFaceSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, FEEvaluationData< dim, Number, true > &fe_eval)
Definition evaluation_kernels_face.h:1890
internal::FEFaceEvaluationImplEvaluateSelector::evaluate_in_face
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)
Definition evaluation_kernels_face.h:1635
internal::FEFaceEvaluationImplEvaluateSelector::project_to_face
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)
Definition evaluation_kernels_face.h:1575
internal::FEFaceEvaluationImplEvaluateSelector::evaluate_tensor_none
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)
Definition evaluation_kernels_face.h:1453
internal::FEFaceEvaluationImplEvaluateSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
Definition evaluation_kernels_face.h:1816
internal::FEFaceEvaluationImplEvaluateSelector::adjust_quadrature_for_face_orientation
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)
Definition evaluation_kernels_face.h:1702
internal::FEFaceEvaluationImplEvaluateSelector::evaluate_tensor
static bool evaluate_tensor(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs_actual, FEEvaluationData< dim, Number, true > &fe_eval)
Definition evaluation_kernels_face.h:1755
internal::FEFaceEvaluationImplGatherEvaluateSelector::Processor::hermite_grad
void hermite_grad(T0 &temp_1, T0 &temp_2, const T1 &src_ptr_1, const T1 &src_ptr_2, const T2 &grad_weight)
Definition evaluation_kernels_face.h:3253
internal::FEFaceEvaluationImplGatherEvaluateSelector::Processor::value_vectorized_indexed
void value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
Definition evaluation_kernels_face.h:3246
internal::FEFaceEvaluationImplGatherEvaluateSelector::Processor::hermite_grad_vectorized
void hermite_grad_vectorized(T0 &temp_1, T0 &temp_2, const T1 src_ptr_1, const T1 src_ptr_2, const T2 &grad_weight)
Definition evaluation_kernels_face.h:3211
internal::FEFaceEvaluationImplGatherEvaluateSelector::Processor::hermite_grad_vectorized_indexed
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)
Definition evaluation_kernels_face.h:3231
internal::FEFaceEvaluationImplGatherEvaluateSelector::run
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)
Definition evaluation_kernels_face.h:3063
internal::FEFaceEvaluationImplGatherEvaluateSelector::supports
static bool supports(const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::ShapeInfo< Number3 > &shape_info, const Number2 *vector_ptr, MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage)
Definition evaluation_kernels_face.h:3176
internal::FEFaceEvaluationImplIntegrateInFaceSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, FEEvaluationData< dim, Number, true > &fe_eval)
Definition evaluation_kernels_face.h:2455
internal::FEFaceEvaluationImplIntegrateScatterSelector::Processor::hermite_grad_vectorized
void hermite_grad_vectorized(const T0 &temp_1, const T1 &temp_2, T2 dst_ptr_1, T3 dst_ptr_2, const T4 &grad_weight)
Definition evaluation_kernels_face.h:3405
internal::FEFaceEvaluationImplIntegrateScatterSelector::Processor::hermite_grad_vectorized_indexed
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)
Definition evaluation_kernels_face.h:3428
internal::FEFaceEvaluationImplIntegrateScatterSelector::Processor::hermite_grad
void hermite_grad(const T0 &temp_1, const T0 &temp_2, T1 &dst_ptr_1, T1 &dst_ptr_2, const T2 &grad_weight)
Definition evaluation_kernels_face.h:3453
internal::FEFaceEvaluationImplIntegrateScatterSelector::Processor::value_vectorized_indexed
void value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
Definition evaluation_kernels_face.h:3445
internal::FEFaceEvaluationImplIntegrateScatterSelector::run
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)
Definition evaluation_kernels_face.h:3283
internal::FEFaceEvaluationImplIntegrateSelector::integrate_tensor
static bool integrate_tensor(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs_actual, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
Definition evaluation_kernels_face.h:2308
internal::FEFaceEvaluationImplIntegrateSelector::integrate_tensor_none
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)
Definition evaluation_kernels_face.h:1921
internal::FEFaceEvaluationImplIntegrateSelector::collect_from_face
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)
Definition evaluation_kernels_face.h:2241
internal::FEFaceEvaluationImplIntegrateSelector::run
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)
Definition evaluation_kernels_face.h:2373
internal::FEFaceEvaluationImplIntegrateSelector::integrate_in_face
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)
Definition evaluation_kernels_face.h:2172
internal::FEFaceEvaluationImplIntegrateSelector::adjust_quadrature_for_face_orientation
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)
Definition evaluation_kernels_face.h:2117
internal::FEFaceEvaluationImplProjectToFaceSelector::run
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
Definition evaluation_kernels_face.h:1843
internal::FEFaceEvaluationImplRaviartThomas::evaluate_or_integrate_in_face
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)
Definition evaluation_kernels_face.h:533
internal::FEFaceEvaluationImplRaviartThomas::Number2
typename FEEvaluationData< dim, Number, true >::shape_info_number_type Number2
Definition evaluation_kernels_face.h:525
internal::FEFaceEvaluationImpl::Eval
EvaluatorTensorProduct< symmetric_evaluate ? evaluate_evenodd :evaluate_general, dim - 1, fe_degree+1, n_q_points_1d, Number, Number2 > Eval
Definition evaluation_kernels_face.h:58
internal::FEFaceEvaluationImpl::evaluate_in_face
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)
Definition evaluation_kernels_face.h:91
internal::FEFaceEvaluationImpl::create_evaluator_tensor_product
static Eval create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const unsigned int subface_index, const unsigned int direction)
Definition evaluation_kernels_face.h:61
internal::FEFaceEvaluationImpl::integrate_in_face
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)
Definition evaluation_kernels_face.h:289
internal::FEFaceEvaluationImpl::Number2
typename FEEvaluationData< dim, Number, true >::shape_info_number_type Number2
Definition evaluation_kernels_face.h:50
internal::FEFaceNormalEvaluationImpl::Number2
typename FEEvaluationData< dim, Number, true >::shape_info_number_type Number2
Definition evaluation_kernels_face.h:857
internal::FEFaceNormalEvaluationImpl::interpolate_quadrature
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)
Definition evaluation_kernels_face.h:897
internal::FEFaceNormalEvaluationImpl::interpolate_generic
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)
Definition evaluation_kernels_face.h:925
internal::FEFaceNormalEvaluationImpl::interpolate
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)
Definition evaluation_kernels_face.h:861
internal::FEFaceNormalEvaluationImpl::interpolate_raviart_thomas
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)
Definition evaluation_kernels_face.h:1026
internal::MatrixFreeFunctions::ShapeInfo::data
std::vector< UnivariateShapeData< Number > > data
Definition shape_info.h:490
internal::MatrixFreeFunctions::ShapeInfo::face_orientations_quad
::Table< 2, unsigned int > face_orientations_quad
Definition shape_info.h:629
tensor_product_kernels.h
Generated by doxygen 1.9.8