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