Reference documentation for deal.II version 9.6.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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_matrix_free_evaluation_kernels_h
17#define dealii_matrix_free_evaluation_kernels_h
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/base/utilities.h>
22#include <deal.II/base/vectorization.h>
23
24#include <deal.II/matrix_free/evaluation_flags.h>
25#include <deal.II/matrix_free/fe_evaluation_data.h>
26#include <deal.II/matrix_free/shape_info.h>
27#include <deal.II/matrix_free/tensor_product_kernels.h>
28
29
30DEAL_II_NAMESPACE_OPEN
31
32
33namespace internal
34{
35 // Select evaluator type from element shape function type
36 template <MatrixFreeFunctions::ElementType element, bool is_long>
37 struct EvaluatorSelector
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
107 template <MatrixFreeFunctions::ElementType type,
108 int dim,
109 int fe_degree,
110 int n_q_points_1d,
111 typename Number>
112 struct FEEvaluationImpl
113 {
114 static const EvaluatorVariant variant =
115 EvaluatorSelector<type, (fe_degree + n_q_points_1d > 4)>::variant;
116 using Number2 =
117 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
118
119 using Eval = EvaluatorTensorProduct<variant,
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,
130 FEEvaluationData<dim, Number, false> &fe_eval);
131
132 static void
133 integrate(const unsigned int n_components,
134 const EvaluationFlags::EvaluationFlags integration_flag,
135 Number *values_dofs_actual,
136 FEEvaluationData<dim, Number, false> &fe_eval,
137 const bool add_into_values_array);
138
139 static Eval
140 create_evaluator_tensor_product(
141 const MatrixFreeFunctions::UnivariateShapeData<Number2>
142 *univariate_shape_data)
143 {
144 if (variant == evaluate_evenodd)
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,
176 FEEvaluationData<dim, Number, false> &fe_eval);
177
178 static void
179 integrate(const unsigned int n_components,
180 const EvaluationFlags::EvaluationFlags integration_flag,
181 Number *values_dofs_actual,
182 FEEvaluationData<dim, Number, false> &fe_eval,
183 const bool add_into_values_array);
184 };
185
186
187
188 template <MatrixFreeFunctions::ElementType type,
189 int dim,
190 int fe_degree,
191 int n_q_points_1d,
192 typename Number>
193 inline void
194 FEEvaluationImpl<type, dim, fe_degree, n_q_points_1d, Number>::evaluate(
195 const unsigned int n_components,
196 const EvaluationFlags::EvaluationFlags evaluation_flag,
197 const Number *values_dofs_actual,
198 FEEvaluationData<dim, Number, false> &fe_eval)
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 =
242 (type == MatrixFreeFunctions::truncated_tensor) ?
243 Utilities::pow(shape_data.front().fe_degree + 1, dim) :
244 fe_eval.get_shape_info().dofs_per_component_on_cell;
245 const Number *values_dofs = values_dofs_actual;
246 if (type == MatrixFreeFunctions::truncated_tensor)
247 {
248 const std::size_t n_dofs_per_comp =
249 fe_eval.get_shape_info().dofs_per_component_on_cell;
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:
366 AssertThrow(false, ExcNotImplemented());
367 }
368
369 // case additional dof for FE_Q_DG0: add values; gradients and second
370 // derivatives evaluate to zero
371 if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 &&
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
385 template <MatrixFreeFunctions::ElementType type,
386 int dim,
387 int fe_degree,
388 int n_q_points_1d,
389 typename Number>
390 inline void
391 FEEvaluationImpl<type, dim, fe_degree, n_q_points_1d, Number>::integrate(
392 const unsigned int n_components,
393 const EvaluationFlags::EvaluationFlags integration_flag,
394 Number *values_dofs_actual,
395 FEEvaluationData<dim, Number, false> &fe_eval,
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 =
436 (type == MatrixFreeFunctions::truncated_tensor) ?
437 Utilities::fixed_power<dim>(shape_data.front().fe_degree + 1) :
438 fe_eval.get_shape_info().dofs_per_component_on_cell;
439 // expand dof_values to tensor product for truncated tensor products
440 Number *values_dofs =
441 (type == MatrixFreeFunctions::truncated_tensor) ?
442 temp1 + 2 * (std::max<std::size_t>(
443 fe_eval.get_shape_info().dofs_per_component_on_cell,
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:
559 AssertThrow(false, ExcNotImplemented());
560 }
561
562 // case FE_Q_DG0: add values, gradients and second derivatives are zero
563 if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
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
584 if (type == MatrixFreeFunctions::truncated_tensor)
585 {
586 const std::size_t n_dofs_per_comp =
587 fe_eval.get_shape_info().dofs_per_component_on_cell;
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
613 FEEvaluationImpl<
614 MatrixFreeFunctions::tensor_none,
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,
621 FEEvaluationData<dim, Number, false> &fe_eval)
622 {
623 Assert(!(evaluation_flag & EvaluationFlags::hessians), ExcNotImplemented());
624
625 const std::size_t n_dofs =
626 fe_eval.get_shape_info().dofs_per_component_on_cell;
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 =
632 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
633
634 if (evaluation_flag & EvaluationFlags::values)
635 {
636 const auto *const shape_values = shape_data.front().shape_values.data();
637 auto *out = fe_eval.begin_values();
638 const auto *in = values_dofs_actual;
639
640 for (unsigned int c = 0; c < n_components; ++c)
641 {
642 apply_matrix_vector_product<evaluate_general,
643 EvaluatorQuantity::value,
644 /* transpose_matrix */ true,
645 /* add */ false,
646 /* consider_strides */ false>(
647 shape_values, in, out, n_dofs, n_q_points, 1, 1);
648
649 out += n_q_points;
650 in += n_dofs;
651 }
652 }
653
654 if (evaluation_flag & EvaluationFlags::gradients)
655 {
656 const auto *const shape_gradients =
657 shape_data.front().shape_gradients.data();
658 auto *out = fe_eval.begin_gradients();
659 const auto *in = values_dofs_actual;
660
661 for (unsigned int c = 0; c < n_components; ++c)
662 {
663 apply_matrix_vector_product<evaluate_general,
664 EvaluatorQuantity::value,
665 /* transpose_matrix */ true,
666 /* add */ false,
667 /* consider_strides */ false>(
668 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
669 out += n_q_points * dim;
670 in += n_dofs;
671 }
672 }
673 }
674
675
676
677 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
678 inline void
679 FEEvaluationImpl<
680 MatrixFreeFunctions::tensor_none,
681 dim,
682 fe_degree,
683 n_q_points_1d,
684 Number>::integrate(const unsigned int n_components,
685 const EvaluationFlags::EvaluationFlags integration_flag,
686 Number *values_dofs_actual,
687 FEEvaluationData<dim, Number, false> &fe_eval,
688 const bool add_into_values_array)
689 {
690 Assert(!(integration_flag & EvaluationFlags::hessians),
691 ExcNotImplemented());
692
693 const std::size_t n_dofs =
694 fe_eval.get_shape_info().dofs_per_component_on_cell;
695 const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
696
697 const auto &shape_data = fe_eval.get_shape_info().data;
698
699 using Number2 =
700 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
701
702 if (integration_flag & EvaluationFlags::values)
703 {
704 const auto *const shape_values = shape_data.front().shape_values.data();
705 auto *in = fe_eval.begin_values();
706 auto *out = values_dofs_actual;
707
708 for (unsigned int c = 0; c < n_components; ++c)
709 {
710 if (add_into_values_array == false)
711 apply_matrix_vector_product<evaluate_general,
712 EvaluatorQuantity::value,
713 /* transpose_matrix */ false,
714 /* add */ false,
715 /* consider_strides */ false>(
716 shape_values, in, out, n_dofs, n_q_points, 1, 1);
717 else
718 apply_matrix_vector_product<evaluate_general,
719 EvaluatorQuantity::value,
720 /* transpose_matrix */ false,
721 /* add */ true,
722 /* consider_strides */ false>(
723 shape_values, in, out, n_dofs, n_q_points, 1, 1);
724
725 in += n_q_points;
726 out += n_dofs;
727 }
728 }
729
730 if (integration_flag & EvaluationFlags::gradients)
731 {
732 const auto *const shape_gradients =
733 shape_data.front().shape_gradients.data();
734 auto *in = fe_eval.begin_gradients();
735 auto *out = values_dofs_actual;
736
737 for (unsigned int c = 0; c < n_components; ++c)
738 {
739 if (add_into_values_array == false &&
740 !(integration_flag & EvaluationFlags::values))
741 apply_matrix_vector_product<evaluate_general,
742 EvaluatorQuantity::value,
743 /* transpose_matrix */ false,
744 /* add */ false,
745 /* consider_strides */ false>(
746 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
747 else
748 apply_matrix_vector_product<evaluate_general,
749 EvaluatorQuantity::value,
750 /* transpose_matrix */ false,
751 /* add */ true,
752 /* consider_strides */ false>(
753 shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
754
755 in += n_q_points * dim;
756 out += 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>
777 struct FEEvaluationImplBasisChange
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
806 DEAL_II_ALWAYS_INLINE
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 =
814 numbers::invalid_unsigned_int,
815 const unsigned int basis_size_2_variable =
816 numbers::invalid_unsigned_int)
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
822 Assert(quantity == EvaluatorQuantity::value, ExcInternalError());
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
833 EvaluatorTensorProduct<variant,
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)
864 FEEvaluationImplBasisChange<variant,
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
927 DEAL_II_ALWAYS_INLINE
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 =
936 numbers::invalid_unsigned_int,
937 const unsigned int basis_size_2_variable =
938 numbers::invalid_unsigned_int)
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,
944 ExcMessage(
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,
950 ExcInternalError());
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;
957 EvaluatorTensorProduct<variant,
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)
1038 FEEvaluationImplBasisChange<variant,
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)
1099 FEEvaluationImplBasisChange<
1100 variant,
1101 EvaluatorQuantity::value,
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)
1138 FEEvaluationImplBasisChange<variant,
1139 EvaluatorQuantity::value,
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
1162 evaluate_gradients_collocation(
1163 const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
1164 const Number *values,
1165 Number *gradients)
1166 {
1167 AssertDimension(shape.shape_gradients_collocation_eo.size(),
1168 (n_points_1d + 1) / 2 * n_points_1d);
1169
1170 EvaluatorTensorProduct<evaluate_evenodd,
1171 dim,
1172 n_points_1d,
1173 n_points_1d,
1174 Number,
1175 Number2>
1176 eval({}, shape.shape_gradients_collocation_eo, {});
1177 EvaluatorTensorProduct<evaluate_evenodd,
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
1216 integrate_gradients_collocation(
1217 const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
1218 Number *values,
1219 const Number *gradients,
1220 const bool add_into_values_array)
1221 {
1222 AssertDimension(shape.shape_gradients_collocation_eo.size(),
1223 (n_points_1d + 1) / 2 * n_points_1d);
1224
1225 EvaluatorTensorProduct<evaluate_evenodd,
1226 dim,
1227 n_points_1d,
1228 n_points_1d,
1229 Number,
1230 Number2>
1231 eval({}, shape.shape_gradients_collocation_eo, {});
1232 EvaluatorTensorProduct<evaluate_evenodd,
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,
1280 FEEvaluationData<dim, Number, false> &fe_eval)
1281 {
1282 using Number2 =
1283 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
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
1288 const MatrixFreeFunctions::UnivariateShapeData<Number2> &data =
1289 fe_eval.get_shape_info().data[0];
1290 AssertDimension(data.shape_gradients_collocation.size(),
1291 data.n_q_points_1d * data.n_q_points_1d);
1292 EvaluatorTensorProduct<evaluate_general,
1293 dim,
1294 n_points_1d,
1295 n_points_1d,
1296 Number,
1297 Number2>
1298 eval({},
1299 data.shape_gradients_collocation.data(),
1300 data.shape_hessians_collocation.data(),
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,
1353 FEEvaluationData<dim, Number, false> &fe_eval,
1354 const bool add_into_values_array)
1355 {
1356 using Number2 =
1357 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
1358
1359 const MatrixFreeFunctions::UnivariateShapeData<Number2> &data =
1360 fe_eval.get_shape_info().data[0];
1361 AssertDimension(data.shape_gradients_collocation.size(),
1362 data.n_q_points_1d * data.n_q_points_1d);
1363 EvaluatorTensorProduct<evaluate_general,
1364 dim,
1365 n_q_points_1d,
1366 n_q_points_1d,
1367 Number,
1368 Number2>
1369 eval({},
1370 data.shape_gradients_collocation.data(),
1371 data.shape_hessians_collocation.data(),
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,
1431 FEEvaluationData<dim, Number, false> &fe_eval)
1432 {
1433 const auto &univariate_shape_data = fe_eval.get_shape_info().data;
1434 using Impl =
1435 FEEvaluationImpl<MatrixFreeFunctions::tensor_general, dim, -1, 0, Number>;
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 +=
1459 fe_eval.get_shape_info().dofs_per_component_on_cell)
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,
1507 ExcNotImplemented(
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,
1524 const FEEvaluationData<dim, Number, false> &fe_eval,
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 =
1530 FEEvaluationImpl<MatrixFreeFunctions::tensor_general, dim, -1, 0, Number>;
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 +=
1554 fe_eval.get_shape_info().dofs_per_component_on_cell)
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,
1617 ExcNotImplemented(
1618 "Only 1d, 2d and 3d implemented for Hessian"));
1619 }
1620 }
1621
1622
1623
1636 template <int dim, int fe_degree, typename Number>
1637 struct FEEvaluationImplCollocation
1638 {
1639 using Number2 =
1640 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
1641 using Eval = EvaluatorTensorProduct<evaluate_evenodd,
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,
1652 FEEvaluationData<dim, Number, false> &fe_eval)
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,
1675 FEEvaluationData<dim, Number, false> &fe_eval,
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>
1718 struct FEEvaluationImplTransformToCollocation
1719 {
1720 static void
1721 evaluate(const unsigned int n_components,
1722 const EvaluationFlags::EvaluationFlags evaluation_flag,
1723 const Number *values_dofs,
1724 FEEvaluationData<dim, Number, false> &fe_eval)
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 {
1738 FEEvaluationImplBasisChange<
1739 evaluate_evenodd,
1740 EvaluatorQuantity::value,
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,
1761 FEEvaluationData<dim, Number, false> &fe_eval,
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
1785 FEEvaluationImplBasisChange<
1786 evaluate_evenodd,
1787 EvaluatorQuantity::value,
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 =
1815 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
1816
1817 template <bool integrate>
1818 static void
1819 evaluate_or_integrate(
1820 const EvaluationFlags::EvaluationFlags evaluation_flag,
1821 Number *values_dofs_actual,
1822 FEEvaluationData<dim, Number, false> &fe_eval,
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
1831 FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
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,
1839 FEEvaluationData<dim, Number, false> &fe_eval,
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 {
1866 EvaluatorTensorProductAnisotropic<dim, fe_degree, n_q_points_1d, false>
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 {
1912 EvaluatorTensorProductAnisotropic<dim, fe_degree, n_q_points_1d, true>
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>
1971 struct FEEvaluationImplSelector
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,
1978 FEEvaluationData<dim, Number, false> &fe_eval,
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,
1993 ExcNotImplemented());
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;
2000 Assert(element_type != MatrixFreeFunctions::tensor_none,
2001 ExcNotImplemented());
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 {
2024 evaluate_or_integrate<
2025 FEEvaluationImplCollocation<dim, fe_degree, Number>>(
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 {
2038 evaluate_or_integrate<
2039 FEEvaluationImplTransformToCollocation<dim,
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_none)
2064 {
2065 evaluate_or_integrate<
2066 FEEvaluationImpl<ElementType::tensor_none, dim, -1, 0, Number>>(
2067 n_components,
2068 actual_flag,
2069 values_dofs,
2070 fe_eval,
2071 sum_into_values_array);
2072 }
2073 else if (element_type == ElementType::tensor_symmetric_plus_dg0)
2074 {
2075 evaluate_or_integrate<
2076 FEEvaluationImpl<ElementType::tensor_symmetric_plus_dg0,
2077 dim,
2078 fe_degree,
2079 n_q_points_1d,
2080 Number>>(n_components,
2081 actual_flag,
2082 values_dofs,
2083 fe_eval,
2084 sum_into_values_array);
2085 }
2086 else if (element_type == ElementType::truncated_tensor)
2087 {
2088 evaluate_or_integrate<FEEvaluationImpl<ElementType::truncated_tensor,
2089 dim,
2090 fe_degree,
2091 n_q_points_1d,
2092 Number>>(
2093 n_components,
2094 actual_flag,
2095 values_dofs,
2096 fe_eval,
2097 sum_into_values_array);
2098 }
2099 else if (element_type == ElementType::tensor_raviart_thomas)
2100 {
2101 if constexpr (fe_degree > 0 && n_q_points_1d > 0 && dim > 1)
2102 {
2103 FEEvaluationImpl<ElementType::tensor_raviart_thomas,
2104 dim,
2105 fe_degree,
2106 n_q_points_1d,
2107 Number>::
2108 template evaluate_or_integrate<do_integrate>(
2109 actual_flag,
2110 const_cast<Number *>(values_dofs),
2111 fe_eval,
2112 sum_into_values_array);
2113 }
2114 else
2115 {
2116 Assert(false,
2117 ExcNotImplemented("Raviart-Thomas currently only possible "
2118 "in 2d/3d and with templated degree"));
2119 }
2120 }
2121 else
2122 {
2123 evaluate_or_integrate<FEEvaluationImpl<ElementType::tensor_general,
2124 dim,
2125 fe_degree,
2126 n_q_points_1d,
2127 Number>>(
2128 n_components,
2129 actual_flag,
2130 values_dofs,
2131 fe_eval,
2132 sum_into_values_array);
2133 }
2134
2135 if ((evaluation_flag & EvaluationFlags::hessians) && !do_integrate)
2136 {
2137 Assert(element_type != MatrixFreeFunctions::tensor_none,
2138 ExcNotImplemented());
2139 if (fe_eval.get_shape_info().data[0].fe_degree <
2140 fe_eval.get_shape_info().data[0].n_q_points_1d)
2141 evaluate_hessians_collocation<n_q_points_1d>(n_components, fe_eval);
2142 else
2143 evaluate_hessians_slow(n_components, values_dofs, fe_eval);
2144 }
2145
2146 return false;
2147 }
2148
2149 private:
2150 template <typename T>
2151 static void
2152 evaluate_or_integrate(
2153 const unsigned int n_components,
2154 const EvaluationFlags::EvaluationFlags evaluation_flag,
2155 const Number *values_dofs,
2156 FEEvaluationData<dim, Number, false> &fe_eval,
2157 const bool sum_into_values_array,
2158 std::bool_constant<false>)
2159 {
2160 (void)sum_into_values_array;
2161
2162 T::evaluate(n_components, evaluation_flag, values_dofs, fe_eval);
2163 }
2164
2165 template <typename T>
2166 static void
2167 evaluate_or_integrate(
2168 const unsigned int n_components,
2169 const EvaluationFlags::EvaluationFlags evaluation_flag,
2170 Number *values_dofs,
2171 FEEvaluationData<dim, Number, false> &fe_eval,
2172 const bool sum_into_values_array,
2173 std::bool_constant<true>)
2174 {
2175 T::integrate(n_components,
2176 evaluation_flag,
2177 values_dofs,
2178 fe_eval,
2179 sum_into_values_array);
2180 }
2181
2182 template <typename T, typename OtherNumber>
2183 static void
2184 evaluate_or_integrate(
2185 const unsigned int n_components,
2186 const EvaluationFlags::EvaluationFlags evaluation_flag,
2187 OtherNumber *values_dofs,
2188 FEEvaluationData<dim, Number, false> &fe_eval,
2189 const bool sum_into_values_array)
2190 {
2191 evaluate_or_integrate<T>(n_components,
2192 evaluation_flag,
2193 values_dofs,
2194 fe_eval,
2195 sum_into_values_array,
2196 std::bool_constant<do_integrate>());
2197 }
2198 };
2199
2200
2201
2206 template <int dim, typename Number>
2207 struct CellwiseInverseMassMatrixImplBasic
2208 {
2209 using Number2 =
2210 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
2211
2212 template <int fe_degree, int = 0>
2213 static bool
2214 run(const unsigned int n_components,
2215 const FEEvaluationData<dim, Number, false> &fe_eval,
2216 const Number *in_array,
2217 Number *out_array)
2218 {
2219 const unsigned int given_degree =
2220 (fe_degree > -1) ? fe_degree :
2221 fe_eval.get_shape_info().data.front().fe_degree;
2222
2223 const unsigned int dofs_per_component =
2224 Utilities::pow(given_degree + 1, dim);
2225
2226 Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
2227 Assert(fe_eval.get_shape_info().element_type <=
2228 MatrixFreeFunctions::tensor_symmetric_no_collocation,
2229 ExcNotImplemented());
2230
2231 EvaluatorTensorProduct<evaluate_evenodd,
2232 dim,
2233 fe_degree + 1,
2234 fe_degree + 1,
2235 Number,
2236 Number2>
2237 evaluator({},
2238 {},
2239 fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
2240 given_degree + 1,
2241 given_degree + 1);
2242
2243 for (unsigned int d = 0; d < n_components; ++d)
2244 {
2245 const Number *in = in_array + d * dofs_per_component;
2246 Number *out = out_array + d * dofs_per_component;
2247 // Need to select 'apply' method with hessian slot because values
2248 // assume symmetries that do not exist in the inverse shapes
2249 evaluator.template hessians<0, true, false>(in, out);
2250 if (dim > 1)
2251 evaluator.template hessians<1, true, false>(out, out);
2252 if (dim > 2)
2253 evaluator.template hessians<2, true, false>(out, out);
2254 }
2255 for (unsigned int q = 0; q < dofs_per_component; ++q)
2256 {
2257 const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
2258 for (unsigned int d = 0; d < n_components; ++d)
2259 out_array[q + d * dofs_per_component] *= inverse_JxW_q;
2260 }
2261 for (unsigned int d = 0; d < n_components; ++d)
2262 {
2263 Number *out = out_array + d * dofs_per_component;
2264 if (dim > 2)
2265 evaluator.template hessians<2, false, false>(out, out);
2266 if (dim > 1)
2267 evaluator.template hessians<1, false, false>(out, out);
2268 evaluator.template hessians<0, false, false>(out, out);
2269 }
2270 return false;
2271 }
2272 };
2273
2274
2275
2282 template <int dim, typename Number>
2283 struct CellwiseInverseMassMatrixImplFlexible
2284 {
2285 using Number2 =
2286 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
2287
2288 template <int fe_degree, int = 0>
2289 static bool
2290 run(const unsigned int n_desired_components,
2291 const FEEvaluationData<dim, Number, false> &fe_eval,
2292 const ArrayView<const Number> &inverse_coefficients,
2293 const bool dyadic_coefficients,
2294 const Number *in_array,
2295 Number *out_array)
2296 {
2297 const unsigned int given_degree =
2298 (fe_degree > -1) ? fe_degree :
2299 fe_eval.get_shape_info().data.front().fe_degree;
2300
2301 const unsigned int dofs_per_component =
2302 Utilities::pow(given_degree + 1, dim);
2303
2304 Assert(inverse_coefficients.size() > 0 &&
2305 inverse_coefficients.size() % dofs_per_component == 0,
2306 ExcMessage(
2307 "Expected diagonal to be a multiple of scalar dof per cells"));
2308
2309 if (!dyadic_coefficients)
2310 {
2311 if (inverse_coefficients.size() != dofs_per_component)
2312 AssertDimension(n_desired_components * dofs_per_component,
2313 inverse_coefficients.size());
2314 }
2315 else
2316 {
2317 AssertDimension(n_desired_components * n_desired_components *
2318 dofs_per_component,
2319 inverse_coefficients.size());
2320 }
2321
2322 Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
2323 Assert(fe_eval.get_shape_info().element_type <=
2324 MatrixFreeFunctions::tensor_symmetric_no_collocation,
2325 ExcNotImplemented());
2326
2327 EvaluatorTensorProduct<evaluate_evenodd,
2328 dim,
2329 fe_degree + 1,
2330 fe_degree + 1,
2331 Number,
2332 Number2>
2333 evaluator({},
2334 {},
2335 fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
2336 given_degree + 1,
2337 given_degree + 1);
2338
2339 const Number *in = in_array;
2340 Number *out = out_array;
2341
2342 const Number *inv_coefficient = inverse_coefficients.data();
2343
2344 const unsigned int shift_coefficient =
2345 inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
2346 0;
2347
2348 const auto n_comp_outer = dyadic_coefficients ? 1 : n_desired_components;
2349 const auto n_comp_inner = dyadic_coefficients ? n_desired_components : 1;
2350
2351 for (unsigned int d = 0; d < n_comp_outer; ++d)
2352 {
2353 for (unsigned int di = 0; di < n_comp_inner; ++di)
2354 {
2355 const Number *in_ = in + di * dofs_per_component;
2356 Number *out_ = out + di * dofs_per_component;
2357 evaluator.template hessians<0, true, false>(in_, out_);
2358 if (dim > 1)
2359 evaluator.template hessians<1, true, false>(out_, out_);
2360 if (dim > 2)
2361 evaluator.template hessians<2, true, false>(out_, out_);
2362 }
2363 if (dyadic_coefficients)
2364 {
2365 const auto n_coeff_components =
2366 n_desired_components * n_desired_components;
2367 if (n_desired_components == dim)
2368 {
2369 for (unsigned int q = 0; q < dofs_per_component; ++q)
2370 vmult<dim>(&inv_coefficient[q * n_coeff_components],
2371 &in[q],
2372 &out[q],
2373 dofs_per_component);
2374 }
2375 else
2376 {
2377 for (unsigned int q = 0; q < dofs_per_component; ++q)
2378 vmult<-1>(&inv_coefficient[q * n_coeff_components],
2379 &in[q],
2380 &out[q],
2381 dofs_per_component,
2382 n_desired_components);
2383 }
2384 }
2385 else
2386 for (unsigned int q = 0; q < dofs_per_component; ++q)
2387 out[q] *= inv_coefficient[q];
2388
2389 for (unsigned int di = 0; di < n_comp_inner; ++di)
2390 {
2391 Number *out_ = out + di * dofs_per_component;
2392 if (dim > 2)
2393 evaluator.template hessians<2, false, false>(out_, out_);
2394 if (dim > 1)
2395 evaluator.template hessians<1, false, false>(out_, out_);
2396 evaluator.template hessians<0, false, false>(out_, out_);
2397 }
2398
2399 in += dofs_per_component;
2400 out += dofs_per_component;
2401 inv_coefficient += shift_coefficient;
2402 }
2403
2404 return false;
2405 }
2406
2407 private:
2408 template <int n_components>
2409 static inline void
2410 vmult(const Number *inverse_coefficients,
2411 const Number *src,
2412 Number *dst,
2413 const unsigned int dofs_per_component,
2414 const unsigned int n_given_components = 0)
2415 {
2416 const unsigned int n_desired_components =
2417 (n_components > -1) ? n_components : n_given_components;
2418
2419 std::array<Number, dim + 2> tmp = {};
2420 Assert(n_desired_components <= dim + 2,
2421 ExcMessage(
2422 "Number of components larger than dim+2 not supported."));
2423
2424 for (unsigned int d = 0; d < n_desired_components; ++d)
2425 tmp[d] = src[d * dofs_per_component];
2426
2427 for (unsigned int d1 = 0; d1 < n_desired_components; ++d1)
2428 {
2429 const Number *inv_coeff_row =
2430 &inverse_coefficients[d1 * n_desired_components];
2431 Number sum = inv_coeff_row[0] * tmp[0];
2432 for (unsigned int d2 = 1; d2 < n_desired_components; ++d2)
2433 sum += inv_coeff_row[d2] * tmp[d2];
2434 dst[d1 * dofs_per_component] = sum;
2435 }
2436 }
2437 };
2438
2439
2440
2447 template <int dim, typename Number>
2448 struct CellwiseInverseMassMatrixImplTransformFromQPoints
2449 {
2450 template <int fe_degree, int n_q_points_1d>
2451 static bool
2452 run(const unsigned int n_desired_components,
2453 const FEEvaluationData<dim, Number, false> &fe_eval,
2454 const Number *in_array,
2455 Number *out_array)
2456 {
2457 static const bool do_inplace =
2458 fe_degree > -1 && (fe_degree + 1 == n_q_points_1d);
2459
2460 Assert(fe_eval.get_shape_info().element_type !=
2461 MatrixFreeFunctions::ElementType::tensor_none,
2462 ExcNotImplemented());
2463
2464 const auto &inverse_shape =
2465 do_inplace ?
2466 fe_eval.get_shape_info().data.front().inverse_shape_values_eo :
2467 fe_eval.get_shape_info().data.front().inverse_shape_values;
2468
2469 const std::size_t dofs_per_component =
2470 do_inplace ? Utilities::pow(fe_degree + 1, dim) :
2471 fe_eval.get_shape_info().dofs_per_component_on_cell;
2472 const std::size_t n_q_points = do_inplace ?
2473 Utilities::pow(fe_degree + 1, dim) :
2474 fe_eval.get_shape_info().n_q_points;
2475
2476 using Number2 =
2477 typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
2478 EvaluatorTensorProduct<do_inplace ? evaluate_evenodd : evaluate_general,
2479 dim,
2480 fe_degree + 1,
2481 n_q_points_1d,
2482 Number,
2483 Number2>
2484 evaluator({},
2485 {},
2486 inverse_shape,
2487 fe_eval.get_shape_info().data.front().fe_degree + 1,
2488 fe_eval.get_shape_info().data.front().n_q_points_1d);
2489
2490 for (unsigned int d = 0; d < n_desired_components; ++d)
2491 {
2492 const Number *in = in_array + d * n_q_points;
2493 Number *out = out_array + d * dofs_per_component;
2494
2495 auto *temp_1 = do_inplace ? out : fe_eval.get_scratch_data().begin();
2496 auto *temp_2 = do_inplace ?
2497 out :
2498 (temp_1 + std::max(n_q_points, dofs_per_component));
2499
2500 if (dim == 3)
2501 {
2502 evaluator.template hessians<2, false, false>(in, temp_1);
2503 evaluator.template hessians<1, false, false>(temp_1, temp_2);
2504 evaluator.template hessians<0, false, false>(temp_2, out);
2505 }
2506 if (dim == 2)
2507 {
2508 evaluator.template hessians<1, false, false>(in, temp_1);
2509 evaluator.template hessians<0, false, false>(temp_1, out);
2510 }
2511 if (dim == 1)
2512 evaluator.template hessians<0, false, false>(in, out);
2513 }
2514 return false;
2515 }
2516 };
2517} // end of namespace internal
2518
2519
2520DEAL_II_NAMESPACE_CLOSE
2521
2522#endif
AlignedVector::data
pointer data()
AlignedVector::size
size_type size() const
ArrayView::begin
iterator begin() const
Definition array_view.h:702
ArrayView::data
value_type * data() const noexcept
Definition array_view.h:661
ArrayView::size
std::size_t size() const
Definition array_view.h:684
FEEvaluationData::shape_info_number_type
ScalarNumber shape_info_number_type
Definition fe_evaluation_data.h:127
FEEvaluationData::get_shape_info
const ShapeInfoType & get_shape_info() const
FEEvaluationData::JxW
Number JxW(const unsigned int q_point) const
FEEvaluationData::begin_gradients
const Number * begin_gradients() const
FEEvaluationData::get_scratch_data
ArrayView< Number > get_scratch_data() const
FEEvaluationData::begin_values
const Number * begin_values() const
FEEvaluationData::begin_hessians
const Number * begin_hessians() const
DEAL_II_ALWAYS_INLINE
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
evaluation_flags.h
fe_evaluation_data.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1638
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1985
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1739
internal::MatrixFreeFunctions::ElementType
ElementType
Definition shape_info.h:58
internal::MatrixFreeFunctions::tensor_symmetric_no_collocation
@ tensor_symmetric_no_collocation
Definition shape_info.h:87
internal::MatrixFreeFunctions::tensor_raviart_thomas
@ tensor_raviart_thomas
Definition shape_info.h:112
internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0
@ tensor_symmetric_plus_dg0
Definition shape_info.h:105
internal::MatrixFreeFunctions::truncated_tensor
@ truncated_tensor
Definition shape_info.h:98
internal::MatrixFreeFunctions::tensor_none
@ tensor_none
Definition shape_info.h:117
internal::MatrixFreeFunctions::tensor_general
@ tensor_general
Definition shape_info.h:92
EvaluationFlags::EvaluationFlags
EvaluationFlags
The EvaluationFlags enum.
Definition evaluation_flags.h:42
Utilities::fixed_power
constexpr T fixed_power(const T t)
Definition utilities.h:942
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
internal::evaluate_hessians_collocation
void evaluate_hessians_collocation(const unsigned int n_components, FEEvaluationData< dim, Number, false > &fe_eval)
Definition evaluation_kernels.h:1279
internal::use_collocation_evaluation
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
Definition tensor_product_kernels.h:2147
internal::integrate_gradients_collocation
void integrate_gradients_collocation(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &shape, Number *values, const Number *gradients, const bool add_into_values_array)
Definition evaluation_kernels.h:1216
internal::evaluate_hessians_slow
void evaluate_hessians_slow(const unsigned int n_components, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
Definition evaluation_kernels.h:1429
internal::apply_matrix_vector_product
std::enable_if_t<(variant==evaluate_general), void > apply_matrix_vector_product(const Number2 *matrix, const Number *in, Number *out)
Definition tensor_product_kernels.h:118
internal::integrate_hessians_collocation
void integrate_hessians_collocation(const unsigned int n_components, FEEvaluationData< dim, Number, false > &fe_eval, const bool add_into_values_array)
Definition evaluation_kernels.h:1352
internal::evaluate_gradients_collocation
void evaluate_gradients_collocation(const MatrixFreeFunctions::UnivariateShapeData< Number2 > &shape, const Number *values, Number *gradients)
Definition evaluation_kernels.h:1162
internal::integrate_hessians_slow
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)
Definition evaluation_kernels.h:1523
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:220
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6960
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6943
shape_info.h
internal::CellwiseInverseMassMatrixImplBasic::Number2
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
Definition evaluation_kernels.h:2209
internal::CellwiseInverseMassMatrixImplBasic::run
static bool run(const unsigned int n_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array)
Definition evaluation_kernels.h:2214
internal::CellwiseInverseMassMatrixImplFlexible::Number2
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
Definition evaluation_kernels.h:2285
internal::CellwiseInverseMassMatrixImplFlexible::run
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)
Definition evaluation_kernels.h:2290
internal::CellwiseInverseMassMatrixImplFlexible::vmult
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)
Definition evaluation_kernels.h:2410
internal::CellwiseInverseMassMatrixImplTransformFromQPoints::run
static bool run(const unsigned int n_desired_components, const FEEvaluationData< dim, Number, false > &fe_eval, const Number *in_array, Number *out_array)
Definition evaluation_kernels.h:2452
internal::FEEvaluationImplBasisChange::do_backward
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)
Definition evaluation_kernels.h:930
internal::FEEvaluationImplBasisChange::do_forward
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)
Definition evaluation_kernels.h:809
internal::FEEvaluationImplBasisChange::do_mass
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)
Definition evaluation_kernels.h:1080
internal::FEEvaluationImplCollocation::evaluate
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
Definition evaluation_kernels.h:1649
internal::FEEvaluationImplCollocation::integrate
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)
Definition evaluation_kernels.h:1672
internal::FEEvaluationImplCollocation::Number2
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
Definition evaluation_kernels.h:1639
internal::FEEvaluationImplSelector::evaluate_or_integrate
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)
Definition evaluation_kernels.h:2184
internal::FEEvaluationImplSelector::evaluate_or_integrate
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 >)
Definition evaluation_kernels.h:2167
internal::FEEvaluationImplSelector::run
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)
Definition evaluation_kernels.h:1975
internal::FEEvaluationImplSelector::evaluate_or_integrate
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 >)
Definition evaluation_kernels.h:2152
internal::FEEvaluationImplTransformToCollocation::integrate
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)
Definition evaluation_kernels.h:1758
internal::FEEvaluationImplTransformToCollocation::evaluate
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
Definition evaluation_kernels.h:1721
internal::FEEvaluationImpl< MatrixFreeFunctions::tensor_raviart_thomas, dim, fe_degree, n_q_points_1d, Number >::Number2
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
Definition evaluation_kernels.h:1814
internal::FEEvaluationImpl::variant
static const EvaluatorVariant variant
Definition evaluation_kernels.h:114
internal::FEEvaluationImpl::integrate
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)
Definition evaluation_kernels.h:391
internal::FEEvaluationImpl::create_evaluator_tensor_product
static Eval create_evaluator_tensor_product(const MatrixFreeFunctions::UnivariateShapeData< Number2 > *univariate_shape_data)
Definition evaluation_kernels.h:140
internal::FEEvaluationImpl::Number2
typename FEEvaluationData< dim, Number, false >::shape_info_number_type Number2
Definition evaluation_kernels.h:116
internal::FEEvaluationImpl::evaluate
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs_actual, FEEvaluationData< dim, Number, false > &fe_eval)
Definition evaluation_kernels.h:194
internal::FEEvaluationImpl::Eval
EvaluatorTensorProduct< variant, dim, fe_degree+1, n_q_points_1d, Number, Number2 > Eval
Definition evaluation_kernels.h:119
internal::MatrixFreeFunctions::ShapeInfo::data
std::vector< UnivariateShapeData< Number > > data
Definition shape_info.h:485
tensor_product_kernels.h
Generated by doxygen 1.11.0