17 #ifndef dealii__matrix_free_evaluation_kernels_h 18 #define dealii__matrix_free_evaluation_kernels_h 20 #include <deal.II/base/config.h> 21 #include <deal.II/base/vectorization.h> 22 #include <deal.II/base/utilities.h> 23 #include <deal.II/matrix_free/tensor_product_kernels.h> 24 #include <deal.II/matrix_free/shape_info.h> 27 DEAL_II_NAMESPACE_OPEN
34 template <MatrixFreeFunctions::ElementType element,
bool is_
long>
35 struct EvaluatorSelector {};
37 template <
bool is_
long>
38 struct EvaluatorSelector<MatrixFreeFunctions::tensor_general,is_long>
44 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric,false>
49 template <>
struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric,true>
54 template <
bool is_
long>
55 struct EvaluatorSelector<MatrixFreeFunctions::truncated_tensor,is_long>
61 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,false>
67 struct EvaluatorSelector<MatrixFreeFunctions::tensor_symmetric_plus_dg0,true>
72 template <
bool is_
long>
73 struct EvaluatorSelector<MatrixFreeFunctions::tensor_gausslobatto,is_long>
91 template <MatrixFreeFunctions::ElementType type,
int dim,
int fe_degree,
92 int n_q_points_1d,
int n_components,
typename Number>
93 struct FEEvaluationImpl
96 void evaluate (
const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
102 const bool evaluate_val,
103 const bool evaluate_grad,
104 const bool evaluate_lapl);
107 void integrate (
const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
112 const bool evaluate_val,
113 const bool evaluate_grad);
117 template <MatrixFreeFunctions::ElementType type,
int dim,
int fe_degree,
118 int n_q_points_1d,
int n_components,
typename Number>
121 FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
122 ::evaluate (
const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
128 const bool evaluate_val,
129 const bool evaluate_grad,
130 const bool evaluate_lapl)
132 if (evaluate_val ==
false && evaluate_grad ==
false && evaluate_lapl ==
false)
136 EvaluatorSelector<type,(fe_degree+n_q_points_1d>4)>::variant;
137 typedef EvaluatorTensorProduct<variant, dim, fe_degree, n_q_points_1d,
140 shape_info.shape_values,
142 shape_info.shape_gradients,
144 shape_info.shape_hessians,
145 shape_info.fe_degree,
146 shape_info.n_q_points_1d);
149 : (Eval::dofs_per_cell > Eval::n_q_points ?
150 Eval::dofs_per_cell : Eval::n_q_points);
151 const unsigned int max_stack_size = 100;
157 temp1 = scratch_data;
158 temp2 = temp1 + std::max(Utilities::fixed_power<dim>(shape_info.fe_degree+1),
159 Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
161 else if (temp_size < max_stack_size)
163 temp1 = &temp_data[0];
164 temp2 = temp1 + temp_size;
168 temp1 = scratch_data;
169 temp2 = temp1 + temp_size;
174 if (type == MatrixFreeFunctions::truncated_tensor)
176 values_dofs = expanded_dof_values;
177 for (
unsigned int c=0; c<n_components; ++c)
178 expanded_dof_values[c] = scratch_data+2*(std::max(shape_info.dofs_per_cell,
179 shape_info.n_q_points)) +
180 c*Utilities::fixed_power<dim>(shape_info.fe_degree+1);
181 const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
182 unsigned int count_p = 0, count_q = 0;
183 for (
int i=0; i<(dim>2?degree+1:1); ++i)
185 for (
int j=0; j<(dim>1?degree+1-i:1); ++j)
187 for (
int k=0; k<degree+1-j-i; ++k, ++count_p, ++count_q)
188 for (
unsigned int c=0; c<n_components; ++c)
189 expanded_dof_values[c][count_q] = values_dofs_actual[c][count_p];
190 for (
int k=degree+1-j-i; k<degree+1; ++k, ++count_q)
191 for (
unsigned int c=0; c<n_components; ++c)
194 for (
int j=degree+1-i; j<degree+1; ++j)
195 for (
int k=0; k<degree+1; ++k, ++count_q)
196 for (
unsigned int c=0; c<n_components; ++c)
199 AssertDimension(count_q, Utilities::fixed_power<dim>(shape_info.fe_degree+1));
205 const unsigned int d1 = dim>1?1:0;
206 const unsigned int d2 = dim>2?2:0;
207 const unsigned int d3 = dim>2?3:0;
208 const unsigned int d4 = dim>2?4:0;
209 const unsigned int d5 = dim>2?5:0;
214 for (
unsigned int c=0; c<n_components; c++)
216 if (evaluate_val ==
true)
217 eval.template values<0,true,false> (values_dofs[c], values_quad[c]);
218 if (evaluate_grad ==
true)
219 eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
220 if (evaluate_lapl ==
true)
221 eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]);
226 for (
unsigned int c=0; c<n_components; c++)
229 if (evaluate_grad ==
true)
231 eval.template gradients<0,true,false> (values_dofs[c], temp1);
232 eval.template values<1,true,false> (temp1, gradients_quad[c][0]);
234 if (evaluate_lapl ==
true)
237 if (evaluate_grad ==
false)
238 eval.template gradients<0,true,false>(values_dofs[c], temp1);
239 eval.template gradients<1,true,false> (temp1, hessians_quad[c][d1+d1]);
242 eval.template hessians<0,true,false>(values_dofs[c], temp1);
243 eval.template values<1,true,false> (temp1, hessians_quad[c][0]);
247 eval.template values<0,true,false> (values_dofs[c], temp1);
248 if (evaluate_grad ==
true)
249 eval.template gradients<1,true,false> (temp1, gradients_quad[c][d1]);
252 if (evaluate_lapl ==
true)
253 eval.template hessians<1,true,false> (temp1, hessians_quad[c][d1]);
256 if (evaluate_val ==
true)
257 eval.template values<1,true,false> (temp1, values_quad[c]);
262 for (
unsigned int c=0; c<n_components; c++)
264 if (evaluate_grad ==
true)
267 eval.template gradients<0,true,false> (values_dofs[c], temp1);
268 eval.template values<1,true,false> (temp1, temp2);
269 eval.template values<2,true,false> (temp2, gradients_quad[c][0]);
272 if (evaluate_lapl ==
true)
275 if (evaluate_grad ==
false)
277 eval.template gradients<0,true,false> (values_dofs[c], temp1);
278 eval.template values<1,true,false> (temp1, temp2);
280 eval.template gradients<2,true,false> (temp2, hessians_quad[c][d4]);
283 eval.template gradients<1,true,false> (temp1, temp2);
284 eval.template values<2,true,false> (temp2, hessians_quad[c][d3]);
287 eval.template hessians<0,true,false>(values_dofs[c], temp1);
288 eval.template values<1,true,false> (temp1, temp2);
289 eval.template values<2,true,false> (temp2, hessians_quad[c][0]);
293 eval.template values<0,true,false> (values_dofs[c], temp1);
294 if (evaluate_grad ==
true)
296 eval.template gradients<1,true,false>(temp1, temp2);
297 eval.template values<2,true,false> (temp2, gradients_quad[c][d1]);
300 if (evaluate_lapl ==
true)
303 if (evaluate_grad ==
false)
304 eval.template gradients<1,true,false>(temp1, temp2);
305 eval.template gradients<2,true,false> (temp2, hessians_quad[c][d5]);
308 eval.template hessians<1,true,false> (temp1, temp2);
309 eval.template values<2,true,false> (temp2, hessians_quad[c][d1]);
313 eval.template values<1,true,false> (temp1, temp2);
314 if (evaluate_grad ==
true)
315 eval.template gradients<2,true,false> (temp2, gradients_quad[c][d2]);
319 if (evaluate_lapl ==
true)
320 eval.template hessians<2,true,false>(temp2, hessians_quad[c][d2]);
323 if (evaluate_val ==
true)
324 eval.template values<2,true,false> (temp2, values_quad[c]);
334 if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 && evaluate_val)
335 for (
unsigned int c=0; c<n_components; ++c)
336 for (
unsigned int q=0; q<shape_info.n_q_points; ++q)
337 values_quad[c][q] += values_dofs[c][shape_info.dofs_per_cell-1];
342 template <MatrixFreeFunctions::ElementType type,
int dim,
int fe_degree,
343 int n_q_points_1d,
int n_components,
typename Number>
346 FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
347 ::integrate (
const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
352 const bool integrate_val,
353 const bool integrate_grad)
356 EvaluatorSelector<type,(fe_degree+n_q_points_1d>4)>::variant;
357 typedef EvaluatorTensorProduct<variant, dim, fe_degree, n_q_points_1d,
360 shape_info.shape_values,
362 shape_info.shape_gradients,
364 shape_info.shape_hessians,
365 shape_info.fe_degree,
366 shape_info.n_q_points_1d);
369 : (Eval::dofs_per_cell > Eval::n_q_points ?
370 Eval::dofs_per_cell : Eval::n_q_points);
371 const unsigned int max_stack_size = 100;
377 temp1 = scratch_data;
378 temp2 = temp1 + std::max(Utilities::fixed_power<dim>(shape_info.fe_degree+1),
379 Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
381 else if (temp_size < max_stack_size)
383 temp1 = &temp_data[0];
384 temp2 = temp1 + temp_size;
388 temp1 = scratch_data;
389 temp2 = temp1 + temp_size;
395 if (type == MatrixFreeFunctions::truncated_tensor)
397 values_dofs = expanded_dof_values;
398 for (
unsigned int c=0; c<n_components; ++c)
399 expanded_dof_values[c] = scratch_data+2*(std::max(shape_info.dofs_per_cell,
400 shape_info.n_q_points)) +
401 c*Utilities::fixed_power<dim>(shape_info.fe_degree+1);
407 const unsigned int d1 = dim>1?1:0;
408 const unsigned int d2 = dim>2?2:0;
413 for (
unsigned int c=0; c<n_components; c++)
415 if (integrate_val ==
true)
416 eval.template values<0,false,false> (values_quad[c], values_dofs[c]);
417 if (integrate_grad ==
true)
419 if (integrate_val ==
true)
420 eval.template gradients<0,false,true> (gradients_quad[c][0], values_dofs[c]);
422 eval.template gradients<0,false,false> (gradients_quad[c][0], values_dofs[c]);
428 for (
unsigned int c=0; c<n_components; c++)
430 if (integrate_val ==
true)
433 eval.template values<0,false,false> (values_quad[c], temp1);
435 if (integrate_grad ==
true)
436 eval.template gradients<0,false,true> (gradients_quad[c][0], temp1);
437 eval.template values<1,false,false>(temp1, values_dofs[c]);
439 if (integrate_grad ==
true)
442 eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
443 if (integrate_val ==
false)
445 eval.template gradients<1,false,false>(temp1, values_dofs[c]);
447 eval.template gradients<0,false,false> (gradients_quad[c][0], temp1);
448 eval.template values<1,false,true> (temp1, values_dofs[c]);
451 eval.template gradients<1,false,true>(temp1, values_dofs[c]);
457 for (
unsigned int c=0; c<n_components; c++)
459 if (integrate_val ==
true)
462 eval.template values<0,false,false> (values_quad[c], temp1);
464 if (integrate_grad ==
true)
465 eval.template gradients<0,false,true> (gradients_quad[c][0], temp1);
466 eval.template values<1,false,false>(temp1, temp2);
467 if (integrate_grad ==
true)
469 eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
470 eval.template gradients<1,false,true>(temp1, temp2);
472 eval.template values<2,false,false> (temp2, values_dofs[c]);
474 else if (integrate_grad ==
true)
476 eval.template gradients<0,false,false>(gradients_quad[c][0], temp1);
477 eval.template values<1,false,false> (temp1, temp2);
478 eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
479 eval.template gradients<1,false,true>(temp1, temp2);
480 eval.template values<2,false,false> (temp2, values_dofs[c]);
482 if (integrate_grad ==
true)
485 eval.template values<0,false,false> (gradients_quad[c][d2], temp1);
486 eval.template values<1,false,false> (temp1, temp2);
487 eval.template gradients<2,false,true> (temp2, values_dofs[c]);
497 if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
500 for (
unsigned int c=0; c<n_components; ++c)
502 values_dofs[c][shape_info.dofs_per_cell-1] = values_quad[c][0];
503 for (
unsigned int q=1; q<shape_info.n_q_points; ++q)
504 values_dofs[c][shape_info.dofs_per_cell-1] += values_quad[c][q];
507 for (
unsigned int c=0; c<n_components; ++c)
511 if (type == MatrixFreeFunctions::truncated_tensor)
513 unsigned int count_p = 0, count_q = 0;
514 const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
515 for (
int i=0; i<(dim>2?degree+1:1); ++i)
517 for (
int j=0; j<(dim>1?degree+1-i:1); ++j)
519 for (
int k=0; k<degree+1-j-i; ++k, ++count_p, ++count_q)
521 for (
unsigned int c=0; c<n_components; ++c)
522 values_dofs_actual[c][count_p] = expanded_dof_values[c][count_q];
526 count_q += i*(degree+1);
528 AssertDimension(count_q, Utilities::fixed_power<dim>(shape_info.fe_degree+1));
534 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
int n_components,
typename Number>
535 struct FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
536 fe_degree, n_q_points_1d, n_components, Number>
539 void evaluate (
const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
545 const bool evaluate_val,
546 const bool evaluate_grad,
547 const bool evaluate_lapl);
550 void integrate (
const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
555 const bool integrate_val,
556 const bool integrate_grad);
559 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
int n_components,
typename Number>
562 FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
563 fe_degree, n_q_points_1d, n_components, Number>
564 ::evaluate (
const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
570 const bool evaluate_val,
571 const bool evaluate_grad,
572 const bool evaluate_lapl)
574 typedef EvaluatorTensorProduct<
evaluate_evenodd, dim, fe_degree, fe_degree+1,
576 Eval eval (shape_info.shape_val_evenodd,
577 shape_info.shape_gra_evenodd,
578 shape_info.shape_hes_evenodd,
579 shape_info.fe_degree,
580 shape_info.n_q_points_1d);
585 const unsigned int d1 = dim>1?1:0;
586 const unsigned int d2 = dim>2?2:0;
587 const unsigned int d3 = dim>2?3:0;
588 const unsigned int d4 = dim>2?4:0;
589 const unsigned int d5 = dim>2?5:0;
594 if (evaluate_val ==
true)
595 std::memcpy (values_quad[0], values_dofs[0],
596 eval.dofs_per_cell * n_components *
597 sizeof (values_dofs[0][0]));
598 for (
unsigned int c=0; c<n_components; c++)
600 if (evaluate_grad ==
true)
601 eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
602 if (evaluate_lapl ==
true)
603 eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]);
608 if (evaluate_val ==
true)
610 std::memcpy (values_quad[0], values_dofs[0],
611 Eval::dofs_per_cell * n_components *
612 sizeof (values_dofs[0][0]));
614 if (evaluate_grad ==
true)
615 for (
unsigned int comp=0; comp<n_components; comp++)
618 eval.template gradients<0,true,false> (values_dofs[comp],
619 gradients_quad[comp][0]);
621 eval.template gradients<1,true,false> (values_dofs[comp],
622 gradients_quad[comp][d1]);
624 if (evaluate_lapl ==
true)
625 for (
unsigned int comp=0; comp<n_components; comp++)
628 eval.template hessians<0,true,false> (values_dofs[comp],
629 hessians_quad[comp][0]);
631 eval.template hessians<1,true,false> (values_dofs[comp],
632 hessians_quad[comp][d1]);
635 eval.template gradients<0,true,false> (values_dofs[comp], scratch_data);
636 eval.template gradients<1,true,false> (scratch_data, hessians_quad[comp][d1+d1]);
641 if (evaluate_val ==
true)
643 std::memcpy (values_quad[0], values_dofs[0],
644 Eval::dofs_per_cell * n_components *
645 sizeof (values_dofs[0][0]));
647 if (evaluate_grad ==
true)
648 for (
unsigned int comp=0; comp<n_components; comp++)
651 eval.template gradients<0,true,false> (values_dofs[comp],
652 gradients_quad[comp][0]);
654 eval.template gradients<1,true,false> (values_dofs[comp],
655 gradients_quad[comp][d1]);
657 eval.template gradients<2,true,false> (values_dofs[comp],
658 gradients_quad[comp][d2]);
660 if (evaluate_lapl ==
true)
661 for (
unsigned int comp=0; comp<n_components; comp++)
664 eval.template hessians<0,true,false> (values_dofs[comp],
665 hessians_quad[comp][0]);
667 eval.template hessians<1,true,false> (values_dofs[comp],
668 hessians_quad[comp][d1]);
670 eval.template hessians<2,true,false> (values_dofs[comp],
671 hessians_quad[comp][d2]);
675 eval.template gradients<0,true,false> (values_dofs[comp], temp1);
676 eval.template gradients<1,true,false> (temp1, hessians_quad[comp][d3]);
678 eval.template gradients<2,true,false> (temp1, hessians_quad[comp][d4]);
680 eval.template gradients<1,true,false> (values_dofs[comp], temp1);
681 eval.template gradients<2,true,false> (temp1, hessians_quad[comp][d5]);
689 template <
int dim,
int fe_degree,
int n_q_po
ints_1d,
int n_components,
typename Number>
692 FEEvaluationImpl<MatrixFreeFunctions::tensor_gausslobatto, dim,
693 fe_degree, n_q_points_1d, n_components, Number>
694 ::integrate (
const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
699 const bool integrate_val,
700 const bool integrate_grad)
702 typedef EvaluatorTensorProduct<
evaluate_evenodd, dim, fe_degree, fe_degree+1,
704 Eval eval (shape_info.shape_val_evenodd,
705 shape_info.shape_gra_evenodd,
706 shape_info.shape_hes_evenodd,
707 shape_info.fe_degree,
708 shape_info.n_q_points_1d);
713 const unsigned int d1 = dim>1?1:0;
714 const unsigned int d2 = dim>2?2:0;
716 if (integrate_val ==
true)
717 std::memcpy (values_dofs[0], values_quad[0],
718 Eval::dofs_per_cell * n_components *
719 sizeof (values_dofs[0][0]));
723 for (
unsigned int c=0; c<n_components; c++)
725 if (integrate_grad ==
true)
727 if (integrate_val ==
true)
728 eval.template gradients<0,false,true> (gradients_quad[c][0],
731 eval.template gradients<0,false,false> (gradients_quad[c][0],
738 if (integrate_grad ==
true)
739 for (
unsigned int comp=0; comp<n_components; comp++)
743 if (integrate_val ==
true)
744 eval.template gradients<0, false, true> (gradients_quad[comp][0],
747 eval.template gradients<0, false, false> (gradients_quad[comp][0],
751 eval.template gradients<1, false, true> (gradients_quad[comp][d1],
757 if (integrate_grad ==
true)
758 for (
unsigned int comp=0; comp<n_components; comp++)
762 if (integrate_val ==
true)
763 eval.template gradients<0, false, true> (gradients_quad[comp][0],
766 eval.template gradients<0, false, false> (gradients_quad[comp][0],
770 eval.template gradients<1, false, true> (gradients_quad[comp][d1],
774 eval.template gradients<2, false, true> (gradients_quad[comp][d2],
787 DEAL_II_NAMESPACE_CLOSE
static const unsigned int invalid_unsigned_int
#define AssertDimension(dim1, dim2)
#define AssertThrow(cond, exc)
static ::ExceptionBase & ExcNotImplemented()