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