Reference documentation for deal.II version 9.5.0
\(\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.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2017 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_matrix_free_evaluation_kernels_h
18#define dealii_matrix_free_evaluation_kernels_h
19
20#include <deal.II/base/config.h>
21
26
32
33
35
36
37namespace internal
38{
39 // Select evaluator type from element shape function type
40 template <MatrixFreeFunctions::ElementType element, bool is_long>
42 {};
43
44 template <bool is_long>
45 struct EvaluatorSelector<MatrixFreeFunctions::tensor_general, is_long>
46 {
47 static const EvaluatorVariant variant = evaluate_general;
48 };
49
50 template <>
51 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, false>
52 {
53 static const EvaluatorVariant variant = evaluate_symmetric;
54 };
55
56 template <>
57 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, true>
58 {
59 static const EvaluatorVariant variant = evaluate_evenodd;
60 };
61
62 template <bool is_long>
63 struct EvaluatorSelector<MatrixFreeFunctions::truncated_tensor, is_long>
64 {
65 static const EvaluatorVariant variant = evaluate_general;
66 };
67
68 template <>
69 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,
70 false>
71 {
72 static const EvaluatorVariant variant = evaluate_general;
73 };
74
75 template <>
76 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0, true>
77 {
78 static const EvaluatorVariant variant = evaluate_evenodd;
79 };
80
81 template <bool is_long>
82 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_collocation,
83 is_long>
84 {
85 static const EvaluatorVariant variant = evaluate_evenodd;
86 };
87
88 template <bool is_long>
89 struct EvaluatorSelector<MatrixFreeFunctions::tensor_raviart_thomas, is_long>
90 {
92 };
93
94
95
113 int dim,
114 int fe_degree,
115 int n_q_points_1d,
116 typename Number>
118 {
120 EvaluatorSelector<type, (fe_degree + n_q_points_1d > 4)>::variant;
121 using Number2 =
123
125 dim,
126 fe_degree + 1,
127 n_q_points_1d,
128 Number,
129 Number2>;
130
131 static void
132 evaluate(const unsigned int n_components,
133 const EvaluationFlags::EvaluationFlags evaluation_flag,
134 const Number * values_dofs_actual,
136
137 static void
138 integrate(const unsigned int n_components,
139 const EvaluationFlags::EvaluationFlags integration_flag,
140 Number * values_dofs_actual,
142 const bool add_into_values_array);
143
144 static Eval
147 *univariate_shape_data)
148 {
150 return Eval(univariate_shape_data->shape_values_eo,
151 univariate_shape_data->shape_gradients_eo,
152 univariate_shape_data->shape_hessians_eo,
153 univariate_shape_data->fe_degree + 1,
154 univariate_shape_data->n_q_points_1d);
155 else
156 return Eval(univariate_shape_data->shape_values,
157 univariate_shape_data->shape_gradients,
158 univariate_shape_data->shape_hessians,
159 univariate_shape_data->fe_degree + 1,
160 univariate_shape_data->n_q_points_1d);
161 }
162 };
163
164
165
170 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
171 struct FEEvaluationImpl<MatrixFreeFunctions::tensor_none,
172 dim,
173 fe_degree,
174 n_q_points_1d,
175 Number>
176 {
177 static void
178 evaluate(const unsigned int n_components,
179 const EvaluationFlags::EvaluationFlags evaluation_flag,
180 const Number * values_dofs_actual,
182
183 static void
184 integrate(const unsigned int n_components,
185 const EvaluationFlags::EvaluationFlags integration_flag,
186 Number * values_dofs_actual,
188 const bool add_into_values_array);
189 };
190
195 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
196 struct FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
197 dim,
198 fe_degree,
199 n_q_points_1d,
200 Number>
201 {
202 template <bool integrate>
203 static void
204 evaluate_or_integrate(
205 const EvaluationFlags::EvaluationFlags evaluation_flag,
206 Number * values_dofs_actual,
208 const bool add_into_values_array = false);
209
210 private:
211 template <typename EvalType, typename Number2>
212 static EvalType
215 {
216 return EvalType(shape_data.shape_values,
217 shape_data.shape_gradients,
218 shape_data.shape_hessians);
219 }
220
221 template <int normal_dir>
222 static void
223 evaluate_tensor_product_per_component(
224 const EvaluationFlags::EvaluationFlags evaluation_flag,
225 Number * values_dofs_actual,
227 const bool add_into_values_array,
228 std::integral_constant<bool, false>);
229
230 template <int normal_dir>
231 static void
232 evaluate_tensor_product_per_component(
233 const EvaluationFlags::EvaluationFlags evaluation_flag,
234 Number * values_dofs_actual,
236 const bool add_into_values_array,
237 std::integral_constant<bool, true>);
238 };
239
240
241
243 int dim,
244 int fe_degree,
245 int n_q_points_1d,
246 typename Number>
247 inline void
249 const unsigned int n_components,
250 const EvaluationFlags::EvaluationFlags evaluation_flag,
251 const Number * values_dofs_actual,
253 {
254 if (evaluation_flag == EvaluationFlags::nothing)
255 return;
256
257 std::array<const MatrixFreeFunctions::UnivariateShapeData<Number2> *, 3>
258 univariate_shape_data;
259
260 const auto &shape_data = fe_eval.get_shape_info().data;
261
262 univariate_shape_data.fill(&shape_data.front());
263
264 if (shape_data.size() == dim)
265 for (int i = 1; i < dim; ++i)
266 univariate_shape_data[i] = &shape_data[i];
267
268 Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
269 Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
270 Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);
271
272 const unsigned int temp_size =
273 Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
274 0 :
275 (Eval::n_rows_of_product > Eval::n_columns_of_product ?
276 Eval::n_rows_of_product :
277 Eval::n_columns_of_product);
278 Number *temp1 = fe_eval.get_scratch_data().begin();
279 Number *temp2;
280 if (temp_size == 0)
281 {
282 temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
283 shape_data.front().fe_degree + 1),
284 Utilities::fixed_power<dim>(
285 shape_data.front().n_q_points_1d));
286 }
287 else
288 {
289 temp2 = temp1 + temp_size;
290 }
291
292 const std::size_t n_q_points = temp_size == 0 ?
293 fe_eval.get_shape_info().n_q_points :
294 Eval::n_columns_of_product;
295 const std::size_t dofs_per_comp =
297 Utilities::pow(shape_data.front().fe_degree + 1, dim) :
299 const Number *values_dofs = values_dofs_actual;
301 {
302 const std::size_t n_dofs_per_comp =
304 Number *values_dofs_tmp =
305 temp1 + 2 * (std::max(n_dofs_per_comp, n_q_points));
306 const int degree =
307 fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
308 for (unsigned int c = 0; c < n_components; ++c)
309 for (int i = 0, count_p = 0, count_q = 0;
310 i < (dim > 2 ? degree + 1 : 1);
311 ++i)
312 {
313 for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
314 {
315 for (int k = 0; k < degree + 1 - j - i;
316 ++k, ++count_p, ++count_q)
317 values_dofs_tmp[c * dofs_per_comp + count_q] =
318 values_dofs_actual[c * n_dofs_per_comp + count_p];
319 for (int k = degree + 1 - j - i; k < degree + 1;
320 ++k, ++count_q)
321 values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
322 }
323 for (int j = degree + 1 - i; j < degree + 1; ++j)
324 for (int k = 0; k < degree + 1; ++k, ++count_q)
325 values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
326 }
327 values_dofs = values_dofs_tmp;
328 }
329
330 Number *values_quad = fe_eval.begin_values();
331 Number *gradients_quad = fe_eval.begin_gradients();
332 Number *hessians_quad = fe_eval.begin_hessians();
333
334 switch (dim)
335 {
336 case 1:
337 for (unsigned int c = 0; c < n_components; ++c)
338 {
339 if (evaluation_flag & EvaluationFlags::values)
340 eval0.template values<0, true, false>(values_dofs, values_quad);
341 if (evaluation_flag & EvaluationFlags::gradients)
342 eval0.template gradients<0, true, false>(values_dofs,
343 gradients_quad);
344 if (evaluation_flag & EvaluationFlags::hessians)
345 eval0.template hessians<0, true, false>(values_dofs,
346 hessians_quad);
347
348 // advance the next component in 1d array
349 values_dofs += dofs_per_comp;
350 values_quad += n_q_points;
351 gradients_quad += n_q_points;
352 hessians_quad += n_q_points;
353 }
354 break;
355
356 case 2:
357 for (unsigned int c = 0; c < n_components; ++c)
358 {
359 // grad x
360 if (evaluation_flag & EvaluationFlags::gradients)
361 {
362 eval0.template gradients<0, true, false>(values_dofs, temp1);
363 eval1.template values<1, true, false>(temp1, gradients_quad);
364 }
365 if (evaluation_flag & EvaluationFlags::hessians)
366 {
367 // grad xy
368 if (!(evaluation_flag & EvaluationFlags::gradients))
369 eval0.template gradients<0, true, false>(values_dofs,
370 temp1);
371 eval1.template gradients<1, true, false>(temp1,
372 hessians_quad +
373 2 * n_q_points);
374
375 // grad xx
376 eval0.template hessians<0, true, false>(values_dofs, temp1);
377 eval1.template values<1, true, false>(temp1, hessians_quad);
378 }
379
380 // grad y
381 eval0.template values<0, true, false>(values_dofs, temp1);
382 if (evaluation_flag & EvaluationFlags::gradients)
383 eval1.template gradients<1, true, false>(temp1,
384 gradients_quad +
385 n_q_points);
386
387 // grad yy
388 if (evaluation_flag & EvaluationFlags::hessians)
389 eval1.template hessians<1, true, false>(temp1,
390 hessians_quad +
391 n_q_points);
392
393 // val: can use values applied in x
394 if (evaluation_flag & EvaluationFlags::values)
395 eval1.template values<1, true, false>(temp1, values_quad);
396
397 // advance to the next component in 1d array
398 values_dofs += dofs_per_comp;
399 values_quad += n_q_points;
400 gradients_quad += 2 * n_q_points;
401 hessians_quad += 3 * n_q_points;
402 }
403 break;
404
405 case 3:
406 for (unsigned int c = 0; c < n_components; ++c)
407 {
408 if (evaluation_flag & EvaluationFlags::gradients)
409 {
410 // grad x
411 eval0.template gradients<0, true, false>(values_dofs, temp1);
412 eval1.template values<1, true, false>(temp1, temp2);
413 eval2.template values<2, true, false>(temp2, gradients_quad);
414 }
415
416 if (evaluation_flag & EvaluationFlags::hessians)
417 {
418 // grad xz
419 if (!(evaluation_flag & EvaluationFlags::gradients))
420 {
421 eval0.template gradients<0, true, false>(values_dofs,
422 temp1);
423 eval1.template values<1, true, false>(temp1, temp2);
424 }
425 eval2.template gradients<2, true, false>(temp2,
426 hessians_quad +
427 4 * n_q_points);
428
429 // grad xy
430 eval1.template gradients<1, true, false>(temp1, temp2);
431 eval2.template values<2, true, false>(temp2,
432 hessians_quad +
433 3 * n_q_points);
434
435 // grad xx
436 eval0.template hessians<0, true, false>(values_dofs, temp1);
437 eval1.template values<1, true, false>(temp1, temp2);
438 eval2.template values<2, true, false>(temp2, hessians_quad);
439 }
440
441 // grad y
442 eval0.template values<0, true, false>(values_dofs, temp1);
443 if (evaluation_flag & EvaluationFlags::gradients)
444 {
445 eval1.template gradients<1, true, false>(temp1, temp2);
446 eval2.template values<2, true, false>(temp2,
447 gradients_quad +
448 n_q_points);
449 }
450
451 if (evaluation_flag & EvaluationFlags::hessians)
452 {
453 // grad yz
454 if (!(evaluation_flag & EvaluationFlags::gradients))
455 eval1.template gradients<1, true, false>(temp1, temp2);
456 eval2.template gradients<2, true, false>(temp2,
457 hessians_quad +
458 5 * n_q_points);
459
460 // grad yy
461 eval1.template hessians<1, true, false>(temp1, temp2);
462 eval2.template values<2, true, false>(temp2,
463 hessians_quad +
464 n_q_points);
465 }
466
467 // grad z: can use the values applied in x direction stored in
468 // temp1
469 eval1.template values<1, true, false>(temp1, temp2);
470 if (evaluation_flag & EvaluationFlags::gradients)
471 eval2.template gradients<2, true, false>(temp2,
472 gradients_quad +
473 2 * n_q_points);
474
475 // grad zz: can use the values applied in x and y direction stored
476 // in temp2
477 if (evaluation_flag & EvaluationFlags::hessians)
478 eval2.template hessians<2, true, false>(temp2,
479 hessians_quad +
480 2 * n_q_points);
481
482 // val: can use the values applied in x & y direction stored in
483 // temp2
484 if (evaluation_flag & EvaluationFlags::values)
485 eval2.template values<2, true, false>(temp2, values_quad);
486
487 // advance to the next component in 1d array
488 values_dofs += dofs_per_comp;
489 values_quad += n_q_points;
490 gradients_quad += 3 * n_q_points;
491 hessians_quad += 6 * n_q_points;
492 }
493 break;
494
495 default:
497 }
498
499 // case additional dof for FE_Q_DG0: add values; gradients and second
500 // derivatives evaluate to zero
502 (evaluation_flag & EvaluationFlags::values))
503 {
504 values_quad -= n_components * n_q_points;
505 values_dofs -= n_components * dofs_per_comp;
506 for (std::size_t c = 0; c < n_components; ++c)
507 for (std::size_t q = 0; q < n_q_points; ++q)
508 values_quad[c * n_q_points + q] +=
509 values_dofs[(c + 1) * dofs_per_comp - 1];
510 }
511 }
512
513
514
516 int dim,
517 int fe_degree,
518 int n_q_points_1d,
519 typename Number>
520 inline void
522 const unsigned int n_components,
523 const EvaluationFlags::EvaluationFlags integration_flag,
524 Number * values_dofs_actual,
526 const bool add_into_values_array)
527 {
528 std::array<const MatrixFreeFunctions::UnivariateShapeData<Number2> *, 3>
529 univariate_shape_data;
530
531 const auto &shape_data = fe_eval.get_shape_info().data;
532 univariate_shape_data.fill(&shape_data.front());
533
534 if (shape_data.size() == dim)
535 for (int i = 1; i < dim; ++i)
536 univariate_shape_data[i] = &shape_data[i];
537
538 Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
539 Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
540 Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);
541
542 const unsigned int temp_size =
543 Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
544 0 :
545 (Eval::n_rows_of_product > Eval::n_columns_of_product ?
546 Eval::n_rows_of_product :
547 Eval::n_columns_of_product);
548 Number *temp1 = fe_eval.get_scratch_data().begin();
549 Number *temp2;
550 if (temp_size == 0)
551 {
552 temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
553 shape_data.front().fe_degree + 1),
554 Utilities::fixed_power<dim>(
555 shape_data.front().n_q_points_1d));
556 }
557 else
558 {
559 temp2 = temp1 + temp_size;
560 }
561
562 const std::size_t n_q_points = temp_size == 0 ?
563 fe_eval.get_shape_info().n_q_points :
564 Eval::n_columns_of_product;
565 const unsigned int dofs_per_comp =
567 Utilities::fixed_power<dim>(shape_data.front().fe_degree + 1) :
569 // expand dof_values to tensor product for truncated tensor products
570 Number *values_dofs =
572 temp1 + 2 * (std::max<std::size_t>(
574 n_q_points)) :
575 values_dofs_actual;
576
577 Number *values_quad = fe_eval.begin_values();
578 Number *gradients_quad = fe_eval.begin_gradients();
579 Number *hessians_quad = fe_eval.begin_hessians();
580
581 switch (dim)
582 {
583 case 1:
584 for (unsigned int c = 0; c < n_components; ++c)
585 {
586 if (integration_flag & EvaluationFlags::values)
587 {
588 if (add_into_values_array == false)
589 eval0.template values<0, false, false>(values_quad,
590 values_dofs);
591 else
592 eval0.template values<0, false, true>(values_quad,
593 values_dofs);
594 }
595 if (integration_flag & EvaluationFlags::gradients)
596 {
597 if (integration_flag & EvaluationFlags::values ||
598 add_into_values_array == true)
599 eval0.template gradients<0, false, true>(gradients_quad,
600 values_dofs);
601 else
602 eval0.template gradients<0, false, false>(gradients_quad,
603 values_dofs);
604 }
605 if ((integration_flag & EvaluationFlags::hessians) != 0u)
606 {
607 if ((integration_flag & EvaluationFlags::values) != 0u ||
608 (integration_flag & EvaluationFlags::gradients) != 0u ||
609 add_into_values_array == true)
610 eval0.template hessians<0, false, true>(hessians_quad,
611 values_dofs);
612 else
613 eval0.template hessians<0, false, false>(hessians_quad,
614 values_dofs);
615 }
616
617 // advance to the next component in 1d array
618 values_dofs += dofs_per_comp;
619 values_quad += n_q_points;
620 gradients_quad += n_q_points;
621 hessians_quad += n_q_points;
622 }
623 break;
624
625 case 2:
626 for (unsigned int c = 0; c < n_components; ++c)
627 {
628 if ((integration_flag & EvaluationFlags::values) &&
629 !(integration_flag & EvaluationFlags::gradients))
630 {
631 eval1.template values<1, false, false>(values_quad, temp1);
632 if (add_into_values_array == false)
633 eval0.template values<0, false, false>(temp1, values_dofs);
634 else
635 eval0.template values<0, false, true>(temp1, values_dofs);
636 }
637 if (integration_flag & EvaluationFlags::gradients)
638 {
639 eval1.template gradients<1, false, false>(gradients_quad +
640 n_q_points,
641 temp1);
642 if (integration_flag & EvaluationFlags::values)
643 eval1.template values<1, false, true>(values_quad, temp1);
644 if (add_into_values_array == false)
645 eval0.template values<0, false, false>(temp1, values_dofs);
646 else
647 eval0.template values<0, false, true>(temp1, values_dofs);
648 eval1.template values<1, false, false>(gradients_quad, temp1);
649 eval0.template gradients<0, false, true>(temp1, values_dofs);
650 }
651 if ((integration_flag & EvaluationFlags::hessians) != 0u)
652 {
653 // grad xx
654 eval1.template values<1, false, false>(hessians_quad, temp1);
655
656 if ((integration_flag & EvaluationFlags::values) != 0u ||
657 (integration_flag & EvaluationFlags::gradients) != 0u ||
658 add_into_values_array == true)
659 eval0.template hessians<0, false, true>(temp1, values_dofs);
660 else
661 eval0.template hessians<0, false, false>(temp1,
662 values_dofs);
663
664 // grad yy
665 eval1.template hessians<1, false, false>(hessians_quad +
666 n_q_points,
667 temp1);
668 eval0.template values<0, false, true>(temp1, values_dofs);
669
670 // grad xy
671 eval1.template gradients<1, false, false>(hessians_quad +
672 2 * n_q_points,
673 temp1);
674 eval0.template gradients<0, false, true>(temp1, values_dofs);
675 }
676
677 // advance to the next component in 1d array
678 values_dofs += dofs_per_comp;
679 values_quad += n_q_points;
680 gradients_quad += 2 * n_q_points;
681 hessians_quad += 3 * n_q_points;
682 }
683 break;
684
685 case 3:
686 for (unsigned int c = 0; c < n_components; ++c)
687 {
688 if ((integration_flag & EvaluationFlags::values) &&
689 !(integration_flag & EvaluationFlags::gradients))
690 {
691 eval2.template values<2, false, false>(values_quad, temp1);
692 eval1.template values<1, false, false>(temp1, temp2);
693 if (add_into_values_array == false)
694 eval0.template values<0, false, false>(temp2, values_dofs);
695 else
696 eval0.template values<0, false, true>(temp2, values_dofs);
697 }
698 if (integration_flag & EvaluationFlags::gradients)
699 {
700 eval2.template gradients<2, false, false>(gradients_quad +
701 2 * n_q_points,
702 temp1);
703 if (integration_flag & EvaluationFlags::values)
704 eval2.template values<2, false, true>(values_quad, temp1);
705 eval1.template values<1, false, false>(temp1, temp2);
706 eval2.template values<2, false, false>(gradients_quad +
707 n_q_points,
708 temp1);
709 eval1.template gradients<1, false, true>(temp1, temp2);
710 if (add_into_values_array == false)
711 eval0.template values<0, false, false>(temp2, values_dofs);
712 else
713 eval0.template values<0, false, true>(temp2, values_dofs);
714 eval2.template values<2, false, false>(gradients_quad, temp1);
715 eval1.template values<1, false, false>(temp1, temp2);
716 eval0.template gradients<0, false, true>(temp2, values_dofs);
717 }
718 if ((integration_flag & EvaluationFlags::hessians) != 0u)
719 {
720 // grad xx
721 eval2.template values<2, false, false>(hessians_quad, temp1);
722 eval1.template values<1, false, false>(temp1, temp2);
723
724 if ((integration_flag & EvaluationFlags::values) != 0u ||
725 (integration_flag & EvaluationFlags::gradients) != 0u ||
726 add_into_values_array == true)
727 eval0.template hessians<0, false, true>(temp2, values_dofs);
728 else
729 eval0.template hessians<0, false, false>(temp2,
730 values_dofs);
731
732 // grad yy
733 eval2.template values<2, false, false>(hessians_quad +
734 n_q_points,
735 temp1);
736 eval1.template hessians<1, false, false>(temp1, temp2);
737 eval0.template values<0, false, true>(temp2, values_dofs);
738
739 // grad zz
740 eval2.template hessians<2, false, false>(hessians_quad +
741 2 * n_q_points,
742 temp1);
743 eval1.template values<1, false, false>(temp1, temp2);
744 eval0.template values<0, false, true>(temp2, values_dofs);
745
746 // grad xy
747 eval2.template values<2, false, false>(hessians_quad +
748 3 * n_q_points,
749 temp1);
750 eval1.template gradients<1, false, false>(temp1, temp2);
751 eval0.template gradients<0, false, true>(temp2, values_dofs);
752
753 // grad xz
754 eval2.template gradients<2, false, false>(hessians_quad +
755 4 * n_q_points,
756 temp1);
757 eval1.template values<1, false, false>(temp1, temp2);
758 eval0.template gradients<0, false, true>(temp2, values_dofs);
759
760 // grad yz
761 eval2.template gradients<2, false, false>(hessians_quad +
762 5 * n_q_points,
763 temp1);
764 eval1.template gradients<1, false, false>(temp1, temp2);
765 eval0.template values<0, false, true>(temp2, values_dofs);
766 }
767
768 // advance to the next component in 1d array
769 values_dofs += dofs_per_comp;
770 values_quad += n_q_points;
771 gradients_quad += 3 * n_q_points;
772 hessians_quad += 6 * n_q_points;
773 }
774 break;
775
776 default:
778 }
779
780 // case FE_Q_DG0: add values, gradients and second derivatives are zero
782 {
783 values_dofs -= n_components * dofs_per_comp - dofs_per_comp + 1;
784 values_quad -= n_components * n_q_points;
785 if (integration_flag & EvaluationFlags::values)
786 for (unsigned int c = 0; c < n_components; ++c)
787 {
788 values_dofs[0] = values_quad[0];
789 for (unsigned int q = 1; q < n_q_points; ++q)
790 values_dofs[0] += values_quad[q];
791 values_dofs += dofs_per_comp;
792 values_quad += n_q_points;
793 }
794 else
795 {
796 for (unsigned int c = 0; c < n_components; ++c)
797 values_dofs[c * dofs_per_comp] = Number();
798 values_dofs += n_components * dofs_per_comp;
799 }
800 }
801
803 {
804 const std::size_t n_dofs_per_comp =
806 values_dofs -= dofs_per_comp * n_components;
807 const int degree =
808 fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
809 for (unsigned int c = 0; c < n_components; ++c)
810 for (int i = 0, count_p = 0, count_q = 0;
811 i < (dim > 2 ? degree + 1 : 1);
812 ++i)
813 {
814 for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
815 {
816 for (int k = 0; k < degree + 1 - j - i;
817 ++k, ++count_p, ++count_q)
818 values_dofs_actual[c * n_dofs_per_comp + count_p] =
819 values_dofs[c * dofs_per_comp + count_q];
820 count_q += j + i;
821 }
822 count_q += i * (degree + 1);
823 }
824 }
825 }
826
827
828
829 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
830 inline void
833 dim,
834 fe_degree,
835 n_q_points_1d,
836 Number>::evaluate(const unsigned int n_components,
837 const EvaluationFlags::EvaluationFlags evaluation_flag,
838 const Number * values_dofs_actual,
840 {
841 Assert(!(evaluation_flag & EvaluationFlags::hessians), ExcNotImplemented());
842
843 const std::size_t n_dofs =
845 const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
846
847 const auto &shape_data = fe_eval.get_shape_info().data;
848
849 using Number2 =
851 using Eval =
853
854 if (evaluation_flag & EvaluationFlags::values)
855 {
856 const auto shape_values = shape_data.front().shape_values.data();
857 auto values_quad_ptr = fe_eval.begin_values();
858 auto values_dofs_actual_ptr = values_dofs_actual;
859
860 Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
861 for (unsigned int c = 0; c < n_components; ++c)
862 {
863 eval.template values<0, true, false>(values_dofs_actual_ptr,
864 values_quad_ptr);
865
866 values_quad_ptr += n_q_points;
867 values_dofs_actual_ptr += n_dofs;
868 }
869 }
870
871 if (evaluation_flag & EvaluationFlags::gradients)
872 {
873 const auto shape_gradients = shape_data.front().shape_gradients.data();
874 auto gradients_quad_ptr = fe_eval.begin_gradients();
875 auto values_dofs_actual_ptr = values_dofs_actual;
876
877 for (unsigned int c = 0; c < n_components; ++c)
878 {
879 for (unsigned int d = 0; d < dim; ++d)
880 {
881 Eval eval(nullptr,
882 shape_gradients + n_q_points * n_dofs * d,
883 nullptr,
884 n_dofs,
885 n_q_points);
886
887 eval.template gradients<0, true, false>(values_dofs_actual_ptr,
888 gradients_quad_ptr);
889
890 gradients_quad_ptr += n_q_points;
891 }
892 values_dofs_actual_ptr += n_dofs;
893 }
894 }
895 }
896
897
898
899 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
900 inline void
903 dim,
904 fe_degree,
905 n_q_points_1d,
906 Number>::integrate(const unsigned int n_components,
907 const EvaluationFlags::EvaluationFlags integration_flag,
908 Number * values_dofs_actual,
910 const bool add_into_values_array)
911 {
912 Assert(!(integration_flag & EvaluationFlags::hessians),
914
915 const std::size_t n_dofs =
917 const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
918
919 const auto &shape_data = fe_eval.get_shape_info().data;
920
921 using Number2 =
923 using Eval =
925
926 if (integration_flag & EvaluationFlags::values)
927 {
928 const auto shape_values = shape_data.front().shape_values.data();
929 auto values_quad_ptr = fe_eval.begin_values();
930 auto values_dofs_actual_ptr = values_dofs_actual;
931
932 Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
933 for (unsigned int c = 0; c < n_components; ++c)
934 {
935 if (add_into_values_array == false)
936 eval.template values<0, false, false>(values_quad_ptr,
937 values_dofs_actual_ptr);
938 else
939 eval.template values<0, false, true>(values_quad_ptr,
940 values_dofs_actual_ptr);
941
942 values_quad_ptr += n_q_points;
943 values_dofs_actual_ptr += n_dofs;
944 }
945 }
946
947 if (integration_flag & EvaluationFlags::gradients)
948 {
949 const auto shape_gradients = shape_data.front().shape_gradients.data();
950 auto gradients_quad_ptr = fe_eval.begin_gradients();
951 auto values_dofs_actual_ptr = values_dofs_actual;
952
953 for (unsigned int c = 0; c < n_components; ++c)
954 {
955 for (unsigned int d = 0; d < dim; ++d)
956 {
957 Eval eval(nullptr,
958 shape_gradients + n_q_points * n_dofs * d,
959 nullptr,
960 n_dofs,
961 n_q_points);
962
963 if ((add_into_values_array == false &&
964 !(integration_flag & EvaluationFlags::values)) &&
965 d == 0)
966 eval.template gradients<0, false, false>(
967 gradients_quad_ptr, values_dofs_actual_ptr);
968 else
969 eval.template gradients<0, false, true>(
970 gradients_quad_ptr, values_dofs_actual_ptr);
971
972 gradients_quad_ptr += n_q_points;
973 }
974 values_dofs_actual_ptr += n_dofs;
975 }
976 }
977 }
978
979
980 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
981 template <bool integrate>
982 inline void
984 dim,
985 fe_degree,
986 n_q_points_1d,
987 Number>::
988 evaluate_or_integrate(
989 const EvaluationFlags::EvaluationFlags evaluation_flag,
990 Number * values_dofs_actual,
992 const bool add_into_values_array)
993 {
994 if (evaluation_flag == EvaluationFlags::nothing)
995 return;
996
997 AssertDimension(fe_eval.get_shape_info().data.size(), 2);
998 // First component:
999 evaluate_tensor_product_per_component<0>(
1000 evaluation_flag,
1001 values_dofs_actual,
1002 fe_eval,
1003 add_into_values_array,
1004 std::integral_constant<bool, integrate>());
1005 // Second component :
1006 evaluate_tensor_product_per_component<1>(
1007 evaluation_flag,
1008 values_dofs_actual,
1009 fe_eval,
1010 add_into_values_array,
1011 std::integral_constant<bool, integrate>());
1012 if (dim == 3)
1013 {
1014 // Third component
1015 evaluate_tensor_product_per_component<2>(
1016 evaluation_flag,
1017 values_dofs_actual,
1018 fe_eval,
1019 add_into_values_array,
1020 std::integral_constant<bool, integrate>());
1021 }
1022 }
1023
1024 // Helper function that applies the 1d evaluation kernels.
1025 // std::integral_constant<bool, false> is the interpolation path, and
1026 // std::integral_constant<bool, true> below is the integration path.
1027 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1028 template <int normal_dir>
1029 inline void
1031 dim,
1032 fe_degree,
1033 n_q_points_1d,
1034 Number>::
1035 evaluate_tensor_product_per_component(
1036 const EvaluationFlags::EvaluationFlags evaluation_flag,
1037 Number * values_dofs_actual,
1039 const bool add_into_values_array,
1040 std::integral_constant<bool, false>)
1041 {
1042 (void)add_into_values_array;
1043
1044 using Number2 =
1046 using EvalNormal =
1048 dim,
1049 (fe_degree == -1) ? 1 : fe_degree + 1,
1050 n_q_points_1d,
1051 normal_dir,
1052 Number,
1053 Number2>;
1054
1055 using EvalTangent =
1057 dim,
1058 (fe_degree == -1) ? 1 : fe_degree,
1059 n_q_points_1d,
1060 normal_dir,
1061 Number,
1062 Number2>;
1063 using Eval0 =
1064 typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
1065 using Eval1 =
1066 typename std::conditional<normal_dir == 1, EvalNormal, EvalTangent>::type;
1067 using Eval2 =
1068 typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
1069
1070 const auto &shape_info = fe_eval.get_shape_info();
1071 Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
1072 ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]));
1073 Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
1074 ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]));
1075 Eval2 eval2 = create_evaluator_tensor_product<Eval2>(
1076 ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1]));
1077
1078 Number *temp1 = fe_eval.get_scratch_data().begin();
1079 Number *temp2;
1080
1081 temp2 =
1082 temp1 +
1083 std::max(Utilities::fixed_power<dim>(shape_info.data[0].fe_degree + 1),
1084 Utilities::fixed_power<dim>(shape_info.data[0].n_q_points_1d));
1085
1086 const std::size_t n_q_points = shape_info.n_q_points;
1087 const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell;
1088
1089 // Initial shift depending on component (normal_dir)
1090 Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
1091 Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
1092 Number *gradients_quad =
1093 fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
1094 Number *hessians_quad =
1095 (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
1096 fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
1097
1098 switch (dim)
1099 {
1100 case 2:
1101 if (evaluation_flag & EvaluationFlags::gradients)
1102 {
1103 eval0.template gradients<0, true, false>(values_dofs, temp1);
1104 eval1.template values<1, true, false>(temp1, gradients_quad);
1105 }
1106 if (evaluation_flag & EvaluationFlags::hessians)
1107 {
1108 // The evaluation/integration here *should* work, however
1109 // the piola transform is not implemented.
1111 // grad xy
1112 if (!(evaluation_flag & EvaluationFlags::gradients))
1113 eval0.template gradients<0, true, false>(values_dofs, temp1);
1114 eval1.template gradients<1, true, false>(temp1,
1115 hessians_quad +
1116 2 * n_q_points);
1117
1118 // grad xx
1119 eval0.template hessians<0, true, false>(values_dofs, temp1);
1120 eval1.template values<1, true, false>(temp1, hessians_quad);
1121 }
1122
1123 // grad y
1124 eval0.template values<0, true, false>(values_dofs, temp1);
1125 if (evaluation_flag & EvaluationFlags::gradients)
1126 eval1.template gradients<1, true, false>(temp1,
1127 gradients_quad +
1128 n_q_points);
1129
1130 // grad yy
1131 if (evaluation_flag & EvaluationFlags::hessians)
1132 eval1.template hessians<1, true, false>(temp1,
1133 hessians_quad + n_q_points);
1134
1135 // val: can use values applied in x
1136 if (evaluation_flag & EvaluationFlags::values)
1137 eval1.template values<1, true, false>(temp1, values_quad);
1138 break;
1139 case 3:
1140 if (evaluation_flag & EvaluationFlags::gradients)
1141 {
1142 // grad x
1143 eval0.template gradients<0, true, false>(values_dofs, temp1);
1144 eval1.template values<1, true, false>(temp1, temp2);
1145 eval2.template values<2, true, false>(temp2, gradients_quad);
1146 }
1147
1148 if (evaluation_flag & EvaluationFlags::hessians)
1149 {
1150 // The evaluation/integration here *should* work, however
1151 // the piola transform is not implemented.
1153 // grad xz
1154 if (!(evaluation_flag & EvaluationFlags::gradients))
1155 {
1156 eval0.template gradients<0, true, false>(values_dofs, temp1);
1157 eval1.template values<1, true, false>(temp1, temp2);
1158 }
1159 eval2.template gradients<2, true, false>(temp2,
1160 hessians_quad +
1161 4 * n_q_points);
1162
1163 // grad xy
1164 eval1.template gradients<1, true, false>(temp1, temp2);
1165 eval2.template values<2, true, false>(temp2,
1166 hessians_quad +
1167 3 * n_q_points);
1168
1169 // grad xx
1170 eval0.template hessians<0, true, false>(values_dofs, temp1);
1171 eval1.template values<1, true, false>(temp1, temp2);
1172 eval2.template values<2, true, false>(temp2, hessians_quad);
1173 }
1174
1175 // grad y
1176 eval0.template values<0, true, false>(values_dofs, temp1);
1177 if (evaluation_flag & EvaluationFlags::gradients)
1178 {
1179 eval1.template gradients<1, true, false>(temp1, temp2);
1180 eval2.template values<2, true, false>(temp2,
1181 gradients_quad +
1182 n_q_points);
1183 }
1184
1185 if (evaluation_flag & EvaluationFlags::hessians)
1186 {
1187 // grad yz
1188 if (!(evaluation_flag & EvaluationFlags::gradients))
1189 eval1.template gradients<1, true, false>(temp1, temp2);
1190 eval2.template gradients<2, true, false>(temp2,
1191 hessians_quad +
1192 5 * n_q_points);
1193
1194 // grad yy
1195 eval1.template hessians<1, true, false>(temp1, temp2);
1196 eval2.template values<2, true, false>(temp2,
1197 hessians_quad + n_q_points);
1198 }
1199
1200 // grad z: can use the values applied in x direction stored in
1201 // temp1
1202 eval1.template values<1, true, false>(temp1, temp2);
1203 if (evaluation_flag & EvaluationFlags::gradients)
1204 eval2.template gradients<2, true, false>(temp2,
1205 gradients_quad +
1206 2 * n_q_points);
1207
1208 // grad zz: can use the values applied in x and y direction stored
1209 // in temp2
1210 if (evaluation_flag & EvaluationFlags::hessians)
1211 eval2.template hessians<2, true, false>(temp2,
1212 hessians_quad +
1213 2 * n_q_points);
1214
1215 // val: can use the values applied in x & y direction stored in
1216 // temp2
1217 if (evaluation_flag & EvaluationFlags::values)
1218 eval2.template values<2, true, false>(temp2, values_quad);
1219 break;
1220 default:
1222 }
1223 }
1224
1225 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1226 template <int normal_dir>
1227 inline void
1229 dim,
1230 fe_degree,
1231 n_q_points_1d,
1232 Number>::
1233 evaluate_tensor_product_per_component(
1234 const EvaluationFlags::EvaluationFlags evaluation_flag,
1235 Number * values_dofs_actual,
1237 const bool add_into_values_array,
1238 std::integral_constant<bool, true>)
1239 {
1240 using Number2 =
1242 using EvalNormal =
1244 dim,
1245 (fe_degree == -1) ? 1 : fe_degree + 1,
1246 n_q_points_1d,
1247 normal_dir,
1248 Number,
1249 Number2>;
1250
1251 using EvalTangent =
1253 dim,
1254 (fe_degree == -1) ? 1 : fe_degree,
1255 n_q_points_1d,
1256 normal_dir,
1257 Number,
1258 Number2>;
1259 using Eval0 =
1260 typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
1261 using Eval1 =
1262 typename std::conditional<normal_dir == 1, EvalNormal, EvalTangent>::type;
1263 using Eval2 =
1264 typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
1265
1266 const auto &shape_info = fe_eval.get_shape_info();
1267 Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
1268 ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]));
1269 Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
1270 ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]));
1271 Eval2 eval2 = create_evaluator_tensor_product<Eval2>(
1272 ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1]));
1273
1274 Number *temp1 = fe_eval.get_scratch_data().begin();
1275 Number *temp2;
1276
1277 temp2 =
1278 temp1 +
1279 std::max(Utilities::fixed_power<dim>(shape_info.data[0].fe_degree + 1),
1280 Utilities::fixed_power<dim>(shape_info.data[0].n_q_points_1d));
1281
1282 const std::size_t n_q_points = shape_info.n_q_points;
1283 const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell;
1284
1285 // Initial shift depending on component (normal_dir)
1286 Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
1287 Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
1288 Number *gradients_quad =
1289 fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
1290 Number *hessians_quad =
1291 (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
1292 fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
1293
1294 // Integrate path
1295 switch (dim)
1296 {
1297 case 2:
1298 if ((evaluation_flag & EvaluationFlags::values) &&
1299 !(evaluation_flag & EvaluationFlags::gradients))
1300 {
1301 eval1.template values<1, false, false>(values_quad, temp1);
1302 if (add_into_values_array == false)
1303 eval0.template values<0, false, false>(temp1, values_dofs);
1304 else
1305 eval0.template values<0, false, true>(temp1, values_dofs);
1306 }
1307 if (evaluation_flag & EvaluationFlags::gradients)
1308 {
1309 eval1.template gradients<1, false, false>(gradients_quad +
1310 n_q_points,
1311 temp1);
1312 if ((evaluation_flag & EvaluationFlags::values))
1313 eval1.template values<1, false, true>(values_quad, temp1);
1314 if (add_into_values_array == false)
1315 eval0.template values<0, false, false>(temp1, values_dofs);
1316 else
1317 eval0.template values<0, false, true>(temp1, values_dofs);
1318 eval1.template values<1, false, false>(gradients_quad, temp1);
1319 eval0.template gradients<0, false, true>(temp1, values_dofs);
1320 }
1321 if (evaluation_flag & EvaluationFlags::hessians)
1322 {
1323 // grad xx
1324 eval1.template values<1, false, false>(hessians_quad, temp1);
1325
1326 if ((evaluation_flag & EvaluationFlags::values) ||
1327 (evaluation_flag & EvaluationFlags::gradients) ||
1328 add_into_values_array == true)
1329 eval0.template hessians<0, false, true>(temp1, values_dofs);
1330 else
1331 eval0.template hessians<0, false, false>(temp1, values_dofs);
1332
1333 // grad yy
1334 eval1.template hessians<1, false, false>(hessians_quad +
1335 n_q_points,
1336 temp1);
1337 eval0.template values<0, false, true>(temp1, values_dofs);
1338
1339 // grad xy
1340 eval1.template gradients<1, false, false>(hessians_quad +
1341 2 * n_q_points,
1342 temp1);
1343 eval0.template gradients<0, false, true>(temp1, values_dofs);
1344 }
1345 break;
1346
1347 case 3:
1348 if ((evaluation_flag & EvaluationFlags::values) &&
1349 !(evaluation_flag & EvaluationFlags::gradients))
1350 {
1351 eval2.template values<2, false, false>(values_quad, temp1);
1352 eval1.template values<1, false, false>(temp1, temp2);
1353 if (add_into_values_array == false)
1354 eval0.template values<0, false, false>(temp2, values_dofs);
1355 else
1356 eval0.template values<0, false, true>(temp2, values_dofs);
1357 }
1358 if (evaluation_flag & EvaluationFlags::gradients)
1359 {
1360 eval2.template gradients<2, false, false>(gradients_quad +
1361 2 * n_q_points,
1362 temp1);
1363 if ((evaluation_flag & EvaluationFlags::values))
1364 eval2.template values<2, false, true>(values_quad, temp1);
1365 eval1.template values<1, false, false>(temp1, temp2);
1366 eval2.template values<2, false, false>(gradients_quad +
1367 n_q_points,
1368 temp1);
1369 eval1.template gradients<1, false, true>(temp1, temp2);
1370 if (add_into_values_array == false)
1371 eval0.template values<0, false, false>(temp2, values_dofs);
1372 else
1373 eval0.template values<0, false, true>(temp2, values_dofs);
1374 eval2.template values<2, false, false>(gradients_quad, temp1);
1375 eval1.template values<1, false, false>(temp1, temp2);
1376 eval0.template gradients<0, false, true>(temp2, values_dofs);
1377 }
1378 if (evaluation_flag & EvaluationFlags::hessians)
1379 {
1380 // grad xx
1381 eval2.template values<2, false, false>(hessians_quad, temp1);
1382 eval1.template values<1, false, false>(temp1, temp2);
1383
1384 if ((evaluation_flag & EvaluationFlags::values) ||
1385 (evaluation_flag & EvaluationFlags::gradients) ||
1386 add_into_values_array == true)
1387 eval0.template hessians<0, false, true>(temp2, values_dofs);
1388 else
1389 eval0.template hessians<0, false, false>(temp2, values_dofs);
1390
1391 // grad yy
1392 eval2.template values<2, false, false>(hessians_quad + n_q_points,
1393 temp1);
1394 eval1.template hessians<1, false, false>(temp1, temp2);
1395 eval0.template values<0, false, true>(temp2, values_dofs);
1396
1397 // grad zz
1398 eval2.template hessians<2, false, false>(hessians_quad +
1399 2 * n_q_points,
1400 temp1);
1401 eval1.template values<1, false, false>(temp1, temp2);
1402 eval0.template values<0, false, true>(temp2, values_dofs);
1403
1404 // grad xy
1405 eval2.template values<2, false, false>(hessians_quad +
1406 3 * n_q_points,
1407 temp1);
1408 eval1.template gradients<1, false, false>(temp1, temp2);
1409 eval0.template gradients<0, false, true>(temp2, values_dofs);
1410
1411 // grad xz
1412 eval2.template gradients<2, false, false>(hessians_quad +
1413 4 * n_q_points,
1414 temp1);
1415 eval1.template values<1, false, false>(temp1, temp2);
1416 eval0.template gradients<0, false, true>(temp2, values_dofs);
1417
1418 // grad yz
1419 eval2.template gradients<2, false, false>(hessians_quad +
1420 5 * n_q_points,
1421 temp1);
1422 eval1.template gradients<1, false, false>(temp1, temp2);
1423 eval0.template values<0, false, true>(temp2, values_dofs);
1424 }
1425
1426 break;
1427 default:
1429 }
1430 }
1431
1449 template <EvaluatorVariant variant,
1450 EvaluatorQuantity quantity,
1451 int dim,
1452 int basis_size_1,
1453 int basis_size_2,
1454 typename = bool,
1455 typename = bool>
1457 {
1458 static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
1459 "The second dimension must not be smaller than the first");
1460
1483 template <typename Number, typename Number2>
1484#ifndef DEBUG
1486#endif
1487 static void
1488 do_forward(const unsigned int n_components,
1489 const AlignedVector<Number2> &transformation_matrix,
1490 const Number * values_in,
1491 Number * values_out,
1492 const unsigned int basis_size_1_variable =
1494 const unsigned int basis_size_2_variable =
1496 {
1497 Assert(
1498 basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
1499 ExcMessage("The second dimension must not be smaller than the first"));
1500
1502
1503 // we do recursion until dim==1 or dim==2 and we have
1504 // basis_size_1==basis_size_2. The latter optimization increases
1505 // optimization possibilities for the compiler but does only work for
1506 // aliased pointers if the sizes are equal.
1507 constexpr int next_dim =
1508 (dim > 2 ||
1509 ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
1510 dim - 1 :
1511 dim;
1512
1513 EvaluatorTensorProduct<variant,
1514 dim,
1515 basis_size_1,
1516 (basis_size_1 == 0 ? 0 : basis_size_2),
1517 Number,
1518 Number2>
1519 eval_val(transformation_matrix,
1520 {},
1521 {},
1522 basis_size_1_variable,
1523 basis_size_2_variable);
1524 const unsigned int np_1 =
1525 basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
1526 const unsigned int np_2 =
1527 basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
1528 Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
1529 ExcMessage("Cannot transform with 0-point basis"));
1530 Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
1531 ExcMessage("Cannot transform with 0-point basis"));
1532
1533 // run loop backwards to ensure correctness if values_in aliases with
1534 // values_out in case with basis_size_1 < basis_size_2
1535 values_in = values_in + n_components * Utilities::fixed_power<dim>(np_1);
1536 values_out =
1537 values_out + n_components * Utilities::fixed_power<dim>(np_2);
1538 for (unsigned int c = n_components; c != 0; --c)
1539 {
1540 values_in -= Utilities::fixed_power<dim>(np_1);
1541 values_out -= Utilities::fixed_power<dim>(np_2);
1542 if (next_dim < dim)
1543 for (unsigned int q = np_1; q != 0; --q)
1545 quantity,
1546 next_dim,
1547 basis_size_1,
1548 basis_size_2>::
1549 do_forward(1,
1550 transformation_matrix,
1551 values_in +
1552 (q - 1) * Utilities::fixed_power<next_dim>(np_1),
1553 values_out +
1554 (q - 1) * Utilities::fixed_power<next_dim>(np_2),
1555 basis_size_1_variable,
1556 basis_size_2_variable);
1557
1558 // the recursion stops if dim==1 or if dim==2 and
1559 // basis_size_1==basis_size_2 (the latter is used because the
1560 // compiler generates nicer code)
1561 if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
1562 {
1563 eval_val.template values<0, true, false>(values_in, values_out);
1564 eval_val.template values<1, true, false>(values_out, values_out);
1565 }
1566 else if (dim == 1)
1567 eval_val.template values<dim - 1, true, false>(values_in,
1568 values_out);
1569 else
1570 eval_val.template values<dim - 1, true, false>(values_out,
1571 values_out);
1572 }
1573 }
1574
1605 template <typename Number, typename Number2>
1606#ifndef DEBUG
1608#endif
1609 static void
1610 do_backward(const unsigned int n_components,
1611 const AlignedVector<Number2> &transformation_matrix,
1612 const bool add_into_result,
1613 Number * values_in,
1614 Number * values_out,
1615 const unsigned int basis_size_1_variable =
1617 const unsigned int basis_size_2_variable =
1619 {
1620 Assert(
1621 basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
1622 ExcMessage("The second dimension must not be smaller than the first"));
1623 Assert(add_into_result == false || values_in != values_out,
1624 ExcMessage(
1625 "Input and output cannot alias with each other when "
1626 "adding the result of the basis change to existing data"));
1627
1628 Assert(quantity == EvaluatorQuantity::value ||
1629 quantity == EvaluatorQuantity::hessian,
1631
1632 constexpr int next_dim =
1633 (dim > 2 ||
1634 ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
1635 dim - 1 :
1636 dim;
1637 EvaluatorTensorProduct<variant,
1638 dim,
1639 basis_size_1,
1640 (basis_size_1 == 0 ? 0 : basis_size_2),
1641 Number,
1642 Number2>
1643 eval_val(transformation_matrix,
1644 transformation_matrix,
1645 transformation_matrix,
1646 basis_size_1_variable,
1647 basis_size_2_variable);
1648 const unsigned int np_1 =
1649 basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
1650 const unsigned int np_2 =
1651 basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
1652 Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
1653 ExcMessage("Cannot transform with 0-point basis"));
1654 Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
1655 ExcMessage("Cannot transform with 0-point basis"));
1656
1657 for (unsigned int c = 0; c < n_components; ++c)
1658 {
1659 if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
1660 {
1661 if (quantity == EvaluatorQuantity::value)
1662 eval_val.template values<1, false, false>(values_in, values_in);
1663 else
1664 eval_val.template hessians<1, false, false>(values_in,
1665 values_in);
1666
1667 if (add_into_result)
1668 {
1669 if (quantity == EvaluatorQuantity::value)
1670 eval_val.template values<0, false, true>(values_in,
1671 values_out);
1672 else
1673 eval_val.template hessians<0, false, true>(values_in,
1674 values_out);
1675 }
1676 else
1677 {
1678 if (quantity == EvaluatorQuantity::value)
1679 eval_val.template values<0, false, false>(values_in,
1680 values_out);
1681 else
1682 eval_val.template hessians<0, false, false>(values_in,
1683 values_out);
1684 }
1685 }
1686 else
1687 {
1688 if (dim == 1 && add_into_result)
1689 {
1690 if (quantity == EvaluatorQuantity::value)
1691 eval_val.template values<0, false, true>(values_in,
1692 values_out);
1693 else
1694 eval_val.template hessians<0, false, true>(values_in,
1695 values_out);
1696 }
1697 else if (dim == 1)
1698 {
1699 if (quantity == EvaluatorQuantity::value)
1700 eval_val.template values<0, false, false>(values_in,
1701 values_out);
1702 else
1703 eval_val.template hessians<0, false, false>(values_in,
1704 values_out);
1705 }
1706 else
1707 {
1708 if (quantity == EvaluatorQuantity::value)
1709 eval_val.template values<dim - 1, false, false>(values_in,
1710 values_in);
1711 else
1712 eval_val.template hessians<dim - 1, false, false>(
1713 values_in, values_in);
1714 }
1715 }
1716 if (next_dim < dim)
1717 for (unsigned int q = 0; q < np_1; ++q)
1719 quantity,
1720 next_dim,
1721 basis_size_1,
1722 basis_size_2>::
1723 do_backward(1,
1724 transformation_matrix,
1725 add_into_result,
1726 values_in +
1727 q * Utilities::fixed_power<next_dim>(np_2),
1728 values_out +
1729 q * Utilities::fixed_power<next_dim>(np_1),
1730 basis_size_1_variable,
1731 basis_size_2_variable);
1732
1733 values_in += Utilities::fixed_power<dim>(np_2);
1734 values_out += Utilities::fixed_power<dim>(np_1);
1735 }
1736 }
1737
1758 template <typename Number, typename Number2>
1759 static void
1760 do_mass(const unsigned int n_components,
1761 const AlignedVector<Number2> &transformation_matrix,
1762 const AlignedVector<Number> & coefficients,
1763 const Number * values_in,
1764 Number * scratch_data,
1765 Number * values_out)
1766 {
1767 constexpr int next_dim = dim > 1 ? dim - 1 : dim;
1768 Number * my_scratch =
1769 basis_size_1 != basis_size_2 ? scratch_data : values_out;
1770
1771 const unsigned int size_per_component = Utilities::pow(basis_size_2, dim);
1772 Assert(coefficients.size() == size_per_component ||
1773 coefficients.size() == n_components * size_per_component,
1774 ExcDimensionMismatch(coefficients.size(), size_per_component));
1775 const unsigned int stride =
1776 coefficients.size() == size_per_component ? 0 : 1;
1777
1778 for (unsigned int q = basis_size_1; q != 0; --q)
1780 variant,
1782 next_dim,
1783 basis_size_1,
1784 basis_size_2>::do_forward(n_components,
1785 transformation_matrix,
1786 values_in +
1787 (q - 1) *
1788 Utilities::pow(basis_size_1, dim - 1),
1789 my_scratch +
1790 (q - 1) *
1791 Utilities::pow(basis_size_2, dim - 1));
1792 EvaluatorTensorProduct<variant,
1793 dim,
1794 basis_size_1,
1795 basis_size_2,
1796 Number,
1797 Number2>
1798 eval_val(transformation_matrix);
1799 const unsigned int n_inner_blocks =
1800 (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
1801 const unsigned int n_blocks = Utilities::pow(basis_size_2, dim - 1);
1802 for (unsigned int ii = 0; ii < n_blocks; ii += n_inner_blocks)
1803 for (unsigned int c = 0; c < n_components; ++c)
1804 {
1805 for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1806 eval_val.template values_one_line<dim - 1, true, false>(
1807 my_scratch + i, my_scratch + i);
1808 for (unsigned int q = 0; q < basis_size_2; ++q)
1809 for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1810 my_scratch[i + q * n_blocks + c * size_per_component] *=
1811 coefficients[i + q * n_blocks +
1812 c * stride * size_per_component];
1813 for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1814 eval_val.template values_one_line<dim - 1, false, false>(
1815 my_scratch + i, my_scratch + i);
1816 }
1817 for (unsigned int q = 0; q < basis_size_1; ++q)
1820 next_dim,
1821 basis_size_1,
1822 basis_size_2>::
1823 do_backward(n_components,
1824 transformation_matrix,
1825 false,
1826 my_scratch + q * Utilities::pow(basis_size_2, dim - 1),
1827 values_out + q * Utilities::pow(basis_size_1, dim - 1));
1828 }
1829 };
1830
1831
1832
1845 template <int dim, int fe_degree, typename Number>
1847 {
1848 using Number2 =
1851 dim,
1852 fe_degree + 1,
1853 fe_degree + 1,
1854 Number,
1855 Number2>;
1856
1857 static void
1858 evaluate(const unsigned int n_components,
1859 const EvaluationFlags::EvaluationFlags evaluation_flag,
1860 const Number * values_dofs,
1862
1863 static void
1865 const EvaluationFlags::EvaluationFlags evaluation_flag,
1866 const Number * values_dofs,
1867 Number * gradients_quad,
1868 Number * hessians_quad);
1869
1870 static void
1871 integrate(const unsigned int n_components,
1872 const EvaluationFlags::EvaluationFlags integration_flag,
1873 Number * values_dofs,
1875 const bool add_into_values_array);
1876
1877 static void
1879 const EvaluationFlags::EvaluationFlags integration_flag,
1880 Number * values_dofs,
1881 Number * gradients_quad,
1882 const Number * hessians_quad,
1883 const bool add_into_values_array);
1884 };
1885
1886
1887
1888 template <int dim, int fe_degree, typename Number>
1889 inline void
1891 const unsigned int n_components,
1892 const EvaluationFlags::EvaluationFlags evaluation_flag,
1893 const Number * values_dofs,
1895 {
1896 constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1897
1898 for (unsigned int c = 0; c < n_components; ++c)
1899 {
1900 if ((evaluation_flag & EvaluationFlags::values) != 0u)
1901 for (unsigned int i = 0; i < n_points; ++i)
1902 fe_eval.begin_values()[n_points * c + i] =
1903 values_dofs[n_points * c + i];
1904
1905 do_evaluate(fe_eval.get_shape_info().data.front(),
1906 evaluation_flag,
1907 values_dofs + c * n_points,
1908 fe_eval.begin_gradients() + c * dim * n_points,
1909 fe_eval.begin_hessians() +
1910 c * dim * (dim + 1) / 2 * n_points);
1911 }
1912 }
1913
1914
1915
1916 template <int dim, int fe_degree, typename Number>
1917 inline void
1920 const EvaluationFlags::EvaluationFlags evaluation_flag,
1921 const Number * values_dofs,
1922 Number * gradients_quad,
1923 Number * hessians_quad)
1924 {
1926 (fe_degree + 2) / 2 * (fe_degree + 1));
1927 constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1928
1929 Eval eval({},
1932 if ((evaluation_flag &
1934 {
1935 eval.template gradients<0, true, false>(values_dofs, gradients_quad);
1936 if (dim > 1)
1937 eval.template gradients<1, true, false>(values_dofs,
1938 gradients_quad + n_points);
1939 if (dim > 2)
1940 eval.template gradients<2, true, false>(values_dofs,
1941 gradients_quad +
1942 2 * n_points);
1943 }
1944 if (evaluation_flag & EvaluationFlags::hessians)
1945 {
1946 eval.template hessians<0, true, false>(values_dofs, hessians_quad);
1947 if (dim > 1)
1948 {
1949 eval.template gradients<1, true, false>(gradients_quad,
1950 hessians_quad +
1951 dim * n_points);
1952 eval.template hessians<1, true, false>(values_dofs,
1953 hessians_quad + n_points);
1954 }
1955 if (dim > 2)
1956 {
1957 eval.template gradients<2, true, false>(gradients_quad,
1958 hessians_quad +
1959 4 * n_points);
1960 eval.template gradients<2, true, false>(gradients_quad + n_points,
1961 hessians_quad +
1962 5 * n_points);
1963 eval.template hessians<2, true, false>(values_dofs,
1964 hessians_quad +
1965 2 * n_points);
1966 }
1967 }
1968 }
1969
1970
1971
1972 template <int dim, int fe_degree, typename Number>
1973 inline void
1975 const unsigned int n_components,
1976 const EvaluationFlags::EvaluationFlags integration_flag,
1977 Number * values_dofs,
1979 const bool add_into_values_array)
1980 {
1981 constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1982
1983 for (unsigned int c = 0; c < n_components; ++c)
1984 {
1985 if ((integration_flag & EvaluationFlags::values) != 0u)
1986 {
1987 if (add_into_values_array)
1988 for (unsigned int i = 0; i < n_points; ++i)
1989 values_dofs[n_points * c + i] +=
1990 fe_eval.begin_values()[n_points * c + i];
1991 else
1992 for (unsigned int i = 0; i < n_points; ++i)
1993 values_dofs[n_points * c + i] =
1994 fe_eval.begin_values()[n_points * c + i];
1995 }
1996
1997 do_integrate(fe_eval.get_shape_info().data.front(),
1998 integration_flag,
1999 values_dofs + c * n_points,
2000 fe_eval.begin_gradients() + c * dim * n_points,
2001 fe_eval.begin_hessians() +
2002 c * dim * (dim + 1) / 2 * n_points,
2003 add_into_values_array ||
2004 ((integration_flag & EvaluationFlags::values) != 0u));
2005 }
2006 }
2007
2008
2009
2010 template <int dim, int fe_degree, typename Number>
2011 inline void
2014 const EvaluationFlags::EvaluationFlags integration_flag,
2015 Number * values_dofs,
2016 Number * gradients_quad,
2017 const Number * hessians_quad,
2018 const bool add_into_values_array)
2019 {
2021 (fe_degree + 2) / 2 * (fe_degree + 1));
2022
2023 Eval eval({},
2026 constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
2027
2028 if ((integration_flag & EvaluationFlags::hessians) != 0u)
2029 {
2030 // diagonal
2031 // grad xx
2032 if (add_into_values_array == true)
2033 eval.template hessians<0, false, true>(hessians_quad, values_dofs);
2034 else
2035 eval.template hessians<0, false, false>(hessians_quad, values_dofs);
2036 // grad yy
2037 if (dim > 1)
2038 eval.template hessians<1, false, true>(hessians_quad + n_points,
2039 values_dofs);
2040 // grad zz
2041 if (dim > 2)
2042 eval.template hessians<2, false, true>(hessians_quad + 2 * n_points,
2043 values_dofs);
2044 // off-diagonal
2045 if (dim == 2)
2046 {
2047 // grad xy, queue into gradient
2048 if (integration_flag & EvaluationFlags::gradients)
2049 eval.template gradients<1, false, true>(hessians_quad +
2050 2 * n_points,
2051 gradients_quad);
2052 else
2053 eval.template gradients<1, false, false>(hessians_quad +
2054 2 * n_points,
2055 gradients_quad);
2056 }
2057 if (dim == 3)
2058 {
2059 // grad xy, queue into gradient
2060 if (integration_flag & EvaluationFlags::gradients)
2061 eval.template gradients<1, false, true>(hessians_quad +
2062 3 * n_points,
2063 gradients_quad);
2064 else
2065 eval.template gradients<1, false, false>(hessians_quad +
2066 3 * n_points,
2067 gradients_quad);
2068
2069 // grad xz
2070 eval.template gradients<2, false, true>(hessians_quad +
2071 4 * n_points,
2072 gradients_quad);
2073
2074 // grad yz
2075 if (integration_flag & EvaluationFlags::gradients)
2076 eval.template gradients<2, false, true>(
2077 hessians_quad + 5 * n_points, gradients_quad + n_points);
2078 else
2079 eval.template gradients<2, false, false>(
2080 hessians_quad + 5 * n_points, gradients_quad + n_points);
2081 }
2082
2083 // if we did not integrate gradients, set the last slot to zero
2084 // which was not touched before, in order to avoid the if
2085 // statement in the gradients loop below
2086 if ((integration_flag & EvaluationFlags::gradients) == 0u)
2087 for (unsigned int q = 0; q < n_points; ++q)
2088 gradients_quad[(dim - 1) * n_points + q] = Number();
2089 }
2090
2091 if ((integration_flag &
2093 {
2094 if (add_into_values_array ||
2095 (integration_flag & EvaluationFlags::hessians) != 0u)
2096 eval.template gradients<0, false, true>(gradients_quad, values_dofs);
2097 else
2098 eval.template gradients<0, false, false>(gradients_quad, values_dofs);
2099 if (dim > 1)
2100 eval.template gradients<1, false, true>(gradients_quad + n_points,
2101 values_dofs);
2102 if (dim > 2)
2103 eval.template gradients<2, false, true>(gradients_quad + 2 * n_points,
2104 values_dofs);
2105 }
2106 }
2107
2108
2109
2120 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2122 {
2123 static void
2124 evaluate(const unsigned int n_components,
2125 const EvaluationFlags::EvaluationFlags evaluation_flag,
2126 const Number * values_dofs,
2128
2129 static void
2130 integrate(const unsigned int n_components,
2131 const EvaluationFlags::EvaluationFlags evaluation_flag,
2132 Number * values_dofs,
2134 const bool add_into_values_array);
2135 };
2136
2137
2138
2139 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2140 inline void
2142 dim,
2143 fe_degree,
2144 n_q_points_1d,
2145 Number>::evaluate(const unsigned int n_components,
2146 const EvaluationFlags::EvaluationFlags evaluation_flag,
2147 const Number * values_dofs,
2149 {
2150 const auto &shape_data = fe_eval.get_shape_info().data.front();
2151
2152 Assert(n_q_points_1d > fe_degree,
2153 ExcMessage("You lose information when going to a collocation space "
2154 "of lower degree, so the evaluation results would be "
2155 "wrong. Thus, this class does not permit the desired "
2156 "operation."));
2157 constexpr std::size_t n_dofs = Utilities::pow(fe_degree + 1, dim);
2158 constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
2159
2160 for (unsigned int c = 0; c < n_components; ++c)
2161 {
2165 dim,
2166 (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
2167 n_q_points_1d>::do_forward(1,
2168 shape_data.shape_values_eo,
2169 values_dofs + c * n_dofs,
2170 fe_eval.begin_values() + c * n_q_points);
2171
2172 // apply derivatives in the collocation space
2173 if (evaluation_flag &
2176 do_evaluate(shape_data,
2177 evaluation_flag & (EvaluationFlags::gradients |
2179 fe_eval.begin_values() + c * n_q_points,
2180 fe_eval.begin_gradients() + c * dim * n_q_points,
2181 fe_eval.begin_hessians() +
2182 c * dim * (dim + 1) / 2 * n_q_points);
2183 }
2184 }
2185
2186
2187
2188 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2189 inline void
2191 dim,
2192 fe_degree,
2193 n_q_points_1d,
2194 Number>::integrate(const unsigned int n_components,
2195 const EvaluationFlags::EvaluationFlags integration_flag,
2196 Number * values_dofs,
2198 const bool add_into_values_array)
2199 {
2200 const auto &shape_data = fe_eval.get_shape_info().data.front();
2201
2202 Assert(n_q_points_1d > fe_degree,
2203 ExcMessage("You lose information when going to a collocation space "
2204 "of lower degree, so the evaluation results would be "
2205 "wrong. Thus, this class does not permit the desired "
2206 "operation."));
2207 constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
2208
2209 for (unsigned int c = 0; c < n_components; ++c)
2210 {
2211 // apply derivatives in collocation space
2212 if (integration_flag &
2215 do_integrate(shape_data,
2216 integration_flag & (EvaluationFlags::gradients |
2218 fe_eval.begin_values() + c * n_q_points,
2219 fe_eval.begin_gradients() + c * dim * n_q_points,
2220 fe_eval.begin_hessians() +
2221 c * dim * (dim + 1) / 2 * n_q_points,
2222 /*add_into_values_array=*/
2223 integration_flag & EvaluationFlags::values);
2224
2225 // transform back to the original space
2229 dim,
2230 (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
2231 n_q_points_1d>::do_backward(1,
2232 shape_data.shape_values_eo,
2233 add_into_values_array,
2234 fe_eval.begin_values() + c * n_q_points,
2235 values_dofs +
2236 c * Utilities::pow(fe_degree + 1, dim));
2237 }
2238 }
2239
2240
2241
2249 constexpr bool
2250 use_collocation_evaluation(const unsigned int fe_degree,
2251 const unsigned int n_q_points_1d)
2252 {
2253 return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
2254 (n_q_points_1d <= 3 * fe_degree / 2 + 1);
2255 }
2256
2257
2273 template <int dim, typename Number, bool do_integrate>
2275 {
2276 template <int fe_degree, int n_q_points_1d, typename OtherNumber>
2277 static bool
2278 run(const unsigned int n_components,
2279 const EvaluationFlags::EvaluationFlags evaluation_flag,
2280 OtherNumber * values_dofs,
2282 const bool sum_into_values_array = false)
2283 {
2284 // `OtherNumber` is either `const Number` (evaluate()) or `Number`
2285 // (integrate())
2286 static_assert(
2287 std::is_same<Number,
2288 typename std::remove_const<OtherNumber>::type>::value,
2289 "Type of Number and of OtherNumber do not match.");
2290
2291 const auto element_type = fe_eval.get_shape_info().element_type;
2292 using ElementType = MatrixFreeFunctions::ElementType;
2293
2294 Assert(fe_eval.get_shape_info().data.size() == 1 ||
2295 (fe_eval.get_shape_info().data.size() == dim &&
2296 element_type == ElementType::tensor_general) ||
2297 element_type == ElementType::tensor_raviart_thomas,
2299
2300 if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
2301 element_type == ElementType::tensor_symmetric_collocation)
2302 {
2305 n_components,
2306 evaluation_flag,
2307 values_dofs,
2308 fe_eval,
2309 sum_into_values_array);
2310 }
2311 // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
2312 // shape_info.h for more details
2313 else if (fe_degree >= 0 &&
2314 use_collocation_evaluation(fe_degree, n_q_points_1d) &&
2315 element_type <= ElementType::tensor_symmetric)
2316 {
2319 fe_degree,
2320 n_q_points_1d,
2321 Number>>(
2322 n_components,
2323 evaluation_flag,
2324 values_dofs,
2325 fe_eval,
2326 sum_into_values_array);
2327 }
2328 else if (fe_degree >= 0 &&
2329 element_type <= ElementType::tensor_symmetric_no_collocation)
2330 {
2331 evaluate_or_integrate<FEEvaluationImpl<ElementType::tensor_symmetric,
2332 dim,
2333 fe_degree,
2334 n_q_points_1d,
2335 Number>>(
2336 n_components,
2337 evaluation_flag,
2338 values_dofs,
2339 fe_eval,
2340 sum_into_values_array);
2341 }
2342 else if (element_type == ElementType::tensor_symmetric_plus_dg0)
2343 {
2345 FEEvaluationImpl<ElementType::tensor_symmetric_plus_dg0,
2346 dim,
2347 fe_degree,
2348 n_q_points_1d,
2349 Number>>(n_components,
2350 evaluation_flag,
2351 values_dofs,
2352 fe_eval,
2353 sum_into_values_array);
2354 }
2355 else if (element_type == ElementType::truncated_tensor)
2356 {
2357 evaluate_or_integrate<FEEvaluationImpl<ElementType::truncated_tensor,
2358 dim,
2359 fe_degree,
2360 n_q_points_1d,
2361 Number>>(
2362 n_components,
2363 evaluation_flag,
2364 values_dofs,
2365 fe_eval,
2366 sum_into_values_array);
2367 }
2368 else if (element_type == ElementType::tensor_none)
2369 {
2370 evaluate_or_integrate<FEEvaluationImpl<ElementType::tensor_none,
2371 dim,
2372 fe_degree,
2373 n_q_points_1d,
2374 Number>>(
2375 n_components,
2376 evaluation_flag,
2377 values_dofs,
2378 fe_eval,
2379 sum_into_values_array);
2380 }
2381 else if (element_type == ElementType::tensor_raviart_thomas)
2382 {
2383 FEEvaluationImpl<ElementType::tensor_raviart_thomas,
2384 dim,
2385 (fe_degree == -1) ? 1 : fe_degree,
2386 (n_q_points_1d < 1) ? 1 : n_q_points_1d,
2387 Number>::
2388 template evaluate_or_integrate<do_integrate>(evaluation_flag,
2389 const_cast<Number *>(
2390 values_dofs),
2391 fe_eval,
2392 sum_into_values_array);
2393 }
2394 else
2395 {
2396 evaluate_or_integrate<FEEvaluationImpl<ElementType::tensor_general,
2397 dim,
2398 fe_degree,
2399 n_q_points_1d,
2400 Number>>(
2401 n_components,
2402 evaluation_flag,
2403 values_dofs,
2404 fe_eval,
2405 sum_into_values_array);
2406 }
2407
2408 return false;
2409 }
2410
2411 private:
2412 template <typename T>
2413 static void
2415 const unsigned int n_components,
2416 const EvaluationFlags::EvaluationFlags evaluation_flag,
2417 const Number * values_dofs,
2419 const bool sum_into_values_array,
2420 std::integral_constant<bool, false>)
2421 {
2422 (void)sum_into_values_array;
2423
2424 T::evaluate(n_components, evaluation_flag, values_dofs, fe_eval);
2425 }
2426
2427 template <typename T>
2428 static void
2430 const unsigned int n_components,
2431 const EvaluationFlags::EvaluationFlags evaluation_flag,
2432 Number * values_dofs,
2434 const bool sum_into_values_array,
2435 std::integral_constant<bool, true>)
2436 {
2437 T::integrate(n_components,
2438 evaluation_flag,
2439 values_dofs,
2440 fe_eval,
2441 sum_into_values_array);
2442 }
2443
2444 template <typename T, typename OtherNumber>
2445 static void
2447 const unsigned int n_components,
2448 const EvaluationFlags::EvaluationFlags evaluation_flag,
2449 OtherNumber * values_dofs,
2451 const bool sum_into_values_array)
2452 {
2453 evaluate_or_integrate<T>(n_components,
2454 evaluation_flag,
2455 values_dofs,
2456 fe_eval,
2457 sum_into_values_array,
2458 std::integral_constant<bool, do_integrate>());
2459 }
2460 };
2461
2462
2463
2464 template <bool symmetric_evaluate,
2465 int dim,
2466 int fe_degree,
2467 int n_q_points_1d,
2468 typename Number>
2470 {
2471 // We enable a transformation to collocation for derivatives if it gives
2472 // correct results (first two conditions), if it is the most efficient
2473 // choice in terms of operation counts (third condition) and if we were
2474 // able to initialize the fields in shape_info.templates.h from the
2475 // polynomials (fourth condition).
2476 using Number2 =
2478
2479 using Eval = EvaluatorTensorProduct<symmetric_evaluate ? evaluate_evenodd :
2481 dim - 1,
2482 fe_degree + 1,
2483 n_q_points_1d,
2484 Number,
2485 Number2>;
2486
2487 static Eval
2490 const unsigned int subface_index,
2491 const unsigned int direction)
2492 {
2493 if (symmetric_evaluate)
2494 return Eval(data.shape_values_eo,
2495 data.shape_gradients_eo,
2496 data.shape_hessians_eo,
2497 data.fe_degree + 1,
2498 data.n_q_points_1d);
2499 else if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
2500 return Eval(data.shape_values,
2501 data.shape_gradients,
2502 data.shape_hessians,
2503 data.fe_degree + 1,
2504 data.n_q_points_1d);
2505 else
2506 {
2507 const unsigned int index =
2508 direction == 0 ? subface_index % 2 : subface_index / 2;
2509 return Eval(data.values_within_subface[index],
2512 data.fe_degree + 1,
2513 data.n_q_points_1d);
2514 }
2515 }
2516
2517 static void
2519 const unsigned int n_components,
2520 const EvaluationFlags::EvaluationFlags evaluation_flag,
2522 Number * values_dofs,
2523 Number * values_quad,
2524 Number * gradients_quad,
2525 Number * hessians_quad,
2526 Number * scratch_data,
2527 const unsigned int subface_index)
2528 {
2529 Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
2530 Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
2531
2532 const std::size_t n_dofs = fe_degree > -1 ?
2533 Utilities::pow(fe_degree + 1, dim - 1) :
2534 Utilities::pow(data.fe_degree + 1, dim - 1);
2535 const std::size_t n_q_points =
2536 fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
2537 Utilities::pow(data.n_q_points_1d, dim - 1);
2538
2539 // keep a copy of the original pointer for the case of the Hessians
2540 Number *values_dofs_ptr = values_dofs;
2541
2542 if ((evaluation_flag & EvaluationFlags::values) != 0u &&
2543 ((evaluation_flag & EvaluationFlags::gradients) == 0u))
2544 for (unsigned int c = 0; c < n_components; ++c)
2545 {
2546 switch (dim)
2547 {
2548 case 3:
2549 eval0.template values<0, true, false>(values_dofs,
2550 values_quad);
2551 eval1.template values<1, true, false>(values_quad,
2552 values_quad);
2553 break;
2554 case 2:
2555 eval0.template values<0, true, false>(values_dofs,
2556 values_quad);
2557 break;
2558 case 1:
2559 values_quad[0] = values_dofs[0];
2560 break;
2561 default:
2562 Assert(false, ExcNotImplemented());
2563 }
2564 // Note: we always keep storage of values, 1st and 2nd derivatives
2565 // in an array
2566 values_dofs += 3 * n_dofs;
2567 values_quad += n_q_points;
2568 }
2569 else if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
2570 for (unsigned int c = 0; c < n_components; ++c)
2571 {
2572 switch (dim)
2573 {
2574 case 3:
2575 if (symmetric_evaluate &&
2576 use_collocation_evaluation(fe_degree, n_q_points_1d))
2577 {
2578 eval0.template values<0, true, false>(values_dofs,
2579 values_quad);
2580 eval0.template values<1, true, false>(values_quad,
2581 values_quad);
2583 dim - 1,
2584 n_q_points_1d,
2585 n_q_points_1d,
2586 Number,
2587 Number2>
2588 eval_grad({}, data.shape_gradients_collocation_eo, {});
2589 eval_grad.template gradients<0, true, false>(
2590 values_quad, gradients_quad);
2591 eval_grad.template gradients<1, true, false>(
2592 values_quad, gradients_quad + n_q_points);
2593 }
2594 else
2595 {
2596 // grad x
2597 eval0.template gradients<0, true, false>(values_dofs,
2598 scratch_data);
2599 eval1.template values<1, true, false>(scratch_data,
2600 gradients_quad);
2601
2602 // grad y
2603 eval0.template values<0, true, false>(values_dofs,
2604 scratch_data);
2605 eval1.template gradients<1, true, false>(scratch_data,
2606 gradients_quad +
2607 n_q_points);
2608
2609 if ((evaluation_flag & EvaluationFlags::values) != 0u)
2610 eval1.template values<1, true, false>(scratch_data,
2611 values_quad);
2612 }
2613 // grad z
2614 eval0.template values<0, true, false>(values_dofs + n_dofs,
2615 scratch_data);
2616 eval1.template values<1, true, false>(
2617 scratch_data, gradients_quad + (dim - 1) * n_q_points);
2618
2619 break;
2620 case 2:
2621 eval0.template values<0, true, false>(values_dofs + n_dofs,
2622 gradients_quad +
2623 n_q_points);
2624 eval0.template gradients<0, true, false>(values_dofs,
2625 gradients_quad);
2626 if ((evaluation_flag & EvaluationFlags::values) != 0u)
2627 eval0.template values<0, true, false>(values_dofs,
2628 values_quad);
2629 break;
2630 case 1:
2631 values_quad[0] = values_dofs[0];
2632 gradients_quad[0] = values_dofs[1];
2633 break;
2634 default:
2636 }
2637 values_dofs += 3 * n_dofs;
2638 values_quad += n_q_points;
2639 gradients_quad += dim * n_q_points;
2640 }
2641
2642 if ((evaluation_flag & EvaluationFlags::hessians) != 0u)
2643 {
2644 values_dofs = values_dofs_ptr;
2645 for (unsigned int c = 0; c < n_components; ++c)
2646 {
2647 switch (dim)
2648 {
2649 case 3:
2650 // grad xx
2651 eval0.template hessians<0, true, false>(values_dofs,
2652 scratch_data);
2653 eval1.template values<1, true, false>(scratch_data,
2654 hessians_quad);
2655
2656 // grad yy
2657 eval0.template values<0, true, false>(values_dofs,
2658 scratch_data);
2659 eval1.template hessians<1, true, false>(scratch_data,
2660 hessians_quad +
2661 n_q_points);
2662
2663 // grad zz
2664 eval0.template values<0, true, false>(values_dofs +
2665 2 * n_dofs,
2666 scratch_data);
2667 eval1.template values<1, true, false>(scratch_data,
2668 hessians_quad +
2669 2 * n_q_points);
2670
2671 // grad xy
2672 eval0.template gradients<0, true, false>(values_dofs,
2673 scratch_data);
2674 eval1.template gradients<1, true, false>(scratch_data,
2675 hessians_quad +
2676 3 * n_q_points);
2677
2678 // grad xz
2679 eval0.template gradients<0, true, false>(values_dofs +
2680 n_dofs,
2681 scratch_data);
2682 eval1.template values<1, true, false>(scratch_data,
2683 hessians_quad +
2684 4 * n_q_points);
2685
2686 // grad yz
2687 eval0.template values<0, true, false>(values_dofs + n_dofs,
2688 scratch_data);
2689 eval1.template gradients<1, true, false>(scratch_data,
2690 hessians_quad +
2691 5 * n_q_points);
2692
2693 break;
2694 case 2:
2695 // grad xx
2696 eval0.template hessians<0, true, false>(values_dofs,
2697 hessians_quad);
2698 // grad yy
2699 eval0.template values<0, true, false>(
2700 values_dofs + 2 * n_dofs, hessians_quad + n_q_points);
2701 // grad xy
2702 eval0.template gradients<0, true, false>(
2703 values_dofs + n_dofs, hessians_quad + 2 * n_q_points);
2704 break;
2705 case 1:
2706 hessians_quad[0] = values_dofs[2];
2707 break;
2708 default:
2710 }
2711 values_dofs += 3 * n_dofs;
2712 hessians_quad += dim * (dim + 1) / 2 * n_q_points;
2713 }
2714 }
2715 }
2716
2717 static void
2719 const unsigned int n_components,
2720 const EvaluationFlags::EvaluationFlags integration_flag,
2722 Number * values_dofs,
2723 Number * values_quad,
2724 Number * gradients_quad,
2725 Number * hessians_quad,
2726 Number * scratch_data,
2727 const unsigned int subface_index)
2728 {
2729 Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
2730 Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
2731
2732 const std::size_t n_dofs =
2733 fe_degree > -1 ?
2734 Utilities::pow(fe_degree + 1, dim - 1) :
2735 (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
2736 const std::size_t n_q_points =
2737 fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) :
2738 Utilities::pow(data.n_q_points_1d, dim - 1);
2739
2740 // keep a copy of the original pointer for the case of the Hessians
2741 Number *values_dofs_ptr = values_dofs;
2742
2743 if ((integration_flag & EvaluationFlags::values) != 0u &&
2744 (integration_flag & EvaluationFlags::gradients) == 0u)
2745 for (unsigned int c = 0; c < n_components; ++c)
2746 {
2747 switch (dim)
2748 {
2749 case 3:
2750 eval1.template values<1, false, false>(values_quad,
2751 values_quad);
2752 eval0.template values<0, false, false>(values_quad,
2753 values_dofs);
2754 break;
2755 case 2:
2756 eval0.template values<0, false, false>(values_quad,
2757 values_dofs);
2758 break;
2759 case 1:
2760 values_dofs[0] = values_quad[0];
2761 break;
2762 default:
2763 Assert(false, ExcNotImplemented());
2764 }
2765 values_dofs += 3 * n_dofs;
2766 values_quad += n_q_points;
2767 }
2768 else if ((integration_flag & EvaluationFlags::gradients) != 0u)
2769 for (unsigned int c = 0; c < n_components; ++c)
2770 {
2771 switch (dim)
2772 {
2773 case 3:
2774 // grad z
2775 eval1.template values<1, false, false>(gradients_quad +
2776 2 * n_q_points,
2777 gradients_quad +
2778 2 * n_q_points);
2779 eval0.template values<0, false, false>(gradients_quad +
2780 2 * n_q_points,
2781 values_dofs + n_dofs);
2782 if (symmetric_evaluate &&
2783 use_collocation_evaluation(fe_degree, n_q_points_1d))
2784 {
2786 dim - 1,
2787 n_q_points_1d,
2788 n_q_points_1d,
2789 Number,
2790 Number2>
2791 eval_grad({}, data.shape_gradients_collocation_eo, {});
2792 if ((integration_flag & EvaluationFlags::values) != 0u)
2793 eval_grad.template gradients<1, false, true>(
2794 gradients_quad + n_q_points, values_quad);
2795 else
2796 eval_grad.template gradients<1, false, false>(
2797 gradients_quad + n_q_points, values_quad);
2798 eval_grad.template gradients<0, false, true>(
2799 gradients_quad, values_quad);
2800 eval0.template values<1, false, false>(values_quad,
2801 values_quad);
2802 eval0.template values<0, false, false>(values_quad,
2803 values_dofs);
2804 }
2805 else
2806 {
2807 if ((integration_flag & EvaluationFlags::values) != 0u)
2808 {
2809 eval1.template values<1, false, false>(values_quad,
2810 scratch_data);
2811 eval1.template gradients<1, false, true>(
2812 gradients_quad + n_q_points, scratch_data);
2813 }
2814 else
2815 eval1.template gradients<1, false, false>(
2816 gradients_quad + n_q_points, scratch_data);
2817
2818 // grad y
2819 eval0.template values<0, false, false>(scratch_data,
2820 values_dofs);
2821
2822 // grad x
2823 eval1.template values<1, false, false>(gradients_quad,
2824 scratch_data);
2825 eval0.template gradients<0, false, true>(scratch_data,
2826 values_dofs);
2827 }
2828 break;
2829 case 2:
2830 eval0.template values<0, false, false>(gradients_quad +
2831 n_q_points,
2832 values_dofs + n_dofs);
2833 eval0.template gradients<0, false, false>(gradients_quad,
2834 values_dofs);
2835 if ((integration_flag & EvaluationFlags::values) != 0u)
2836 eval0.template values<0, false, true>(values_quad,
2837 values_dofs);
2838 break;
2839 case 1:
2840 values_dofs[0] = values_quad[0];
2841 values_dofs[1] = gradients_quad[0];
2842 break;
2843 default:
2845 }
2846 values_dofs += 3 * n_dofs;
2847 values_quad += n_q_points;
2848 gradients_quad += dim * n_q_points;
2849 }
2850
2851 if ((integration_flag & EvaluationFlags::hessians) != 0u)
2852 {
2853 values_dofs = values_dofs_ptr;
2854 for (unsigned int c = 0; c < n_components; ++c)
2855 {
2856 switch (dim)
2857 {
2858 case 3:
2859 // grad xx
2860 eval1.template values<1, false, false>(hessians_quad,
2861 scratch_data);
2862 if ((integration_flag & (EvaluationFlags::values |
2864 eval0.template hessians<0, false, true>(scratch_data,
2865 values_dofs);
2866 else
2867 eval0.template hessians<0, false, false>(scratch_data,
2868 values_dofs);
2869
2870 // grad yy
2871 eval1.template hessians<1, false, false>(hessians_quad +
2872 n_q_points,
2873 scratch_data);
2874 eval0.template values<0, false, true>(scratch_data,
2875 values_dofs);
2876
2877 // grad zz
2878 eval1.template values<1, false, false>(hessians_quad +
2879 2 * n_q_points,
2880 scratch_data);
2881 eval0.template values<0, false, false>(scratch_data,
2882 values_dofs +
2883 2 * n_dofs);
2884
2885 // grad xy
2886 eval1.template gradients<1, false, false>(hessians_quad +
2887 3 * n_q_points,
2888 scratch_data);
2889 eval0.template gradients<0, false, true>(scratch_data,
2890 values_dofs);
2891
2892 // grad xz
2893 eval1.template values<1, false, false>(hessians_quad +
2894 4 * n_q_points,
2895 scratch_data);
2896 if ((integration_flag & EvaluationFlags::gradients) != 0u)
2897 eval0.template gradients<0, false, true>(scratch_data,
2898 values_dofs +
2899 n_dofs);
2900 else
2901 eval0.template gradients<0, false, false>(scratch_data,
2902 values_dofs +
2903 n_dofs);
2904
2905 // grad yz
2906 eval1.template gradients<1, false, false>(hessians_quad +
2907 5 * n_q_points,
2908 scratch_data);
2909 eval0.template values<0, false, true>(scratch_data,
2910 values_dofs + n_dofs);
2911
2912 break;
2913 case 2:
2914 // grad xx
2915 if ((integration_flag & (EvaluationFlags::values |
2917 eval0.template hessians<0, false, true>(hessians_quad,
2918 values_dofs);
2919 else
2920 eval0.template hessians<0, false, false>(hessians_quad,
2921 values_dofs);
2922
2923 // grad yy
2924 eval0.template values<0, false, false>(
2925 hessians_quad + n_q_points, values_dofs + 2 * n_dofs);
2926 // grad xy
2927 if ((integration_flag & EvaluationFlags::gradients) != 0u)
2928 eval0.template gradients<0, false, true>(
2929 hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
2930 else
2931 eval0.template gradients<0, false, false>(
2932 hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
2933 break;
2934 case 1:
2935 values_dofs[2] = hessians_quad[0];
2936 if ((integration_flag & EvaluationFlags::values) == 0u)
2937 values_dofs[0] = 0;
2938 if ((integration_flag & EvaluationFlags::gradients) == 0u)
2939 values_dofs[1] = 0;
2940 break;
2941 default:
2943 }
2944 values_dofs += 3 * n_dofs;
2945 hessians_quad += dim * (dim + 1) / 2 * n_q_points;
2946 }
2947 }
2948 }
2949 };
2950
2951 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
2953 {
2954 using Number2 =
2957 dim - 1,
2958 fe_degree,
2959 n_q_points_1d,
2960 Number,
2961 Number2>;
2962
2963 template <typename EvalType>
2964 static EvalType
2967 const unsigned int subface_index,
2968 const unsigned int direction)
2969 {
2970 if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
2971 return EvalType(data.shape_values,
2972 data.shape_gradients,
2973 data.shape_hessians);
2974 else
2975 {
2976 const unsigned int index =
2977 direction == 0 ? subface_index % 2 : subface_index / 2;
2978 return EvalType(data.values_within_subface[index],
2981 }
2982 }
2983
2984 template <bool integrate>
2985 static void
2987 const EvaluationFlags::EvaluationFlags evaluation_flag,
2988 Number * values_dofs,
2990 Number * scratch_data,
2991 const unsigned int subface_index,
2992 const unsigned int face_no)
2993 {
2994 const unsigned int face_direction = face_no / 2;
2995
2996 // We first evaluate the anisotropic faces, i.e the faces where
2997 // face_direction != component. Note that the call order here is not
2998 // important, since the pointers are shifted accordingly within the
2999 // function. However, this is the order in which the components will be in
3000 // the quadrature points. Furthermore, the isotropic faces have no "normal
3001 // direction" but we still pass in normal_dir = 2 since this is used for
3002 // the pointers.
3003 // -----------------------------------------------------------------------------------
3004 // | | Anisotropic faces | Isotropic faces|
3005 // | Face dir | comp, coords, normal_dir | comp, coords, normal_dir | comp, coords |
3006 // | --------------------------------------------------------------------------------|
3007 // | 0 | 1, y, 0 | - | 0, y |
3008 // | 1 | 0, x, 0 | - | 1, x |
3009 // | --------------------------------------------------------------------------------|
3010 // | 0 | 1, yz, 0 | 2, yz, 1 | 0, yz |
3011 // | 1 | 2, zx, 0 | 0, zx, 1 | 1, zx |
3012 // | 2 | 0, xy, 0 | 1, xy, 1 | 2, xy |
3013 // -----------------------------------------------------------------------------------
3014 evaluate_in_face_apply<0>(values_dofs,
3015 fe_eval,
3016 scratch_data,
3017 evaluation_flag,
3018 face_direction,
3019 subface_index,
3020 std::integral_constant<bool, integrate>());
3021
3022 if (dim == 3)
3023 evaluate_in_face_apply<1>(values_dofs,
3024 fe_eval,
3025 scratch_data,
3026 evaluation_flag,
3027 face_direction,
3028 subface_index,
3029 std::integral_constant<bool, integrate>());
3030
3031 evaluate_in_face_apply<2>(values_dofs,
3032 fe_eval,
3033 scratch_data,
3034 evaluation_flag,
3035 face_direction,
3036 subface_index,
3037 std::integral_constant<bool, integrate>());
3038 }
3039
3040 /*
3041 * Helper function which applies the 1d kernels for on one
3042 * component in a face. normal_dir indicates the direction of the continuous
3043 * component of the RT space. std::integral_constant<bool, false> is the
3044 * evaluation path, and std::integral_constant<bool, true> below is the
3045 * integration path. These two functions can be fused together since all
3046 * offsets and pointers are the exact same.
3047 */
3048 template <int normal_dir>
3049 static inline void
3051 Number * values_dofs,
3053 Number * scratch_data,
3054 const EvaluationFlags::EvaluationFlags evaluation_flag,
3055 const unsigned int face_direction,
3056 const unsigned int subface_index,
3057 std::integral_constant<bool, false>)
3058 {
3059 using EvalNormal =
3061 dim - 1,
3062 (fe_degree == -1) ? 1 : fe_degree + 1,
3063 n_q_points_1d,
3064 normal_dir,
3065 Number,
3066 Number2>;
3067 using EvalTangent =
3069 dim - 1,
3070 (fe_degree == -1) ? 1 : fe_degree,
3071 n_q_points_1d,
3072 normal_dir,
3073 Number,
3074 Number2>;
3075
3076 using TempEval0 = typename std::
3077 conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
3078 using TempEval1 = typename std::
3079 conditional<normal_dir == 0, EvalTangent, EvalNormal>::type;
3080 using Eval0 = typename std::
3081 conditional<normal_dir == 2, EvalGeneral, TempEval0>::type;
3082 using Eval1 = typename std::
3083 conditional<normal_dir == 2, EvalGeneral, TempEval1>::type;
3084
3085 const auto &shape_info = fe_eval.get_shape_info();
3086 Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
3087 ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]),
3088 subface_index,
3089 0);
3090 Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
3091 ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]),
3092 subface_index,
3093 1);
3094
3095 constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
3096 const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face;
3097 const std::size_t n_dofs_normal =
3098 n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
3099 const std::size_t dofs_stride =
3100 (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
3101 n_dofs_tangent;
3102
3103 static constexpr ::ndarray<unsigned int, 3, 3> component_table = {
3104 {{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
3105 const unsigned int component =
3106 (dim == 2 && normal_dir == 0 && face_direction == 1) ?
3107 0 :
3108 component_table[face_direction][normal_dir];
3109
3110 // Initial offsets
3111 values_dofs +=
3112 3 * ((component == 0) ?
3113 0 :
3114 ((component == 1) ?
3115 ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) :
3116 ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent :
3117 n_dofs_normal + n_dofs_tangent)));
3118 const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir;
3119 Number *values_quad = fe_eval.begin_values() + n_q_points * shift;
3120 Number *gradients_quad =
3121 fe_eval.begin_gradients() + dim * n_q_points * shift;
3122 Number *hessians_quad =
3123 fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * shift;
3124
3125 // Evaluation path
3126 if ((evaluation_flag & EvaluationFlags::values) &&
3127 !(evaluation_flag & EvaluationFlags::gradients))
3128 {
3129 switch (dim)
3130 {
3131 case 3:
3132 eval0.template values<0, true, false>(values_dofs, values_quad);
3133 eval1.template values<1, true, false>(values_quad, values_quad);
3134 break;
3135 case 2:
3136 eval0.template values<0, true, false>(values_dofs, values_quad);
3137 break;
3138 default:
3139 Assert(false, ExcNotImplemented());
3140 }
3141 }
3142 else if (evaluation_flag & EvaluationFlags::gradients)
3143 {
3144 switch (dim)
3145 {
3146 case 3:
3147 // grad x
3148 eval0.template gradients<0, true, false>(values_dofs,
3149 scratch_data);
3150 eval1.template values<1, true, false>(scratch_data,
3151 gradients_quad);
3152
3153 // grad y
3154 eval0.template values<0, true, false>(values_dofs,
3155 scratch_data);
3156 eval1.template gradients<1, true, false>(scratch_data,
3157 gradients_quad +
3158 n_q_points);
3159
3160 if (evaluation_flag & EvaluationFlags::values)
3161 eval1.template values<1, true, false>(scratch_data,
3162 values_quad);
3163
3164 // grad z
3165 eval0.template values<0, true, false>(values_dofs + dofs_stride,
3166 scratch_data);
3167 eval1.template values<1, true, false>(scratch_data,
3168 gradients_quad +
3169 2 * n_q_points);
3170
3171 break;
3172 case 2:
3173 eval0.template values<0, true, false>(values_dofs + dofs_stride,
3174 gradients_quad +
3175 n_q_points);
3176 eval0.template gradients<0, true, false>(values_dofs,
3177 gradients_quad);
3178 if ((evaluation_flag & EvaluationFlags::values))
3179 eval0.template values<0, true, false>(values_dofs,
3180 values_quad);
3181 break;
3182 default:
3184 }
3185 }
3186
3187 if (evaluation_flag & EvaluationFlags::hessians)
3188 {
3189 switch (dim)
3190 {
3191 case 3:
3192 // grad xx
3193 eval0.template hessians<0, true, false>(values_dofs,
3194 scratch_data);
3195 eval1.template values<1, true, false>(scratch_data,
3196 hessians_quad);
3197
3198 // grad yy
3199 eval0.template values<0, true, false>(values_dofs,
3200 scratch_data);
3201 eval1.template hessians<1, true, false>(scratch_data,
3202 hessians_quad +
3203 n_q_points);
3204
3205 // grad zz
3206 eval0.template values<0, true, false>(values_dofs +
3207 2 * dofs_stride,
3208 scratch_data);
3209 eval1.template values<1, true, false>(scratch_data,
3210 hessians_quad +
3211 2 * n_q_points);
3212
3213 // grad xy
3214 eval0.template gradients<0, true, false>(values_dofs,
3215 scratch_data);
3216 eval1.template gradients<1, true, false>(scratch_data,
3217 hessians_quad +
3218 3 * n_q_points);
3219
3220 // grad xz
3221 eval0.template gradients<0, true, false>(values_dofs +
3222 dofs_stride,
3223 scratch_data);
3224 eval1.template values<1, true, false>(scratch_data,
3225 hessians_quad +
3226 4 * n_q_points);
3227
3228 // grad yz
3229 eval0.template values<0, true, false>(values_dofs + dofs_stride,
3230 scratch_data);
3231 eval1.template gradients<1, true, false>(scratch_data,
3232 hessians_quad +
3233 5 * n_q_points);
3234
3235 break;
3236 case 2:
3237 // grad xx
3238 eval0.template hessians<0, true, false>(values_dofs,
3239 hessians_quad);
3240 // grad yy
3241 eval0.template values<0, true, false>(
3242 values_dofs + 2 * dofs_stride, hessians_quad + n_q_points);
3243 // grad xy
3244 eval0.template gradients<0, true, false>(
3245 values_dofs + dofs_stride, hessians_quad + 2 * n_q_points);
3246 break;
3247 default:
3249 }
3250 }
3251 }
3252
3253 template <int normal_dir>
3254 static inline void
3256 Number * values_dofs,
3258 Number * scratch_data,
3259 const EvaluationFlags::EvaluationFlags evaluation_flag,
3260 const unsigned int face_direction,
3261 const unsigned int subface_index,
3262 std::integral_constant<bool, true>)
3263 {
3264 using EvalNormal =
3266 dim - 1,
3267 (fe_degree == -1) ? 1 : fe_degree + 1,
3268 n_q_points_1d,
3269 normal_dir,
3270 Number,
3271 Number2>;
3272 using EvalTangent =
3274 dim - 1,
3275 (fe_degree == -1) ? 1 : fe_degree,
3276 n_q_points_1d,
3277 normal_dir,
3278 Number,
3279 Number2>;
3280
3281 using TempEval0 = typename std::
3282 conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
3283 using TempEval1 = typename std::
3284 conditional<normal_dir == 0, EvalTangent, EvalNormal>::type;
3285 using Eval0 = typename std::
3286 conditional<normal_dir == 2, EvalGeneral, TempEval0>::type;
3287 using Eval1 = typename std::
3288 conditional<normal_dir == 2, EvalGeneral, TempEval1>::type;
3289
3290 const auto &shape_info = fe_eval.get_shape_info();
3291 Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
3292 ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]),
3293 subface_index,
3294 0);
3295 Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
3296 ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]),
3297 subface_index,
3298 1);
3299
3300 constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
3301 const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face;
3302 const std::size_t n_dofs_normal =
3303 n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
3304 const std::size_t dofs_stride =
3305 (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
3306 n_dofs_tangent;
3307
3308 static constexpr ::ndarray<unsigned int, 3, 3> component_table = {
3309 {{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
3310 const unsigned int component =
3311 (dim == 2 && normal_dir == 0 && face_direction == 1) ?
3312 0 :
3313 component_table[face_direction][normal_dir];
3314
3315 // Initial offsets
3316 values_dofs +=
3317 3 * ((component == 0) ?
3318 0 :
3319 ((component == 1) ?
3320 ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) :
3321 ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent :
3322 n_dofs_normal + n_dofs_tangent)));
3323 const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir;
3324 Number *values_quad = fe_eval.begin_values() + n_q_points * shift;
3325 Number *gradients_quad =
3326 fe_eval.begin_gradients() + dim * n_q_points * shift;
3327 Number *hessians_quad =
3328 fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * shift;
3329
3330 // Integration path
3331 if ((evaluation_flag & EvaluationFlags::values) &&
3332 !(evaluation_flag & EvaluationFlags::gradients))
3333 {
3334 switch (dim)
3335 {
3336 case 3:
3337 eval1.template values<1, false, false>(values_quad,
3338 values_quad);
3339 eval0.template values<0, false, false>(values_quad,
3340 values_dofs);
3341 break;
3342 case 2:
3343 eval0.template values<0, false, false>(values_quad,
3344 values_dofs);
3345 break;
3346 default:
3347 Assert(false, ExcNotImplemented());
3348 }
3349 }
3350 else if (evaluation_flag & EvaluationFlags::gradients)
3351 {
3352 switch (dim)
3353 {
3354 case 3:
3355 // grad z
3356 eval1.template values<1, false, false>(gradients_quad +
3357 2 * n_q_points,
3358 gradients_quad +
3359 2 * n_q_points);
3360 eval0.template values<0, false, false>(
3361 gradients_quad + 2 * n_q_points, values_dofs + dofs_stride);
3362
3363 if (evaluation_flag & EvaluationFlags::values)
3364 {
3365 eval1.template values<1, false, false>(values_quad,
3366 scratch_data);
3367 eval1.template gradients<1, false, true>(gradients_quad +
3368 n_q_points,
3369 scratch_data);
3370 }
3371 else
3372 eval1.template gradients<1, false, false>(gradients_quad +
3373 n_q_points,
3374 scratch_data);
3375
3376 // grad y
3377 eval0.template values<0, false, false>(scratch_data,
3378 values_dofs);
3379
3380 // grad x
3381 eval1.template values<1, false, false>(gradients_quad,
3382 scratch_data);
3383 eval0.template gradients<0, false, true>(scratch_data,
3384 values_dofs);
3385
3386 break;
3387 case 2:
3388 eval0.template values<0, false, false>(
3389 gradients_quad + n_q_points, values_dofs + dofs_stride);
3390 eval0.template gradients<0, false, false>(gradients_quad,
3391 values_dofs);
3392 if (evaluation_flag & EvaluationFlags::values)
3393 eval0.template values<0, false, true>(values_quad,
3394 values_dofs);
3395 break;
3396 default:
3398 }
3399 }
3400
3401 if (evaluation_flag & EvaluationFlags::hessians)
3402 {
3403 switch (dim)
3404 {
3405 case 3:
3406 // grad xx
3407 eval1.template values<1, false, false>(hessians_quad,
3408 scratch_data);
3409 if ((evaluation_flag &
3411 eval0.template hessians<0, false, true>(scratch_data,
3412 values_dofs);
3413 else
3414 eval0.template hessians<0, false, false>(scratch_data,
3415 values_dofs);
3416
3417 // grad yy
3418 eval1.template hessians<1, false, false>(hessians_quad +
3419 n_q_points,
3420 scratch_data);
3421 eval0.template values<0, false, true>(scratch_data,
3422 values_dofs);
3423
3424 // grad zz
3425 eval1.template values<1, false, false>(hessians_quad +
3426 2 * n_q_points,
3427 scratch_data);
3428 eval0.template values<0, false, false>(scratch_data,
3429 values_dofs +
3430 2 * dofs_stride);
3431
3432 // grad xy
3433 eval1.template gradients<1, false, false>(hessians_quad +
3434 3 * n_q_points,
3435 scratch_data);
3436 eval0.template gradients<0, false, true>(scratch_data,
3437 values_dofs);
3438
3439 // grad xz
3440 eval1.template values<1, false, false>(hessians_quad +
3441 4 * n_q_points,
3442 scratch_data);
3443 if ((evaluation_flag & EvaluationFlags::gradients))
3444 eval0.template gradients<0, false, true>(scratch_data,
3445 values_dofs +
3446 dofs_stride);
3447 else
3448 eval0.template gradients<0, false, false>(scratch_data,
3449 values_dofs +
3450 dofs_stride);
3451
3452 // grad yz
3453 eval1.template gradients<1, false, false>(hessians_quad +
3454 5 * n_q_points,
3455 scratch_data);
3456 eval0.template values<0, false, true>(scratch_data,
3457 values_dofs +
3458 dofs_stride);
3459
3460 break;
3461 case 2:
3462 // grad xx
3463 if (evaluation_flag &
3465 eval0.template hessians<0, false, true>(hessians_quad,
3466 values_dofs);
3467 else
3468 eval0.template hessians<0, false, false>(hessians_quad,
3469 values_dofs);
3470
3471 // grad yy
3472 eval0.template values<0, false, false>(
3473 hessians_quad + n_q_points, values_dofs + 2 * dofs_stride);
3474 // grad xy
3475 if ((evaluation_flag & EvaluationFlags::gradients))
3476 eval0.template gradients<0, false, true>(
3477 hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
3478 else
3479 eval0.template gradients<0, false, false>(
3480 hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
3481 break;
3482 default:
3484 }
3485 }
3486 }
3487 };
3488
3489
3490 template <int dim, int fe_degree, typename Number>
3492 {
3493 using Number2 =
3495
3496 template <bool do_evaluate, bool add_into_output>
3497 static void
3498 interpolate(const unsigned int n_components,
3501 const Number * input,
3502 Number * output,
3503 const unsigned int face_no)
3504 {
3505 Assert(static_cast<unsigned int>(fe_degree) ==
3506 shape_info.data.front().fe_degree ||
3507 fe_degree == -1,
3510 interpolate_generic_raviart_thomas<do_evaluate, add_into_output>(
3511 n_components, input, output, flags, face_no, shape_info);
3512 else
3513 interpolate_generic<do_evaluate, add_into_output>(
3514 n_components,
3515 input,
3516 output,
3517 flags,
3518 face_no,
3519 shape_info.data.front().fe_degree + 1,
3520 shape_info.data.front().shape_data_on_face,
3521 shape_info.dofs_per_component_on_cell,
3522 3 * shape_info.dofs_per_component_on_face);
3523 }
3524
3528 template <bool do_evaluate, bool add_into_output>
3529 static void
3531 const unsigned int n_components,
3534 const Number * input,
3535 Number * output,
3536 const unsigned int face_no)
3537 {
3538 Assert(static_cast<unsigned int>(fe_degree + 1) ==
3539 shape_info.data.front().n_q_points_1d ||
3540 fe_degree == -1,
3542
3543 interpolate_generic<do_evaluate, add_into_output>(
3544 n_components,
3545 input,
3546 output,
3547 flags,
3548 face_no,
3549 shape_info.data.front().quadrature.size(),
3550 shape_info.data.front().quadrature_data_on_face,
3551 shape_info.n_q_points,
3552 shape_info.n_q_points_face);
3553 }
3554
3555 private:
3556 template <bool do_evaluate, bool add_into_output, int face_direction = 0>
3557 static void
3558 interpolate_generic(const unsigned int n_components,
3559 const Number * input,
3560 Number * output,
3562 const unsigned int face_no,
3563 const unsigned int n_points_1d,
3564 const std::array<AlignedVector<Number2>, 2> &shape_data,
3565 const unsigned int dofs_per_component_on_cell,
3566 const unsigned int dofs_per_component_on_face)
3567 {
3568 if (face_direction == face_no / 2)
3569 {
3571 dim,
3572 fe_degree + 1,
3573 0,
3574 Number,
3575 Number2>
3576 evalf(shape_data[face_no % 2].begin(),
3577 nullptr,
3578 nullptr,
3579 n_points_1d,
3580 0);
3581
3582 const unsigned int in_stride = do_evaluate ?
3583 dofs_per_component_on_cell :
3584 dofs_per_component_on_face;
3585 const unsigned int out_stride = do_evaluate ?
3586 dofs_per_component_on_face :
3587 dofs_per_component_on_cell;
3588
3589 for (unsigned int c = 0; c < n_components; ++c)
3590 {
3591 if (flag & EvaluationFlags::hessians)
3592 evalf.template apply_face<face_direction,
3593 do_evaluate,
3594 add_into_output,
3595 2>(input, output);
3596 else if (flag & EvaluationFlags::gradients)
3597 evalf.template apply_face<face_direction,
3598 do_evaluate,
3599 add_into_output,
3600 1>(input, output);
3601 else
3602 evalf.template apply_face<face_direction,
3603 do_evaluate,
3604 add_into_output,
3605 0>(input, output);
3606 input += in_stride;
3607 output += out_stride;
3608 }
3609 }
3610 else if (face_direction < dim)
3611 {
3612 interpolate_generic<do_evaluate,
3613 add_into_output,
3614 std::min(face_direction + 1, dim - 1)>(
3615 n_components,
3616 input,
3617 output,
3618 flag,
3619 face_no,
3620 n_points_1d,
3621 shape_data,
3622 dofs_per_component_on_cell,
3623 dofs_per_component_on_face);
3624 }
3625 }
3626
3627 template <typename EvalType>
3628 static EvalType
3631 const unsigned int face_no)
3632 {
3633 return EvalType(data.shape_data_on_face[face_no % 2], {}, {});
3634 }
3635
3636 template <bool do_evaluate,
3637 bool add_into_output,
3638 int face_direction = 0,
3639 int max_derivative = 0>
3640 static void
3642 const unsigned int n_components,
3643 const Number * input,
3644 Number * output,
3646 const unsigned int face_no,
3648 {
3649 if (dim == 1)
3650 {
3651 // This should never happen since the FE_RaviartThomasNodal is not
3652 // defined for dim = 1. It prevents compiler warnings of infinite
3653 // recursion.
3654 Assert(false, ExcInternalError());
3655 return;
3656 }
3657
3658 bool increase_max_der = false;
3659 if ((flag & EvaluationFlags::hessians && max_derivative < 2) ||
3660 (flag & EvaluationFlags::gradients && max_derivative < 1))
3661 increase_max_der = true;
3662
3663 if (face_direction == face_no / 2 && !increase_max_der)
3664 {
3666 add_into_output,
3667 face_direction,
3668 max_derivative>(
3669 shape_info, face_no, input, output);
3670 }
3671 else if (face_direction == face_no / 2)
3672 {
3673 // Only increase max_derivative
3675 add_into_output,
3676 face_direction,
3677 std::min(max_derivative + 1, 2)>(
3678 n_components, input, output, flag, face_no, shape_info);
3679 }
3680 else if (face_direction < dim)
3681 {
3682 if (increase_max_der)
3683 {
3685 do_evaluate,
3686 add_into_output,
3687 std::min(face_direction + 1, dim - 1),
3688 std::min(max_derivative + 1, 2)>(
3689 n_components, input, output, flag, face_no, shape_info);
3690 }
3691 else
3692 {
3694 add_into_output,
3695 std::min(face_direction + 1,
3696 dim - 1),
3697 max_derivative>(
3698 n_components, input, output, flag, face_no, shape_info);
3699 }
3700 }
3701 }
3702
3703 /* Help function for interpolate_generic_raviart_thomas */
3704 template <bool do_evaluate,
3705 bool add_into_output,
3706 int face_direction,
3707 int max_derivative>
3708 static inline void
3711 const unsigned int face_no,
3712 const Number * input,
3713 Number * output)
3714 {
3715 // These types are evaluators in either normal or tangential direction
3716 // depending on the face direction, with different normal directions for
3717 // the different components.
3718 using Evalf0 = typename std::conditional<
3719 face_direction == 0,
3721 dim,
3722 (fe_degree == -1) ? 1 : fe_degree + 1,
3723 0,
3724 0,
3725 Number,
3726 Number2>,
3728 dim,
3729 (fe_degree == -1) ? 1 : fe_degree,
3730 0,
3731 0,
3732 Number,
3733 Number2>>::type;
3734 using Evalf1 = typename std::conditional<
3735 face_direction == 1,
3737 dim,
3738 (fe_degree == -1) ? 1 : fe_degree + 1,
3739 0,
3740 1,
3741 Number,
3742 Number2>,
3744 dim,
3745 (fe_degree == -1) ? 1 : fe_degree,
3746 0,
3747 1,
3748 Number,
3749 Number2>>::type;
3750 using Evalf2 = typename std::conditional<
3751 face_direction == 2,
3753 dim,
3754 (fe_degree == -1) ? 1 : fe_degree + 1,
3755 0,
3756 2,
3757 Number,
3758 Number2>,
3760 dim,
3761 (fe_degree == -1) ? 1 : fe_degree,
3762 0,
3763 2,
3764 Number,
3765 Number2>>::type;
3766
3767 Evalf0 evalf0 =
3768 create_evaluator_tensor_product<Evalf0>((face_direction == 0) ?
3769 shape_info.data[0] :
3770 shape_info.data[1],
3771 face_no);
3772 Evalf1 evalf1 =
3773 create_evaluator_tensor_product<Evalf1>((face_direction == 1) ?
3774 shape_info.data[0] :
3775 shape_info.data[1],
3776 face_no);
3777 Evalf2 evalf2 =
3778 create_evaluator_tensor_product<Evalf2>((face_direction == 2) ?
3779 shape_info.data[0] :
3780 shape_info.data[1],
3781 face_no);
3782
3783 const unsigned int dofs_per_component_on_cell =
3784 shape_info.dofs_per_component_on_cell;
3785 const unsigned int dofs_per_component_on_face =
3786 3 * shape_info.dofs_per_component_on_face;
3787
3788 // NOTE! dofs_per_component_on_face is in the tangent direction,
3789 // i.e (fe.degree+1)*fe.degree. Normal faces are only
3790 // fe.degree*fe.degree
3791 const unsigned int in_stride =
3792 do_evaluate ? dofs_per_component_on_cell : dofs_per_component_on_face;
3793 const unsigned int out_stride =
3794 do_evaluate ? dofs_per_component_on_face : dofs_per_component_on_cell;
3795
3796 const unsigned int in_stride_after_normal =
3797 do_evaluate ?
3798 dofs_per_component_on_cell :
3799 dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2);
3800 const unsigned int out_stride_after_normal =
3801 do_evaluate ?
3802 dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2) :
3803 dofs_per_component_on_cell;
3804
3805 evalf0.template apply_face<face_direction,
3806 do_evaluate,
3807 add_into_output,
3808 max_derivative>(input, output);
3809 // stride to next component
3810 input += (face_direction == 0) ? in_stride_after_normal : in_stride;
3811 output += (face_direction == 0) ? out_stride_after_normal : out_stride;
3812
3813 evalf1.template apply_face<face_direction,
3814 do_evaluate,
3815 add_into_output,
3816 max_derivative>(input, output);
3817
3818 if (dim == 3)
3819 {
3820 // stride to next component
3821 input += (face_direction == 1) ? in_stride_after_normal : in_stride;
3822 output +=
3823 (face_direction == 1) ? out_stride_after_normal : out_stride;
3824
3825 evalf2.template apply_face<face_direction,
3826 do_evaluate,
3827 add_into_output,
3828 max_derivative>(input, output);
3829 }
3830 }
3831 };
3832
3833
3834
3835 // internal helper function for reading data; base version of different types
3836 template <typename VectorizedArrayType, typename Number2>
3837 void
3838 do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
3839 {
3840 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3841 dst[v] = src_ptr[v];
3842 }
3843
3844
3845
3846 // internal helper function for reading data; specialized version where we
3847 // can use a dedicated load function
3848 template <typename Number, std::size_t width>
3849 void
3851 {
3852 dst.load(src_ptr);
3853 }
3854
3855
3856
3857 // internal helper function for reading data; base version of different types
3858 template <typename VectorizedArrayType, typename Number2>
3859 void
3860 do_vectorized_gather(const Number2 * src_ptr,
3861 const unsigned int * indices,
3862 VectorizedArrayType &dst)
3863 {
3864 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3865 dst[v] = src_ptr[indices[v]];
3866 }
3867
3868
3869
3870 // internal helper function for reading data; specialized version where we
3871 // can use a dedicated gather function
3872 template <typename Number, std::size_t width>
3873 void
3874 do_vectorized_gather(const Number * src_ptr,
3875 const unsigned int * indices,
3877 {
3878 dst.gather(src_ptr, indices);
3879 }
3880
3881
3882
3883 // internal helper function for reading data; base version of different types
3884 template <typename VectorizedArrayType, typename Number2>
3885 void
3886 do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
3887 {
3888 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3889 dst_ptr[v] += src[v];
3890 }
3891
3892
3893
3894 // internal helper function for reading data; specialized version where we
3895 // can use a dedicated load function
3896 template <typename Number, std::size_t width>
3897 void
3899 {
3901 tmp.load(dst_ptr);
3902 (tmp + src).store(dst_ptr);
3903 }
3904
3905
3906
3907 // internal helper function for reading data; base version of different types
3908 template <typename VectorizedArrayType, typename Number2>
3909 void
3910 do_vectorized_scatter_add(const VectorizedArrayType src,
3911 const unsigned int * indices,
3912 Number2 * dst_ptr)
3913 {
3914 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3915 dst_ptr[indices[v]] += src[v];
3916 }
3917
3918
3919
3920 // internal helper function for reading data; specialized version where we
3921 // can use a dedicated gather function
3922 template <typename Number, std::size_t width>
3923 void
3925 const unsigned int * indices,
3926 Number * dst_ptr)
3927 {
3928#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
3929 for (unsigned int v = 0; v < width; ++v)
3930 dst_ptr[indices[v]] += src[v];
3931#else
3933 tmp.gather(dst_ptr, indices);
3934 (tmp + src).scatter(indices, dst_ptr);
3935#endif
3936 }
3937
3938
3939
3940 template <typename Number>
3941 void
3942 adjust_for_face_orientation(const unsigned int dim,
3943 const unsigned int n_components,
3945 const unsigned int *orientation,
3946 const bool integrate,
3947 const std::size_t n_q_points,
3948 Number * tmp_values,
3949 Number * values_quad,
3950 Number * gradients_quad,
3951 Number * hessians_quad)
3952 {
3953 for (unsigned int c = 0; c < n_components; ++c)
3954 {
3955 if (flag & EvaluationFlags::values)
3956 {
3957 if (integrate)
3958 for (unsigned int q = 0; q < n_q_points; ++q)
3959 tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
3960 else
3961 for (unsigned int q = 0; q < n_q_points; ++q)
3962 tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
3963 for (unsigned int q = 0; q < n_q_points; ++q)
3964 values_quad[c * n_q_points + q] = tmp_values[q];
3965 }
3966 if (flag & EvaluationFlags::gradients)
3967 for (unsigned int d = 0; d < dim; ++d)
3968 {
3969 if (integrate)
3970 for (unsigned int q = 0; q < n_q_points; ++q)
3971 tmp_values[q] =
3972 gradients_quad[(c * dim + d) * n_q_points + orientation[q]];
3973 else
3974 for (unsigned int q = 0; q < n_q_points; ++q)
3975 tmp_values[orientation[q]] =
3976 gradients_quad[(c * dim + d) * n_q_points + q];
3977 for (unsigned int q = 0; q < n_q_points; ++q)
3978 gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q];
3979 }
3980 if (flag & EvaluationFlags::hessians)
3981 {
3982 const unsigned int hdim = (dim * (dim + 1)) / 2;
3983 for (unsigned int d = 0; d < hdim; ++d)
3984 {
3985 if (integrate)
3986 for (unsigned int q = 0; q < n_q_points; ++q)
3987 tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
3988 orientation[q]];
3989 else
3990 for (unsigned int q = 0; q < n_q_points; ++q)
3991 tmp_values[orientation[q]] =
3992 hessians_quad[(c * hdim + d) * n_q_points + q];
3993 for (unsigned int q = 0; q < n_q_points; ++q)
3994 hessians_quad[(c * hdim + d) * n_q_points + q] =
3995 tmp_values[q];
3996 }
3997 }
3998 }
3999 }
4000
4001
4002
4003 template <typename Number, typename VectorizedArrayType>
4004 void
4006 const unsigned int dim,
4007 const unsigned int n_components,
4008 const unsigned int v,
4010 const unsigned int * orientation,
4011 const bool integrate,
4012 const std::size_t n_q_points,
4013 Number * tmp_values,
4014 VectorizedArrayType * values_quad,
4015 VectorizedArrayType * gradients_quad = nullptr,
4016 VectorizedArrayType * hessians_quad = nullptr)
4017 {
4018 for (unsigned int c = 0; c < n_components; ++c)
4019 {
4020 if (flag & EvaluationFlags::values)
4021 {
4022 if (integrate)
4023 for (unsigned int q = 0; q < n_q_points; ++q)
4024 tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v];
4025 else
4026 for (unsigned int q = 0; q < n_q_points; ++q)
4027 tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v];
4028 for (unsigned int q = 0; q < n_q_points; ++q)
4029 values_quad[c * n_q_points + q][v] = tmp_values[q];
4030 }
4031 if (flag & EvaluationFlags::gradients)
4032 for (unsigned int d = 0; d < dim; ++d)
4033 {
4034 Assert(gradients_quad != nullptr, ExcInternalError());
4035 if (integrate)
4036 for (unsigned int q = 0; q < n_q_points; ++q)
4037 tmp_values[q] = gradients_quad[(c * dim + d) * n_q_points +
4038 orientation[q]][v];
4039 else
4040 for (unsigned int q = 0; q < n_q_points; ++q)
4041 tmp_values[orientation[q]] =
4042 gradients_quad[(c * dim + d) * n_q_points + q][v];
4043 for (unsigned int q = 0; q < n_q_points; ++q)
4044 gradients_quad[(c * dim + d) * n_q_points + q][v] =
4045 tmp_values[q];
4046 }
4047 if (flag & EvaluationFlags::hessians)
4048 {
4049 Assert(hessians_quad != nullptr, ExcInternalError());
4050 const unsigned int hdim = (dim * (dim + 1)) / 2;
4051 for (unsigned int d = 0; d < hdim; ++d)
4052 {
4053 if (integrate)
4054 for (unsigned int q = 0; q < n_q_points; ++q)
4055 tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
4056 orientation[q]][v];
4057 else
4058 for (unsigned int q = 0; q < n_q_points; ++q)
4059 tmp_values[orientation[q]] =
4060 hessians_quad[(c * hdim + d) * n_q_points + q][v];
4061 for (unsigned int q = 0; q < n_q_points; ++q)
4062 hessians_quad[(c * hdim + d) * n_q_points + q][v] =
4063 tmp_values[q];
4064 }
4065 }
4066 }
4067 }
4068
4069
4070
4071 template <int dim, typename Number>
4073 {
4074 template <int fe_degree, int n_q_points_1d>
4075 static bool
4076 run(const unsigned int n_components,
4077 const EvaluationFlags::EvaluationFlags evaluation_flag,
4078 const Number * values_dofs,
4080 {
4081 const auto &shape_info = fe_eval.get_shape_info();
4082 const auto &shape_data = shape_info.data.front();
4083 using Number2 =
4085
4086 if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
4087 {
4088 Assert((fe_eval.get_dof_access_index() ==
4090 fe_eval.is_interior_face() == false) == false,
4092
4093 const unsigned int face_no = fe_eval.get_face_no();
4094 const unsigned int face_orientation = fe_eval.get_face_orientation();
4095 const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
4096 const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
4097
4098 using Eval =
4100
4101 if (evaluation_flag & EvaluationFlags::values)
4102 {
4103 const auto shape_values =
4104 &shape_data.shape_values_face(face_no, face_orientation, 0);
4105
4106 auto values_quad_ptr = fe_eval.begin_values();
4107 auto values_dofs_actual_ptr = values_dofs;
4108
4109 Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
4110 for (unsigned int c = 0; c < n_components; ++c)
4111 {
4112 eval.template values<0, true, false>(values_dofs_actual_ptr,
4113 values_quad_ptr);
4114
4115 values_quad_ptr += n_q_points;
4116 values_dofs_actual_ptr += n_dofs;
4117 }
4118 }
4119
4120 if (evaluation_flag & EvaluationFlags::gradients)
4121 {
4122 auto gradients_quad_ptr = fe_eval.begin_gradients();
4123 auto values_dofs_actual_ptr = values_dofs;
4124
4125 std::array<const Number2 *, dim> shape_gradients;
4126 for (unsigned int d = 0; d < dim; ++d)
4127 shape_gradients[d] = &shape_data.shape_gradients_face(
4128 face_no, face_orientation, d, 0);
4129
4130 for (unsigned int c = 0; c < n_components; ++c)
4131 {
4132 for (unsigned int d = 0; d < dim; ++d)
4133 {
4134 Eval eval(nullptr,
4135 shape_gradients[d],
4136 nullptr,
4137 n_dofs,
4138 n_q_points);
4139
4140 eval.template gradients<0, true, false>(
4141 values_dofs_actual_ptr, gradients_quad_ptr);
4142
4143 gradients_quad_ptr += n_q_points;
4144 }
4145 values_dofs_actual_ptr += n_dofs;
4146 }
4147 }
4148
4149 Assert(!(evaluation_flag & EvaluationFlags::hessians),
4151
4152 return true;
4153 }
4154
4155 const unsigned int dofs_per_face =
4156 fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
4157 Utilities::pow(shape_data.fe_degree + 1, dim - 1);
4158
4159 // Note: we always keep storage of values, 1st and 2nd derivatives in an
4160 // array, so reserve space for all three here
4161 Number *temp = fe_eval.get_scratch_data().begin();
4162 Number *scratch_data = temp + 3 * n_components * dofs_per_face;
4163
4164 bool use_vectorization = true;
4165
4166 if (fe_eval.get_dof_access_index() ==
4168 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
4169 for (unsigned int v = 0; v < Number::size(); ++v)
4170 if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
4171 fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
4172 use_vectorization = false;
4173
4174 if (use_vectorization == false)
4175 {
4176 for (unsigned int v = 0; v < Number::size(); ++v)
4177 {
4178 // the loop breaks once an invalid_unsigned_int is hit for
4179 // all cases except the exterior faces in the ECL loop (where
4180 // some faces might be at the boundaries but others not)
4181 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4182 {
4183 for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4184 ++i)
4185 temp[i][v] = 0;
4186 continue;
4187 }
4188
4190 template interpolate<true, false>(n_components,
4191 evaluation_flag,
4192 shape_info,
4193 values_dofs,
4194 scratch_data,
4195 fe_eval.get_face_no(v));
4196
4197 for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4198 ++i)
4199 temp[i][v] = scratch_data[i][v];
4200 }
4201 }
4202 else
4204 template interpolate<true, false>(n_components,
4205 evaluation_flag,
4206 shape_info,
4207 values_dofs,
4208 temp,
4209 fe_eval.get_face_no());
4210
4211 const unsigned int subface_index = fe_eval.get_subface_index();
4212 constexpr unsigned int n_q_points_1d_actual =
4213 fe_degree > -1 ? n_q_points_1d : 0;
4214
4215 if (fe_degree >= 1 &&
4216 shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
4217 {
4219 (fe_degree == -1) ? 1 : fe_degree,
4220 (n_q_points_1d < 1) ? 1 :
4221 n_q_points_1d,
4222 Number>::
4223 template evaluate_or_integrate_in_face<false>(
4224 evaluation_flag,
4225 temp,
4226 fe_eval,
4227 scratch_data,
4228 subface_index,
4229 fe_eval.get_face_no());
4230 }
4231 else if (fe_degree > -1 &&
4233 shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
4235 dim,
4236 fe_degree,
4237 n_q_points_1d_actual,
4238 Number>::evaluate_in_face(n_components,
4239 evaluation_flag,
4240 shape_data,
4241 temp,
4242 fe_eval.begin_values(),
4243 fe_eval
4244 .begin_gradients(),
4245 fe_eval.begin_hessians(),
4246 scratch_data,
4247 subface_index);
4248 else
4250 dim,
4251 fe_degree,
4252 n_q_points_1d_actual,
4253 Number>::evaluate_in_face(n_components,
4254 evaluation_flag,
4255 shape_data,
4256 temp,
4257 fe_eval.begin_values(),
4258 fe_eval
4259 .begin_gradients(),
4260 fe_eval.begin_hessians(),
4261 scratch_data,
4262 subface_index);
4263
4264 if (use_vectorization == false)
4265 {
4266 for (unsigned int v = 0; v < Number::size(); ++v)
4267 {
4268 // the loop breaks once an invalid_unsigned_int is hit for
4269 // all cases except the exterior faces in the ECL loop (where
4270 // some faces might be at the boundaries but others not)
4271 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4272 continue;
4273
4274 if (fe_eval.get_face_orientation(v) != 0)
4276 dim,
4277 n_components,
4278 v,
4279 evaluation_flag,
4281 fe_eval.get_face_orientation(v), 0),
4282 false,
4283 shape_info.n_q_points_face,
4284 &temp[0][0],
4285 fe_eval.begin_values(),
4286 fe_eval.begin_gradients(),
4287 fe_eval.begin_hessians());
4288 }
4289 }
4290 else if (fe_eval.get_face_orientation() != 0)
4292 dim,
4293 n_components,
4294 evaluation_flag,
4296 fe_eval.get_face_orientation(), 0),
4297 false,
4298 shape_info.n_q_points_face,
4299 temp,
4300 fe_eval.begin_values(),
4301 fe_eval.begin_gradients(),
4302 fe_eval.begin_hessians());
4303
4304 return false;
4305 }
4306 };
4307
4308
4309
4310 template <int dim, typename Number>
4312 {
4313 template <int fe_degree, int n_q_points_1d>
4314 static bool
4315 run(const unsigned int n_components,
4316 const EvaluationFlags::EvaluationFlags integration_flag,
4317 Number * values_dofs,
4319 {
4320 const auto &shape_info = fe_eval.get_shape_info();
4321 const auto &shape_data = shape_info.data.front();
4322 using Number2 =
4324
4325 if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
4326 {
4327 Assert((fe_eval.get_dof_access_index() ==
4329 fe_eval.is_interior_face() == false) == false,
4331
4332 const unsigned int face_no = fe_eval.get_face_no();
4333 const unsigned int face_orientation = fe_eval.get_face_orientation();
4334 const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
4335 const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
4336
4337 using Eval =
4339
4340 if (integration_flag & EvaluationFlags::values)
4341 {
4342 const auto shape_values =
4343 &shape_data.shape_values_face(face_no, face_orientation, 0);
4344
4345 auto values_quad_ptr = fe_eval.begin_values();
4346 auto values_dofs_actual_ptr = values_dofs;
4347
4348 Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
4349 for (unsigned int c = 0; c < n_components; ++c)
4350 {
4351 eval.template values<0, false, false>(values_quad_ptr,
4352 values_dofs_actual_ptr);
4353
4354 values_quad_ptr += n_q_points;
4355 values_dofs_actual_ptr += n_dofs;
4356 }
4357 }
4358
4359 if (integration_flag & EvaluationFlags::gradients)
4360 {
4361 auto gradients_quad_ptr = fe_eval.begin_gradients();
4362 auto values_dofs_actual_ptr = values_dofs;
4363
4364 std::array<const Number2 *, dim> shape_gradients;
4365 for (unsigned int d = 0; d < dim; ++d)
4366 shape_gradients[d] = &shape_data.shape_gradients_face(
4367 face_no, face_orientation, d, 0);
4368
4369 for (unsigned int c = 0; c < n_components; ++c)
4370 {
4371 for (unsigned int d = 0; d < dim; ++d)
4372 {
4373 Eval eval(nullptr,
4374 shape_gradients[d],
4375 nullptr,
4376 n_dofs,
4377 n_q_points);
4378
4379 if (!(integration_flag & EvaluationFlags::values) &&
4380 d == 0)
4381 eval.template gradients<0, false, false>(
4382 gradients_quad_ptr, values_dofs_actual_ptr);
4383 else
4384 eval.template gradients<0, false, true>(
4385 gradients_quad_ptr, values_dofs_actual_ptr);
4386
4387 gradients_quad_ptr += n_q_points;
4388 }
4389 values_dofs_actual_ptr += n_dofs;
4390 }
4391 }
4392
4393 Assert(!(integration_flag & EvaluationFlags::hessians),
4395
4396 return true;
4397 }
4398
4399 const unsigned int dofs_per_face =
4400 fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
4401 Utilities::pow(shape_data.fe_degree + 1, dim - 1);
4402
4403 Number *temp = fe_eval.get_scratch_data().begin();
4404 Number *scratch_data = temp + 3 * n_components * dofs_per_face;
4405
4406 bool use_vectorization = true;
4407
4408 if (fe_eval.get_dof_access_index() ==
4410 fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
4411 use_vectorization =
4413 std::all_of(fe_eval.get_cell_ids().begin() + 1,
4414 fe_eval.get_cell_ids().end(),
4415 [&](const auto &v) {
4416 return v == fe_eval.get_cell_ids()[0] ||
4417 v == numbers::invalid_unsigned_int;
4418 });
4419
4420 if (use_vectorization == false)
4421 {
4422 for (unsigned int v = 0; v < Number::size(); ++v)
4423 {
4424 // the loop breaks once an invalid_unsigned_int is hit for
4425 // all cases except the exterior faces in the ECL loop (where
4426 // some faces might be at the boundaries but others not)
4427 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4428 continue;
4429
4430 if (fe_eval.get_face_orientation(v) != 0)
4432 dim,
4433 n_components,
4434 v,
4435 integration_flag,
4437 fe_eval.get_face_orientation(v), 0),
4438 true,
4439 shape_info.n_q_points_face,
4440 &temp[0][0],
4441 fe_eval.begin_values(),
4442 fe_eval.begin_gradients(),
4443 fe_eval.begin_hessians());
4444 }
4445 }
4446 else if (fe_eval.get_face_orientation() != 0)
4448 dim,
4449 n_components,
4450 integration_flag,
4452 fe_eval.get_face_orientation(), 0),
4453 true,
4454 shape_info.n_q_points_face,
4455 temp,
4456 fe_eval.begin_values(),
4457 fe_eval.begin_gradients(),
4458 fe_eval.begin_hessians());
4459
4460 const unsigned int n_q_points_1d_actual =
4461 fe_degree > -1 ? n_q_points_1d : 0;
4462 const unsigned int subface_index = fe_eval.get_subface_index();
4463
4464 if (fe_degree >= 1 &&
4465 shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
4466 {
4468 (fe_degree == -1) ? 1 : fe_degree,
4469 (n_q_points_1d < 1) ? 1 :
4470 n_q_points_1d,
4471 Number>::
4472 template evaluate_or_integrate_in_face<true>(integration_flag,
4473 temp,
4474 fe_eval,
4475 scratch_data,
4476 subface_index,
4477 fe_eval.get_face_no());
4478 }
4479 else if (fe_degree > -1 &&
4480 fe_eval.get_subface_index() >=
4482 shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
4484 true,
4485 dim,
4486 fe_degree,
4487 n_q_points_1d_actual,
4488 Number>::integrate_in_face(n_components,
4489 integration_flag,
4490 shape_data,
4491 temp,
4492 fe_eval.begin_values(),
4493 fe_eval.begin_gradients(),
4494 fe_eval.begin_hessians(),
4495 scratch_data,
4496 subface_index);
4497 else
4499 false,
4500 dim,
4501 fe_degree,
4502 n_q_points_1d_actual,
4503 Number>::integrate_in_face(n_components,
4504 integration_flag,
4505 shape_data,
4506 temp,
4507 fe_eval.begin_values(),
4508 fe_eval.begin_gradients(),
4509 fe_eval.begin_hessians(),
4510 scratch_data,
4511 subface_index);
4512
4513 if (use_vectorization == false)
4514 {
4515 for (unsigned int v = 0; v < Number::size(); ++v)
4516 {
4517 // the loop breaks once an invalid_unsigned_int is hit for
4518 // all cases except the exterior faces in the ECL loop (where
4519 // some faces might be at the boundaries but others not)
4520 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4521 continue;
4522
4524 template interpolate<false, false>(n_components,
4525 integration_flag,
4526 shape_info,
4527 values_dofs,
4528 scratch_data,
4529 fe_eval.get_face_no(v));
4530
4531 for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
4532 ++i)
4533 temp[i][v] = scratch_data[i][v];
4534 }
4535 }
4536 else
4538 template interpolate<false, false>(n_components,
4539 integration_flag,
4540 shape_info,
4541 temp,
4542 values_dofs,
4543 fe_eval.get_face_no());
4544 return false;
4545 }
4546 };
4547
4548
4549
4550 template <int n_face_orientations,
4551 typename Processor,
4552 typename EvaluationData,
4553 const bool check_face_orientations = false>
4554 void
4556 Processor & proc,
4557 const unsigned int n_components,
4558 const EvaluationFlags::EvaluationFlags evaluation_flag,
4559 typename Processor::Number2_ * global_vector_ptr,
4560 const std::vector<ArrayView<const typename Processor::Number2_>> *sm_ptr,
4561 const EvaluationData & fe_eval,
4562 typename Processor::VectorizedArrayType_ * temp1)
4563 {
4564 constexpr int dim = Processor::dim_;
4565 constexpr int fe_degree = Processor::fe_degree_;
4566 using VectorizedArrayType = typename Processor::VectorizedArrayType_;
4567 constexpr int n_lanes = VectorizedArrayType::size();
4568
4569 using Number = typename Processor::Number_;
4570 using Number2_ = typename Processor::Number2_;
4571
4572 const auto & shape_data = fe_eval.get_shape_info().data.front();
4573 constexpr bool integrate = Processor::do_integrate;
4574 const unsigned int face_no = fe_eval.get_face_no();
4575 const auto & dof_info = fe_eval.get_dof_info();
4576 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
4577 const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index =
4578 fe_eval.get_dof_access_index();
4579 AssertIndexRange(cell,
4580 dof_info.index_storage_variants[dof_access_index].size());
4581 constexpr unsigned int dofs_per_face =
4582 Utilities::pow(fe_degree + 1, dim - 1);
4583 const unsigned int subface_index = fe_eval.get_subface_index();
4584
4585 const unsigned int n_filled_lanes =
4586 dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
4587
4588 bool all_faces_are_same = n_filled_lanes == n_lanes;
4589 if (n_face_orientations == n_lanes)
4590 for (unsigned int v = 1; v < n_lanes; ++v)
4591 if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) ||
4592 fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0))
4593 {
4594 all_faces_are_same = false;
4595 break;
4596 }
4597
4598 // check for re-orientation ...
4599 std::array<const unsigned int *, n_face_orientations> orientation = {};
4600
4601 if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same &&
4602 fe_eval.is_interior_face() == 0)
4603 for (unsigned int v = 0; v < n_lanes; ++v)
4604 {
4605 // the loop breaks once an invalid_unsigned_int is hit for
4606 // all cases except the exterior faces in the ECL loop (where
4607 // some faces might be at the boundaries but others not)
4608 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4609 continue;
4610
4611 if (shape_data.nodal_at_cell_boundaries &&
4612 fe_eval.get_face_orientation(v) != 0)
4613 {
4614 // ... and in case we detect a re-orientation, go to the other
4615 // version of this function that actually allows for this
4616 if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
4617 check_face_orientations == false)
4618 {
4619 fe_face_evaluation_process_and_io<n_face_orientations,
4620 Processor,
4621 EvaluationData,
4622 true>(proc,
4623 n_components,
4624 evaluation_flag,
4625 global_vector_ptr,
4626 sm_ptr,
4627 fe_eval,
4628 temp1);
4629 return;
4630 }
4631 orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
4632 fe_eval.get_face_orientation(v), 0);
4633 }
4634 }
4635 else if (dim == 3 && fe_eval.get_face_orientation() != 0)
4636 {
4637 // go to the other version of this function
4638 if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
4639 check_face_orientations == false)
4640 {
4641 fe_face_evaluation_process_and_io<n_face_orientations,
4642 Processor,
4643 EvaluationData,
4644 true>(proc,
4645 n_components,
4646 evaluation_flag,
4647 global_vector_ptr,
4648 sm_ptr,
4649 fe_eval,
4650 temp1);
4651 return;
4652 }
4653 for (unsigned int v = 0; v < n_face_orientations; ++v)
4654 orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
4655 fe_eval.get_face_orientation(), 0);
4656 }
4657
4658 // we know that the gradient weights for the Hermite case on the
4659 // right (side==1) are the negative from the value at the left
4660 // (side==0), so we only read out one of them.
4661 VectorizedArrayType grad_weight =
4662 shape_data
4663 .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) :
4664 (1 + face_no % 2))];
4665
4666 // face_to_cell_index_hermite
4667 std::array<const unsigned int *, n_face_orientations> index_array_hermite =
4668 {};
4669 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4670 {
4671 if (n_face_orientations == 1)
4672 index_array_hermite[0] =
4673 &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0);
4674 else
4675 {
4676 for (unsigned int v = 0; v < n_lanes; ++v)
4677 {
4678 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4679 continue;
4680
4681 const auto face_no = fe_eval.get_face_no(v);
4682
4683 grad_weight[v] =
4684 shape_data.shape_data_on_face[0][fe_degree +
4685 (integrate ?
4686 (2 - (face_no % 2)) :
4687 (1 + (face_no % 2)))][0];
4688
4689 index_array_hermite[v] =
4690 &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no,
4691 0);
4692 }
4693 }
4694 }
4695
4696 // face_to_cell_index_nodal
4697 std::array<const unsigned int *, n_face_orientations> index_array_nodal =
4698 {};
4699 if (shape_data.nodal_at_cell_boundaries == true)
4700 {
4701 if (n_face_orientations == 1)
4702 index_array_nodal[0] =
4703 &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0);
4704 else
4705 {
4706 for (unsigned int v = 0; v < n_lanes; ++v)
4707 {
4708 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
4709 continue;
4710
4711 const auto face_no = fe_eval.get_face_no(v);
4712
4713 index_array_nodal[v] =
4714 &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no,
4715 0);
4716 }
4717 }
4718 }
4719
4720
4721 const auto reorientate = [&](const unsigned int v, const unsigned int i) {
4722 return (!check_face_orientations || orientation[v] == nullptr) ?
4723 i :
4724 orientation[v][i];
4725 };
4726
4727 const unsigned int cell_index =
4729 fe_eval.get_cell_ids()[0] :
4730 cell * n_lanes;
4731 const unsigned int *dof_indices =
4732 &dof_info.dof_indices_contiguous[dof_access_index][cell_index];
4733
4734 for (unsigned int comp = 0; comp < n_components; ++comp)
4735 {
4736 const std::size_t index_offset =
4737 dof_info.component_dof_indices_offset
4738 [fe_eval.get_active_fe_index()]
4739 [fe_eval.get_first_selected_component()] +
4740 comp * Utilities::pow(fe_degree + 1, dim);
4741
4742 // case 1: contiguous and interleaved indices
4743 if (n_face_orientations == 1 &&
4744 dof_info.index_storage_variants[dof_access_index][cell] ==
4746 interleaved_contiguous)
4747 {
4749 dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
4750 n_lanes);
4751 Number2_ *vector_ptr =
4752 global_vector_ptr + dof_indices[0] + index_offset * n_lanes;
4753
4754 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4755 {
4756 for (unsigned int i = 0; i < dofs_per_face; ++i)
4757 {
4758 Assert(n_face_orientations == 1, ExcNotImplemented());
4759
4760 const unsigned int ind1 = index_array_hermite[0][2 * i];
4761 const unsigned int ind2 = index_array_hermite[0][2 * i + 1];
4762 const unsigned int i_ = reorientate(0, i);
4763 proc.hermite_grad_vectorized(temp1[i_],
4764 temp1[i_ + dofs_per_face],
4765 vector_ptr + ind1 * n_lanes,
4766 vector_ptr + ind2 * n_lanes,
4767 grad_weight);
4768 }
4769 }
4770 else
4771 {
4772 for (unsigned int i = 0; i < dofs_per_face; ++i)
4773 {
4774 Assert(n_face_orientations == 1, ExcNotImplemented());
4775
4776 const unsigned int i_ = reorientate(0, i);
4777 const unsigned int ind = index_array_nodal[0][i];
4778 proc.value_vectorized(temp1[i_],
4779 vector_ptr + ind * n_lanes);
4780 }
4781 }
4782 }
4783
4784 // case 2: contiguous and interleaved indices with fixed stride
4785 else if (n_face_orientations == 1 &&
4786 dof_info.index_storage_variants[dof_access_index][cell] ==
4788 interleaved_contiguous_strided)
4789 {
4791 dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
4792 n_lanes);
4793 Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes;
4794 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4795 {
4796 for (unsigned int i = 0; i < dofs_per_face; ++i)
4797 {
4798 Assert(n_face_orientations == 1, ExcNotImplemented());
4799
4800 const unsigned int i_ = reorientate(0, i);
4801 const unsigned int ind1 =
4802 index_array_hermite[0][2 * i] * n_lanes;
4803 const unsigned int ind2 =
4804 index_array_hermite[0][2 * i + 1] * n_lanes;
4805 proc.hermite_grad_vectorized_indexed(
4806 temp1[i_],
4807 temp1[i_ + dofs_per_face],
4808 vector_ptr + ind1,
4809 vector_ptr + ind2,
4810 grad_weight,
4811 dof_indices,
4812 dof_indices);
4813 }
4814 }
4815 else
4816 {
4817 for (unsigned int i = 0; i < dofs_per_face; ++i)
4818 {
4819 Assert(n_face_orientations == 1, ExcNotImplemented());
4820
4821 const unsigned int i_ = reorientate(0, i);
4822 const unsigned int ind = index_array_nodal[0][i] * n_lanes;
4823 proc.value_vectorized_indexed(temp1[i_],
4824 vector_ptr + ind,
4825 dof_indices);
4826 }
4827 }
4828 }
4829
4830 // case 3: contiguous and interleaved indices with mixed stride
4831 else if (n_face_orientations == 1 &&
4832 dof_info.index_storage_variants[dof_access_index][cell] ==
4834 interleaved_contiguous_mixed_strides)
4835 {
4836 const unsigned int *strides =
4837 &dof_info.dof_indices_interleave_strides[dof_access_index]
4838 [cell * n_lanes];
4839 unsigned int indices[n_lanes];
4840 for (unsigned int v = 0; v < n_lanes; ++v)
4841 indices[v] = dof_indices[v] + index_offset * strides[v];
4842 const unsigned int n_filled_lanes =
4843 dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
4844
4845 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
4846 {
4847 if (n_filled_lanes == n_lanes)
4848 for (unsigned int i = 0; i < dofs_per_face; ++i)
4849 {
4850 Assert(n_face_orientations == 1, ExcNotImplemented());
4851
4852 const unsigned int i_ = reorientate(0, i);
4853 unsigned int ind1[n_lanes];
4855 for (unsigned int v = 0; v < n_lanes; ++v)
4856 ind1[v] = indices[v] +
4857 index_array_hermite[0][2 * i] * strides[v];
4858 unsigned int ind2[n_lanes];
4860 for (unsigned int v = 0; v < n_lanes; ++v)
4861 ind2[v] =
4862 indices[v] +
4863 // TODO
4864 index_array_hermite[0][2 * i + 1] * strides[v];
4865 proc.hermite_grad_vectorized_indexed(
4866 temp1[i_],
4867 temp1[i_ + dofs_per_face],
4868 global_vector_ptr,
4869 global_vector_ptr,
4870 grad_weight,
4871 ind1,
4872 ind2);
4873 }
4874 else
4875 {
4876 if (integrate == false)
4877 for (unsigned int i = 0; i < 2 * dofs_per_face; ++i)
4878 temp1[i] = VectorizedArrayType();
4879
4880 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4881 for (unsigned int i = 0; i < dofs_per_face; ++i)
4882 {
4883 const unsigned int i_ =
4884 reorientate(n_face_orientations == 1 ? 0 : v, i);
4885 proc.hermite_grad(
4886 temp1[i_][v],
4887 temp1[i_ + dofs_per_face][v],
4888 global_vector_ptr
4889 [indices[v] +
4890 index_array_hermite
4891 [n_face_orientations == 1 ? 0 : v][2 * i] *
4892 strides[v]],
4893 global_vector_ptr
4894 [indices[v] +
4895 index_array_hermite[n_face_orientations == 1 ?
4896 0 :
4897 v][2 * i + 1] *
4898 strides[v]],
4899 grad_weight[n_face_orientations == 1 ? 0 : v]);
4900 }
4901 }
4902 }
4903 else
4904 {
4905 if (n_filled_lanes == n_lanes)
4906 for (unsigned int i = 0; i < dofs_per_face; ++i)
4907 {
4908 Assert(n_face_orientations == 1, ExcInternalError());
4909 unsigned int ind[n_lanes];
4911 for (unsigned int v = 0; v < n_lanes; ++v)
4912 ind[v] =
4913 indices[v] + index_array_nodal[0][i] * strides[v];
4914 const unsigned int i_ = reorientate(0, i);
4915 proc.value_vectorized_indexed(temp1[i_],
4916 global_vector_ptr,
4917 ind);
4918 }
4919 else
4920 {
4921 if (integrate == false)
4922 for (unsigned int i = 0; i < dofs_per_face; ++i)
4923 temp1[i] = VectorizedArrayType();
4924
4925 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4926 for (unsigned int i = 0; i < dofs_per_face; ++i)
4927 proc.value(
4928 temp1[reorientate(n_face_orientations == 1 ? 0 : v,
4929 i)][v],
4930 global_vector_ptr
4931 [indices[v] +
4932 index_array_nodal[n_face_orientations == 1 ? 0 : v]
4933 [i] *
4934 strides[v]]);
4935 }
4936 }
4937 }
4938
4939 // case 4: contiguous indices without interleaving
4940 else if (n_face_orientations > 1 ||
4941 dof_info.index_storage_variants[dof_access_index][cell] ==
4943 contiguous)
4944 {
4945 Number2_ *vector_ptr = global_vector_ptr + index_offset;
4946
4947 const bool vectorization_possible =
4948 all_faces_are_same && (sm_ptr == nullptr);
4949
4950 std::array<Number2_ *, n_lanes> vector_ptrs;
4951 std::array<unsigned int, n_lanes> reordered_indices;
4952
4953 if (vectorization_possible == false)
4954 {
4955 vector_ptrs = {};
4956 if (n_face_orientations == 1)
4957 {
4958 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4959 if (sm_ptr == nullptr)
4960 {
4961 vector_ptrs[v] = vector_ptr + dof_indices[v];
4962 }
4963 else
4964 {
4965 const auto &temp =
4966 dof_info
4967 .dof_indices_contiguous_sm[dof_access_index]
4968 [cell * n_lanes + v];
4969 vector_ptrs[v] = const_cast<Number2_ *>(
4970 sm_ptr->operator[](temp.first).data() +
4971 temp.second + index_offset);
4972 }
4973 }
4974 else if (n_face_orientations == n_lanes)
4975 {
4976 const auto &cells = fe_eval.get_cell_ids();
4977 for (unsigned int v = 0; v < n_lanes; ++v)
4978 if (cells[v] != numbers::invalid_unsigned_int)
4979 {
4980 if (sm_ptr == nullptr)
4981 {
4982 vector_ptrs[v] =
4983 vector_ptr +
4984 dof_info
4985 .dof_indices_contiguous[dof_access_index]
4986 [cells[v]];
4987 }
4988 else
4989 {
4990 const auto &temp =
4991 dof_info
4992 .dof_indices_contiguous_sm[dof_access_index]
4993 [cells[v]];
4994 vector_ptrs[v] = const_cast<Number2_ *>(
4995 sm_ptr->operator[](temp.first).data() +
4996 temp.second + index_offset);
4997 }
4998 }
4999 }
5000 else
5001 {
5002 Assert(false, ExcNotImplemented());
5003 }
5004 }
5005 else if (n_face_orientations == n_lanes)
5006 {
5007 for (unsigned int v = 0; v < n_lanes; ++v)
5008 reordered_indices[v] =
5009 dof_info.dof_indices_contiguous[dof_access_index]
5010 [fe_eval.get_cell_ids()[v]];
5011 dof_indices = reordered_indices.data();
5012 }
5013
5014 if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
5015 {
5016 if (vectorization_possible)
5017 for (unsigned int i = 0; i < dofs_per_face; ++i)
5018 {
5019 const unsigned int ind1 = index_array_hermite[0][2 * i];
5020 const unsigned int ind2 =
5021 index_array_hermite[0][2 * i + 1];
5022 const unsigned int i_ = reorientate(0, i);
5023
5024 proc.hermite_grad_vectorized_indexed(
5025 temp1[i_],
5026 temp1[i_ + dofs_per_face],
5027 vector_ptr + ind1,
5028 vector_ptr + ind2,
5029 grad_weight,
5030 dof_indices,
5031 dof_indices);
5032 }
5033 else if (n_face_orientations == 1)
5034 for (unsigned int i = 0; i < dofs_per_face; ++i)
5035 {
5036 const unsigned int ind1 = index_array_hermite[0][2 * i];
5037 const unsigned int ind2 =
5038 index_array_hermite[0][2 * i + 1];
5039 const unsigned int i_ = reorientate(0, i);
5040
5041 for (unsigned int v = 0; v < n_filled_lanes; ++v)
5042 proc.hermite_grad(temp1[i_][v],
5043 temp1[i_ + dofs_per_face][v],
5044 vector_ptrs[v][ind1],
5045 vector_ptrs[v][ind2],
5046 grad_weight[v]);
5047
5048 if (integrate == false)
5049 for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
5050 {
5051 temp1[i][v] = 0.0;
5052 temp1[i + dofs_per_face][v] = 0.0;
5053 }
5054 }
5055 else
5056 {
5057 if (integrate == false && n_filled_lanes < n_lanes)
5058 for (unsigned int i = 0; i < dofs_per_face; ++i)
5059 temp1[i] = temp1[i + dofs_per_face] = Number();
5060
5061 for (unsigned int v = 0; v < n_filled_lanes; ++v)
5062 for (unsigned int i = 0; i < dofs_per_face; ++i)
5063 proc.hermite_grad(
5064 temp1[reorientate(v, i)][v],
5065 temp1[reorientate(v, i) + dofs_per_face][v],
5066 vector_ptrs[v][index_array_hermite[v][2 * i]],
5067 vector_ptrs[v][index_array_hermite[v][2 * i + 1]],
5068 grad_weight[v]);
5069 }
5070 }
5071 else
5072 {
5073 if (vectorization_possible)
5074 for (unsigned int i = 0; i < dofs_per_face; ++i)
5075 {
5076 const unsigned int ind = index_array_nodal[0][i];
5077 const unsigned int i_ = reorientate(0, i);
5078
5079 proc.value_vectorized_indexed(temp1[i_],
5080 vector_ptr + ind,
5081 dof_indices);
5082 }
5083 else if (n_face_orientations == 1)
5084 for (unsigned int i = 0; i < dofs_per_face; ++i)
5085 {
5086 const unsigned int ind = index_array_nodal[0][i];
5087 const unsigned int i_ = reorientate(0, i);
5088
5089 for (unsigned int v = 0; v < n_filled_lanes; ++v)
5090 proc.value(temp1[i_][v], vector_ptrs[v][ind]);
5091
5092 if (integrate == false)
5093 for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
5094 temp1[i_][v] = 0.0;
5095 }
5096 else
5097 {
5098 if (integrate == false && n_filled_lanes < n_lanes)
5099 for (unsigned int i = 0; i < dofs_per_face; ++i)
5100 temp1[i] = Number();
5101
5102 for (unsigned int v = 0; v < n_filled_lanes; ++v)
5103 for (unsigned int i = 0; i < dofs_per_face; ++i)
5104 proc.value(temp1[reorientate(v, i)][v],
5105 vector_ptrs[v][index_array_nodal[v][i]]);
5106 }
5107 }
5108 }
5109 else
5110 {
5111 // We should not end up here, this should be caught by
5112 // FEFaceEvaluationImplGatherEvaluateSelector::supports()
5113 Assert(false, ExcInternalError());
5114 }
5115 temp1 += 3 * dofs_per_face;
5116 }
5117 }
5118
5119
5120
5121 template <int dim, typename Number2, typename VectorizedArrayType>
5123 {
5124 using Number = typename VectorizedArrayType::value_type;
5125
5126 template <int fe_degree, int n_q_points_1d>
5127 static bool
5128 run(const unsigned int n_components,
5129 const EvaluationFlags::EvaluationFlags evaluation_flag,
5130 const Number2 * src_ptr,
5131 const std::vector<ArrayView<const Number2>> * sm_ptr,
5133 {
5134 Assert(fe_degree > -1, ExcInternalError());
5138
5139 const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
5140
5141 VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
5142 VectorizedArrayType *scratch_data =
5143 temp + 3 * n_components * dofs_per_face;
5144
5146
5147 if (fe_eval.get_dof_access_index() ==
5149 fe_eval.is_interior_face() == false)
5150 fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
5151 p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
5152 else
5153 fe_face_evaluation_process_and_io<1>(
5154 p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
5155
5156 const unsigned int subface_index = fe_eval.get_subface_index();
5157
5158 if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
5160 dim,
5161 fe_degree,
5162 n_q_points_1d,
5163 VectorizedArrayType>::
5164 evaluate_in_face(n_components,
5165 evaluation_flag,
5166 fe_eval.get_shape_info().data.front(),
5167 temp,
5168 fe_eval.begin_values(),
5169 fe_eval.begin_gradients(),
5170 fe_eval.begin_hessians(),
5171 scratch_data,
5172 subface_index);
5173 else
5175 dim,
5176 fe_degree,
5177 n_q_points_1d,
5178 VectorizedArrayType>::
5179 evaluate_in_face(n_components,
5180 evaluation_flag,
5181 fe_eval.get_shape_info().data.front(),
5182 temp,
5183 fe_eval.begin_values(),
5184 fe_eval.begin_gradients(),
5185 fe_eval.begin_hessians(),
5186 scratch_data,
5187 subface_index);
5188
5189 // re-orientation for cases not possible with above algorithm
5190 if (subface_index < GeometryInfo<dim>::max_children_per_cell)
5191 {
5192 if (fe_eval.get_dof_access_index() ==
5194 fe_eval.is_interior_face() == false)
5195 {
5196 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5197 {
5198 // the loop breaks once an invalid_unsigned_int is hit for
5199 // all cases except the exterior faces in the ECL loop (where
5200 // some faces might be at the boundaries but others not)
5201 if (fe_eval.get_cell_ids()[v] ==
5203 continue;
5204
5205 if (fe_eval.get_face_orientation(v) != 0)
5207 dim,
5208 n_components,
5209 v,
5210 evaluation_flag,
5212 fe_eval.get_face_orientation(v), 0),
5213 false,
5214 Utilities::pow(n_q_points_1d, dim - 1),
5215 &temp[0][0],
5216 fe_eval.begin_values(),
5217 fe_eval.begin_gradients(),
5218 fe_eval.begin_hessians());
5219 }
5220 }
5221 else if (fe_eval.get_face_orientation() != 0)
5223 dim,
5224 n_components,
5225 evaluation_flag,
5227 fe_eval.get_face_orientation(), 0),
5228 false,
5229 Utilities::pow(n_q_points_1d, dim - 1),
5230 temp,
5231 fe_eval.begin_values(),
5232 fe_eval.begin_gradients(),
5233 fe_eval.begin_hessians());
5234 }
5235
5236 return false;
5237 }
5238
5239 template <typename Number3>
5240 static bool
5242 const MatrixFreeFunctions::ShapeInfo<Number3> & shape_info,
5243 const Number2 * vector_ptr,
5245 {
5246 const unsigned int fe_degree = shape_info.data.front().fe_degree;
5247 if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries ||
5248 (evaluation_flag & EvaluationFlags::gradients &&
5249 (fe_degree < 2 ||
5250 shape_info.data.front().element_type !=
5252 (evaluation_flag & EvaluationFlags::hessians) ||
5253 vector_ptr == nullptr ||
5254 shape_info.data.front().element_type >
5256 storage <
5258 return false;
5259 else
5260 return true;
5261 }
5262
5263 private:
5264 template <int fe_degree>
5266 {
5267 static const bool do_integrate = false;
5268 static const int dim_ = dim;
5269 static const int fe_degree_ = fe_degree;
5270 using VectorizedArrayType_ = VectorizedArrayType;
5272 using Number2_ = const Number2;
5273
5274 template <typename T0, typename T1, typename T2>
5275 void
5277 T0 & temp_2,
5278 const T1 src_ptr_1,
5279 const T1 src_ptr_2,
5280 const T2 &grad_weight)
5281 {
5282 do_vectorized_read(src_ptr_1, temp_1);
5283 do_vectorized_read(src_ptr_2, temp_2);
5284 temp_2 = grad_weight * (temp_1 - temp_2);
5285 }
5286
5287 template <typename T1, typename T2>
5288 void
5289 value_vectorized(T1 &temp, const T2 src_ptr)
5290 {
5291 do_vectorized_read(src_ptr, temp);
5292 }
5293
5294 template <typename T0, typename T1, typename T2, typename T3>
5295 void
5297 T0 & temp_2,
5298 const T1 src_ptr_1,
5299 const T1 src_ptr_2,
5300 const T2 &grad_weight,
5301 const T3 &indices_1,
5302 const T3 &indices_2)
5303 {
5304 do_vectorized_gather(src_ptr_1, indices_1, temp_1);
5305 do_vectorized_gather(src_ptr_2, indices_2, temp_2);
5306 temp_2 = grad_weight * (temp_1 - temp_2);
5307 }
5308
5309 template <typename T0, typename T1, typename T2>
5310 void
5311 value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
5312 {
5313 do_vectorized_gather(src_ptr, indices, temp);
5314 }
5315
5316 template <typename T0, typename T1, typename T2>
5317 void
5318 hermite_grad(T0 & temp_1,
5319 T0 & temp_2,
5320 const T1 &src_ptr_1,
5321 const T1 &src_ptr_2,
5322 const T2 &grad_weight)
5323 {
5324 // case 3a)
5325 temp_1 = src_ptr_1;
5326 temp_2 = grad_weight * (temp_1 - src_ptr_2);
5327 }
5328
5329 template <typename T1, typename T2>
5330 void
5331 value(T1 &temp, const T2 &src_ptr)
5332 {
5333 // case 3b)
5334 temp = src_ptr;
5335 }
5336 };
5337 };
5338
5339
5340
5341 template <int dim, typename Number2, typename VectorizedArrayType>
5343 {
5344 using Number = typename VectorizedArrayType::value_type;
5345
5346 template <int fe_degree, int n_q_points_1d>
5347 static bool
5348 run(const unsigned int n_components,
5349 const EvaluationFlags::EvaluationFlags integration_flag,
5350 Number2 * dst_ptr,
5351 const std::vector<ArrayView<const Number2>> * sm_ptr,
5353 {
5354 Assert(fe_degree > -1, ExcInternalError());
5358
5359 const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
5360
5361 VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
5362 VectorizedArrayType *scratch_data =
5363 temp + 3 * n_components * dofs_per_face;
5364
5365 const unsigned int subface_index = fe_eval.get_subface_index();
5366
5367 // re-orientation for cases not possible with the io function below
5368 if (subface_index < GeometryInfo<dim>::max_children_per_cell)
5369 {
5370 if (fe_eval.get_dof_access_index() ==
5372 fe_eval.is_interior_face() == false)
5373 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5374 {
5375 // the loop breaks once an invalid_unsigned_int is hit for
5376 // all cases except the exterior faces in the ECL loop (where
5377 // some faces might be at the boundaries but others not)
5378 if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
5379 continue;
5380
5381 if (fe_eval.get_face_orientation(v) != 0)
5383 dim,
5384 n_components,
5385 v,
5386 integration_flag,
5388 fe_eval.get_face_orientation(v), 0),
5389 true,
5390 Utilities::pow(n_q_points_1d, dim - 1),
5391 &temp[0][0],
5392 fe_eval.begin_values(),
5393 fe_eval.begin_gradients(),
5394 fe_eval.begin_hessians());
5395 }
5396 else if (fe_eval.get_face_orientation() != 0)
5398 dim,
5399 n_components,
5400 integration_flag,
5402 fe_eval.get_face_orientation(), 0),
5403 true,
5404 Utilities::pow(n_q_points_1d, dim - 1),
5405 temp,
5406 fe_eval.begin_values(),
5407 fe_eval.begin_gradients(),
5408 fe_eval.begin_hessians());
5409 }
5410
5411 if (fe_degree > -1 && fe_eval.get_subface_index() >=
5412 GeometryInfo<dim - 1>::max_children_per_cell)
5414 dim,
5415 fe_degree,
5416 n_q_points_1d,
5417 VectorizedArrayType>::
5418 integrate_in_face(n_components,
5419 integration_flag,
5420 fe_eval.get_shape_info().data.front(),
5421 temp,
5422 fe_eval.begin_values(),
5423 fe_eval.begin_gradients(),
5424 fe_eval.begin_hessians(),
5425 scratch_data,
5426 subface_index);
5427 else
5429 dim,
5430 fe_degree,
5431 n_q_points_1d,
5432 VectorizedArrayType>::
5433 integrate_in_face(n_components,
5434 integration_flag,
5435 fe_eval.get_shape_info().data.front(),
5436 temp,
5437 fe_eval.begin_values(),
5438 fe_eval.begin_gradients(),
5439 fe_eval.begin_hessians(),
5440 scratch_data,
5441 subface_index);
5442
5444
5445 if (fe_eval.get_dof_access_index() ==
5447 fe_eval.is_interior_face() == false)
5448 fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
5449 p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
5450 else
5451 fe_face_evaluation_process_and_io<1>(
5452 p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
5453
5454 return false;
5455 }
5456
5457 private:
5458 template <int fe_degree>
5460 {
5461 static const bool do_integrate = true;
5462 static const int dim_ = dim;
5463 static const int fe_degree_ = fe_degree;
5464 using VectorizedArrayType_ = VectorizedArrayType;
5466 using Number2_ = Number2;
5467
5468 template <typename T0, typename T1, typename T2, typename T3, typename T4>
5469 void
5470 hermite_grad_vectorized(const T0 &temp_1,
5471 const T1 &temp_2,
5472 T2 dst_ptr_1,
5473 T3 dst_ptr_2,
5474 const T4 &grad_weight)
5475 {
5476 // case 1a)
5477 const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
5478 const VectorizedArrayType grad = grad_weight * temp_2;
5479 do_vectorized_add(val, dst_ptr_1);
5480 do_vectorized_add(grad, dst_ptr_2);
5481 }
5482
5483 template <typename T0, typename T1>
5484 void
5485 value_vectorized(const T0 &temp, T1 dst_ptr)
5486 {
5487 // case 1b)
5488 do_vectorized_add(temp, dst_ptr);
5489 }
5490
5491 template <typename T0, typename T1, typename T2, typename T3>
5492 void
5494 const T0 &temp_2,
5495 T1 dst_ptr_1,
5496 T1 dst_ptr_2,
5497 const T2 &grad_weight,
5498 const T3 &indices_1,
5499 const T3 &indices_2)
5500 {
5501 // case 2a)
5502 const VectorizedArrayType val = temp_1 - grad_weight * temp_2;
5503 const VectorizedArrayType grad = grad_weight * temp_2;
5504 do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
5505 do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
5506 }
5507
5508 template <typename T0, typename T1, typename T2>
5509 void
5510 value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
5511 {
5512 // case 2b)
5513 do_vectorized_scatter_add(temp, indices, dst_ptr);
5514 }
5515
5516 template <typename T0, typename T1, typename T2>
5517 void
5518 hermite_grad(const T0 &temp_1,
5519 const T0 &temp_2,
5520 T1 & dst_ptr_1,
5521 T1 & dst_ptr_2,
5522 const T2 &grad_weight)
5523 {
5524 // case 3a)
5525 const Number val = temp_1 - grad_weight * temp_2;
5526 const Number grad = grad_weight * temp_2;
5527 dst_ptr_1 += val;
5528 dst_ptr_2 += grad;
5529 }
5530
5531 template <typename T0, typename T1>
5532 void
5533 value(const T0 &temp, T1 &dst_ptr)
5534 {
5535 // case 3b)
5536 dst_ptr += temp;
5537 }
5538 };
5539 };
5540
5541
5542
5547 template <int dim, typename Number>
5549 {
5550 using Number2 =
5552
5553 template <int fe_degree, int = 0>
5554 static bool
5555 run(const unsigned int n_components,
5557 const Number * in_array,
5558 Number * out_array)
5559 {
5560 const unsigned int given_degree =
5561 (fe_degree > -1) ? fe_degree :
5562 fe_eval.get_shape_info().data.front().fe_degree;
5563
5564 const unsigned int dofs_per_component =
5565 Utilities::pow(given_degree + 1, dim);
5566
5567 Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
5571
5573 dim,
5574 fe_degree + 1,
5575 fe_degree + 1,
5576 Number,
5577 Number2>
5578 evaluator({},
5579 {},
5580 fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
5581 given_degree + 1,
5582 given_degree + 1);
5583
5584 for (unsigned int d = 0; d < n_components; ++d)
5585 {
5586 const Number *in = in_array + d * dofs_per_component;
5587 Number * out = out_array + d * dofs_per_component;
5588 // Need to select 'apply' method with hessian slot because values
5589 // assume symmetries that do not exist in the inverse shapes
5590 evaluator.template hessians<0, true, false>(in, out);
5591 if (dim > 1)
5592 evaluator.template hessians<1, true, false>(out, out);
5593 if (dim > 2)
5594 evaluator.template hessians<2, true, false>(out, out);
5595 }
5596 for (unsigned int q = 0; q < dofs_per_component; ++q)
5597 {
5598 const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
5599 for (unsigned int d = 0; d < n_components; ++d)
5600 out_array[q + d * dofs_per_component] *= inverse_JxW_q;
5601 }
5602 for (unsigned int d = 0; d < n_components; ++d)
5603 {
5604 Number *out = out_array + d * dofs_per_component;
5605 if (dim > 2)
5606 evaluator.template hessians<2, false, false>(out, out);
5607 if (dim > 1)
5608 evaluator.template hessians<1, false, false>(out, out);
5609 evaluator.template hessians<0, false, false>(out, out);
5610 }
5611 return false;
5612 }
5613 };
5614
5615
5616
5623 template <int dim, typename Number>
5625 {
5626 using Number2 =
5628
5629 template <int fe_degree, int = 0>
5630 static bool
5631 run(const unsigned int n_desired_components,
5633 const ArrayView<const Number> & inverse_coefficients,
5634 const bool dyadic_coefficients,
5635 const Number * in_array,
5636 Number * out_array)
5637 {
5638 const unsigned int given_degree =
5639 (fe_degree > -1) ? fe_degree :
5640 fe_eval.get_shape_info().data.front().fe_degree;
5641
5642 const unsigned int dofs_per_component =
5643 Utilities::pow(given_degree + 1, dim);
5644
5645 Assert(inverse_coefficients.size() > 0 &&
5646 inverse_coefficients.size() % dofs_per_component == 0,
5647 ExcMessage(
5648 "Expected diagonal to be a multiple of scalar dof per cells"));
5649
5650 if (!dyadic_coefficients)
5651 {
5652 if (inverse_coefficients.size() != dofs_per_component)
5653 AssertDimension(n_desired_components * dofs_per_component,
5654 inverse_coefficients.size());
5655 }
5656 else
5657 {
5658 AssertDimension(n_desired_components * n_desired_components *
5659 dofs_per_component,
5660 inverse_coefficients.size());
5661 }
5662
5663 Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
5667
5669 dim,
5670 fe_degree + 1,
5671 fe_degree + 1,
5672 Number,
5673 Number2>
5674 evaluator({},
5675 {},
5676 fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
5677 given_degree + 1,
5678 given_degree + 1);
5679
5680 const Number *in = in_array;
5681 Number * out = out_array;
5682
5683 const Number *inv_coefficient = inverse_coefficients.data();
5684
5685 const unsigned int shift_coefficient =
5686 inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
5687 0;
5688
5689 const auto n_comp_outer = dyadic_coefficients ? 1 : n_desired_components;
5690 const auto n_comp_inner = dyadic_coefficients ? n_desired_components : 1;
5691
5692 for (unsigned int d = 0; d < n_comp_outer; ++d)
5693 {
5694 for (unsigned int di = 0; di < n_comp_inner; ++di)
5695 {
5696 const Number *in_ = in + di * dofs_per_component;
5697 Number * out_ = out + di * dofs_per_component;
5698 evaluator.template hessians<0, true, false>(in_, out_);
5699 if (dim > 1)
5700 evaluator.template hessians<1, true, false>(out_, out_);
5701 if (dim > 2)
5702 evaluator.template hessians<2, true, false>(out_, out_);
5703 }
5704 if (dyadic_coefficients)
5705 {
5706 const auto n_coeff_components =
5707 n_desired_components * n_desired_components;
5708 if (n_desired_components == dim)
5709 {
5710 for (unsigned int q = 0; q < dofs_per_component; ++q)
5711 vmult<dim>(&inv_coefficient[q * n_coeff_components],
5712 &in[q],
5713 &out[q],
5714 dofs_per_component);
5715 }
5716 else
5717 {
5718 for (unsigned int q = 0; q < dofs_per_component; ++q)
5719 vmult<-1>(&inv_coefficient[q * n_coeff_components],
5720 &in[q],
5721 &out[q],
5722 dofs_per_component,
5723 n_desired_components);
5724 }
5725 }
5726 else
5727 for (unsigned int q = 0; q < dofs_per_component; ++q)
5728 out[q] *= inv_coefficient[q];
5729
5730 for (unsigned int di = 0; di < n_comp_inner; ++di)
5731 {
5732 Number *out_ = out + di * dofs_per_component;
5733 if (dim > 2)
5734 evaluator.template hessians<2, false, false>(out_, out_);
5735 if (dim > 1)
5736 evaluator.template hessians<1, false, false>(out_, out_);
5737 evaluator.template hessians<0, false, false>(out_, out_);
5738 }
5739
5740 in += dofs_per_component;
5741 out += dofs_per_component;
5742 inv_coefficient += shift_coefficient;
5743 }
5744
5745 return false;
5746 }
5747
5748 private:
5749 template <int n_components>
5750 static inline void
5751 vmult(const Number * inverse_coefficients,
5752 const Number * src,
5753 Number * dst,
5754 const unsigned int dofs_per_component,
5755 const unsigned int n_given_components = 0)
5756 {
5757 const unsigned int n_desired_components =
5758 (n_components > -1) ? n_components : n_given_components;
5759
5760 std::array<Number, dim + 2> tmp = {};
5761 Assert(n_desired_components <= dim + 2,
5762 ExcMessage(
5763 "Number of components larger than dim+2 not supported."));
5764
5765 for (unsigned int d = 0; d < n_desired_components; ++d)
5766 tmp[d] = src[d * dofs_per_component];
5767
5768 for (unsigned int d1 = 0; d1 < n_desired_components; ++d1)
5769 {
5770 const Number *inv_coeff_row =
5771 &inverse_coefficients[d1 * n_desired_components];
5772 Number sum = inv_coeff_row[0] * tmp[0];
5773 for (unsigned int d2 = 1; d2 < n_desired_components; ++d2)
5774 sum += inv_coeff_row[d2] * tmp[d2];
5775 dst[d1 * dofs_per_component] = sum;
5776 }
5777 }
5778 };
5779
5780
5781
5788 template <int dim, typename Number>
5790 {
5791 template <int fe_degree, int n_q_points_1d>
5792 static bool
5793 run(const unsigned int n_desired_components,
5795 const Number * in_array,
5796 Number * out_array)
5797 {
5798 static const bool do_inplace =
5799 fe_degree > -1 && (fe_degree + 1 == n_q_points_1d);
5800
5804
5805 const auto &inverse_shape =
5806 do_inplace ?
5807 fe_eval.get_shape_info().data.front().inverse_shape_values_eo :
5808 fe_eval.get_shape_info().data.front().inverse_shape_values;
5809
5810 const std::size_t dofs_per_component =
5811 do_inplace ? Utilities::pow(fe_degree + 1, dim) :
5813 const std::size_t n_q_points = do_inplace ?
5814 Utilities::pow(fe_degree + 1, dim) :
5815 fe_eval.get_shape_info().n_q_points;
5816
5817 using Number2 =
5820 dim,
5821 fe_degree + 1,
5822 n_q_points_1d,
5823 Number,
5824 Number2>
5825 evaluator({},
5826 {},
5827 inverse_shape,
5828 fe_eval.get_shape_info().data.front().fe_degree + 1,
5829 fe_eval.get_shape_info().data.front().n_q_points_1d);
5830
5831 for (unsigned int d = 0; d < n_desired_components; ++d)
5832 {
5833 const Number *in = in_array + d * n_q_points;
5834 Number * out = out_array + d * dofs_per_component;
5835
5836 auto temp_1 = do_inplace ? out : fe_eval.get_scratch_data().begin();
5837 auto temp_2 = do_inplace ?
5838 out :
5839 (temp_1 + std::max(n_q_points, dofs_per_component));
5840
5841 if (dim == 3)
5842 {
5843 evaluator.template hessians<2, false, false>(in, temp_1);
5844 evaluator.template hessians<1, false, false>(temp_1, temp_2);
5845 evaluator.template hessians<0, false, false>(temp_2, out);
5846 }
5847 if (dim == 2)
5848 {
5849 evaluator.template hessians<1, false, false>(in, temp_1);
5850 evaluator.template hessians<0, false, false>(temp_1, out);
5851 }
5852 if (dim == 1)
5853 evaluator.template hessians<0, false, false>(in, out);
5854 }
5855 return false;
5856 }
5857 };
5858
5865 {
5866 template <int fe_degree, int n_q_points_1d>
5867 static bool
5869 {
5870 return fe_degree != -1;
5871 }
5872 };
5873} // end of namespace internal
5874
5875
5877
5878#endif
size_type size() const
iterator begin() const
Definition array_view.h:594
value_type * data() const noexcept
Definition array_view.h:553
std::size_t size() const
Definition array_view.h:576
std::uint8_t get_face_no(const unsigned int v=0) const
internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex get_dof_access_index() const
const ShapeInfoType & get_shape_info() const
Number JxW(const unsigned int q_point) const
const std::array< unsigned int, n_lanes > & get_cell_ids() const
const Number * begin_gradients() const
unsigned int get_subface_index() const
bool is_interior_face() const
ArrayView< Number > get_scratch_data() const
const Number * begin_values() const
std::uint8_t get_face_orientation(const unsigned int v=0) const
const Number * begin_hessians() const
void gather(const Number *base_ptr, const unsigned int *offsets)
void load(const OtherNumber *ptr)
#define DEAL_II_ALWAYS_INLINE
Definition config.h:106
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:138
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
unsigned int cell_index
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
EvaluationFlags
The EvaluationFlags enum.
constexpr T pow(const T base, const int iexp)
Definition utilities.h:447
T fixed_power(const T t)
Definition utilities.h:983
void do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
void adjust_for_face_orientation_per_lane(const unsigned int dim, const unsigned int n_components, const unsigned int v, const EvaluationFlags::EvaluationFlags flag, const unsigned int *orientation, const bool integrate, const std::size_t n_q_points, Number *tmp_values, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad=nullptr, VectorizedArrayType *hessians_quad=nullptr)
void fe_face_evaluation_process_and_io(Processor &proc, const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, typename Processor::Number2_ *global_vector_ptr, const std::vector< ArrayView< const typename Processor::Number2_ > > *sm_ptr, const EvaluationData &fe_eval, typename Processor::VectorizedArrayType_ *temp1)
void do_vectorized_scatter_add(const VectorizedArrayType src, const unsigned int *indices, Number2 *dst_ptr)
void do_vectorized_gather(const Number2 *src_ptr, const unsigned int *indices, VectorizedArrayType &dst)
void do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
void adjust_for_face_orientation(const unsigned int dim, const unsigned int n_components, const EvaluationFlags::EvaluationFlags flag, const unsigned int *orientation, const bool integrate, const std::size_t n_q_points, Number *tmp_values, Number *values_quad, Number *gradients_quad, Number *hessians_quad)
static const unsigned int invalid_unsigned_int
Definition types.h:213
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
static bool run(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array)
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
static bool run(const unsigned int n_desired_components, const FEEvaluationData< dim, Number, false > &fe_eval, const ArrayView< const Number > &inverse_coefficients, const bool dyadic_coefficients, const Number *in_array, Number *out_array)
static void vmult(const Number *inverse_coefficients, const Number *src, Number *dst, const unsigned int dofs_per_component, const unsigned int n_given_components=0)
static bool run(const unsigned int n_desired_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array)
static void do_mass(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const AlignedVector< Number > &coefficients, const Number *values_in, Number *scratch_data, Number *values_out)
static void do_forward(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const Number *values_in, Number *values_out, const unsigned int basis_size_1_variable=numbers::invalid_unsigned_int, const unsigned int basis_size_2_variable=numbers::invalid_unsigned_int)
static void do_backward(const unsigned int n_components, const AlignedVector< Number2 > &transformation_matrix, const bool add_into_result, Number *values_in, Number *values_out, const unsigned int basis_size_1_variable=numbers::invalid_unsigned_int, const unsigned int basis_size_2_variable=numbers::invalid_unsigned_int)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
static void do_evaluate(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &shape, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, Number *gradients_quad, Number *hessians_quad)
static void do_integrate(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &shape, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, Number *gradients_quad, const Number *hessians_quad, const bool add_into_values_array)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
static void evaluate_or_integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, OtherNumber *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, OtherNumber *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array=false)
static void evaluate_or_integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array, std::integral_constant< bool, false >)
static void evaluate_or_integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array, std::integral_constant< bool, true >)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static EvalType create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &shape_data)
static const EvaluatorVariant variant
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
EvaluatorTensorProduct< variant, dim, fe_degree+1, n_q_points_1d, Number, Number2 > Eval
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs_actual, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
static Eval create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number2 > *univariate_shape_data)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs_actual, FEEvaluationData< dim, Number, false > &fe_eval)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
void hermite_grad(T0 &temp_1, T0 &temp_2, const T1 &src_ptr_1, const T1 &src_ptr_2, const T2 &grad_weight)
void value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
void hermite_grad_vectorized(T0 &temp_1, T0 &temp_2, const T1 src_ptr_1, const T1 src_ptr_2, const T2 &grad_weight)
void hermite_grad_vectorized_indexed(T0 &temp_1, T0 &temp_2, const T1 src_ptr_1, const T1 src_ptr_2, const T2 &grad_weight, const T3 &indices_1, const T3 &indices_2)
typename VectorizedArrayType::value_type Number
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number2 *src_ptr, const std::vector< ArrayView< const Number2 > > *sm_ptr, FEEvaluationData< dim, VectorizedArrayType, true > &fe_eval)
static bool supports(const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::ShapeInfo< Number3 > &shape_info, const Number2 *vector_ptr, MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage)
void hermite_grad_vectorized(const T0 &temp_1, const T1 &temp_2, T2 dst_ptr_1, T3 dst_ptr_2, const T4 &grad_weight)
void hermite_grad_vectorized_indexed(const T0 &temp_1, const T0 &temp_2, T1 dst_ptr_1, T1 dst_ptr_2, const T2 &grad_weight, const T3 &indices_1, const T3 &indices_2)
void hermite_grad(const T0 &temp_1, const T0 &temp_2, T1 &dst_ptr_1, T1 &dst_ptr_2, const T2 &grad_weight)
void value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
typename VectorizedArrayType::value_type Number
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number2 *dst_ptr, const std::vector< ArrayView< const Number2 > > *sm_ptr, FEEvaluationData< dim, VectorizedArrayType, true > &fe_eval)
static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static void evaluate_or_integrate_in_face(const EvaluationFlags::EvaluationFlags evaluation_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, Number *scratch_data, const unsigned int subface_index, const unsigned int face_no)
static void evaluate_in_face_apply(Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, Number *scratch_data, const EvaluationFlags::EvaluationFlags evaluation_flag, const unsigned int face_direction, const unsigned int subface_index, std::integral_constant< bool, false >)
static EvalType create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const unsigned int subface_index, const unsigned int direction)
static void evaluate_in_face_apply(Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, Number *scratch_data, const EvaluationFlags::EvaluationFlags evaluation_flag, const unsigned int face_direction, const unsigned int subface_index, std::integral_constant< bool, true >)
typename FEEvaluationData< dim, Number, true >::shape_info_number_type Number2
EvaluatorTensorProduct< symmetric_evaluate ? evaluate_evenodd :evaluate_general, dim - 1, fe_degree+1, n_q_points_1d, Number, Number2 > Eval
static void evaluate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data, const unsigned int subface_index)
static Eval create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const unsigned int subface_index, const unsigned int direction)
static void integrate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, Number *values_dofs, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data, const unsigned int subface_index)
typename FEEvaluationData< dim, Number, true >::shape_info_number_type Number2
static EvalType create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &data, const unsigned int face_no)
static void interpolate_generic_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)
typename FEEvaluationData< dim, Number, true >::shape_info_number_type Number2
static void interpolate_generic_raviart_thomas_apply_face(const MatrixFreeFunctions::ShapeInfo< Number2 > &shape_info, const unsigned int face_no, const Number *input, Number *output)
static void interpolate_quadrature(const unsigned int n_components, const EvaluationFlags::EvaluationFlags flags, const MatrixFreeFunctions::ShapeInfo< Number2 > &shape_info, const Number *input, Number *output, const unsigned int face_no)
static void interpolate_generic(const unsigned int n_components, const Number *input, Number *output, const EvaluationFlags::EvaluationFlags flag, const unsigned int face_no, const unsigned int n_points_1d, const std::array< AlignedVector< Number2 >, 2 > &shape_data, const unsigned int dofs_per_component_on_cell, const unsigned int dofs_per_component_on_face)
static void interpolate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags flags, const MatrixFreeFunctions::ShapeInfo< Number2 > &shape_info, const Number *input, Number *output, const unsigned int face_no)
std::vector< UnivariateShapeData< Number > > data
Definition shape_info.h:441
::Table< 2, unsigned int > face_orientations_quad
Definition shape_info.h:580
std::array< AlignedVector< Number >, 2 > shape_data_on_face
Definition shape_info.h:253
std::array< AlignedVector< Number >, 2 > hessians_within_subface
Definition shape_info.h:282
std::array< AlignedVector< Number >, 2 > values_within_subface
Definition shape_info.h:270
std::array< AlignedVector< Number >, 2 > gradients_within_subface
Definition shape_info.h:276