Reference documentation for deal.II version 9.1.0-pre
evaluation_kernels.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2018 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>
23 #include <deal.II/base/vectorization.h>
24 
25 #include <deal.II/matrix_free/shape_info.h>
26 #include <deal.II/matrix_free/tensor_product_kernels.h>
27 
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 
32 
33 namespace 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 
104  template <MatrixFreeFunctions::ElementType type,
105  int dim,
106  int fe_degree,
107  int n_q_points_1d,
108  int n_components,
109  typename Number>
110  struct FEEvaluationImpl
111  {
112  static void
113  evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
114  const Number * values_dofs_actual,
115  Number * values_quad,
116  Number * gradients_quad,
117  Number * hessians_quad,
118  Number * scratch_data,
119  const bool evaluate_values,
120  const bool evaluate_gradients,
121  const bool evaluate_hessians);
122 
123  static void
124  integrate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
125  Number * values_dofs_actual,
126  Number * values_quad,
127  Number * gradients_quad,
128  Number * scratch_data,
129  const bool integrate_values,
130  const bool integrate_gradients,
131  const bool add_into_values_array);
132  };
133 
134 
135 
136  template <MatrixFreeFunctions::ElementType type,
137  int dim,
138  int fe_degree,
139  int n_q_points_1d,
140  int n_components,
141  typename Number>
142  inline void
143  FEEvaluationImpl<type, dim, fe_degree, n_q_points_1d, n_components, Number>::
144  evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
145  const Number * values_dofs_actual,
146  Number * values_quad,
147  Number * gradients_quad,
148  Number * hessians_quad,
149  Number * scratch_data,
150  const bool evaluate_values,
151  const bool evaluate_gradients,
152  const bool evaluate_hessians)
153  {
154  if (evaluate_values == false && evaluate_gradients == false &&
155  evaluate_hessians == false)
156  return;
157 
158  const EvaluatorVariant variant =
159  EvaluatorSelector<type, (fe_degree + n_q_points_1d > 4)>::variant;
160  using Eval = EvaluatorTensorProduct<variant,
161  dim,
162  fe_degree + 1,
163  n_q_points_1d,
164  Number>;
165  Eval eval(variant == evaluate_evenodd ? shape_info.shape_values_eo :
166  shape_info.shape_values,
167  variant == evaluate_evenodd ? shape_info.shape_gradients_eo :
168  shape_info.shape_gradients,
169  variant == evaluate_evenodd ? shape_info.shape_hessians_eo :
170  shape_info.shape_hessians,
171  shape_info.fe_degree + 1,
172  shape_info.n_q_points_1d);
173 
174  const unsigned int temp_size =
175  Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
176  0 :
177  (Eval::n_rows_of_product > Eval::n_columns_of_product ?
178  Eval::n_rows_of_product :
179  Eval::n_columns_of_product);
180  Number *temp1;
181  Number *temp2;
182  if (temp_size == 0)
183  {
184  temp1 = scratch_data;
185  temp2 = temp1 +
186  std::max(Utilities::fixed_power<dim>(shape_info.fe_degree + 1),
187  Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
188  }
189  else
190  {
191  temp1 = scratch_data;
192  temp2 = temp1 + temp_size;
193  }
194 
195  const unsigned int n_q_points =
196  temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
197  const unsigned int dofs_per_comp =
198  (type == MatrixFreeFunctions::truncated_tensor) ?
199  Utilities::fixed_power<dim>(shape_info.fe_degree + 1) :
200  shape_info.dofs_per_component_on_cell;
201  const Number *values_dofs = values_dofs_actual;
202  if (type == MatrixFreeFunctions::truncated_tensor)
203  {
204  Number *values_dofs_tmp =
205  scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
206  shape_info.n_q_points));
207  const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
208  unsigned int count_p = 0, count_q = 0;
209  for (int i = 0; i < (dim > 2 ? degree + 1 : 1); ++i)
210  {
211  for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
212  {
213  for (int k = 0; k < degree + 1 - j - i;
214  ++k, ++count_p, ++count_q)
215  for (unsigned int c = 0; c < n_components; ++c)
216  values_dofs_tmp[c * dofs_per_comp + count_q] =
217  values_dofs_actual
218  [c * shape_info.dofs_per_component_on_cell + count_p];
219  for (int k = degree + 1 - j - i; k < degree + 1; ++k, ++count_q)
220  for (unsigned int c = 0; c < n_components; ++c)
221  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
222  }
223  for (int j = degree + 1 - i; j < degree + 1; ++j)
224  for (int k = 0; k < degree + 1; ++k, ++count_q)
225  for (unsigned int c = 0; c < n_components; ++c)
226  values_dofs_tmp[c * dofs_per_comp + count_q] = Number();
227  }
228  AssertDimension(count_q, dofs_per_comp);
229  values_dofs = values_dofs_tmp;
230  }
231 
232  switch (dim)
233  {
234  case 1:
235  for (unsigned int c = 0; c < n_components; c++)
236  {
237  if (evaluate_values == true)
238  eval.template values<0, true, false>(values_dofs, values_quad);
239  if (evaluate_gradients == true)
240  eval.template gradients<0, true, false>(values_dofs,
241  gradients_quad);
242  if (evaluate_hessians == true)
243  eval.template hessians<0, true, false>(values_dofs,
244  hessians_quad);
245 
246  // advance the next component in 1D array
247  values_dofs += dofs_per_comp;
248  values_quad += n_q_points;
249  gradients_quad += n_q_points;
250  hessians_quad += n_q_points;
251  }
252  break;
253 
254  case 2:
255  for (unsigned int c = 0; c < n_components; c++)
256  {
257  // grad x
258  if (evaluate_gradients == true)
259  {
260  eval.template gradients<0, true, false>(values_dofs, temp1);
261  eval.template values<1, true, false>(temp1, gradients_quad);
262  }
263  if (evaluate_hessians == true)
264  {
265  // grad xy
266  if (evaluate_gradients == false)
267  eval.template gradients<0, true, false>(values_dofs, temp1);
268  eval.template gradients<1, true, false>(temp1,
269  hessians_quad +
270  2 * n_q_points);
271 
272  // grad xx
273  eval.template hessians<0, true, false>(values_dofs, temp1);
274  eval.template values<1, true, false>(temp1, hessians_quad);
275  }
276 
277  // grad y
278  eval.template values<0, true, false>(values_dofs, temp1);
279  if (evaluate_gradients == true)
280  eval.template gradients<1, true, false>(temp1,
281  gradients_quad +
282  n_q_points);
283 
284  // grad yy
285  if (evaluate_hessians == true)
286  eval.template hessians<1, true, false>(temp1,
287  hessians_quad +
288  n_q_points);
289 
290  // val: can use values applied in x
291  if (evaluate_values == true)
292  eval.template values<1, true, false>(temp1, values_quad);
293 
294  // advance to the next component in 1D array
295  values_dofs += dofs_per_comp;
296  values_quad += n_q_points;
297  gradients_quad += 2 * n_q_points;
298  hessians_quad += 3 * n_q_points;
299  }
300  break;
301 
302  case 3:
303  for (unsigned int c = 0; c < n_components; c++)
304  {
305  if (evaluate_gradients == true)
306  {
307  // grad x
308  eval.template gradients<0, true, false>(values_dofs, temp1);
309  eval.template values<1, true, false>(temp1, temp2);
310  eval.template values<2, true, false>(temp2, gradients_quad);
311  }
312 
313  if (evaluate_hessians == true)
314  {
315  // grad xz
316  if (evaluate_gradients == false)
317  {
318  eval.template gradients<0, true, false>(values_dofs,
319  temp1);
320  eval.template values<1, true, false>(temp1, temp2);
321  }
322  eval.template gradients<2, true, false>(temp2,
323  hessians_quad +
324  4 * n_q_points);
325 
326  // grad xy
327  eval.template gradients<1, true, false>(temp1, temp2);
328  eval.template values<2, true, false>(temp2,
329  hessians_quad +
330  3 * n_q_points);
331 
332  // grad xx
333  eval.template hessians<0, true, false>(values_dofs, temp1);
334  eval.template values<1, true, false>(temp1, temp2);
335  eval.template values<2, true, false>(temp2, hessians_quad);
336  }
337 
338  // grad y
339  eval.template values<0, true, false>(values_dofs, temp1);
340  if (evaluate_gradients == true)
341  {
342  eval.template gradients<1, true, false>(temp1, temp2);
343  eval.template values<2, true, false>(temp2,
344  gradients_quad +
345  n_q_points);
346  }
347 
348  if (evaluate_hessians == true)
349  {
350  // grad yz
351  if (evaluate_gradients == false)
352  eval.template gradients<1, true, false>(temp1, temp2);
353  eval.template gradients<2, true, false>(temp2,
354  hessians_quad +
355  5 * n_q_points);
356 
357  // grad yy
358  eval.template hessians<1, true, false>(temp1, temp2);
359  eval.template values<2, true, false>(temp2,
360  hessians_quad +
361  n_q_points);
362  }
363 
364  // grad z: can use the values applied in x direction stored in
365  // temp1
366  eval.template values<1, true, false>(temp1, temp2);
367  if (evaluate_gradients == true)
368  eval.template gradients<2, true, false>(temp2,
369  gradients_quad +
370  2 * n_q_points);
371 
372  // grad zz: can use the values applied in x and y direction stored
373  // in temp2
374  if (evaluate_hessians == true)
375  eval.template hessians<2, true, false>(temp2,
376  hessians_quad +
377  2 * n_q_points);
378 
379  // val: can use the values applied in x & y direction stored in
380  // temp2
381  if (evaluate_values == true)
382  eval.template values<2, true, false>(temp2, values_quad);
383 
384  // advance to the next component in 1D array
385  values_dofs += dofs_per_comp;
386  values_quad += n_q_points;
387  gradients_quad += 3 * n_q_points;
388  hessians_quad += 6 * n_q_points;
389  }
390  break;
391 
392  default:
393  AssertThrow(false, ExcNotImplemented());
394  }
395 
396  // case additional dof for FE_Q_DG0: add values; gradients and second
397  // derivatives evaluate to zero
398  if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 &&
399  evaluate_values)
400  {
401  values_quad -= n_components * n_q_points;
402  values_dofs -= n_components * dofs_per_comp;
403  for (unsigned int c = 0; c < n_components; ++c)
404  for (unsigned int q = 0; q < shape_info.n_q_points; ++q)
405  values_quad[c * shape_info.n_q_points + q] +=
406  values_dofs[(c + 1) * shape_info.dofs_per_component_on_cell - 1];
407  }
408  }
409 
410 
411 
412  template <MatrixFreeFunctions::ElementType type,
413  int dim,
414  int fe_degree,
415  int n_q_points_1d,
416  int n_components,
417  typename Number>
418  inline void
419  FEEvaluationImpl<type, dim, fe_degree, n_q_points_1d, n_components, Number>::
420  integrate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
421  Number * values_dofs_actual,
422  Number * values_quad,
423  Number * gradients_quad,
424  Number * scratch_data,
425  const bool integrate_values,
426  const bool integrate_gradients,
427  const bool add_into_values_array)
428  {
429  const EvaluatorVariant variant =
430  EvaluatorSelector<type, (fe_degree + n_q_points_1d > 4)>::variant;
431  using Eval = EvaluatorTensorProduct<variant,
432  dim,
433  fe_degree + 1,
434  n_q_points_1d,
435  Number>;
436  Eval eval(variant == evaluate_evenodd ? shape_info.shape_values_eo :
437  shape_info.shape_values,
438  variant == evaluate_evenodd ? shape_info.shape_gradients_eo :
439  shape_info.shape_gradients,
440  variant == evaluate_evenodd ? shape_info.shape_hessians_eo :
441  shape_info.shape_hessians,
442  shape_info.fe_degree + 1,
443  shape_info.n_q_points_1d);
444 
445  const unsigned int temp_size =
446  Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
447  0 :
448  (Eval::n_rows_of_product > Eval::n_columns_of_product ?
449  Eval::n_rows_of_product :
450  Eval::n_columns_of_product);
451  Number *temp1;
452  Number *temp2;
453  if (temp_size == 0)
454  {
455  temp1 = scratch_data;
456  temp2 = temp1 +
457  std::max(Utilities::fixed_power<dim>(shape_info.fe_degree + 1),
458  Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
459  }
460  else
461  {
462  temp1 = scratch_data;
463  temp2 = temp1 + temp_size;
464  }
465 
466  const unsigned int n_q_points =
467  temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
468  const unsigned int dofs_per_comp =
469  (type == MatrixFreeFunctions::truncated_tensor) ?
470  Utilities::fixed_power<dim>(shape_info.fe_degree + 1) :
471  shape_info.dofs_per_component_on_cell;
472  // expand dof_values to tensor product for truncated tensor products
473  Number *values_dofs =
474  (type == MatrixFreeFunctions::truncated_tensor) ?
475  scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
476  shape_info.n_q_points)) :
477  values_dofs_actual;
478 
479  switch (dim)
480  {
481  case 1:
482  for (unsigned int c = 0; c < n_components; c++)
483  {
484  if (integrate_values == true)
485  {
486  if (add_into_values_array == false)
487  eval.template values<0, false, false>(values_quad,
488  values_dofs);
489  else
490  eval.template values<0, false, true>(values_quad,
491  values_dofs);
492  }
493  if (integrate_gradients == true)
494  {
495  if (integrate_values == true || add_into_values_array == true)
496  eval.template gradients<0, false, true>(gradients_quad,
497  values_dofs);
498  else
499  eval.template gradients<0, false, false>(gradients_quad,
500  values_dofs);
501  }
502 
503  // advance to the next component in 1D array
504  values_dofs += dofs_per_comp;
505  values_quad += n_q_points;
506  gradients_quad += n_q_points;
507  }
508  break;
509 
510  case 2:
511  for (unsigned int c = 0; c < n_components; c++)
512  {
513  if (integrate_values == true && integrate_gradients == false)
514  {
515  eval.template values<1, false, false>(values_quad, temp1);
516  if (add_into_values_array == false)
517  eval.template values<0, false, false>(temp1, values_dofs);
518  else
519  eval.template values<0, false, true>(temp1, values_dofs);
520  }
521  if (integrate_gradients == true)
522  {
523  eval.template gradients<1, false, false>(gradients_quad +
524  n_q_points,
525  temp1);
526  if (integrate_values)
527  eval.template values<1, false, true>(values_quad, temp1);
528  if (add_into_values_array == false)
529  eval.template values<0, false, false>(temp1, values_dofs);
530  else
531  eval.template values<0, false, true>(temp1, values_dofs);
532  eval.template values<1, false, false>(gradients_quad, temp1);
533  eval.template gradients<0, false, true>(temp1, values_dofs);
534  }
535 
536  // advance to the next component in 1D array
537  values_dofs += dofs_per_comp;
538  values_quad += n_q_points;
539  gradients_quad += 2 * n_q_points;
540  }
541  break;
542 
543  case 3:
544  for (unsigned int c = 0; c < n_components; c++)
545  {
546  if (integrate_values == true && integrate_gradients == false)
547  {
548  eval.template values<2, false, false>(values_quad, temp1);
549  eval.template values<1, false, false>(temp1, temp2);
550  if (add_into_values_array == false)
551  eval.template values<0, false, false>(temp2, values_dofs);
552  else
553  eval.template values<0, false, true>(temp2, values_dofs);
554  }
555  if (integrate_gradients == true)
556  {
557  eval.template gradients<2, false, false>(gradients_quad +
558  2 * n_q_points,
559  temp1);
560  if (integrate_values)
561  eval.template values<2, false, true>(values_quad, temp1);
562  eval.template values<1, false, false>(temp1, temp2);
563  eval.template values<2, false, false>(gradients_quad +
564  n_q_points,
565  temp1);
566  eval.template gradients<1, false, true>(temp1, temp2);
567  if (add_into_values_array == false)
568  eval.template values<0, false, false>(temp2, values_dofs);
569  else
570  eval.template values<0, false, true>(temp2, values_dofs);
571  eval.template values<2, false, false>(gradients_quad, temp1);
572  eval.template values<1, false, false>(temp1, temp2);
573  eval.template gradients<0, false, true>(temp2, values_dofs);
574  }
575 
576  // advance to the next component in 1D array
577  values_dofs += dofs_per_comp;
578  values_quad += n_q_points;
579  gradients_quad += 3 * n_q_points;
580  }
581  break;
582 
583  default:
584  AssertThrow(false, ExcNotImplemented());
585  }
586 
587  // case FE_Q_DG0: add values, gradients and second derivatives are zero
588  if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
589  {
590  values_dofs -= n_components * dofs_per_comp -
591  shape_info.dofs_per_component_on_cell + 1;
592  values_quad -= n_components * n_q_points;
593  if (integrate_values)
594  for (unsigned int c = 0; c < n_components; ++c)
595  {
596  values_dofs[0] = values_quad[0];
597  for (unsigned int q = 1; q < shape_info.n_q_points; ++q)
598  values_dofs[0] += values_quad[q];
599  values_dofs += dofs_per_comp;
600  values_quad += n_q_points;
601  }
602  else
603  {
604  for (unsigned int c = 0; c < n_components; ++c)
605  values_dofs[c * shape_info.dofs_per_component_on_cell] = Number();
606  values_dofs += n_components * shape_info.dofs_per_component_on_cell;
607  }
608  }
609 
610  if (type == MatrixFreeFunctions::truncated_tensor)
611  {
612  values_dofs -= dofs_per_comp * n_components;
613  unsigned int count_p = 0, count_q = 0;
614  const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
615  for (int i = 0; i < (dim > 2 ? degree + 1 : 1); ++i)
616  {
617  for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
618  {
619  for (int k = 0; k < degree + 1 - j - i;
620  ++k, ++count_p, ++count_q)
621  {
622  for (unsigned int c = 0; c < n_components; ++c)
623  values_dofs_actual
624  [c * shape_info.dofs_per_component_on_cell + count_p] =
625  values_dofs[c * dofs_per_comp + count_q];
626  }
627  count_q += j + i;
628  }
629  count_q += i * (degree + 1);
630  }
631  AssertDimension(count_q,
632  Utilities::fixed_power<dim>(shape_info.fe_degree + 1));
633  }
634  }
635 
636 
637 
649  template <EvaluatorVariant variant,
650  int dim,
651  int basis_size_1,
652  int basis_size_2,
653  int n_components,
654  typename Number,
655  typename Number2>
656  struct FEEvaluationImplBasisChange
657  {
658  static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
659  "The second dimension must not be smaller than the first");
660 
682 #ifndef DEBUG
683  DEAL_II_ALWAYS_INLINE
684 #endif
685  static void
686  do_forward(
687  const AlignedVector<Number2> &transformation_matrix,
688  const Number * values_in,
689  Number * values_out,
690  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
691  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
692  {
693  Assert(
694  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
695  ExcMessage("The second dimension must not be smaller than the first"));
696 
697  // we do recursion until dim==1 or dim==2 and we have
698  // basis_size_1==basis_size_2. The latter optimization increases
699  // optimization possibilities for the compiler but does only work for
700  // aliased pointers if the sizes are equal.
701  constexpr int next_dim =
702  (dim > 2 ||
703  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
704  dim - 1 :
705  dim;
706 
707  EvaluatorTensorProduct<variant,
708  dim,
709  basis_size_1,
710  (basis_size_1 == 0 ? 0 : basis_size_2),
711  Number,
712  Number2>
713  eval_val(transformation_matrix,
714  AlignedVector<Number2>(),
715  AlignedVector<Number2>(),
716  basis_size_1_variable,
717  basis_size_2_variable);
718  const unsigned int np_1 =
719  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
720  const unsigned int np_2 =
721  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
722  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
723  ExcMessage("Cannot transform with 0-point basis"));
724  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
725  ExcMessage("Cannot transform with 0-point basis"));
726 
727  // run loop backwards to ensure correctness if values_in aliases with
728  // values_out in case with basis_size_1 < basis_size_2
729  values_in = values_in + n_components * Utilities::fixed_power<dim>(np_1);
730  values_out =
731  values_out + n_components * Utilities::fixed_power<dim>(np_2);
732  for (unsigned int c = n_components; c != 0; --c)
733  {
734  values_in -= Utilities::fixed_power<dim>(np_1);
735  values_out -= Utilities::fixed_power<dim>(np_2);
736  if (next_dim < dim)
737  for (unsigned int q = np_1; q != 0; --q)
738  FEEvaluationImplBasisChange<
739  variant,
740  next_dim,
741  basis_size_1,
742  basis_size_2,
743  1,
744  Number,
745  Number2>::do_forward(transformation_matrix,
746  values_in +
747  (q - 1) *
748  Utilities::fixed_power<next_dim>(np_1),
749  values_out +
750  (q - 1) *
751  Utilities::fixed_power<next_dim>(np_2),
752  basis_size_1_variable,
753  basis_size_2_variable);
754 
755  // the recursion stops if dim==1 or if dim==2 and
756  // basis_size_1==basis_size_2 (the latter is used because the
757  // compiler generates nicer code)
758  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
759  {
760  eval_val.template values<0, true, false>(values_in, values_out);
761  eval_val.template values<1, true, false>(values_out, values_out);
762  }
763  else if (dim == 1)
764  eval_val.template values<dim - 1, true, false>(values_in,
765  values_out);
766  else
767  eval_val.template values<dim - 1, true, false>(values_out,
768  values_out);
769  }
770  }
771 
801 #ifndef DEBUG
802  DEAL_II_ALWAYS_INLINE
803 #endif
804  static void
805  do_backward(
806  const AlignedVector<Number2> &transformation_matrix,
807  const bool add_into_result,
808  Number * values_in,
809  Number * values_out,
810  const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
811  const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
812  {
813  Assert(
814  basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable,
815  ExcMessage("The second dimension must not be smaller than the first"));
816  Assert(add_into_result == false || values_in != values_out,
817  ExcMessage(
818  "Input and output cannot alias with each other when "
819  "adding the result of the basis change to existing data"));
820 
821  constexpr int next_dim =
822  (dim > 2 ||
823  ((basis_size_1 == 0 || basis_size_2 > basis_size_1) && dim > 1)) ?
824  dim - 1 :
825  dim;
826  EvaluatorTensorProduct<variant,
827  dim,
828  basis_size_1,
829  (basis_size_1 == 0 ? 0 : basis_size_2),
830  Number,
831  Number2>
832  eval_val(transformation_matrix,
833  AlignedVector<Number2>(),
834  AlignedVector<Number2>(),
835  basis_size_1_variable,
836  basis_size_2_variable);
837  const unsigned int np_1 =
838  basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
839  const unsigned int np_2 =
840  basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
841  Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
842  ExcMessage("Cannot transform with 0-point basis"));
843  Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
844  ExcMessage("Cannot transform with 0-point basis"));
845 
846  for (unsigned int c = 0; c < n_components; ++c)
847  {
848  if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
849  {
850  eval_val.template values<1, false, false>(values_in, values_in);
851  if (add_into_result)
852  eval_val.template values<0, false, true>(values_in, values_out);
853  else
854  eval_val.template values<0, false, false>(values_in,
855  values_out);
856  }
857  else
858  {
859  if (dim == 1 && add_into_result)
860  eval_val.template values<0, false, true>(values_in, values_out);
861  else if (dim == 1)
862  eval_val.template values<0, false, false>(values_in,
863  values_out);
864  else
865  eval_val.template values<dim - 1, false, false>(values_in,
866  values_in);
867  }
868  if (next_dim < dim)
869  for (unsigned int q = 0; q < np_1; ++q)
870  FEEvaluationImplBasisChange<variant,
871  next_dim,
872  basis_size_1,
873  basis_size_2,
874  1,
875  Number,
876  Number2>::
877  do_backward(transformation_matrix,
878  add_into_result,
879  values_in +
880  q * Utilities::fixed_power<next_dim>(np_2),
881  values_out +
882  q * Utilities::fixed_power<next_dim>(np_1),
883  basis_size_1_variable,
884  basis_size_2_variable);
885 
886  values_in += Utilities::fixed_power<dim>(np_2);
887  values_out += Utilities::fixed_power<dim>(np_1);
888  }
889  }
890 
910  static void
911  do_mass(const AlignedVector<Number2> &transformation_matrix,
912  const AlignedVector<Number> & coefficients,
913  const Number * values_in,
914  Number * scratch_data,
915  Number * values_out)
916  {
917  constexpr int next_dim = dim > 1 ? dim - 1 : dim;
918  Number * my_scratch =
919  basis_size_1 != basis_size_2 ? scratch_data : values_out;
920 
921  const unsigned int size_per_component = Utilities::pow(basis_size_2, dim);
922  Assert(coefficients.size() == size_per_component ||
923  coefficients.size() == n_components * size_per_component,
924  ExcDimensionMismatch(coefficients.size(), size_per_component));
925  const unsigned int stride =
926  coefficients.size() == size_per_component ? 0 : 1;
927 
928  for (unsigned int q = basis_size_1; q != 0; --q)
929  FEEvaluationImplBasisChange<
930  variant,
931  next_dim,
932  basis_size_1,
933  basis_size_2,
934  n_components,
935  Number,
936  Number2>::do_forward(transformation_matrix,
937  values_in +
938  (q - 1) *
939  Utilities::pow(basis_size_1, dim - 1),
940  my_scratch +
941  (q - 1) *
942  Utilities::pow(basis_size_2, dim - 1));
943  EvaluatorTensorProduct<variant,
944  dim,
945  basis_size_1,
946  basis_size_2,
947  Number,
948  Number2>
949  eval_val(transformation_matrix);
950  const unsigned int n_inner_blocks =
951  (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
952  const unsigned int n_blocks = Utilities::pow(basis_size_2, dim - 1);
953  for (unsigned int ii = 0; ii < n_blocks; ii += n_inner_blocks)
954  for (unsigned int c = 0; c < n_components; ++c)
955  {
956  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
957  eval_val.template values_one_line<dim - 1, true, false>(
958  my_scratch + i, my_scratch + i);
959  for (unsigned int q = 0; q < basis_size_2; ++q)
960  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
961  my_scratch[i + q * n_blocks + c * size_per_component] *=
962  coefficients[i + q * n_blocks +
963  c * stride * size_per_component];
964  for (unsigned int i = ii; i < ii + n_inner_blocks; ++i)
965  eval_val.template values_one_line<dim - 1, false, false>(
966  my_scratch + i, my_scratch + i);
967  }
968  for (unsigned int q = 0; q < basis_size_1; ++q)
969  FEEvaluationImplBasisChange<
970  variant,
971  next_dim,
972  basis_size_1,
973  basis_size_2,
974  n_components,
975  Number,
976  Number2>::do_backward(transformation_matrix,
977  false,
978  my_scratch +
979  q * Utilities::pow(basis_size_2, dim - 1),
980  values_out +
981  q * Utilities::pow(basis_size_1, dim - 1));
982  }
983  };
984 
985 
986 
1001  template <int dim, int fe_degree, int n_components, typename Number>
1002  struct FEEvaluationImplCollocation
1003  {
1004  static void
1005  evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1006  const Number * values_dofs,
1007  Number * values_quad,
1008  Number * gradients_quad,
1009  Number * hessians_quad,
1010  Number * scratch_data,
1011  const bool evaluate_values,
1012  const bool evaluate_gradients,
1013  const bool evaluate_hessians);
1014 
1015  static void
1016  integrate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1017  Number * values_dofs,
1018  Number * values_quad,
1019  Number * gradients_quad,
1020  Number * scratch_data,
1021  const bool integrate_values,
1022  const bool integrate_gradients,
1023  const bool add_into_values_array);
1024  };
1025 
1026 
1027 
1028  template <int dim, int fe_degree, int n_components, typename Number>
1029  inline void
1030  FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>::evaluate(
1031  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1032  const Number * values_dofs,
1033  Number * values_quad,
1034  Number * gradients_quad,
1035  Number * hessians_quad,
1036  Number *,
1037  const bool evaluate_values,
1038  const bool evaluate_gradients,
1039  const bool evaluate_hessians)
1040  {
1041  AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
1042  (fe_degree + 2) / 2 * (fe_degree + 1));
1043 
1044  EvaluatorTensorProduct<evaluate_evenodd,
1045  dim,
1046  fe_degree + 1,
1047  fe_degree + 1,
1048  Number>
1049  eval(AlignedVector<Number>(),
1050  shape_info.shape_gradients_collocation_eo,
1051  shape_info.shape_hessians_collocation_eo);
1052  constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
1053 
1054  for (unsigned int c = 0; c < n_components; c++)
1055  {
1056  if (evaluate_values == true)
1057  for (unsigned int i = 0; i < n_q_points; ++i)
1058  values_quad[i] = values_dofs[i];
1059  if (evaluate_gradients == true || evaluate_hessians == true)
1060  {
1061  eval.template gradients<0, true, false>(values_dofs,
1062  gradients_quad);
1063  if (dim > 1)
1064  eval.template gradients<1, true, false>(values_dofs,
1065  gradients_quad +
1066  n_q_points);
1067  if (dim > 2)
1068  eval.template gradients<2, true, false>(values_dofs,
1069  gradients_quad +
1070  2 * n_q_points);
1071  }
1072  if (evaluate_hessians == true)
1073  {
1074  eval.template hessians<0, true, false>(values_dofs, hessians_quad);
1075  if (dim > 1)
1076  {
1077  eval.template gradients<1, true, false>(gradients_quad,
1078  hessians_quad +
1079  dim * n_q_points);
1080  eval.template hessians<1, true, false>(values_dofs,
1081  hessians_quad +
1082  n_q_points);
1083  }
1084  if (dim > 2)
1085  {
1086  eval.template gradients<2, true, false>(gradients_quad,
1087  hessians_quad +
1088  4 * n_q_points);
1089  eval.template gradients<2, true, false>(
1090  gradients_quad + n_q_points, hessians_quad + 5 * n_q_points);
1091  eval.template hessians<2, true, false>(values_dofs,
1092  hessians_quad +
1093  2 * n_q_points);
1094  }
1095  hessians_quad += (dim * (dim + 1)) / 2 * n_q_points;
1096  }
1097  gradients_quad += dim * n_q_points;
1098  values_quad += n_q_points;
1099  values_dofs += n_q_points;
1100  }
1101  }
1102 
1103 
1104 
1105  template <int dim, int fe_degree, int n_components, typename Number>
1106  inline void
1107  FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>::integrate(
1108  const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1109  Number * values_dofs,
1110  Number * values_quad,
1111  Number * gradients_quad,
1112  Number *,
1113  const bool integrate_values,
1114  const bool integrate_gradients,
1115  const bool add_into_values_array)
1116  {
1117  AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
1118  (fe_degree + 2) / 2 * (fe_degree + 1));
1119 
1120  EvaluatorTensorProduct<evaluate_evenodd,
1121  dim,
1122  fe_degree + 1,
1123  fe_degree + 1,
1124  Number>
1125  eval(AlignedVector<Number>(),
1126  shape_info.shape_gradients_collocation_eo,
1127  shape_info.shape_hessians_collocation_eo);
1128  constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
1129 
1130  for (unsigned int c = 0; c < n_components; c++)
1131  {
1132  if (integrate_values == true && add_into_values_array == false)
1133  for (unsigned int i = 0; i < n_q_points; ++i)
1134  values_dofs[i] = values_quad[i];
1135  else if (integrate_values == true)
1136  for (unsigned int i = 0; i < n_q_points; ++i)
1137  values_dofs[i] += values_quad[i];
1138  if (integrate_gradients == true)
1139  {
1140  if (integrate_values == true || add_into_values_array == true)
1141  eval.template gradients<0, false, true>(gradients_quad,
1142  values_dofs);
1143  else
1144  eval.template gradients<0, false, false>(gradients_quad,
1145  values_dofs);
1146  if (dim > 1)
1147  eval.template gradients<1, false, true>(gradients_quad +
1148  n_q_points,
1149  values_dofs);
1150  if (dim > 2)
1151  eval.template gradients<2, false, true>(gradients_quad +
1152  2 * n_q_points,
1153  values_dofs);
1154  }
1155  gradients_quad += dim * n_q_points;
1156  values_quad += n_q_points;
1157  values_dofs += n_q_points;
1158  }
1159  }
1160 
1161 
1162 
1175  template <int dim,
1176  int fe_degree,
1177  int n_q_points_1d,
1178  int n_components,
1179  typename Number>
1180  struct FEEvaluationImplTransformToCollocation
1181  {
1182  static void
1183  evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1184  const Number * values_dofs,
1185  Number * values_quad,
1186  Number * gradients_quad,
1187  Number * hessians_quad,
1188  Number * scratch_data,
1189  const bool evaluate_values,
1190  const bool evaluate_gradients,
1191  const bool evaluate_hessians);
1192 
1193  static void
1194  integrate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1195  Number * values_dofs,
1196  Number * values_quad,
1197  Number * gradients_quad,
1198  Number * scratch_data,
1199  const bool integrate_values,
1200  const bool integrate_gradients,
1201  const bool add_into_values_array);
1202  };
1203 
1204 
1205 
1206  template <int dim,
1207  int fe_degree,
1208  int n_q_points_1d,
1209  int n_components,
1210  typename Number>
1211  inline void
1212  FEEvaluationImplTransformToCollocation<
1213  dim,
1214  fe_degree,
1215  n_q_points_1d,
1216  n_components,
1217  Number>::evaluate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1218  const Number * values_dofs,
1219  Number * values_quad,
1220  Number *gradients_quad,
1221  Number *hessians_quad,
1222  Number *,
1223  const bool,
1224  const bool evaluate_gradients,
1225  const bool evaluate_hessians)
1226  {
1227  Assert(n_q_points_1d > fe_degree,
1228  ExcMessage("You lose information when going to a collocation space "
1229  "of lower degree, so the evaluation results would be "
1230  "wrong. Thus, this class does not permit the desired "
1231  "operation."));
1232  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
1233 
1234  for (unsigned int c = 0; c < n_components; c++)
1235  {
1236  FEEvaluationImplBasisChange<
1237  evaluate_evenodd,
1238  dim,
1239  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1240  n_q_points_1d,
1241  1,
1242  Number,
1243  Number>::do_forward(shape_info.shape_values_eo,
1244  values_dofs,
1245  values_quad);
1246 
1247  // apply derivatives in the collocation space
1248  if (evaluate_gradients == true || evaluate_hessians == true)
1249  FEEvaluationImplCollocation<dim, n_q_points_1d - 1, 1, Number>::
1250  evaluate(shape_info,
1251  values_quad,
1252  nullptr,
1253  gradients_quad,
1254  hessians_quad,
1255  nullptr,
1256  false,
1257  evaluate_gradients,
1258  evaluate_hessians);
1259 
1260  values_dofs += shape_info.dofs_per_component_on_cell;
1261  values_quad += n_q_points;
1262  gradients_quad += dim * n_q_points;
1263  hessians_quad += (dim * (dim + 1)) / 2 * n_q_points;
1264  }
1265  }
1266 
1267 
1268 
1269  template <int dim,
1270  int fe_degree,
1271  int n_q_points_1d,
1272  int n_components,
1273  typename Number>
1274  inline void
1275  FEEvaluationImplTransformToCollocation<
1276  dim,
1277  fe_degree,
1278  n_q_points_1d,
1279  n_components,
1280  Number>::integrate(const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
1281  Number *values_dofs,
1282  Number *values_quad,
1283  Number *gradients_quad,
1284  Number *,
1285  const bool integrate_values,
1286  const bool integrate_gradients,
1287  const bool add_into_values_array)
1288  {
1289  Assert(n_q_points_1d > fe_degree,
1290  ExcMessage("You lose information when going to a collocation space "
1291  "of lower degree, so the evaluation results would be "
1292  "wrong. Thus, this class does not permit the desired "
1293  "operation."));
1294  AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
1295  (n_q_points_1d + 1) / 2 * n_q_points_1d);
1296  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
1297 
1298  for (unsigned int c = 0; c < n_components; c++)
1299  {
1300  // apply derivatives in collocation space
1301  if (integrate_gradients == true)
1302  FEEvaluationImplCollocation<dim, n_q_points_1d - 1, 1, Number>::
1303  integrate(shape_info,
1304  values_quad,
1305  nullptr,
1306  gradients_quad,
1307  nullptr,
1308  false,
1309  integrate_gradients,
1310  /*add_into_values_array=*/integrate_values);
1311 
1312  // transform back to the original space
1313  FEEvaluationImplBasisChange<
1314  evaluate_evenodd,
1315  dim,
1316  (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
1317  n_q_points_1d,
1318  1,
1319  Number,
1320  Number>::do_backward(shape_info.shape_values_eo,
1321  add_into_values_array,
1322  values_quad,
1323  values_dofs);
1324  gradients_quad += dim * n_q_points;
1325  values_quad += n_q_points;
1326  values_dofs += shape_info.dofs_per_component_on_cell;
1327  }
1328  }
1329 
1330 
1331 
1332  template <bool symmetric_evaluate,
1333  int dim,
1334  int fe_degree,
1335  int n_q_points_1d,
1336  int n_components,
1337  typename Number>
1338  struct FEFaceEvaluationImpl
1339  {
1340  static void
1341  evaluate_in_face(const MatrixFreeFunctions::ShapeInfo<Number> &data,
1342  Number * values_dofs,
1343  Number * values_quad,
1344  Number * gradients_quad,
1345  Number * scratch_data,
1346  const bool evaluate_val,
1347  const bool evaluate_grad,
1348  const unsigned int subface_index)
1349  {
1350  const AlignedVector<Number> &val1 =
1351  symmetric_evaluate ?
1352  data.shape_values_eo :
1353  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1354  data.shape_values :
1355  data.values_within_subface[subface_index % 2]);
1356  const AlignedVector<Number> &val2 =
1357  symmetric_evaluate ?
1358  data.shape_values_eo :
1359  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1360  data.shape_values :
1361  data.values_within_subface[subface_index / 2]);
1362 
1363  const AlignedVector<Number> &grad1 =
1364  symmetric_evaluate ?
1365  data.shape_gradients_eo :
1366  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1367  data.shape_gradients :
1368  data.gradients_within_subface[subface_index % 2]);
1369  const AlignedVector<Number> &grad2 =
1370  symmetric_evaluate ?
1371  data.shape_gradients_eo :
1372  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1373  data.shape_gradients :
1374  data.gradients_within_subface[subface_index / 2]);
1375 
1376  using Eval =
1377  internal::EvaluatorTensorProduct<symmetric_evaluate ?
1378  internal::evaluate_evenodd :
1379  internal::evaluate_general,
1380  dim - 1,
1381  fe_degree + 1,
1382  n_q_points_1d,
1383  Number>;
1384  Eval eval1(val1,
1385  grad1,
1386  AlignedVector<Number>(),
1387  data.fe_degree + 1,
1388  data.n_q_points_1d);
1389  Eval eval2(val2,
1390  grad2,
1391  AlignedVector<Number>(),
1392  data.fe_degree + 1,
1393  data.n_q_points_1d);
1394 
1395  const unsigned int size_deg =
1396  fe_degree > -1 ?
1397  Utilities::pow(fe_degree + 1, dim - 1) :
1398  (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
1399 
1400  const unsigned int n_q_points = fe_degree > -1 ?
1401  Utilities::pow(n_q_points_1d, dim - 1) :
1402  data.n_q_points_face;
1403 
1404  if (evaluate_grad == false)
1405  for (unsigned int c = 0; c < n_components; ++c)
1406  {
1407  switch (dim)
1408  {
1409  case 3:
1410  eval1.template values<0, true, false>(values_dofs,
1411  values_quad);
1412  eval2.template values<1, true, false>(values_quad,
1413  values_quad);
1414  break;
1415  case 2:
1416  eval1.template values<0, true, false>(values_dofs,
1417  values_quad);
1418  break;
1419  case 1:
1420  values_quad[0] = values_dofs[0];
1421  break;
1422  default:
1423  Assert(false, ExcNotImplemented());
1424  }
1425  values_dofs += 2 * size_deg;
1426  values_quad += n_q_points;
1427  }
1428  else
1429  for (unsigned int c = 0; c < n_components; ++c)
1430  {
1431  switch (dim)
1432  {
1433  case 3:
1434  if (symmetric_evaluate && n_q_points_1d > fe_degree)
1435  {
1436  eval1.template values<0, true, false>(values_dofs,
1437  values_quad);
1438  eval1.template values<1, true, false>(values_quad,
1439  values_quad);
1440  internal::EvaluatorTensorProduct<
1441  internal::evaluate_evenodd,
1442  dim - 1,
1443  n_q_points_1d,
1444  n_q_points_1d,
1445  Number>
1446  eval_grad(AlignedVector<Number>(),
1447  data.shape_gradients_collocation_eo,
1448  AlignedVector<Number>());
1449  eval_grad.template gradients<0, true, false>(
1450  values_quad, gradients_quad);
1451  eval_grad.template gradients<1, true, false>(
1452  values_quad, gradients_quad + n_q_points);
1453  }
1454  else
1455  {
1456  eval1.template gradients<0, true, false>(values_dofs,
1457  scratch_data);
1458  eval2.template values<1, true, false>(scratch_data,
1459  gradients_quad);
1460 
1461  eval1.template values<0, true, false>(values_dofs,
1462  scratch_data);
1463  eval2.template gradients<1, true, false>(scratch_data,
1464  gradients_quad +
1465  n_q_points);
1466 
1467  if (evaluate_val == true)
1468  eval2.template values<1, true, false>(scratch_data,
1469  values_quad);
1470  }
1471  eval1.template values<0, true, false>(values_dofs + size_deg,
1472  scratch_data);
1473  eval2.template values<1, true, false>(
1474  scratch_data, gradients_quad + (dim - 1) * n_q_points);
1475 
1476  break;
1477  case 2:
1478  eval1.template values<0, true, false>(values_dofs + size_deg,
1479  gradients_quad +
1480  (dim - 1) *
1481  n_q_points);
1482  eval1.template gradients<0, true, false>(values_dofs,
1483  gradients_quad);
1484  if (evaluate_val == true)
1485  eval1.template values<0, true, false>(values_dofs,
1486  values_quad);
1487  break;
1488  case 1:
1489  values_quad[0] = values_dofs[0];
1490  gradients_quad[0] = values_dofs[1];
1491  break;
1492  default:
1493  AssertThrow(false, ExcNotImplemented());
1494  }
1495  values_dofs += 2 * size_deg;
1496  values_quad += n_q_points;
1497  gradients_quad += dim * n_q_points;
1498  }
1499  }
1500 
1501  static void
1502  integrate_in_face(const MatrixFreeFunctions::ShapeInfo<Number> &data,
1503  Number * values_dofs,
1504  Number * values_quad,
1505  Number * gradients_quad,
1506  Number * scratch_data,
1507  const bool integrate_val,
1508  const bool integrate_grad,
1509  const unsigned int subface_index)
1510  {
1511  const AlignedVector<Number> &val1 =
1512  symmetric_evaluate ?
1513  data.shape_values_eo :
1514  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1515  data.shape_values :
1516  data.values_within_subface[subface_index % 2]);
1517  const AlignedVector<Number> &val2 =
1518  symmetric_evaluate ?
1519  data.shape_values_eo :
1520  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1521  data.shape_values :
1522  data.values_within_subface[subface_index / 2]);
1523 
1524  const AlignedVector<Number> &grad1 =
1525  symmetric_evaluate ?
1526  data.shape_gradients_eo :
1527  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1528  data.shape_gradients :
1529  data.gradients_within_subface[subface_index % 2]);
1530  const AlignedVector<Number> &grad2 =
1531  symmetric_evaluate ?
1532  data.shape_gradients_eo :
1533  (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
1534  data.shape_gradients :
1535  data.gradients_within_subface[subface_index / 2]);
1536 
1537  using Eval =
1538  internal::EvaluatorTensorProduct<symmetric_evaluate ?
1539  internal::evaluate_evenodd :
1540  internal::evaluate_general,
1541  dim - 1,
1542  fe_degree + 1,
1543  n_q_points_1d,
1544  Number>;
1545  Eval eval1(val1, grad1, val1, data.fe_degree + 1, data.n_q_points_1d);
1546  Eval eval2(val2, grad2, val1, data.fe_degree + 1, data.n_q_points_1d);
1547 
1548  const unsigned int size_deg =
1549  fe_degree > -1 ?
1550  Utilities::pow(fe_degree + 1, dim - 1) :
1551  (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
1552 
1553  const unsigned int n_q_points =
1554  fe_degree > -1 ?
1555  Utilities::fixed_int_power<n_q_points_1d, dim - 1>::value :
1556  data.n_q_points_face;
1557 
1558  if (integrate_grad == false)
1559  for (unsigned int c = 0; c < n_components; ++c)
1560  {
1561  switch (dim)
1562  {
1563  case 3:
1564  eval2.template values<1, false, false>(values_quad,
1565  values_quad);
1566  eval1.template values<0, false, false>(values_quad,
1567  values_dofs);
1568  break;
1569  case 2:
1570  eval1.template values<0, false, false>(values_quad,
1571  values_dofs);
1572  break;
1573  case 1:
1574  values_dofs[0] = values_quad[0];
1575  break;
1576  default:
1577  Assert(false, ExcNotImplemented());
1578  }
1579  values_dofs += 2 * size_deg;
1580  values_quad += n_q_points;
1581  }
1582  else
1583  for (unsigned int c = 0; c < n_components; ++c)
1584  {
1585  switch (dim)
1586  {
1587  case 3:
1588  eval2.template values<1, false, false>(gradients_quad +
1589  2 * n_q_points,
1590  gradients_quad +
1591  2 * n_q_points);
1592  eval1.template values<0, false, false>(
1593  gradients_quad + 2 * n_q_points, values_dofs + size_deg);
1594  if (symmetric_evaluate && n_q_points_1d > fe_degree)
1595  {
1596  internal::EvaluatorTensorProduct<
1597  internal::evaluate_evenodd,
1598  dim - 1,
1599  n_q_points_1d,
1600  n_q_points_1d,
1601  Number>
1602  eval_grad(AlignedVector<Number>(),
1603  data.shape_gradients_collocation_eo,
1604  AlignedVector<Number>());
1605  if (integrate_val)
1606  eval_grad.template gradients<1, false, true>(
1607  gradients_quad + n_q_points, values_quad);
1608  else
1609  eval_grad.template gradients<1, false, false>(
1610  gradients_quad + n_q_points, values_quad);
1611  eval_grad.template gradients<0, false, true>(
1612  gradients_quad, values_quad);
1613  eval1.template values<1, false, false>(values_quad,
1614  values_quad);
1615  eval1.template values<0, false, false>(values_quad,
1616  values_dofs);
1617  }
1618  else
1619  {
1620  if (integrate_val)
1621  {
1622  eval2.template values<1, false, false>(values_quad,
1623  scratch_data);
1624  eval2.template gradients<1, false, true>(
1625  gradients_quad + n_q_points, scratch_data);
1626  }
1627  else
1628  eval2.template gradients<1, false, false>(
1629  gradients_quad + n_q_points, scratch_data);
1630 
1631  eval1.template values<0, false, false>(scratch_data,
1632  values_dofs);
1633  eval2.template values<1, false, false>(gradients_quad,
1634  scratch_data);
1635  eval1.template gradients<0, false, true>(scratch_data,
1636  values_dofs);
1637  }
1638  break;
1639  case 2:
1640  eval1.template values<0, false, false>(
1641  gradients_quad + n_q_points, values_dofs + size_deg);
1642  eval1.template gradients<0, false, false>(gradients_quad,
1643  values_dofs);
1644  if (integrate_val == true)
1645  eval1.template values<0, false, true>(values_quad,
1646  values_dofs);
1647  break;
1648  case 1:
1649  values_dofs[0] = values_quad[0];
1650  values_dofs[1] = gradients_quad[0];
1651  break;
1652  default:
1653  AssertThrow(false, ExcNotImplemented());
1654  }
1655  values_dofs += 2 * size_deg;
1656  values_quad += n_q_points;
1657  gradients_quad += dim * n_q_points;
1658  }
1659  }
1660  };
1661 
1662 
1663 
1664  template <int dim, int fe_degree, int n_components, typename Number>
1665  struct FEFaceNormalEvaluationImpl
1666  {
1667  template <bool do_evaluate, bool add_into_output>
1668  static void
1669  interpolate(const MatrixFreeFunctions::ShapeInfo<Number> &data,
1670  const Number * input,
1671  Number * output,
1672  const bool do_gradients,
1673  const unsigned int face_no)
1674  {
1675  internal::EvaluatorTensorProduct<internal::evaluate_general,
1676  dim,
1677  fe_degree + 1,
1678  0,
1679  Number>
1680  evalf(data.shape_data_on_face[face_no % 2],
1681  AlignedVector<Number>(),
1682  AlignedVector<Number>(),
1683  data.fe_degree + 1,
1684  0);
1685 
1686  const unsigned int in_stride = do_evaluate ?
1687  data.dofs_per_component_on_cell :
1688  2 * data.dofs_per_component_on_face;
1689  const unsigned int out_stride = do_evaluate ?
1690  2 * data.dofs_per_component_on_face :
1691  data.dofs_per_component_on_cell;
1692  const unsigned int face_direction = face_no / 2;
1693  for (unsigned int c = 0; c < n_components; c++)
1694  {
1695  if (do_gradients)
1696  {
1697  if (face_direction == 0)
1698  evalf.template apply_face<0, do_evaluate, add_into_output, 1>(
1699  input, output);
1700  else if (face_direction == 1)
1701  evalf.template apply_face<1, do_evaluate, add_into_output, 1>(
1702  input, output);
1703  else
1704  evalf.template apply_face<2, do_evaluate, add_into_output, 1>(
1705  input, output);
1706  }
1707  else
1708  {
1709  if (face_direction == 0)
1710  evalf.template apply_face<0, do_evaluate, add_into_output, 0>(
1711  input, output);
1712  else if (face_direction == 1)
1713  evalf.template apply_face<1, do_evaluate, add_into_output, 0>(
1714  input, output);
1715  else
1716  evalf.template apply_face<2, do_evaluate, add_into_output, 0>(
1717  input, output);
1718  }
1719  input += in_stride;
1720  output += out_stride;
1721  }
1722  }
1723  };
1724 } // end of namespace internal
1725 
1726 
1727 DEAL_II_NAMESPACE_CLOSE
1728 
1729 #endif
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:173
Utilities::fixed_int_power
Definition: utilities.h:330
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1366
Utilities::pow
constexpr unsigned int pow(const unsigned int base, const unsigned int iexp)
Definition: utilities.h:353
internal::MatrixFreeFunctions::ShapeInfo::shape_values_eo
AlignedVector< Number > shape_values_eo
Definition: shape_info.h:164
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1329
internal::MatrixFreeFunctions::ShapeInfo::shape_hessians
AlignedVector< Number > shape_hessians
Definition: shape_info.h:157
internal::FEEvaluationImplBasisChange::do_forward
static void do_forward(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:686
internal::MatrixFreeFunctions::ShapeInfo::shape_values
AlignedVector< Number > shape_values
Definition: shape_info.h:139
internal::FEEvaluationImplBasisChange::do_mass
static void do_mass(const AlignedVector< Number2 > &transformation_matrix, const AlignedVector< Number > &coefficients, const Number *values_in, Number *scratch_data, Number *values_out)
Definition: evaluation_kernels.h:911
AlignedVector::size
size_type size() const
Utilities::fixed_power
T fixed_power(const T t)
Definition: utilities.h:912
StandardExceptions::ExcMessage
static::ExceptionBase & ExcMessage(std::string arg1)
internal::FEEvaluationImplBasisChange::do_backward
static void do_backward(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:805
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1227
internal::MatrixFreeFunctions::ShapeInfo::shape_gradients_collocation_eo
AlignedVector< Number > shape_gradients_collocation_eo
Definition: shape_info.h:185
StandardExceptions::ExcDimensionMismatch
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
internal::MatrixFreeFunctions::ShapeInfo::shape_data_on_face
AlignedVector< Number > shape_data_on_face[2]
Definition: shape_info.h:199
internal::MatrixFreeFunctions::ShapeInfo::shape_hessians_eo
AlignedVector< Number > shape_hessians_eo
Definition: shape_info.h:178
internal::MatrixFreeFunctions::ShapeInfo::gradients_within_subface
AlignedVector< Number > gradients_within_subface[2]
Definition: shape_info.h:211
internal::MatrixFreeFunctions::ShapeInfo::shape_gradients_eo
AlignedVector< Number > shape_gradients_eo
Definition: shape_info.h:171
internal::MatrixFreeFunctions::ShapeInfo::shape_gradients
AlignedVector< Number > shape_gradients
Definition: shape_info.h:148
StandardExceptions::ExcNotImplemented
static::ExceptionBase & ExcNotImplemented()
internal::MatrixFreeFunctions::ShapeInfo::values_within_subface
AlignedVector< Number > values_within_subface[2]
Definition: shape_info.h:205
internal::MatrixFreeFunctions::ShapeInfo::shape_hessians_collocation_eo
AlignedVector< Number > shape_hessians_collocation_eo
Definition: shape_info.h:192
DoFTools::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
Generated by   doxygen 1.8.11