Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
evaluation_selector.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2019 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_selector_h
18 #define dealii_matrix_free_evaluation_selector_h
19 
20 #include <deal.II/matrix_free/evaluation_kernels.h>
21 
22 DEAL_II_NAMESPACE_OPEN
23 
24 #ifndef DOXYGEN
25 namespace internal
26 {
27  namespace EvaluationSelectorImplementation
28  {
29  // The following classes serve the purpose of choosing the correct template
30  // specialization of the FEEvaluationImpl* classes in case fe_degree
31  // and n_q_points_1d are only given as runtime parameters.
32  // The logic is the following:
33  // 1. Start with fe_degree=0, n_q_points_1d=0 and DEPTH=0.
34  // 2. If the current assumption on fe_degree doesn't match the runtime
35  // parameter, increase fe_degree by one and try again.
36  // If fe_degree==10 use the class Default which serves as a fallback.
37  // 3. After fixing the fe_degree, DEPTH is increased (DEPTH=1) and we start
38  // with
39  // n_q_points=fe_degree+1.
40  // 4. If the current assumption on n_q_points_1d doesn't match the runtime
41  // parameter, increase n_q_points_1d by one and try again.
42  // If n_q_points_1d==degree+3 use the class Default which serves as a
43  // fallback.
44 
49  template <int dim, int n_components, typename Number>
50  struct Default
51  {
52  static inline void
53  evaluate(
55  Number * values_dofs_actual,
56  Number * values_quad,
57  Number * gradients_quad,
58  Number * hessians_quad,
59  Number * scratch_data,
60  const bool evaluate_values,
61  const bool evaluate_gradients,
62  const bool evaluate_hessians)
63  {
66  dim,
67  -1,
68  0,
70  Number>::evaluate(shape_info,
71  values_dofs_actual,
72  values_quad,
73  gradients_quad,
74  hessians_quad,
75  scratch_data,
76  evaluate_values,
77  evaluate_gradients,
78  evaluate_hessians);
79  }
80 
81  static inline void
82  integrate(
84  Number * values_dofs_actual,
85  Number * values_quad,
86  Number * gradients_quad,
87  Number * scratch_data,
88  const bool integrate_values,
89  const bool integrate_gradients,
90  const bool sum_into_values_array = false)
91  {
94  dim,
95  -1,
96  0,
98  Number>::integrate(shape_info,
99  values_dofs_actual,
100  values_quad,
101  gradients_quad,
102  scratch_data,
103  integrate_values,
104  integrate_gradients,
105  sum_into_values_array);
106  }
107  };
108 
109 
113  template <int dim,
114  int n_components,
115  typename Number,
116  int DEPTH = 0,
117  int degree = 0,
118  int n_q_points_1d = 0,
119  class Enable = void>
120  struct Factory : Default<dim, n_components, Number>
121  {};
122 
128  template <int n_q_points_1d, int dim, int n_components, typename Number>
129  struct Factory<dim, n_components, Number, 0, 10, n_q_points_1d>
130  : Default<dim, n_components, Number>
131  {};
132 
138  template <int degree,
139  int n_q_points_1d,
140  int dim,
141  int n_components,
142  typename Number>
143  struct Factory<dim,
144  n_components,
145  Number,
146  1,
147  degree,
148  n_q_points_1d,
149  typename std::enable_if<n_q_points_1d == degree + 3>::type>
150  : Default<dim, n_components, Number>
151  {};
152 
156  template <int degree,
157  int n_q_points_1d,
158  int dim,
159  int n_components,
160  typename Number>
161  struct Factory<dim, n_components, Number, 0, degree, n_q_points_1d>
162  {
163  static inline void
164  evaluate(
166  Number * values_dofs_actual,
167  Number * values_quad,
168  Number * gradients_quad,
169  Number * hessians_quad,
170  Number * scratch_data,
171  const bool evaluate_values,
172  const bool evaluate_gradients,
173  const bool evaluate_hessians)
174  {
175  const unsigned int runtime_degree = shape_info.fe_degree;
176  constexpr unsigned int start_n_q_points = degree + 1;
177  if (runtime_degree == degree)
178  Factory<dim, n_components, Number, 1, degree, start_n_q_points>::
179  evaluate(shape_info,
180  values_dofs_actual,
181  values_quad,
182  gradients_quad,
183  hessians_quad,
184  scratch_data,
185  evaluate_values,
186  evaluate_gradients,
187  evaluate_hessians);
188  else
189  Factory<dim, n_components, Number, 0, degree + 1, n_q_points_1d>::
190  evaluate(shape_info,
191  values_dofs_actual,
192  values_quad,
193  gradients_quad,
194  hessians_quad,
195  scratch_data,
196  evaluate_values,
197  evaluate_gradients,
198  evaluate_hessians);
199  }
200 
201  static inline void
202  integrate(
204  Number * values_dofs_actual,
205  Number * values_quad,
206  Number * gradients_quad,
207  Number * scratch_data,
208  const bool integrate_values,
209  const bool integrate_gradients,
210  const bool sum_into_values_array = false)
211  {
212  const int runtime_degree = shape_info.fe_degree;
213  constexpr unsigned int start_n_q_points = degree + 1;
214  if (runtime_degree == degree)
215  Factory<dim, n_components, Number, 1, degree, start_n_q_points>::
216  integrate(shape_info,
217  values_dofs_actual,
218  values_quad,
219  gradients_quad,
220  scratch_data,
221  integrate_values,
222  integrate_gradients,
223  sum_into_values_array);
224  else
225  Factory<dim, n_components, Number, 0, degree + 1, n_q_points_1d>::
226  integrate(shape_info,
227  values_dofs_actual,
228  values_quad,
229  gradients_quad,
230  scratch_data,
231  integrate_values,
232  integrate_gradients,
233  sum_into_values_array);
234  }
235  };
236 
241  template <int degree,
242  int n_q_points_1d,
243  int dim,
244  int n_components,
245  typename Number>
246  struct Factory<dim,
247  n_components,
248  Number,
249  1,
250  degree,
251  n_q_points_1d,
252  typename std::enable_if<(n_q_points_1d < degree + 3)>::type>
253  {
261  static constexpr bool use_collocation =
262  n_q_points_1d > degree &&n_q_points_1d <= 3 * degree / 2 + 1 &&
263  n_q_points_1d < 200;
264 
265  static inline void
266  evaluate(
268  Number * values_dofs_actual,
269  Number * values_quad,
270  Number * gradients_quad,
271  Number * hessians_quad,
272  Number * scratch_data,
273  const bool evaluate_values,
274  const bool evaluate_gradients,
275  const bool evaluate_hessians)
276  {
277  const int runtime_n_q_points_1d = shape_info.n_q_points_1d;
278  if (runtime_n_q_points_1d == n_q_points_1d)
279  {
280  if (n_q_points_1d == degree + 1 &&
281  shape_info.element_type ==
283  internal::
284  FEEvaluationImplCollocation<dim, degree, n_components, Number>::
285  evaluate(shape_info,
286  values_dofs_actual,
287  values_quad,
288  gradients_quad,
289  hessians_quad,
290  scratch_data,
291  evaluate_values,
292  evaluate_gradients,
293  evaluate_hessians);
294  else if (use_collocation)
296  dim,
297  degree,
298  n_q_points_1d,
299  n_components,
300  Number>::evaluate(shape_info,
301  values_dofs_actual,
302  values_quad,
303  gradients_quad,
304  hessians_quad,
305  scratch_data,
306  evaluate_values,
307  evaluate_gradients,
308  evaluate_hessians);
309  else
312  dim,
313  degree,
314  n_q_points_1d,
315  n_components,
316  Number>::evaluate(shape_info,
317  values_dofs_actual,
318  values_quad,
319  gradients_quad,
320  hessians_quad,
321  scratch_data,
322  evaluate_values,
323  evaluate_gradients,
324  evaluate_hessians);
325  }
326  else
327  Factory<dim, n_components, Number, 1, degree, n_q_points_1d + 1>::
328  evaluate(shape_info,
329  values_dofs_actual,
330  values_quad,
331  gradients_quad,
332  hessians_quad,
333  scratch_data,
334  evaluate_values,
335  evaluate_gradients,
336  evaluate_hessians);
337  }
338 
339  static inline void
340  integrate(
342  Number * values_dofs_actual,
343  Number * values_quad,
344  Number * gradients_quad,
345  Number * scratch_data,
346  const bool integrate_values,
347  const bool integrate_gradients,
348  const bool sum_into_values_array)
349  {
350  const int runtime_n_q_points_1d = shape_info.n_q_points_1d;
351  if (runtime_n_q_points_1d == n_q_points_1d)
352  {
353  if (n_q_points_1d == degree + 1 &&
354  shape_info.element_type ==
356  internal::
357  FEEvaluationImplCollocation<dim, degree, n_components, Number>::
358  integrate(shape_info,
359  values_dofs_actual,
360  values_quad,
361  gradients_quad,
362  scratch_data,
363  integrate_values,
364  integrate_gradients,
365  sum_into_values_array);
366  else if (use_collocation)
368  dim,
369  degree,
370  n_q_points_1d,
371  n_components,
372  Number>::integrate(shape_info,
373  values_dofs_actual,
374  values_quad,
375  gradients_quad,
376  scratch_data,
377  integrate_values,
378  integrate_gradients,
379  sum_into_values_array);
380  else
383  dim,
384  degree,
385  n_q_points_1d,
386  n_components,
387  Number>::integrate(shape_info,
388  values_dofs_actual,
389  values_quad,
390  gradients_quad,
391  scratch_data,
392  integrate_values,
393  integrate_gradients,
394  sum_into_values_array);
395  }
396  else
397  Factory<dim, n_components, Number, 1, degree, n_q_points_1d + 1>::
398  integrate(shape_info,
399  values_dofs_actual,
400  values_quad,
401  gradients_quad,
402  scratch_data,
403  integrate_values,
404  integrate_gradients,
405  sum_into_values_array);
406  }
407  };
408 
409 
410 
415  template <int dim, int n_components, typename Number>
416  void
417  symmetric_selector_evaluate(
419  Number * values_dofs_actual,
420  Number * values_quad,
421  Number * gradients_quad,
422  Number * hessians_quad,
423  Number * scratch_data,
424  const bool evaluate_values,
425  const bool evaluate_gradients,
426  const bool evaluate_hessians)
427  {
428  Assert(shape_info.element_type <=
430  ExcInternalError());
431  Factory<dim, n_components, Number>::evaluate(shape_info,
432  values_dofs_actual,
433  values_quad,
434  gradients_quad,
435  hessians_quad,
436  scratch_data,
437  evaluate_values,
438  evaluate_gradients,
439  evaluate_hessians);
440  }
441 
442 
443 
448  template <int dim, int n_components, typename Number>
449  void
450  symmetric_selector_integrate(
452  Number * values_dofs_actual,
453  Number * values_quad,
454  Number * gradients_quad,
455  Number * scratch_data,
456  const bool integrate_values,
457  const bool integrate_gradients,
458  const bool sum_into_values_array = false)
459  {
460  Assert(shape_info.element_type <=
462  ExcInternalError());
463  Factory<dim, n_components, Number>::integrate(shape_info,
464  values_dofs_actual,
465  values_quad,
466  gradients_quad,
467  scratch_data,
468  integrate_values,
469  integrate_gradients,
470  sum_into_values_array);
471  }
472  } // namespace EvaluationSelectorImplementation
473 } // namespace internal
474 #endif
475 
476 
488 template <int dim,
489  int fe_degree,
490  int n_q_points_1d,
491  int n_components,
492  typename Number>
494 {
502  static constexpr bool use_collocation =
503  n_q_points_1d > fe_degree &&n_q_points_1d <= 3 * fe_degree / 2 + 1 &&
504  n_q_points_1d < 200;
505 
513  static void
515  Number * values_dofs_actual,
516  Number * values_quad,
517  Number * gradients_quad,
518  Number * hessians_quad,
519  Number * scratch_data,
520  const bool evaluate_values,
521  const bool evaluate_gradients,
522  const bool evaluate_hessians);
523 
531  static void
533  Number * values_dofs_actual,
534  Number * values_quad,
535  Number * gradients_quad,
536  Number * scratch_data,
537  const bool integrate_values,
538  const bool integrate_gradients,
539  const bool sum_into_values_array = false);
540 };
541 
552 template <int dim, int n_q_points_1d, int n_components, typename Number>
553 struct SelectEvaluator<dim, -1, n_q_points_1d, n_components, Number>
554 {
563  static void
565  Number * values_dofs_actual,
566  Number * values_quad,
567  Number * gradients_quad,
568  Number * hessians_quad,
569  Number * scratch_data,
570  const bool evaluate_values,
571  const bool evaluate_gradients,
572  const bool evaluate_hessians);
573 
582  static void
584  Number * values_dofs_actual,
585  Number * values_quad,
586  Number * gradients_quad,
587  Number * scratch_data,
588  const bool integrate_values,
589  const bool integrate_gradients,
590  const bool sum_into_values_array = false);
591 };
592 
593 //----------------------Implementation for SelectEvaluator---------------------
594 #ifndef DOXYGEN
595 
596 template <int dim,
597  int fe_degree,
598  int n_q_points_1d,
599  int n_components,
600  typename Number>
601 inline void
604  Number * values_dofs_actual,
605  Number * values_quad,
606  Number * gradients_quad,
607  Number * hessians_quad,
608  Number * scratch_data,
609  const bool evaluate_values,
610  const bool evaluate_gradients,
611  const bool evaluate_hessians)
612 {
613  Assert(fe_degree >= 0 && n_q_points_1d > 0, ExcInternalError());
614 
615  if (fe_degree + 1 == n_q_points_1d &&
616  shape_info.element_type ==
618  {
619  internal::
620  FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>::
621  evaluate(shape_info,
622  values_dofs_actual,
623  values_quad,
624  gradients_quad,
625  hessians_quad,
626  scratch_data,
627  evaluate_values,
628  evaluate_gradients,
629  evaluate_hessians);
630  }
631  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
632  // shape_info.h for more details
633  else if (use_collocation && shape_info.element_type <=
635  {
637  dim,
638  fe_degree,
639  n_q_points_1d,
640  n_components,
641  Number>::evaluate(shape_info,
642  values_dofs_actual,
643  values_quad,
644  gradients_quad,
645  hessians_quad,
646  scratch_data,
647  evaluate_values,
648  evaluate_gradients,
649  evaluate_hessians);
650  }
651  else if (shape_info.element_type <=
653  {
656  dim,
657  fe_degree,
658  n_q_points_1d,
659  n_components,
660  Number>::evaluate(shape_info,
661  values_dofs_actual,
662  values_quad,
663  gradients_quad,
664  hessians_quad,
665  scratch_data,
666  evaluate_values,
667  evaluate_gradients,
668  evaluate_hessians);
669  }
670  else if (shape_info.element_type ==
672  {
675  dim,
676  fe_degree,
677  n_q_points_1d,
678  n_components,
679  Number>::evaluate(shape_info,
680  values_dofs_actual,
681  values_quad,
682  gradients_quad,
683  hessians_quad,
684  scratch_data,
685  evaluate_values,
686  evaluate_gradients,
687  evaluate_hessians);
688  }
689  else if (shape_info.element_type ==
691  {
694  dim,
695  fe_degree,
696  n_q_points_1d,
697  n_components,
698  Number>::evaluate(shape_info,
699  values_dofs_actual,
700  values_quad,
701  gradients_quad,
702  hessians_quad,
703  scratch_data,
704  evaluate_values,
705  evaluate_gradients,
706  evaluate_hessians);
707  }
708  else if (shape_info.element_type ==
710  {
712  dim,
713  fe_degree,
714  n_q_points_1d,
715  n_components,
716  Number>::evaluate(shape_info,
717  values_dofs_actual,
718  values_quad,
719  gradients_quad,
720  hessians_quad,
721  scratch_data,
722  evaluate_values,
723  evaluate_gradients,
724  evaluate_hessians);
725  }
726  else
727  AssertThrow(false, ExcNotImplemented());
728 }
729 
730 
731 
732 template <int dim,
733  int fe_degree,
734  int n_q_points_1d,
735  int n_components,
736  typename Number>
737 inline void
740  Number * values_dofs_actual,
741  Number * values_quad,
742  Number * gradients_quad,
743  Number * scratch_data,
744  const bool integrate_values,
745  const bool integrate_gradients,
746  const bool sum_into_values_array)
747 {
748  Assert(fe_degree >= 0 && n_q_points_1d > 0, ExcInternalError());
749 
750  if (fe_degree + 1 == n_q_points_1d &&
751  shape_info.element_type ==
753  {
754  internal::
755  FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>::
756  integrate(shape_info,
757  values_dofs_actual,
758  values_quad,
759  gradients_quad,
760  scratch_data,
761  integrate_values,
762  integrate_gradients,
763  sum_into_values_array);
764  }
765  // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
766  // shape_info.h for more details
767  else if (use_collocation && shape_info.element_type <=
769  {
771  dim,
772  fe_degree,
773  n_q_points_1d,
774  n_components,
775  Number>::integrate(shape_info,
776  values_dofs_actual,
777  values_quad,
778  gradients_quad,
779  scratch_data,
780  integrate_values,
781  integrate_gradients,
782  sum_into_values_array);
783  }
784  else if (shape_info.element_type <=
786  {
789  dim,
790  fe_degree,
791  n_q_points_1d,
792  n_components,
793  Number>::integrate(shape_info,
794  values_dofs_actual,
795  values_quad,
796  gradients_quad,
797  scratch_data,
798  integrate_values,
799  integrate_gradients,
800  sum_into_values_array);
801  }
802  else if (shape_info.element_type ==
804  {
807  dim,
808  fe_degree,
809  n_q_points_1d,
810  n_components,
811  Number>::integrate(shape_info,
812  values_dofs_actual,
813  values_quad,
814  gradients_quad,
815  scratch_data,
816  integrate_values,
817  integrate_gradients,
818  sum_into_values_array);
819  }
820  else if (shape_info.element_type ==
822  {
825  dim,
826  fe_degree,
827  n_q_points_1d,
828  n_components,
829  Number>::integrate(shape_info,
830  values_dofs_actual,
831  values_quad,
832  gradients_quad,
833  scratch_data,
834  integrate_values,
835  integrate_gradients,
836  sum_into_values_array);
837  }
838  else if (shape_info.element_type ==
840  {
842  dim,
843  fe_degree,
844  n_q_points_1d,
845  n_components,
846  Number>::integrate(shape_info,
847  values_dofs_actual,
848  values_quad,
849  gradients_quad,
850  scratch_data,
851  integrate_values,
852  integrate_gradients,
853  sum_into_values_array);
854  }
855  else
856  AssertThrow(false, ExcNotImplemented());
857 }
858 
859 
860 
861 template <int dim, int dummy, int n_components, typename Number>
862 inline void
865  Number * values_dofs_actual,
866  Number * values_quad,
867  Number * gradients_quad,
868  Number * hessians_quad,
869  Number * scratch_data,
870  const bool evaluate_values,
871  const bool evaluate_gradients,
872  const bool evaluate_hessians)
873 {
874  if (shape_info.element_type ==
876  {
879  dim,
880  -1,
881  0,
882  n_components,
883  Number>::evaluate(shape_info,
884  values_dofs_actual,
885  values_quad,
886  gradients_quad,
887  hessians_quad,
888  scratch_data,
889  evaluate_values,
890  evaluate_gradients,
891  evaluate_hessians);
892  }
893  else if (shape_info.element_type ==
895  {
898  dim,
899  -1,
900  0,
901  n_components,
902  Number>::evaluate(shape_info,
903  values_dofs_actual,
904  values_quad,
905  gradients_quad,
906  hessians_quad,
907  scratch_data,
908  evaluate_values,
909  evaluate_gradients,
910  evaluate_hessians);
911  }
912  else if (shape_info.element_type ==
915  dim,
916  -1,
917  0,
918  n_components,
919  Number>::evaluate(shape_info,
920  values_dofs_actual,
921  values_quad,
922  gradients_quad,
923  hessians_quad,
924  scratch_data,
925  evaluate_values,
926  evaluate_gradients,
927  evaluate_hessians);
928  else
929  internal::EvaluationSelectorImplementation::
930  symmetric_selector_evaluate<dim, n_components, Number>(shape_info,
931  values_dofs_actual,
932  values_quad,
933  gradients_quad,
934  hessians_quad,
935  scratch_data,
936  evaluate_values,
937  evaluate_gradients,
938  evaluate_hessians);
939 }
940 
941 
942 
943 template <int dim, int dummy, int n_components, typename Number>
944 inline void
947  Number * values_dofs_actual,
948  Number * values_quad,
949  Number * gradients_quad,
950  Number * scratch_data,
951  const bool integrate_values,
952  const bool integrate_gradients,
953  const bool sum_into_values_array)
954 {
955  if (shape_info.element_type ==
957  {
960  dim,
961  -1,
962  0,
963  n_components,
964  Number>::integrate(shape_info,
965  values_dofs_actual,
966  values_quad,
967  gradients_quad,
968  scratch_data,
969  integrate_values,
970  integrate_gradients,
971  sum_into_values_array);
972  }
973  else if (shape_info.element_type ==
975  {
978  dim,
979  -1,
980  0,
981  n_components,
982  Number>::integrate(shape_info,
983  values_dofs_actual,
984  values_quad,
985  gradients_quad,
986  scratch_data,
987  integrate_values,
988  integrate_gradients,
989  sum_into_values_array);
990  }
991  else if (shape_info.element_type ==
994  dim,
995  -1,
996  0,
997  n_components,
998  Number>::integrate(shape_info,
999  values_dofs_actual,
1000  values_quad,
1001  gradients_quad,
1002  scratch_data,
1003  integrate_values,
1004  integrate_gradients,
1005  sum_into_values_array);
1006  else
1007  internal::EvaluationSelectorImplementation::
1008  symmetric_selector_integrate<dim, n_components, Number>(
1009  shape_info,
1010  values_dofs_actual,
1011  values_quad,
1012  gradients_quad,
1013  scratch_data,
1014  integrate_values,
1015  integrate_gradients,
1016  sum_into_values_array);
1017 }
1018 #endif // DOXYGEN
1019 
1020 DEAL_II_NAMESPACE_CLOSE
1021 
1022 #endif
static void integrate(const internal::MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs_actual, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool integrate_values, const bool integrate_gradients, const bool sum_into_values_array=false)
static constexpr bool use_collocation
STL namespace.
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
static void evaluate(const internal::MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs_actual, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians)
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcNotImplemented()
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
static ::ExceptionBase & ExcInternalError()