Reference documentation for deal.II version GIT relicensing-478-g3275795167 2024-04-24 07:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
evaluation_kernels.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_matrix_free_evaluation_kernels_h
17#define dealii_matrix_free_evaluation_kernels_h
18
19#include <deal.II/base/config.h>
20
23
28
29
31
32
33namespace internal
34{
35 // Select evaluator type from element shape function type
36 template <MatrixFreeFunctions::ElementType element, bool is_long>
38 {};
39
40 template <bool is_long>
41 struct EvaluatorSelector<MatrixFreeFunctions::tensor_general, is_long>
42 {
43 static const EvaluatorVariant variant = evaluate_general;
44 };
45
46 template <>
47 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, false>
48 {
49 static const EvaluatorVariant variant = evaluate_symmetric;
50 };
51
52 template <>
53 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric, true>
54 {
55 static const EvaluatorVariant variant = evaluate_evenodd;
56 };
57
58 template <bool is_long>
59 struct EvaluatorSelector<MatrixFreeFunctions::truncated_tensor, is_long>
60 {
61 static const EvaluatorVariant variant = evaluate_general;
62 };
63
64 template <>
65 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,
66 false>
67 {
68 static const EvaluatorVariant variant = evaluate_general;
69 };
70
71 template <>
72 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0, true>
73 {
74 static const EvaluatorVariant variant = evaluate_evenodd;
75 };
76
77 template <bool is_long>
78 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_collocation,
79 is_long>
80 {
81 static const EvaluatorVariant variant = evaluate_evenodd;
82 };
83
84
85
108 int dim,
109 int fe_degree,
110 int n_q_points_1d,
111 typename Number>
113 {
115 EvaluatorSelector<type, (fe_degree + n_q_points_1d > 4)>::variant;
116 using Number2 =
118
120 dim,
121 fe_degree + 1,
122 n_q_points_1d,
123 Number,
124 Number2>;
125
126 static void
127 evaluate(const unsigned int n_components,
128 const EvaluationFlags::EvaluationFlags evaluation_flag,
129 const Number *values_dofs_actual,
131
132 static void
133 integrate(const unsigned int n_components,
134 const EvaluationFlags::EvaluationFlags integration_flag,
135 Number *values_dofs_actual,
137 const bool add_into_values_array);
138
139 static Eval
142 *univariate_shape_data)
143 {
145 return Eval(univariate_shape_data->shape_values_eo,
146 univariate_shape_data->shape_gradients_eo,
147 univariate_shape_data->shape_hessians_eo,
148 univariate_shape_data->fe_degree + 1,
149 univariate_shape_data->n_q_points_1d);
150 else
151 return Eval(univariate_shape_data->shape_values,
152 univariate_shape_data->shape_gradients,
153 univariate_shape_data->shape_hessians,
154 univariate_shape_data->fe_degree + 1,
155 univariate_shape_data->n_q_points_1d);
156 }
157 };
158
159
160
165 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
166 struct FEEvaluationImpl<MatrixFreeFunctions::tensor_none,
167 dim,
168 fe_degree,
169 n_q_points_1d,
170 Number>
171 {
172 static void
173 evaluate(const unsigned int n_components,
174 const EvaluationFlags::EvaluationFlags evaluation_flag,
175 const Number *values_dofs_actual,
177
178 static void
179 integrate(const unsigned int n_components,
180 const EvaluationFlags::EvaluationFlags integration_flag,
181 Number *values_dofs_actual,
183 const bool add_into_values_array);
184 };
185
186
187
189 int dim,
190 int fe_degree,
191 int n_q_points_1d,
192 typename Number>
193 inline void
195 const unsigned int n_components,
196 const EvaluationFlags::EvaluationFlags evaluation_flag,
197 const Number *values_dofs_actual,
199 {
200 if (evaluation_flag == EvaluationFlags::nothing)
201 return;
202
203 std::array<const MatrixFreeFunctions::UnivariateShapeData<Number2> *, 3>
204 univariate_shape_data;
205
206 const auto &shape_data = fe_eval.get_shape_info().data;
207
208 univariate_shape_data.fill(&shape_data.front());
209
210 if (shape_data.size() == dim)
211 for (int i = 1; i < dim; ++i)
212 univariate_shape_data[i] = &shape_data[i];
213
214 Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
215 Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
216 Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);
217
218 const unsigned int temp_size =
219 Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
220 0 :
221 (Eval::n_rows_of_product > Eval::n_columns_of_product ?
222 Eval::n_rows_of_product :
223 Eval::n_columns_of_product);
224 Number *temp1 = fe_eval.get_scratch_data().begin();
225 Number *temp2;
226 if (temp_size == 0)
227 {
228 temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
229 shape_data.front().fe_degree + 1),
230 Utilities::fixed_power<dim>(
231 shape_data.front().n_q_points_1d));
232 }
233 else
234 {
235 temp2 = temp1 + temp_size;
236 }
237
238 const std::size_t n_q_points = temp_size == 0 ?
239 fe_eval.get_shape_info().n_q_points :
240 Eval::n_columns_of_product;
241 const std::size_t dofs_per_comp =
243 Utilities::pow(shape_data.front().fe_degree + 1, dim) :
245 const Number *values_dofs = values_dofs_actual;
247 {
248 const std::size_t n_dofs_per_comp =
250 Number *values_dofs_tmp =
251 temp1 + 2 * (std::max(n_dofs_per_comp, n_q_points));
252 const int degree =
253 fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
254 for (unsigned int c = 0; c < n_components; ++c)
255 for (int i = 0, count_p = 0, count_q = 0;
256 i < (dim > 2 ? degree + 1 : 1);
257 ++i)
258 {
259 for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
260 {
261 for (int k = 0; k < degree + 1 - j - i;
262 ++k, ++count_p, ++count_q)
263 values_dofs_tmp[c * dofs_per_comp + count_q] =
264 values_dofs_actual[c * n_dofs_per_comp + count_p];
265 for (int k = degree + 1 - j - i; k < degree + 1;
266 ++k, ++count_q)
267 values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
268 }
269 for (int j = degree + 1 - i; j < degree + 1; ++j)
270 for (int k = 0; k < degree + 1; ++k, ++count_q)
271 values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
272 }
273 values_dofs = values_dofs_tmp;
274 }
275
276 Number *values_quad = fe_eval.begin_values();
277 Number *gradients_quad = fe_eval.begin_gradients();
278
279 switch (dim)
280 {
281 case 1:
282 for (unsigned int c = 0; c < n_components; ++c)
283 {
284 if (evaluation_flag & EvaluationFlags::values)
285 eval0.template values<0, true, false>(values_dofs, values_quad);
286 if (evaluation_flag & EvaluationFlags::gradients)
287 eval0.template gradients<0, true, false>(values_dofs,
288 gradients_quad);
289
290 // advance the next component in 1d array
291 values_dofs += dofs_per_comp;
292 values_quad += n_q_points;
293 gradients_quad += n_q_points;
294 }
295 break;
296
297 case 2:
298 for (unsigned int c = 0; c < n_components; ++c)
299 {
300 // grad x
301 if (evaluation_flag & EvaluationFlags::gradients)
302 {
303 eval0.template gradients<0, true, false>(values_dofs, temp1);
304 eval1.template values<1, true, false, 2>(temp1,
305 gradients_quad);
306 }
307
308 // grad y
309 eval0.template values<0, true, false>(values_dofs, temp1);
310 if (evaluation_flag & EvaluationFlags::gradients)
311 eval1.template gradients<1, true, false, 2>(temp1,
312 gradients_quad + 1);
313
314 // val: can use values applied in x
315 if (evaluation_flag & EvaluationFlags::values)
316 eval1.template values<1, true, false>(temp1, values_quad);
317
318 // advance to the next component in 1d array
319 values_dofs += dofs_per_comp;
320 values_quad += n_q_points;
321 gradients_quad += 2 * n_q_points;
322 }
323 break;
324
325 case 3:
326 for (unsigned int c = 0; c < n_components; ++c)
327 {
328 if (evaluation_flag & EvaluationFlags::gradients)
329 {
330 // grad x
331 eval0.template gradients<0, true, false>(values_dofs, temp1);
332 eval1.template values<1, true, false>(temp1, temp2);
333 eval2.template values<2, true, false, 3>(temp2,
334 gradients_quad);
335 }
336
337 // grad y
338 eval0.template values<0, true, false>(values_dofs, temp1);
339 if (evaluation_flag & EvaluationFlags::gradients)
340 {
341 eval1.template gradients<1, true, false>(temp1, temp2);
342 eval2.template values<2, true, false, 3>(temp2,
343 gradients_quad + 1);
344 }
345
346 // grad z: can use the values applied in x direction stored in
347 // temp1
348 eval1.template values<1, true, false>(temp1, temp2);
349 if (evaluation_flag & EvaluationFlags::gradients)
350 eval2.template gradients<2, true, false, 3>(temp2,
351 gradients_quad + 2);
352
353 // val: can use the values applied in x & y direction stored in
354 // temp2
355 if (evaluation_flag & EvaluationFlags::values)
356 eval2.template values<2, true, false>(temp2, values_quad);
357
358 // advance to the next component in 1d array
359 values_dofs += dofs_per_comp;
360 values_quad += n_q_points;
361 gradients_quad += 3 * n_q_points;
362 }
363 break;
364
365 default:
367 }
368
369 // case additional dof for FE_Q_DG0: add values; gradients and second
370 // derivatives evaluate to zero
372 (evaluation_flag & EvaluationFlags::values))
373 {
374 values_quad -= n_components * n_q_points;
375 values_dofs -= n_components * dofs_per_comp;
376 for (std::size_t c = 0; c < n_components; ++c)
377 for (std::size_t q = 0; q < n_q_points; ++q)
378 values_quad[c * n_q_points + q] +=
379 values_dofs[(c + 1) * dofs_per_comp - 1];
380 }
381 }
382
383
384
386 int dim,
387 int fe_degree,
388 int n_q_points_1d,
389 typename Number>
390 inline void
392 const unsigned int n_components,
393 const EvaluationFlags::EvaluationFlags integration_flag,
394 Number *values_dofs_actual,
396 const bool add_into_values_array)
397 {
398 std::array<const MatrixFreeFunctions::UnivariateShapeData<Number2> *, 3>
399 univariate_shape_data;
400
401 const auto &shape_data = fe_eval.get_shape_info().data;
402 univariate_shape_data.fill(&shape_data.front());
403
404 if (shape_data.size() == dim)
405 for (int i = 1; i < dim; ++i)
406 univariate_shape_data[i] = &shape_data[i];
407
408 Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
409 Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
410 Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);
411
412 const unsigned int temp_size =
413 Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
414 0 :
415 (Eval::n_rows_of_product > Eval::n_columns_of_product ?
416 Eval::n_rows_of_product :
417 Eval::n_columns_of_product);
418 Number *temp1 = fe_eval.get_scratch_data().begin();
419 Number *temp2;
420 if (temp_size == 0)
421 {
422 temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
423 shape_data.front().fe_degree + 1),
424 Utilities::fixed_power<dim>(
425 shape_data.front().n_q_points_1d));
426 }
427 else
428 {
429 temp2 = temp1 + temp_size;
430 }
431
432 const std::size_t n_q_points = temp_size == 0 ?
433 fe_eval.get_shape_info().n_q_points :
434 Eval::n_columns_of_product;
435 const unsigned int dofs_per_comp =
437 Utilities::fixed_power<dim>(shape_data.front().fe_degree + 1) :
439 // expand dof_values to tensor product for truncated tensor products
440 Number *values_dofs =
442 temp1 + 2 * (std::max<std::size_t>(
444 n_q_points)) :
445 values_dofs_actual;
446
447 Number *values_quad = fe_eval.begin_values();
448 Number *gradients_quad = fe_eval.begin_gradients();
449
450 switch (dim)
451 {
452 case 1:
453 for (unsigned int c = 0; c < n_components; ++c)
454 {
455 if (integration_flag & EvaluationFlags::values)
456 {
457 if (add_into_values_array == false)
458 eval0.template values<0, false, false>(values_quad,
459 values_dofs);
460 else
461 eval0.template values<0, false, true>(values_quad,
462 values_dofs);
463 }
464 if (integration_flag & EvaluationFlags::gradients)
465 {
466 if (integration_flag & EvaluationFlags::values ||
467 add_into_values_array == true)
468 eval0.template gradients<0, false, true>(gradients_quad,
469 values_dofs);
470 else
471 eval0.template gradients<0, false, false>(gradients_quad,
472 values_dofs);
473 }
474
475 // advance to the next component in 1d array
476 values_dofs += dofs_per_comp;
477 values_quad += n_q_points;
478 gradients_quad += n_q_points;
479 }
480 break;
481
482 case 2:
483 for (unsigned int c = 0; c < n_components; ++c)
484 {
485 if ((integration_flag & EvaluationFlags::values) &&
486 !(integration_flag & EvaluationFlags::gradients))
487 {
488 eval1.template values<1, false, false>(values_quad, temp1);
489 if (add_into_values_array == false)
490 eval0.template values<0, false, false>(temp1, values_dofs);
491 else
492 eval0.template values<0, false, true>(temp1, values_dofs);
493 }
494 if (integration_flag & EvaluationFlags::gradients)
495 {
496 eval1.template gradients<1, false, false, 2>(gradients_quad +
497 1,
498 temp1);
499 if (integration_flag & EvaluationFlags::values)
500 eval1.template values<1, false, true>(values_quad, temp1);
501 if (add_into_values_array == false)
502 eval0.template values<0, false, false>(temp1, values_dofs);
503 else
504 eval0.template values<0, false, true>(temp1, values_dofs);
505 eval1.template values<1, false, false, 2>(gradients_quad,
506 temp1);
507 eval0.template gradients<0, false, true>(temp1, values_dofs);
508 }
509
510 // advance to the next component in 1d array
511 values_dofs += dofs_per_comp;
512 values_quad += n_q_points;
513 gradients_quad += 2 * n_q_points;
514 }
515 break;
516
517 case 3:
518 for (unsigned int c = 0; c < n_components; ++c)
519 {
520 if ((integration_flag & EvaluationFlags::values) &&
521 !(integration_flag & EvaluationFlags::gradients))
522 {
523 eval2.template values<2, false, false>(values_quad, temp1);
524 eval1.template values<1, false, false>(temp1, temp2);
525 if (add_into_values_array == false)
526 eval0.template values<0, false, false>(temp2, values_dofs);
527 else
528 eval0.template values<0, false, true>(temp2, values_dofs);
529 }
530 if (integration_flag & EvaluationFlags::gradients)
531 {
532 eval2.template gradients<2, false, false, 3>(gradients_quad +
533 2,
534 temp1);
535 if (integration_flag & EvaluationFlags::values)
536 eval2.template values<2, false, true>(values_quad, temp1);
537 eval1.template values<1, false, false>(temp1, temp2);
538 eval2.template values<2, false, false, 3>(gradients_quad + 1,
539 temp1);
540 eval1.template gradients<1, false, true>(temp1, temp2);
541 if (add_into_values_array == false)
542 eval0.template values<0, false, false>(temp2, values_dofs);
543 else
544 eval0.template values<0, false, true>(temp2, values_dofs);
545 eval2.template values<2, false, false, 3>(gradients_quad,
546 temp1);
547 eval1.template values<1, false, false>(temp1, temp2);
548 eval0.template gradients<0, false, true>(temp2, values_dofs);
549 }
550
551 // advance to the next component in 1d array
552 values_dofs += dofs_per_comp;
553 values_quad += n_q_points;
554 gradients_quad += 3 * n_q_points;
555 }
556 break;
557
558 default:
560 }
561
562 // case FE_Q_DG0: add values, gradients and second derivatives are zero
564 {
565 values_dofs -= n_components * dofs_per_comp - dofs_per_comp + 1;
566 values_quad -= n_components * n_q_points;
567 if (integration_flag & EvaluationFlags::values)
568 for (unsigned int c = 0; c < n_components; ++c)
569 {
570 values_dofs[0] = values_quad[0];
571 for (unsigned int q = 1; q < n_q_points; ++q)
572 values_dofs[0] += values_quad[q];
573 values_dofs += dofs_per_comp;
574 values_quad += n_q_points;
575 }
576 else
577 {
578 for (unsigned int c = 0; c < n_components; ++c)
579 values_dofs[c * dofs_per_comp] = Number();
580 values_dofs += n_components * dofs_per_comp;
581 }
582 }
583
585 {
586 const std::size_t n_dofs_per_comp =
588 values_dofs -= dofs_per_comp * n_components;
589 const int degree =
590 fe_degree != -1 ? fe_degree : shape_data.front().fe_degree;
591 for (unsigned int c = 0; c < n_components; ++c)
592 for (int i = 0, count_p = 0, count_q = 0;
593 i < (dim > 2 ? degree + 1 : 1);
594 ++i)
595 {
596 for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
597 {
598 for (int k = 0; k < degree + 1 - j - i;
599 ++k, ++count_p, ++count_q)
600 values_dofs_actual[c * n_dofs_per_comp + count_p] =
601 values_dofs[c * dofs_per_comp + count_q];
602 count_q += j + i;
603 }
604 count_q += i * (degree + 1);
605 }
606 }
607 }
608
609
610
611 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
612 inline void
615 dim,
616 fe_degree,
617 n_q_points_1d,
618 Number>::evaluate(const unsigned int n_components,
619 const EvaluationFlags::EvaluationFlags evaluation_flag,
620 const Number *values_dofs_actual,
622 {
623 Assert(!(evaluation_flag & EvaluationFlags::hessians), ExcNotImplemented());
624
625 const std::size_t n_dofs =
627 const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
628
629 const auto &shape_data = fe_eval.get_shape_info().data;
630
631 using Number2 =
633 using Eval =
635
636 if (evaluation_flag & EvaluationFlags::values)
637 {
638 const auto *const shape_values = shape_data.front().shape_values.data();
639 auto *values_quad_ptr = fe_eval.begin_values();
640 const auto *values_dofs_actual_ptr = values_dofs_actual;
641
642 Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
643 for (unsigned int c = 0; c < n_components; ++c)
644 {
645 eval.template values<0, true, false>(values_dofs_actual_ptr,
646 values_quad_ptr);
647
648 values_quad_ptr += n_q_points;
649 values_dofs_actual_ptr += n_dofs;
650 }
651 }
652
653 if (evaluation_flag & EvaluationFlags::gradients)
654 {
655 const auto *const shape_gradients =
656 shape_data.front().shape_gradients.data();
657 auto *gradients_quad_ptr = fe_eval.begin_gradients();
658 const auto *values_dofs_actual_ptr = values_dofs_actual;
659
660 for (unsigned int c = 0; c < n_components; ++c)
661 {
662 for (unsigned int d = 0; d < dim; ++d)
663 {
664 Eval eval(nullptr,
665 shape_gradients + n_q_points * n_dofs * d,
666 nullptr,
667 n_dofs,
668 n_q_points);
669
670 eval.template gradients<0, true, false, dim>(
671 values_dofs_actual_ptr, gradients_quad_ptr + d);
672 }
673 gradients_quad_ptr += n_q_points * dim;
674 values_dofs_actual_ptr += n_dofs;
675 }
676 }
677 }
678
679
680
681 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
682 inline void
685 dim,
686 fe_degree,
687 n_q_points_1d,
688 Number>::integrate(const unsigned int n_components,
689 const EvaluationFlags::EvaluationFlags integration_flag,
690 Number *values_dofs_actual,
692 const bool add_into_values_array)
693 {
694 Assert(!(integration_flag & EvaluationFlags::hessians),
696
697 const std::size_t n_dofs =
699 const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
700
701 const auto &shape_data = fe_eval.get_shape_info().data;
702
703 using Number2 =
705 using Eval =
707
708 if (integration_flag & EvaluationFlags::values)
709 {
710 const auto *const shape_values = shape_data.front().shape_values.data();
711 auto *values_quad_ptr = fe_eval.begin_values();
712 auto *values_dofs_actual_ptr = values_dofs_actual;
713
714 Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
715 for (unsigned int c = 0; c < n_components; ++c)
716 {
717 if (add_into_values_array == false)
718 eval.template values<0, false, false>(values_quad_ptr,
719 values_dofs_actual_ptr);
720 else
721 eval.template values<0, false, true>(values_quad_ptr,
722 values_dofs_actual_ptr);
723
724 values_quad_ptr += n_q_points;
725 values_dofs_actual_ptr += n_dofs;
726 }
727 }
728
729 if (integration_flag & EvaluationFlags::gradients)
730 {
731 const auto *const shape_gradients =
732 shape_data.front().shape_gradients.data();
733 auto *gradients_quad_ptr = fe_eval.begin_gradients();
734 auto *values_dofs_actual_ptr = values_dofs_actual;
735
736 for (unsigned int c = 0; c < n_components; ++c)
737 {
738 for (unsigned int d = 0; d < dim; ++d)
739 {
740 Eval eval(nullptr,
741 shape_gradients + n_q_points * n_dofs * d,
742 nullptr,
743 n_dofs,
744 n_q_points);
745
746 if ((add_into_values_array == false &&
747 !(integration_flag & EvaluationFlags::values)) &&
748 d == 0)
749 eval.template gradients<0, false, false, dim>(
750 gradients_quad_ptr + d, values_dofs_actual_ptr);
751 else
752 eval.template gradients<0, false, true, dim>(
753 gradients_quad_ptr + d, values_dofs_actual_ptr);
754 }
755 gradients_quad_ptr += n_q_points * dim;
756 values_dofs_actual_ptr += n_dofs;
757 }
758 }
759 }
760
761
762
772 template <EvaluatorVariant variant,
773 EvaluatorQuantity quantity,
774 int dim,
775 int basis_size_1,
776 int basis_size_2>
778 {
779 static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
780 "The second dimension must not be smaller than the first");
781
804 template <typename Number, typename Number2>
805#ifndef DEBUG
807#endif
808 static void
809 do_forward(const unsigned int n_components,
810 const AlignedVector<Number2> &transformation_matrix,
811 const Number *values_in,
812 Number *values_out,
813 const unsigned int basis_size_1_variable =
815 const unsigned int basis_size_2_variable =
817 {
818 Assert(
819 basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
820 ExcMessage("The second dimension must not be smaller than the first"));
821
823
824 // we do recursion until dim==1 or dim==2 and we have
825 // basis_size_1==basis_size_2. The latter optimization increases
826 // optimization possibilities for the compiler but does only work for
827 // aliased pointers if the sizes are equal.
828 constexpr int next_dim = (dim == 1 || (dim == 2 && basis_size_1 > 0 &&
829 basis_size_1 == basis_size_2)) ?
830 dim :
831 dim - 1;
832
834 dim,
835 basis_size_1,
836 (basis_size_1 == 0 ? 0 : basis_size_2),
837 Number,
838 Number2>
839 eval_val(transformation_matrix,
840 {},
841 {},
842 basis_size_1_variable,
843 basis_size_2_variable);
844 const unsigned int np_1 =
845 basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
846 const unsigned int np_2 =
847 basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
848 Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
849 ExcMessage("Cannot transform with 0-point basis"));
850 Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
851 ExcMessage("Cannot transform with 0-point basis"));
852
853 // run loop backwards to ensure correctness if values_in aliases with
854 // values_out in case with basis_size_1 < basis_size_2
855 values_in = values_in + n_components * Utilities::fixed_power<dim>(np_1);
856 values_out =
857 values_out + n_components * Utilities::fixed_power<dim>(np_2);
858 for (unsigned int c = n_components; c != 0; --c)
859 {
860 values_in -= Utilities::fixed_power<dim>(np_1);
861 values_out -= Utilities::fixed_power<dim>(np_2);
862 if (next_dim < dim)
863 for (unsigned int q = np_1; q != 0; --q)
865 quantity,
866 next_dim,
867 basis_size_1,
868 basis_size_2>::
869 do_forward(1,
870 transformation_matrix,
871 values_in +
872 (q - 1) * Utilities::fixed_power<next_dim>(np_1),
873 values_out +
874 (q - 1) * Utilities::fixed_power<next_dim>(np_2),
875 basis_size_1_variable,
876 basis_size_2_variable);
877
878 // the recursion stops if dim==1 or if dim==2 and
879 // basis_size_1==basis_size_2 (the latter is used because the
880 // compiler generates nicer code)
881 if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
882 {
883 eval_val.template values<0, true, false>(values_in, values_out);
884 eval_val.template values<1, true, false>(values_out, values_out);
885 }
886 else if (dim == 1)
887 eval_val.template values<dim - 1, true, false>(values_in,
888 values_out);
889 else
890 eval_val.template values<dim - 1, true, false>(values_out,
891 values_out);
892 }
893 }
894
925 template <typename Number, typename Number2>
926#ifndef DEBUG
928#endif
929 static void
930 do_backward(const unsigned int n_components,
931 const AlignedVector<Number2> &transformation_matrix,
932 const bool add_into_result,
933 Number *values_in,
934 Number *values_out,
935 const unsigned int basis_size_1_variable =
937 const unsigned int basis_size_2_variable =
939 {
940 Assert(
941 basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
942 ExcMessage("The second dimension must not be smaller than the first"));
943 Assert(add_into_result == false || values_in != values_out,
945 "Input and output cannot alias with each other when "
946 "adding the result of the basis change to existing data"));
947
948 Assert(quantity == EvaluatorQuantity::value ||
949 quantity == EvaluatorQuantity::hessian,
951
952 constexpr int next_dim =
953 (dim > 2 ||
954 ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
955 dim - 1 :
956 dim;
958 dim,
959 basis_size_1,
960 (basis_size_1 == 0 ? 0 : basis_size_2),
961 Number,
962 Number2>
963 eval_val(transformation_matrix,
964 transformation_matrix,
965 transformation_matrix,
966 basis_size_1_variable,
967 basis_size_2_variable);
968 const unsigned int np_1 =
969 basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
970 const unsigned int np_2 =
971 basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
972 Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
973 ExcMessage("Cannot transform with 0-point basis"));
974 Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
975 ExcMessage("Cannot transform with 0-point basis"));
976
977 for (unsigned int c = 0; c < n_components; ++c)
978 {
979 if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
980 {
981 if (quantity == EvaluatorQuantity::value)
982 eval_val.template values<1, false, false>(values_in, values_in);
983 else
984 eval_val.template hessians<1, false, false>(values_in,
985 values_in);
986
987 if (add_into_result)
988 {
989 if (quantity == EvaluatorQuantity::value)
990 eval_val.template values<0, false, true>(values_in,
991 values_out);
992 else
993 eval_val.template hessians<0, false, true>(values_in,
994 values_out);
995 }
996 else
997 {
998 if (quantity == EvaluatorQuantity::value)
999 eval_val.template values<0, false, false>(values_in,
1000 values_out);
1001 else
1002 eval_val.template hessians<0, false, false>(values_in,
1003 values_out);
1004 }
1005 }
1006 else
1007 {
1008 if (dim == 1 && add_into_result)
1009 {
1010 if (quantity == EvaluatorQuantity::value)
1011 eval_val.template values<0, false, true>(values_in,
1012 values_out);
1013 else
1014 eval_val.template hessians<0, false, true>(values_in,
1015 values_out);
1016 }
1017 else if (dim == 1)
1018 {
1019 if (quantity == EvaluatorQuantity::value)
1020 eval_val.template values<0, false, false>(values_in,
1021 values_out);
1022 else
1023 eval_val.template hessians<0, false, false>(values_in,
1024 values_out);
1025 }
1026 else
1027 {
1028 if (quantity == EvaluatorQuantity::value)
1029 eval_val.template values<dim - 1, false, false>(values_in,
1030 values_in);
1031 else
1032 eval_val.template hessians<dim - 1, false, false>(
1033 values_in, values_in);
1034 }
1035 }
1036 if (next_dim < dim)
1037 for (unsigned int q = 0; q < np_1; ++q)
1039 quantity,
1040 next_dim,
1041 basis_size_1,
1042 basis_size_2>::
1043 do_backward(1,
1044 transformation_matrix,
1045 add_into_result,
1046 values_in +
1047 q * Utilities::fixed_power<next_dim>(np_2),
1048 values_out +
1049 q * Utilities::fixed_power<next_dim>(np_1),
1050 basis_size_1_variable,
1051 basis_size_2_variable);
1052
1053 values_in += Utilities::fixed_power<dim>(np_2);
1054 values_out += Utilities::fixed_power<dim>(np_1);
1055 }
1056 }
1057
1078 template <typename Number, typename Number2>
1079 static void
1080 do_mass(const unsigned int n_components,
1081 const AlignedVector<Number2> &transformation_matrix,
1082 const AlignedVector<Number> &coefficients,
1083 const Number *values_in,
1084 Number *scratch_data,
1085 Number *values_out)
1086 {
1087 constexpr int next_dim = dim > 1 ? dim - 1 : dim;
1088 Number *my_scratch =
1089 basis_size_1 != basis_size_2 ? scratch_data : values_out;
1090
1091 const unsigned int size_per_component = Utilities::pow(basis_size_2, dim);
1092 Assert(coefficients.size() == size_per_component ||
1093 coefficients.size() == n_components * size_per_component,
1094 ExcDimensionMismatch(coefficients.size(), size_per_component));
1095 const unsigned int stride =
1096 coefficients.size() == size_per_component ? 0 : 1;
1097
1098 for (unsigned int q = basis_size_1; q != 0; --q)
1100 variant,
1102 next_dim,
1103 basis_size_1,
1104 basis_size_2>::do_forward(n_components,
1105 transformation_matrix,
1106 values_in +
1107 (q - 1) *
1108 Utilities::pow(basis_size_1, dim - 1),
1109 my_scratch +
1110 (q - 1) *
1111 Utilities::pow(basis_size_2, dim - 1));
1112 EvaluatorTensorProduct<variant,
1113 dim,
1114 basis_size_1,
1115 basis_size_2,
1116 Number,
1117 Number2>
1118 eval_val(transformation_matrix);
1119 const unsigned int n_inner_blocks =
1120 (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
1121 const unsigned int n_blocks = Utilities::pow(basis_size_2, dim - 1);
1122 for (unsigned int ii = 0; ii < n_blocks; ii += n_inner_blocks)
1123 for (unsigned int c = 0; c < n_components; ++c)
1124 {
1125 for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1126 eval_val.template values_one_line<dim - 1, true, false>(
1127 my_scratch + i, my_scratch + i);
1128 for (unsigned int q = 0; q < basis_size_2; ++q)
1129 for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1130 my_scratch[i + q * n_blocks + c * size_per_component] *=
1131 coefficients[i + q * n_blocks +
1132 c * stride * size_per_component];
1133 for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
1134 eval_val.template values_one_line<dim - 1, false, false>(
1135 my_scratch + i, my_scratch + i);
1136 }
1137 for (unsigned int q = 0; q < basis_size_1; ++q)
1140 next_dim,
1141 basis_size_1,
1142 basis_size_2>::
1143 do_backward(n_components,
1144 transformation_matrix,
1145 false,
1146 my_scratch + q * Utilities::pow(basis_size_2, dim - 1),
1147 values_out + q * Utilities::pow(basis_size_1, dim - 1));
1148 }
1149 };
1150
1151
1152
1160 template <int n_points_1d, int dim, typename Number, typename Number2>
1161 inline void
1164 const Number *values,
1165 Number *gradients)
1166 {
1168 (n_points_1d + 1) / 2 * n_points_1d);
1169
1171 dim,
1172 n_points_1d,
1173 n_points_1d,
1174 Number,
1175 Number2>
1176 eval({}, shape.shape_gradients_collocation_eo, {});
1178 2,
1179 n_points_1d,
1180 n_points_1d,
1181 Number,
1182 Number2>
1183 eval_2d({}, shape.shape_gradients_collocation_eo, {});
1184
1185 if (dim == 1)
1186 eval.template gradients<0, true, false>(values, gradients);
1187 else
1188 {
1189 if (dim > 2)
1190 eval.template gradients<2, true, false, dim>(values, gradients + 2);
1191 constexpr unsigned int loop_bound = (dim > 2 ? n_points_1d : 1);
1192 constexpr unsigned int n_points_2d = n_points_1d * n_points_1d;
1193 const Number *in = values + (loop_bound - 1) * n_points_2d;
1194 Number *out = gradients + (loop_bound - 1) * dim * n_points_2d;
1195 for (unsigned int l = 0; l < loop_bound; ++l)
1196 {
1197 eval_2d.template gradients<0, true, false, dim>(in, out);
1198 eval_2d.template gradients<1, true, false, dim>(in, out + 1);
1199 in -= n_points_2d;
1200 out -= dim * n_points_2d;
1201 }
1202 }
1203 }
1204
1205
1206
1214 template <int n_points_1d, int dim, typename Number, typename Number2>
1215 inline void
1218 Number *values,
1219 const Number *gradients,
1220 const bool add_into_values_array)
1221 {
1223 (n_points_1d + 1) / 2 * n_points_1d);
1224
1226 dim,
1227 n_points_1d,
1228 n_points_1d,
1229 Number,
1230 Number2>
1231 eval({}, shape.shape_gradients_collocation_eo, {});
1233 2,
1234 n_points_1d,
1235 n_points_1d,
1236 Number,
1237 Number2>
1238 eval_2d({}, shape.shape_gradients_collocation_eo, {});
1239
1240 if (dim == 1)
1241 {
1242 if (add_into_values_array)
1243 eval.template gradients<0, false, true>(gradients, values);
1244 else
1245 eval.template gradients<0, false, false>(gradients, values);
1246 }
1247 else
1248 {
1249 constexpr unsigned int loop_bound = (dim > 2 ? n_points_1d : 1);
1250 constexpr unsigned int n_points_2d = n_points_1d * n_points_1d;
1251
1252 const Number *in = gradients + (loop_bound - 1) * dim * n_points_2d;
1253 Number *out = values + (loop_bound - 1) * n_points_2d;
1254 for (unsigned int l = 0; l < loop_bound; ++l)
1255 {
1256 if (add_into_values_array)
1257 eval_2d.template gradients<0, false, true, dim>(in, out);
1258 else
1259 eval_2d.template gradients<0, false, false, dim>(in, out);
1260 eval_2d.template gradients<1, false, true, dim>(in + 1, out);
1261 in -= dim * n_points_2d;
1262 out -= n_points_2d;
1263 }
1264 }
1265 if (dim > 2)
1266 eval.template gradients<2, false, true, dim>(gradients + 2, values);
1267 }
1268
1269
1270
1277 template <int n_points_1d, int dim, typename Number>
1278 inline void
1279 evaluate_hessians_collocation(const unsigned int n_components,
1281 {
1282 using Number2 =
1284
1285 // might have non-symmetric quadrature formula, so use the more
1286 // conservative 'evaluate_general' scheme rather than 'even_odd' as the
1287 // Hessians are not used very often
1289 fe_eval.get_shape_info().data[0];
1291 data.n_q_points_1d * data.n_q_points_1d);
1293 dim,
1294 n_points_1d,
1295 n_points_1d,
1296 Number,
1297 Number2>
1298 eval({},
1301 data.n_q_points_1d,
1302 data.n_q_points_1d);
1303
1304 const Number *values = fe_eval.begin_values();
1305 Number *hessians = fe_eval.begin_hessians();
1306 Number *scratch = fe_eval.get_scratch_data().begin();
1307 const std::size_t n_points = fe_eval.get_shape_info().n_q_points;
1308 for (unsigned int comp = 0; comp < n_components; ++comp)
1309 {
1310 // xx derivative
1311 eval.template hessians<0, true, false>(values, hessians);
1312 if (dim > 1)
1313 {
1314 // xy derivative: we might or might not have the gradients already
1315 // computed elsewhere, but we recompute them here since it adds
1316 // only moderate extra work (at most 25%)
1317 eval.template gradients<0, true, false>(values, scratch);
1318 eval.template gradients<1, true, false>(scratch,
1319 hessians + dim * n_points);
1320 // yy derivative
1321 eval.template hessians<1, true, false>(values, hessians + n_points);
1322 }
1323 if (dim > 2)
1324 {
1325 // xz derivative
1326 eval.template gradients<2, true, false>(scratch,
1327 hessians + 4 * n_points);
1328 // yz derivative
1329 eval.template gradients<1, true, false>(values, scratch);
1330 eval.template gradients<2, true, false>(scratch,
1331 hessians + 5 * n_points);
1332 // zz derivative
1333 eval.template hessians<2, true, false>(values,
1334 hessians + 2 * n_points);
1335 }
1336
1337 values += n_points;
1338 hessians += (dim * (dim + 1)) / 2 * n_points;
1339 }
1340 }
1341
1342
1343
1350 template <int n_q_points_1d, int dim, typename Number>
1351 inline void
1352 integrate_hessians_collocation(const unsigned int n_components,
1354 const bool add_into_values_array)
1355 {
1356 using Number2 =
1358
1360 fe_eval.get_shape_info().data[0];
1362 data.n_q_points_1d * data.n_q_points_1d);
1364 dim,
1365 n_q_points_1d,
1366 n_q_points_1d,
1367 Number,
1368 Number2>
1369 eval({},
1372 data.n_q_points_1d,
1373 data.n_q_points_1d);
1374 Number *values = fe_eval.begin_values();
1375 const Number *hessians = fe_eval.begin_hessians();
1376 Number *scratch = fe_eval.get_scratch_data().begin();
1377 const std::size_t n_points = fe_eval.get_shape_info().n_q_points;
1378
1379 for (unsigned int comp = 0; comp < n_components; ++comp)
1380 {
1381 // xx derivative
1382 if (add_into_values_array == true)
1383 eval.template hessians<0, false, true>(hessians, values);
1384 else
1385 eval.template hessians<0, false, false>(hessians, values);
1386
1387 // yy derivative
1388 if (dim > 1)
1389 eval.template hessians<1, false, true>(hessians + n_points, values);
1390 if (dim > 2)
1391 {
1392 // zz derivative
1393 eval.template hessians<2, false, true>(hessians + 2 * n_points,
1394 values);
1395 // yz derivative
1396 eval.template gradients<2, false, false>(hessians + 5 * n_points,
1397 scratch);
1398 eval.template gradients<1, false, true>(scratch, values);
1399
1400 // xz derivative
1401 eval.template gradients<2, false, false>(hessians + 4 * n_points,
1402 scratch);
1403 }
1404
1405 if (dim > 1)
1406 {
1407 // xy derivative, combined with xz in 3d
1408 eval.template gradients<1, false, (dim > 2)>(hessians +
1409 dim * n_points,
1410 scratch);
1411 eval.template gradients<0, false, true>(scratch, values);
1412 }
1413
1414 values += n_points;
1415 hessians += (dim * (dim + 1)) / 2 * n_points;
1416 }
1417 }
1418
1419
1420
1427 template <int dim, typename Number>
1428 void
1429 evaluate_hessians_slow(const unsigned int n_components,
1430 const Number *values_dofs,
1432 {
1433 const auto &univariate_shape_data = fe_eval.get_shape_info().data;
1434 using Impl =
1436 using Eval = typename Impl::Eval;
1437 Eval eval0 =
1438 Impl::create_evaluator_tensor_product(&univariate_shape_data[0]);
1439 Eval eval1 = Impl::create_evaluator_tensor_product(
1440 &univariate_shape_data[std::min<int>(1,
1441 univariate_shape_data.size() - 1)]);
1442 Eval eval2 = Impl::create_evaluator_tensor_product(
1443 &univariate_shape_data[std::min<int>(2,
1444 univariate_shape_data.size() - 1)]);
1445
1446 const unsigned int n_points = fe_eval.get_shape_info().n_q_points;
1447 Number *tmp1 = fe_eval.get_scratch_data().begin();
1448 Number *tmp2 =
1449 tmp1 + std::max(Utilities::fixed_power<dim>(
1450 univariate_shape_data.front().fe_degree + 1),
1451 Utilities::fixed_power<dim>(
1452 univariate_shape_data.front().n_q_points_1d));
1453 Number *hessians = fe_eval.begin_hessians();
1454
1455 for (unsigned int comp = 0; comp < n_components;
1456 ++comp,
1457 hessians += n_points * dim * (dim + 1) / 2,
1458 values_dofs +=
1460 switch (dim)
1461 {
1462 case 1:
1463 eval0.template hessians<0, true, false>(values_dofs, hessians);
1464 break;
1465 case 2:
1466 // xx derivative
1467 eval0.template hessians<0, true, false>(values_dofs, tmp1);
1468 eval1.template values<1, true, false>(tmp1, hessians);
1469 // xy derivative
1470 eval0.template gradients<0, true, false>(values_dofs, tmp1);
1471 eval1.template gradients<1, true, false>(tmp1,
1472 hessians + 2 * n_points);
1473 // yy derivative
1474 eval0.template values<0, true, false>(values_dofs, tmp1);
1475 eval1.template hessians<1, true, false>(tmp1, hessians + n_points);
1476 break;
1477 case 3:
1478 // xx derivative
1479 eval0.template hessians<0, true, false>(values_dofs, tmp1);
1480 eval1.template values<1, true, false>(tmp1, tmp2);
1481 eval2.template values<2, true, false>(tmp2, hessians);
1482 // xy derivative
1483 eval0.template gradients<0, true, false>(values_dofs, tmp1);
1484 eval1.template gradients<1, true, false>(tmp1, tmp2);
1485 eval2.template values<2, true, false>(tmp2,
1486 hessians + 3 * n_points);
1487 // xz derivative
1488 eval1.template values<1, true, false>(tmp1, tmp2);
1489 eval2.template gradients<2, true, false>(tmp2,
1490 hessians + 4 * n_points);
1491 // yy derivative
1492 eval0.template values<0, true, false>(values_dofs, tmp1);
1493 eval1.template hessians<1, true, false>(tmp1, tmp2);
1494 eval2.template values<2, true, false>(tmp2, hessians + n_points);
1495 // yz derivative
1496 eval1.template gradients<1, true, false>(tmp1, tmp2);
1497 eval2.template gradients<2, true, false>(tmp2,
1498 hessians + 5 * n_points);
1499 // zz derivative
1500 eval1.template values<1, true, false>(tmp1, tmp2);
1501 eval2.template hessians<2, true, false>(tmp2,
1502 hessians + 2 * n_points);
1503 break;
1504
1505 default:
1506 Assert(false,
1508 "Only 1d, 2d and 3d implemented for Hessian"));
1509 }
1510 }
1511
1512
1513
1521 template <int dim, typename Number>
1522 void
1523 integrate_hessians_slow(const unsigned int n_components,
1525 Number *values_dofs,
1526 const bool add_into_values_array)
1527 {
1528 const auto &univariate_shape_data = fe_eval.get_shape_info().data;
1529 using Impl =
1531 using Eval = typename Impl::Eval;
1532 Eval eval0 =
1533 Impl::create_evaluator_tensor_product(&univariate_shape_data[0]);
1534 Eval eval1 = Impl::create_evaluator_tensor_product(
1535 &univariate_shape_data[std::min<int>(1,
1536 univariate_shape_data.size() - 1)]);
1537 Eval eval2 = Impl::create_evaluator_tensor_product(
1538 &univariate_shape_data[std::min<int>(2,
1539 univariate_shape_data.size() - 1)]);
1540
1541 const unsigned int n_points = fe_eval.get_shape_info().n_q_points;
1542 Number *tmp1 = fe_eval.get_scratch_data().begin();
1543 Number *tmp2 =
1544 tmp1 + std::max(Utilities::fixed_power<dim>(
1545 univariate_shape_data.front().fe_degree + 1),
1546 Utilities::fixed_power<dim>(
1547 univariate_shape_data.front().n_q_points_1d));
1548 const Number *hessians = fe_eval.begin_hessians();
1549
1550 for (unsigned int comp = 0; comp < n_components;
1551 ++comp,
1552 hessians += n_points * dim * (dim + 1) / 2,
1553 values_dofs +=
1555 switch (dim)
1556 {
1557 case 1:
1558 if (add_into_values_array)
1559 eval0.template hessians<0, false, true>(hessians, values_dofs);
1560 else
1561 eval0.template hessians<0, false, false>(hessians, values_dofs);
1562 break;
1563 case 2:
1564 // xx derivative
1565 eval1.template values<1, false, false>(hessians, tmp1);
1566 if (add_into_values_array)
1567 eval0.template hessians<0, false, true>(tmp1, values_dofs);
1568 else
1569 eval0.template hessians<0, false, false>(tmp1, values_dofs);
1570
1571 // xy derivative
1572 eval1.template gradients<1, false, false>(hessians + 2 * n_points,
1573 tmp1);
1574 eval0.template gradients<0, false, true>(tmp1, values_dofs);
1575 // yy derivative
1576 eval1.template hessians<1, false, false>(hessians + n_points, tmp1);
1577 eval0.template values<0, false, true>(tmp1, values_dofs);
1578 break;
1579 case 3:
1580 // xx derivative
1581 eval2.template values<2, false, false>(hessians, tmp1);
1582 eval1.template values<1, false, false>(tmp1, tmp2);
1583
1584 if (add_into_values_array)
1585 eval0.template hessians<0, false, true>(tmp2, values_dofs);
1586 else
1587 eval0.template hessians<0, false, false>(tmp2, values_dofs);
1588
1589 // xy derivative
1590 eval2.template values<2, false, false>(hessians + 3 * n_points,
1591 tmp1);
1592 eval1.template gradients<1, false, false>(tmp1, tmp2);
1593 // xz derivative
1594 eval2.template gradients<2, false, false>(hessians + 4 * n_points,
1595 tmp1);
1596 eval1.template values<1, false, true>(tmp1, tmp2);
1597 eval1.template values<0, false, true>(tmp2, values_dofs);
1598
1599 // yy derivative
1600 eval2.template values<2, false, false>(hessians + n_points, tmp1);
1601 eval1.template hessians<1, false, false>(tmp1, tmp2);
1602
1603 // yz derivative
1604 eval2.template gradients<2, false, false>(hessians + 5 * n_points,
1605 tmp1);
1606 eval1.template gradients<1, false, true>(tmp1, tmp2);
1607
1608 // zz derivative
1609 eval2.template hessians<2, false, false>(hessians + 2 * n_points,
1610 tmp1);
1611 eval1.template values<1, false, true>(tmp1, tmp2);
1612 eval0.template values<0, false, true>(tmp2, values_dofs);
1613 break;
1614
1615 default:
1616 Assert(false,
1618 "Only 1d, 2d and 3d implemented for Hessian"));
1619 }
1620 }
1621
1622
1623
1636 template <int dim, int fe_degree, typename Number>
1638 {
1639 using Number2 =
1642 dim,
1643 fe_degree + 1,
1644 fe_degree + 1,
1645 Number,
1646 Number2>;
1647
1648 static void
1649 evaluate(const unsigned int n_components,
1650 const EvaluationFlags::EvaluationFlags evaluation_flag,
1651 const Number *values_dofs,
1653 {
1654 constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1655
1656 for (unsigned int c = 0; c < n_components; ++c)
1657 {
1658 if ((evaluation_flag & EvaluationFlags::values) != 0u)
1659 for (unsigned int i = 0; i < n_points; ++i)
1660 fe_eval.begin_values()[n_points * c + i] =
1661 values_dofs[n_points * c + i];
1662
1663 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
1664 evaluate_gradients_collocation<fe_degree + 1, dim>(
1665 fe_eval.get_shape_info().data.front(),
1666 values_dofs + c * n_points,
1667 fe_eval.begin_gradients() + c * dim * n_points);
1668 }
1669 }
1670
1671 static void
1672 integrate(const unsigned int n_components,
1673 const EvaluationFlags::EvaluationFlags integration_flag,
1674 Number *values_dofs,
1676 const bool add_into_values_array)
1677 {
1678 constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
1679
1680 for (unsigned int c = 0; c < n_components; ++c)
1681 {
1682 if ((integration_flag & EvaluationFlags::values) != 0u)
1683 {
1684 if (add_into_values_array)
1685 for (unsigned int i = 0; i < n_points; ++i)
1686 values_dofs[n_points * c + i] +=
1687 fe_eval.begin_values()[n_points * c + i];
1688 else
1689 for (unsigned int i = 0; i < n_points; ++i)
1690 values_dofs[n_points * c + i] =
1691 fe_eval.begin_values()[n_points * c + i];
1692 }
1693
1694 if ((integration_flag & EvaluationFlags::gradients) != 0u)
1695 integrate_gradients_collocation<fe_degree + 1, dim>(
1696 fe_eval.get_shape_info().data.front(),
1697 values_dofs + c * n_points,
1698 fe_eval.begin_gradients() + c * dim * n_points,
1699 add_into_values_array ||
1700 ((integration_flag & EvaluationFlags::values) != 0u));
1701 }
1702 }
1703 };
1704
1705
1706
1717 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1719 {
1720 static void
1721 evaluate(const unsigned int n_components,
1722 const EvaluationFlags::EvaluationFlags evaluation_flag,
1723 const Number *values_dofs,
1725 {
1726 const auto &shape_data = fe_eval.get_shape_info().data.front();
1727
1728 Assert(n_q_points_1d > fe_degree,
1729 ExcMessage("You lose information when going to a collocation "
1730 "space of lower degree, so the evaluation results "
1731 "would be wrong. Thus, this class does not permit "
1732 "the chosen operation."));
1733 constexpr std::size_t n_dofs = Utilities::pow(fe_degree + 1, dim);
1734 constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
1735
1736 for (unsigned int c = 0; c < n_components; ++c)
1737 {
1741 dim,
1742 (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1743 n_q_points_1d>::do_forward(1,
1744 shape_data.shape_values_eo,
1745 values_dofs + c * n_dofs,
1746 fe_eval.begin_values() + c * n_q_points);
1747
1748 // apply derivatives in the collocation space
1749 if (evaluation_flag & EvaluationFlags::gradients)
1750 evaluate_gradients_collocation<n_q_points_1d, dim>(
1751 shape_data,
1752 fe_eval.begin_values() + c * n_q_points,
1753 fe_eval.begin_gradients() + c * dim * n_q_points);
1754 }
1755 }
1756
1757 static void
1758 integrate(const unsigned int n_components,
1759 const EvaluationFlags::EvaluationFlags integration_flag,
1760 Number *values_dofs,
1762 const bool add_into_values_array)
1763 {
1764 const auto &shape_data = fe_eval.get_shape_info().data.front();
1765
1766 Assert(n_q_points_1d > fe_degree,
1767 ExcMessage("You lose information when going to a collocation "
1768 "space of lower degree, so the evaluation results "
1769 "would be wrong. Thus, this class does not permit "
1770 "the chosen operation."));
1771 constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
1772
1773 for (unsigned int c = 0; c < n_components; ++c)
1774 {
1775 // apply derivatives in collocation space
1776 if (integration_flag & EvaluationFlags::gradients)
1777 integrate_gradients_collocation<n_q_points_1d, dim>(
1778 shape_data,
1779 fe_eval.begin_values() + c * n_q_points,
1780 fe_eval.begin_gradients() + c * dim * n_q_points,
1781 /*add_into_values_array=*/
1782 integration_flag & EvaluationFlags::values);
1783
1784 // transform back to the original space
1788 dim,
1789 (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1790 n_q_points_1d>::do_backward(1,
1791 shape_data.shape_values_eo,
1792 add_into_values_array,
1793 fe_eval.begin_values() + c * n_q_points,
1794 values_dofs +
1795 c *
1796 Utilities::pow(fe_degree + 1, dim));
1797 }
1798 }
1799 };
1800
1801
1802
1807 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1808 struct FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
1809 dim,
1810 fe_degree,
1811 n_q_points_1d,
1812 Number>
1813 {
1814 using Number2 =
1816
1817 template <bool integrate>
1818 static void
1819 evaluate_or_integrate(
1820 const EvaluationFlags::EvaluationFlags evaluation_flag,
1821 Number *values_dofs_actual,
1823 const bool add_into_values_array = false);
1824 };
1825
1826
1827
1828 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
1829 template <bool integrate>
1830 inline void
1832 dim,
1833 fe_degree,
1834 n_q_points_1d,
1835 Number>::
1836 evaluate_or_integrate(
1837 const EvaluationFlags::EvaluationFlags evaluation_flag,
1838 Number *values_dofs,
1840 const bool add)
1841 {
1842 Assert(dim == 2 || dim == 3,
1843 ExcMessage("Only dim = 2,3 implemented for Raviart-Thomas "
1844 "evaluation/integration"));
1845
1846 if (evaluation_flag == EvaluationFlags::nothing)
1847 return;
1848
1849 AssertDimension(fe_eval.get_shape_info().data.size(), 2);
1850 AssertDimension(n_q_points_1d,
1851 fe_eval.get_shape_info().data[0].n_q_points_1d);
1852 AssertDimension(n_q_points_1d,
1853 fe_eval.get_shape_info().data[1].n_q_points_1d);
1854 AssertDimension(fe_degree, fe_eval.get_shape_info().data[0].fe_degree);
1855 AssertDimension(fe_degree, fe_eval.get_shape_info().data[1].fe_degree + 1);
1856
1857 const auto &shape_data = fe_eval.get_shape_info().data;
1858 const unsigned int dofs_per_component =
1859 Utilities::pow(fe_degree, dim - 1) * (fe_degree + 1);
1860 const unsigned int n_points = Utilities::pow(n_q_points_1d, dim);
1861 Number *gradients = fe_eval.begin_gradients();
1862 Number *values = fe_eval.begin_values();
1863
1864 if (integrate)
1865 {
1867 eval;
1868
1869 const bool do_values = evaluation_flag & EvaluationFlags::values;
1870 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
1871 integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
1872 values,
1873 gradients,
1874 do_values);
1875 if constexpr (dim > 2)
1876 eval.template tangential<2, 0>(shape_data[1], values, values);
1877 eval.template tangential<1, 0>(shape_data[1], values, values);
1878 eval.template normal<0>(shape_data[0], values, values_dofs, add);
1879
1880 values += n_points;
1881 gradients += n_points * dim;
1882 values_dofs += dofs_per_component;
1883
1884 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
1885 integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
1886 values,
1887 gradients,
1888 do_values);
1889 if constexpr (dim > 2)
1890 eval.template tangential<2, 1>(shape_data[1], values, values);
1891 eval.template tangential<0, 1>(shape_data[1], values, values);
1892 eval.template normal<1>(shape_data[0], values, values_dofs, add);
1893
1894 if constexpr (dim > 2)
1895 {
1896 values += n_points;
1897 gradients += n_points * dim;
1898 values_dofs += dofs_per_component;
1899
1900 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
1901 integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
1902 values,
1903 gradients,
1904 do_values);
1905 eval.template tangential<1, 2>(shape_data[1], values, values);
1906 eval.template tangential<0, 2>(shape_data[1], values, values);
1907 eval.template normal<0>(shape_data[0], values, values_dofs, add);
1908 }
1909 }
1910 else
1911 {
1913 eval;
1914 eval.template normal<0>(shape_data[0], values_dofs, values);
1915 eval.template tangential<1, 0>(shape_data[1], values, values);
1916 if constexpr (dim > 2)
1917 eval.template tangential<2, 0>(shape_data[1], values, values);
1918 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
1919 evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
1920 values,
1921 gradients);
1922
1923 values += n_points;
1924 gradients += n_points * dim;
1925 values_dofs += dofs_per_component;
1926
1927 eval.template normal<1>(shape_data[0], values_dofs, values);
1928 eval.template tangential<0, 1>(shape_data[1], values, values);
1929 if constexpr (dim > 2)
1930 eval.template tangential<2, 1>(shape_data[1], values, values);
1931 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
1932 evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
1933 values,
1934 gradients);
1935
1936 if constexpr (dim > 2)
1937 {
1938 values += n_points;
1939 gradients += n_points * dim;
1940 values_dofs += dofs_per_component;
1941
1942 eval.template normal<2>(shape_data[0], values_dofs, values);
1943 eval.template tangential<0, 2>(shape_data[1], values, values);
1944 eval.template tangential<1, 2>(shape_data[1], values, values);
1945 if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
1946 evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
1947 values,
1948 gradients);
1949 }
1950 }
1951 }
1952
1953
1954
1970 template <int dim, typename Number, bool do_integrate>
1972 {
1973 template <int fe_degree, int n_q_points_1d, typename OtherNumber>
1974 static bool
1975 run(const unsigned int n_components,
1976 const EvaluationFlags::EvaluationFlags evaluation_flag,
1977 OtherNumber *values_dofs,
1979 const bool sum_into_values_array_in = false)
1980 {
1981 // `OtherNumber` is either `const Number` (evaluate()) or `Number`
1982 // (integrate())
1983 static_assert(std::is_same_v<Number, std::remove_const_t<OtherNumber>>,
1984 "Type of Number and of OtherNumber do not match.");
1985
1986 const auto element_type = fe_eval.get_shape_info().element_type;
1987 using ElementType = MatrixFreeFunctions::ElementType;
1988
1989 Assert(fe_eval.get_shape_info().data.size() == 1 ||
1990 (fe_eval.get_shape_info().data.size() == dim &&
1991 element_type == ElementType::tensor_general) ||
1992 element_type == ElementType::tensor_raviart_thomas,
1994
1995 EvaluationFlags::EvaluationFlags actual_flag = evaluation_flag;
1996 bool sum_into_values_array = sum_into_values_array_in;
1997 if (evaluation_flag & EvaluationFlags::hessians)
1998 {
1999 actual_flag |= EvaluationFlags::values;
2002 if constexpr (do_integrate)
2003 {
2004 if (fe_eval.get_shape_info().data[0].fe_degree <
2005 fe_eval.get_shape_info().data[0].n_q_points_1d)
2006 integrate_hessians_collocation<n_q_points_1d>(
2007 n_components,
2008 fe_eval,
2009 evaluation_flag & EvaluationFlags::values);
2010 else
2011 {
2012 integrate_hessians_slow(n_components,
2013 fe_eval,
2014 values_dofs,
2015 sum_into_values_array);
2016 sum_into_values_array = true;
2017 }
2018 }
2019 }
2020
2021 if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
2022 element_type == ElementType::tensor_symmetric_collocation)
2023 {
2026 n_components,
2027 actual_flag,
2028 values_dofs,
2029 fe_eval,
2030 sum_into_values_array);
2031 }
2032 // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
2033 // shape_info.h for more details
2034 else if (fe_degree >= 0 &&
2035 use_collocation_evaluation(fe_degree, n_q_points_1d) &&
2036 element_type <= ElementType::tensor_symmetric)
2037 {
2040 fe_degree,
2041 n_q_points_1d,
2042 Number>>(
2043 n_components,
2044 actual_flag,
2045 values_dofs,
2046 fe_eval,
2047 sum_into_values_array);
2048 }
2049 else if (fe_degree >= 0 &&
2050 element_type <= ElementType::tensor_symmetric_no_collocation)
2051 {
2052 evaluate_or_integrate<FEEvaluationImpl<ElementType::tensor_symmetric,
2053 dim,
2054 fe_degree,
2055 n_q_points_1d,
2056 Number>>(
2057 n_components,
2058 actual_flag,
2059 values_dofs,
2060 fe_eval,
2061 sum_into_values_array);
2062 }
2063 else if (element_type == ElementType::tensor_symmetric_plus_dg0)
2064 {
2066 FEEvaluationImpl<ElementType::tensor_symmetric_plus_dg0,
2067 dim,
2068 fe_degree,
2069 n_q_points_1d,
2070 Number>>(n_components,
2071 actual_flag,
2072 values_dofs,
2073 fe_eval,
2074 sum_into_values_array);
2075 }
2076 else if (element_type == ElementType::truncated_tensor)
2077 {
2078 evaluate_or_integrate<FEEvaluationImpl<ElementType::truncated_tensor,
2079 dim,
2080 fe_degree,
2081 n_q_points_1d,
2082 Number>>(
2083 n_components,
2084 actual_flag,
2085 values_dofs,
2086 fe_eval,
2087 sum_into_values_array);
2088 }
2089 else if (element_type == ElementType::tensor_none)
2090 {
2091 evaluate_or_integrate<FEEvaluationImpl<ElementType::tensor_none,
2092 dim,
2093 fe_degree,
2094 n_q_points_1d,
2095 Number>>(
2096 n_components,
2097 actual_flag,
2098 values_dofs,
2099 fe_eval,
2100 sum_into_values_array);
2101 }
2102 else if (element_type == ElementType::tensor_raviart_thomas)
2103 {
2104 if constexpr (fe_degree > 0 && n_q_points_1d > 0 && dim > 1)
2105 {
2106 FEEvaluationImpl<ElementType::tensor_raviart_thomas,
2107 dim,
2108 fe_degree,
2109 n_q_points_1d,
2110 Number>::
2111 template evaluate_or_integrate<do_integrate>(
2112 actual_flag,
2113 const_cast<Number *>(values_dofs),
2114 fe_eval,
2115 sum_into_values_array);
2116 }
2117 else
2118 {
2119 Assert(false,
2120 ExcNotImplemented("Raviart-Thomas currently only possible "
2121 "in 2d/3d and with templated degree"));
2122 }
2123 }
2124 else
2125 {
2126 evaluate_or_integrate<FEEvaluationImpl<ElementType::tensor_general,
2127 dim,
2128 fe_degree,
2129 n_q_points_1d,
2130 Number>>(
2131 n_components,
2132 actual_flag,
2133 values_dofs,
2134 fe_eval,
2135 sum_into_values_array);
2136 }
2137
2138 if ((evaluation_flag & EvaluationFlags::hessians) && !do_integrate)
2139 {
2142 if (fe_eval.get_shape_info().data[0].fe_degree <
2143 fe_eval.get_shape_info().data[0].n_q_points_1d)
2144 evaluate_hessians_collocation<n_q_points_1d>(n_components, fe_eval);
2145 else
2146 evaluate_hessians_slow(n_components, values_dofs, fe_eval);
2147 }
2148
2149 return false;
2150 }
2151
2152 private:
2153 template <typename T>
2154 static void
2156 const unsigned int n_components,
2157 const EvaluationFlags::EvaluationFlags evaluation_flag,
2158 const Number *values_dofs,
2160 const bool sum_into_values_array,
2161 std::bool_constant<false>)
2162 {
2163 (void)sum_into_values_array;
2164
2165 T::evaluate(n_components, evaluation_flag, values_dofs, fe_eval);
2166 }
2167
2168 template <typename T>
2169 static void
2171 const unsigned int n_components,
2172 const EvaluationFlags::EvaluationFlags evaluation_flag,
2173 Number *values_dofs,
2175 const bool sum_into_values_array,
2176 std::bool_constant<true>)
2177 {
2178 T::integrate(n_components,
2179 evaluation_flag,
2180 values_dofs,
2181 fe_eval,
2182 sum_into_values_array);
2183 }
2184
2185 template <typename T, typename OtherNumber>
2186 static void
2188 const unsigned int n_components,
2189 const EvaluationFlags::EvaluationFlags evaluation_flag,
2190 OtherNumber *values_dofs,
2192 const bool sum_into_values_array)
2193 {
2194 evaluate_or_integrate<T>(n_components,
2195 evaluation_flag,
2196 values_dofs,
2197 fe_eval,
2198 sum_into_values_array,
2199 std::bool_constant<do_integrate>());
2200 }
2201 };
2202
2203
2204
2209 template <int dim, typename Number>
2211 {
2212 using Number2 =
2214
2215 template <int fe_degree, int = 0>
2216 static bool
2217 run(const unsigned int n_components,
2219 const Number *in_array,
2220 Number *out_array)
2221 {
2222 const unsigned int given_degree =
2223 (fe_degree > -1) ? fe_degree :
2224 fe_eval.get_shape_info().data.front().fe_degree;
2225
2226 const unsigned int dofs_per_component =
2227 Utilities::pow(given_degree + 1, dim);
2228
2229 Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
2233
2235 dim,
2236 fe_degree + 1,
2237 fe_degree + 1,
2238 Number,
2239 Number2>
2240 evaluator({},
2241 {},
2242 fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
2243 given_degree + 1,
2244 given_degree + 1);
2245
2246 for (unsigned int d = 0; d < n_components; ++d)
2247 {
2248 const Number *in = in_array + d * dofs_per_component;
2249 Number *out = out_array + d * dofs_per_component;
2250 // Need to select 'apply' method with hessian slot because values
2251 // assume symmetries that do not exist in the inverse shapes
2252 evaluator.template hessians<0, true, false>(in, out);
2253 if (dim > 1)
2254 evaluator.template hessians<1, true, false>(out, out);
2255 if (dim > 2)
2256 evaluator.template hessians<2, true, false>(out, out);
2257 }
2258 for (unsigned int q = 0; q < dofs_per_component; ++q)
2259 {
2260 const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
2261 for (unsigned int d = 0; d < n_components; ++d)
2262 out_array[q + d * dofs_per_component] *= inverse_JxW_q;
2263 }
2264 for (unsigned int d = 0; d < n_components; ++d)
2265 {
2266 Number *out = out_array + d * dofs_per_component;
2267 if (dim > 2)
2268 evaluator.template hessians<2, false, false>(out, out);
2269 if (dim > 1)
2270 evaluator.template hessians<1, false, false>(out, out);
2271 evaluator.template hessians<0, false, false>(out, out);
2272 }
2273 return false;
2274 }
2275 };
2276
2277
2278
2285 template <int dim, typename Number>
2287 {
2288 using Number2 =
2290
2291 template <int fe_degree, int = 0>
2292 static bool
2293 run(const unsigned int n_desired_components,
2295 const ArrayView<const Number> &inverse_coefficients,
2296 const bool dyadic_coefficients,
2297 const Number *in_array,
2298 Number *out_array)
2299 {
2300 const unsigned int given_degree =
2301 (fe_degree > -1) ? fe_degree :
2302 fe_eval.get_shape_info().data.front().fe_degree;
2303
2304 const unsigned int dofs_per_component =
2305 Utilities::pow(given_degree + 1, dim);
2306
2307 Assert(inverse_coefficients.size() > 0 &&
2308 inverse_coefficients.size() % dofs_per_component == 0,
2309 ExcMessage(
2310 "Expected diagonal to be a multiple of scalar dof per cells"));
2311
2312 if (!dyadic_coefficients)
2313 {
2314 if (inverse_coefficients.size() != dofs_per_component)
2315 AssertDimension(n_desired_components * dofs_per_component,
2316 inverse_coefficients.size());
2317 }
2318 else
2319 {
2320 AssertDimension(n_desired_components * n_desired_components *
2321 dofs_per_component,
2322 inverse_coefficients.size());
2323 }
2324
2325 Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
2329
2331 dim,
2332 fe_degree + 1,
2333 fe_degree + 1,
2334 Number,
2335 Number2>
2336 evaluator({},
2337 {},
2338 fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
2339 given_degree + 1,
2340 given_degree + 1);
2341
2342 const Number *in = in_array;
2343 Number *out = out_array;
2344
2345 const Number *inv_coefficient = inverse_coefficients.data();
2346
2347 const unsigned int shift_coefficient =
2348 inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
2349 0;
2350
2351 const auto n_comp_outer = dyadic_coefficients ? 1 : n_desired_components;
2352 const auto n_comp_inner = dyadic_coefficients ? n_desired_components : 1;
2353
2354 for (unsigned int d = 0; d < n_comp_outer; ++d)
2355 {
2356 for (unsigned int di = 0; di < n_comp_inner; ++di)
2357 {
2358 const Number *in_ = in + di * dofs_per_component;
2359 Number *out_ = out + di * dofs_per_component;
2360 evaluator.template hessians<0, true, false>(in_, out_);
2361 if (dim > 1)
2362 evaluator.template hessians<1, true, false>(out_, out_);
2363 if (dim > 2)
2364 evaluator.template hessians<2, true, false>(out_, out_);
2365 }
2366 if (dyadic_coefficients)
2367 {
2368 const auto n_coeff_components =
2369 n_desired_components * n_desired_components;
2370 if (n_desired_components == dim)
2371 {
2372 for (unsigned int q = 0; q < dofs_per_component; ++q)
2373 vmult<dim>(&inv_coefficient[q * n_coeff_components],
2374 &in[q],
2375 &out[q],
2376 dofs_per_component);
2377 }
2378 else
2379 {
2380 for (unsigned int q = 0; q < dofs_per_component; ++q)
2381 vmult<-1>(&inv_coefficient[q * n_coeff_components],
2382 &in[q],
2383 &out[q],
2384 dofs_per_component,
2385 n_desired_components);
2386 }
2387 }
2388 else
2389 for (unsigned int q = 0; q < dofs_per_component; ++q)
2390 out[q] *= inv_coefficient[q];
2391
2392 for (unsigned int di = 0; di < n_comp_inner; ++di)
2393 {
2394 Number *out_ = out + di * dofs_per_component;
2395 if (dim > 2)
2396 evaluator.template hessians<2, false, false>(out_, out_);
2397 if (dim > 1)
2398 evaluator.template hessians<1, false, false>(out_, out_);
2399 evaluator.template hessians<0, false, false>(out_, out_);
2400 }
2401
2402 in += dofs_per_component;
2403 out += dofs_per_component;
2404 inv_coefficient += shift_coefficient;
2405 }
2406
2407 return false;
2408 }
2409
2410 private:
2411 template <int n_components>
2412 static inline void
2413 vmult(const Number *inverse_coefficients,
2414 const Number *src,
2415 Number *dst,
2416 const unsigned int dofs_per_component,
2417 const unsigned int n_given_components = 0)
2418 {
2419 const unsigned int n_desired_components =
2420 (n_components > -1) ? n_components : n_given_components;
2421
2422 std::array<Number, dim + 2> tmp = {};
2423 Assert(n_desired_components <= dim + 2,
2424 ExcMessage(
2425 "Number of components larger than dim+2 not supported."));
2426
2427 for (unsigned int d = 0; d < n_desired_components; ++d)
2428 tmp[d] = src[d * dofs_per_component];
2429
2430 for (unsigned int d1 = 0; d1 < n_desired_components; ++d1)
2431 {
2432 const Number *inv_coeff_row =
2433 &inverse_coefficients[d1 * n_desired_components];
2434 Number sum = inv_coeff_row[0] * tmp[0];
2435 for (unsigned int d2 = 1; d2 < n_desired_components; ++d2)
2436 sum += inv_coeff_row[d2] * tmp[d2];
2437 dst[d1 * dofs_per_component] = sum;
2438 }
2439 }
2440 };
2441
2442
2443
2450 template <int dim, typename Number>
2452 {
2453 template <int fe_degree, int n_q_points_1d>
2454 static bool
2455 run(const unsigned int n_desired_components,
2457 const Number *in_array,
2458 Number *out_array)
2459 {
2460 static const bool do_inplace =
2461 fe_degree > -1 && (fe_degree + 1 == n_q_points_1d);
2462
2466
2467 const auto &inverse_shape =
2468 do_inplace ?
2469 fe_eval.get_shape_info().data.front().inverse_shape_values_eo :
2470 fe_eval.get_shape_info().data.front().inverse_shape_values;
2471
2472 const std::size_t dofs_per_component =
2473 do_inplace ? Utilities::pow(fe_degree + 1, dim) :
2475 const std::size_t n_q_points = do_inplace ?
2476 Utilities::pow(fe_degree + 1, dim) :
2477 fe_eval.get_shape_info().n_q_points;
2478
2479 using Number2 =
2482 dim,
2483 fe_degree + 1,
2484 n_q_points_1d,
2485 Number,
2486 Number2>
2487 evaluator({},
2488 {},
2489 inverse_shape,
2490 fe_eval.get_shape_info().data.front().fe_degree + 1,
2491 fe_eval.get_shape_info().data.front().n_q_points_1d);
2492
2493 for (unsigned int d = 0; d < n_desired_components; ++d)
2494 {
2495 const Number *in = in_array + d * n_q_points;
2496 Number *out = out_array + d * dofs_per_component;
2497
2498 auto *temp_1 = do_inplace ? out : fe_eval.get_scratch_data().begin();
2499 auto *temp_2 = do_inplace ?
2500 out :
2501 (temp_1 + std::max(n_q_points, dofs_per_component));
2502
2503 if (dim == 3)
2504 {
2505 evaluator.template hessians<2, false, false>(in, temp_1);
2506 evaluator.template hessians<1, false, false>(temp_1, temp_2);
2507 evaluator.template hessians<0, false, false>(temp_2, out);
2508 }
2509 if (dim == 2)
2510 {
2511 evaluator.template hessians<1, false, false>(in, temp_1);
2512 evaluator.template hessians<0, false, false>(temp_1, out);
2513 }
2514 if (dim == 1)
2515 evaluator.template hessians<0, false, false>(in, out);
2516 }
2517 return false;
2518 }
2519 };
2520} // end of namespace internal
2521
2522
2524
2525#endif
pointer data()
size_type size() const
iterator begin() const
Definition array_view.h:702
value_type * data() const noexcept
Definition array_view.h:661
std::size_t size() const
Definition array_view.h:684
const ShapeInfoType & get_shape_info() const
Number JxW(const unsigned int q_point) const
const Number * begin_gradients() const
ArrayView< Number > get_scratch_data() const
const Number * begin_values() const
const Number * begin_hessians() const
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
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:966
void evaluate_hessians_collocation(const unsigned int n_components, FEEvaluationData< dim, Number, false > &fe_eval)
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
void integrate_gradients_collocation(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &shape, Number *values, const Number *gradients, const bool add_into_values_array)
void evaluate_hessians_slow(const unsigned int n_components, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
void integrate_hessians_collocation(const unsigned int n_components, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
void evaluate_gradients_collocation(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &shape, const Number *values, Number *gradients)
void integrate_hessians_slow(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, Number *values_dofs, const bool add_into_values_array)
static const unsigned int invalid_unsigned_int
Definition types.h:220
::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_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 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_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 evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
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)
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 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::bool_constant< true >)
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_in=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::bool_constant< false >)
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 evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
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)
std::vector< UnivariateShapeData< Number > > data
Definition shape_info.h:485