From: Martin Kronbichler Date: Wed, 29 Mar 2017 19:40:40 +0000 (+0200) Subject: Improve documentation, make names of data structures uniform X-Git-Tag: v9.0.0-rc1~1768^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=be84c1f2e57120731b62aab18165671b7e21253d;p=dealii.git Improve documentation, make names of data structures uniform --- diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 337b5a9618..c5f5af51a3 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -70,7 +70,7 @@ namespace internal }; template - struct EvaluatorSelector + struct EvaluatorSelector { static const EvaluatorVariant variant = evaluate_evenodd; }; @@ -86,11 +86,12 @@ namespace internal * the template classes EvaluatorTensorProduct which in turn are selected * from the MatrixFreeFunctions::ElementType template argument. * - * There are two specialized implementation classes FEEvaluationImplSpectral - * (for Gauss-Lobatto elements where the 'values' operation is identity) and - * FEEvaluationImplTransformSpectral (which can be transformed to a spectral - * evaluation and uses the identity in these contexts), which both allow for - * shorter code. + * There are two specialized implementation classes + * FEEvaluationImplCollocation (for Gauss-Lobatto elements where the nodal + * points and the quadrature points coincide and the 'values' operation is + * identity) and FEEvaluationImplTransformToCollocation (which can be + * transformed to a collocation space and can then use the identity in these + * spaces), which both allow for shorter code. * * @author Katharina Kormann, Martin Kronbichler, 2012, 2014, 2017 */ @@ -105,9 +106,9 @@ namespace internal VectorizedArray *gradients_quad[][dim], VectorizedArray *hessians_quad[][(dim*(dim+1))/2], VectorizedArray *scratch_data, - const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl); + const bool evaluate_values, + const bool evaluate_gradients, + const bool evaluate_hessians); static void integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, @@ -115,8 +116,8 @@ namespace internal VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], VectorizedArray *scratch_data, - const bool evaluate_val, - const bool evaluate_grad); + const bool evaluate_values, + const bool evaluate_gradients); }; @@ -131,22 +132,22 @@ namespace internal VectorizedArray *gradients_quad[][dim], VectorizedArray *hessians_quad[][(dim*(dim+1))/2], VectorizedArray *scratch_data, - const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl) + const bool evaluate_values, + const bool evaluate_gradients, + const bool evaluate_hessians) { - if (evaluate_val == false && evaluate_grad == false && evaluate_lapl == false) + if (evaluate_values == false && evaluate_gradients == false && evaluate_hessians == false) return; const EvaluatorVariant variant = EvaluatorSelector4)>::variant; typedef EvaluatorTensorProduct > Eval; - Eval eval (variant == evaluate_evenodd ? shape_info.shape_val_evenodd : + Eval eval (variant == evaluate_evenodd ? shape_info.shape_values_eo : shape_info.shape_values, - variant == evaluate_evenodd ? shape_info.shape_gra_evenodd : + variant == evaluate_evenodd ? shape_info.shape_gradients_eo : shape_info.shape_gradients, - variant == evaluate_evenodd ? shape_info.shape_hes_evenodd : + variant == evaluate_evenodd ? shape_info.shape_hessians_eo : shape_info.shape_hessians, shape_info.fe_degree, shape_info.n_q_points_1d); @@ -218,11 +219,11 @@ namespace internal case 1: for (unsigned int c=0; c (values_dofs[c], values_quad[c]); - if (evaluate_grad == true) + if (evaluate_gradients == true) eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]); - if (evaluate_lapl == true) + if (evaluate_hessians == true) eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]); } break; @@ -231,15 +232,15 @@ namespace internal for (unsigned int c=0; c (values_dofs[c], temp1); eval.template values<1,true,false> (temp1, gradients_quad[c][0]); } - if (evaluate_lapl == true) + if (evaluate_hessians == true) { // grad xy - if (evaluate_grad == false) + if (evaluate_gradients == false) eval.template gradients<0,true,false>(values_dofs[c], temp1); eval.template gradients<1,true,false> (temp1, hessians_quad[c][d1+d1]); @@ -250,15 +251,15 @@ namespace internal // grad y eval.template values<0,true,false> (values_dofs[c], temp1); - if (evaluate_grad == true) + if (evaluate_gradients == true) eval.template gradients<1,true,false> (temp1, gradients_quad[c][d1]); // grad yy - if (evaluate_lapl == true) + if (evaluate_hessians == true) eval.template hessians<1,true,false> (temp1, hessians_quad[c][d1]); // val: can use values applied in x - if (evaluate_val == true) + if (evaluate_values == true) eval.template values<1,true,false> (temp1, values_quad[c]); } break; @@ -266,7 +267,7 @@ namespace internal case 3: for (unsigned int c=0; c (values_dofs[c], temp1); @@ -274,10 +275,10 @@ namespace internal eval.template values<2,true,false> (temp2, gradients_quad[c][0]); } - if (evaluate_lapl == true) + if (evaluate_hessians == true) { // grad xz - if (evaluate_grad == false) + if (evaluate_gradients == false) { eval.template gradients<0,true,false> (values_dofs[c], temp1); eval.template values<1,true,false> (temp1, temp2); @@ -296,16 +297,16 @@ namespace internal // grad y eval.template values<0,true,false> (values_dofs[c], temp1); - if (evaluate_grad == true) + if (evaluate_gradients == true) { eval.template gradients<1,true,false>(temp1, temp2); eval.template values<2,true,false> (temp2, gradients_quad[c][d1]); } - if (evaluate_lapl == true) + if (evaluate_hessians == true) { // grad yz - if (evaluate_grad == false) + if (evaluate_gradients == false) eval.template gradients<1,true,false>(temp1, temp2); eval.template gradients<2,true,false> (temp2, hessians_quad[c][d5]); @@ -316,16 +317,16 @@ namespace internal // grad z: can use the values applied in x direction stored in temp1 eval.template values<1,true,false> (temp1, temp2); - if (evaluate_grad == true) + if (evaluate_gradients == true) eval.template gradients<2,true,false> (temp2, gradients_quad[c][d2]); // grad zz: can use the values applied in x and y direction stored // in temp2 - if (evaluate_lapl == true) + if (evaluate_hessians == true) eval.template hessians<2,true,false>(temp2, hessians_quad[c][d2]); // val: can use the values applied in x & y direction stored in temp2 - if (evaluate_val == true) + if (evaluate_values == true) eval.template values<2,true,false> (temp2, values_quad[c]); } break; @@ -336,7 +337,7 @@ namespace internal // case additional dof for FE_Q_DG0: add values; gradients and second // derivatives evaluate to zero - if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 && evaluate_val) + if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 && evaluate_values) for (unsigned int c=0; c *values_quad[], VectorizedArray *gradients_quad[][dim], VectorizedArray *scratch_data, - const bool integrate_val, - const bool integrate_grad) + const bool integrate_values, + const bool integrate_gradients) { const EvaluatorVariant variant = EvaluatorSelector4)>::variant; typedef EvaluatorTensorProduct > Eval; - Eval eval (variant == evaluate_evenodd ? shape_info.shape_val_evenodd : + Eval eval (variant == evaluate_evenodd ? shape_info.shape_values_eo : shape_info.shape_values, - variant == evaluate_evenodd ? shape_info.shape_gra_evenodd : + variant == evaluate_evenodd ? shape_info.shape_gradients_eo : shape_info.shape_gradients, - variant == evaluate_evenodd ? shape_info.shape_hes_evenodd : + variant == evaluate_evenodd ? shape_info.shape_hessians_eo : shape_info.shape_hessians, shape_info.fe_degree, shape_info.n_q_points_1d); @@ -416,11 +417,11 @@ namespace internal case 1: for (unsigned int c=0; c (values_quad[c], values_dofs[c]); - if (integrate_grad == true) + if (integrate_gradients == true) { - if (integrate_val == true) + if (integrate_values == true) eval.template gradients<0,false,true> (gradients_quad[c][0], values_dofs[c]); else eval.template gradients<0,false,false> (gradients_quad[c][0], values_dofs[c]); @@ -431,20 +432,20 @@ namespace internal case 2: for (unsigned int c=0; c (values_quad[c], temp1); //grad x - if (integrate_grad == true) + if (integrate_gradients == true) eval.template gradients<0,false,true> (gradients_quad[c][0], temp1); eval.template values<1,false,false>(temp1, values_dofs[c]); } - if (integrate_grad == true) + if (integrate_gradients == true) { // grad y eval.template values<0,false,false> (gradients_quad[c][d1], temp1); - if (integrate_val == false) + if (integrate_values == false) { eval.template gradients<1,false,false>(temp1, values_dofs[c]); //grad x @@ -460,22 +461,22 @@ namespace internal case 3: for (unsigned int c=0; c (values_quad[c], temp1); //grad x: can sum to temporary value in temp1 - if (integrate_grad == true) + if (integrate_gradients == true) eval.template gradients<0,false,true> (gradients_quad[c][0], temp1); eval.template values<1,false,false>(temp1, temp2); - if (integrate_grad == true) + if (integrate_gradients == true) { eval.template values<0,false,false> (gradients_quad[c][d1], temp1); eval.template gradients<1,false,true>(temp1, temp2); } eval.template values<2,false,false> (temp2, values_dofs[c]); } - else if (integrate_grad == true) + else if (integrate_gradients == true) { eval.template gradients<0,false,false>(gradients_quad[c][0], temp1); eval.template values<1,false,false> (temp1, temp2); @@ -483,7 +484,7 @@ namespace internal eval.template gradients<1,false,true>(temp1, temp2); eval.template values<2,false,false> (temp2, values_dofs[c]); } - if (integrate_grad == true) + if (integrate_gradients == true) { // grad z: can sum to temporary x and y value in output eval.template values<0,false,false> (gradients_quad[c][d2], temp1); @@ -500,7 +501,7 @@ namespace internal // case FE_Q_DG0: add values, gradients and second derivatives are zero if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0) { - if (integrate_val) + if (integrate_values) for (unsigned int c=0; c - struct FEEvaluationImplTransformSpectral + struct FEEvaluationImplTransformToCollocation { static void evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, @@ -555,9 +558,9 @@ namespace internal VectorizedArray *gradients_quad[][dim], VectorizedArray *hessians_quad[][(dim*(dim+1))/2], VectorizedArray *scratch_data, - const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl); + const bool evaluate_values, + const bool evaluate_gradients, + const bool evaluate_hessians); static void integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, @@ -565,14 +568,14 @@ namespace internal VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], VectorizedArray *scratch_data, - const bool integrate_val, - const bool integrate_grad); + const bool integrate_values, + const bool integrate_gradients); }; template inline void - FEEvaluationImplTransformSpectral + FEEvaluationImplTransformToCollocation ::evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, VectorizedArray *values_dofs[], VectorizedArray *values_quad[], @@ -580,19 +583,19 @@ namespace internal VectorizedArray *hessians_quad[][(dim*(dim+1))/2], VectorizedArray *, const bool , - const bool evaluate_grad, - const bool evaluate_lapl) + const bool evaluate_gradients, + const bool evaluate_hessians) { typedef EvaluatorTensorProduct > Eval; - Eval eval_val (shape_info.shape_val_evenodd, - shape_info.shape_gra_evenodd, - shape_info.shape_hes_evenodd, + Eval eval_val (shape_info.shape_values_eo, + AlignedVector >(), + AlignedVector >(), shape_info.fe_degree, shape_info.n_q_points_1d); - Eval eval(shape_info.shape_values, - shape_info.shape_grad_spectral_eo, - shape_info.shape_hessian_spectral_eo, + Eval eval(AlignedVector >(), + shape_info.shape_gradients_collocation_eo, + shape_info.shape_hessians_collocation_eo, shape_info.fe_degree, shape_info.n_q_points_1d); @@ -607,7 +610,9 @@ namespace internal for (unsigned int c=0; c(values_dofs[c], values_quad[c]); else if (dim == 2) @@ -622,31 +627,29 @@ namespace internal eval_val.template values<2,true,false>(gradients_quad[c][0], values_quad[c]); } - // apply derivatives in spectral space - if (evaluate_grad == true) + // apply derivatives in the collocation space + if (evaluate_gradients == true || evaluate_hessians == true) { eval.template gradients<0,true,false>(values_quad[c], gradients_quad[c][0]); - if (dim >= 2) + if (dim > 1) eval.template gradients<1,true,false>(values_quad[c], gradients_quad[c][d1]); - if (dim >= 3) + if (dim > 2) eval.template gradients<2,true,false>(values_quad[c], gradients_quad[c][d2]); } - if (evaluate_lapl == true) + if (evaluate_hessians == true) { eval.template hessians<0,true,false> (values_quad[c], hessians_quad[c][0]); if (dim > 1) { - eval.template gradients<0,true,false> (values_quad[c], hessians_quad[c][d2]); - eval.template gradients<1,true,false> (hessians_quad[c][d2], hessians_quad[c][d3]); + // re-use grad_x already in gradients + eval.template gradients<1,true,false> (gradients_quad[c][0], hessians_quad[c][d3]); eval.template hessians<1,true,false> (values_quad[c], hessians_quad[c][d1]); } if (dim > 2) { - // note that grad x is already in hessians_quad[c][d2] - eval.template gradients<2,true,false> (hessians_quad[c][d2], hessians_quad[c][d4]); - - eval.template gradients<1,true,false> (values_quad[c], hessians_quad[c][d2]); - eval.template gradients<2,true,false> (hessians_quad[c][d2], hessians_quad[c][d5]); + // re-use grad_x and grad_y already in gradients + eval.template gradients<2,true,false> (gradients_quad[c][0], hessians_quad[c][d4]); + eval.template gradients<2,true,false> (gradients_quad[c][d1], hessians_quad[c][d5]); eval.template hessians<2,true,false> (values_quad[c], hessians_quad[c][d2]); } } @@ -656,25 +659,25 @@ namespace internal template inline void - FEEvaluationImplTransformSpectral + FEEvaluationImplTransformToCollocation ::integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, VectorizedArray *values_dofs[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], VectorizedArray *, - const bool integrate_val, - const bool integrate_grad) + const bool integrate_values, + const bool integrate_gradients) { typedef EvaluatorTensorProduct > Eval; - Eval eval_val (shape_info.shape_val_evenodd, - shape_info.shape_gra_evenodd, - shape_info.shape_hes_evenodd, + Eval eval_val (shape_info.shape_values_eo, + AlignedVector >(), + AlignedVector >(), shape_info.fe_degree, shape_info.n_q_points_1d); - Eval eval(shape_info.shape_values, - shape_info.shape_grad_spectral_eo, - shape_info.shape_hessian_spectral_eo, + Eval eval(AlignedVector >(), + shape_info.shape_gradients_collocation_eo, + shape_info.shape_hessians_collocation_eo, shape_info.fe_degree, shape_info.n_q_points_1d); @@ -686,20 +689,20 @@ namespace internal for (unsigned int c=0; c(gradients_quad[c][0], values_quad[c]); else eval.template gradients<0,false,false>(gradients_quad[c][0], values_quad[c]); - if (dim >= 2) + if (dim > 1) eval.template gradients<1,false,true>(gradients_quad[c][d1], values_quad[c]); - if (dim >= 3) + if (dim > 2) eval.template gradients<2,false,true>(gradients_quad[c][d2], values_quad[c]); } - // transform back to the usual space + // transform back to the original space if (dim == 1) eval_val.template values<0,false,false>(values_quad[c], values_dofs[c]); else if (dim == 2) @@ -720,14 +723,20 @@ namespace internal /** * This struct performs the evaluation of function values, gradients and - * Hessians for tensor-product finite elements. This a specialization for - * "spectral" elements like Gauss-Lobatto elements where the 'values' - * operation is identity, which allows us to write shorter code. + * Hessians for tensor-product finite elements. This a specialization for + * elements where the nodal points coincide with the quadrature points like + * FE_Q shape functions on Gauss-Lobatto elements integrated with + * Gauss-Lobatto quadrature. The assumption of this class is that the shape + * 'values' operation is identity, which allows us to write shorter code. + * + * In literature, this form of evaluation is often called spectral + * evaluation, spectral collocation or simply collocation, meaning the same + * location for shape functions and evaluation space (quadrature points). * * @author Katharina Kormann, 2012 */ template - struct FEEvaluationImplSpectral + struct FEEvaluationImplCollocation { static void evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, @@ -736,9 +745,9 @@ namespace internal VectorizedArray *gradients_quad[][dim], VectorizedArray *hessians_quad[][(dim*(dim+1))/2], VectorizedArray *scratch_data, - const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl); + const bool evaluate_values, + const bool evaluate_gradients, + const bool evaluate_hessians); static void integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, @@ -746,31 +755,31 @@ namespace internal VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], VectorizedArray *scratch_data, - const bool integrate_val, - const bool integrate_grad); + const bool integrate_values, + const bool integrate_gradients); }; template inline void - FEEvaluationImplSpectral + FEEvaluationImplCollocation ::evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, VectorizedArray *values_dofs[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], VectorizedArray *hessians_quad[][(dim*(dim+1))/2], VectorizedArray *, - const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl) + const bool evaluate_values, + const bool evaluate_gradients, + const bool evaluate_hessians) { typedef EvaluatorTensorProduct > Eval; - Eval eval (shape_info.shape_val_evenodd, - shape_info.shape_gra_evenodd, - shape_info.shape_hes_evenodd, - shape_info.fe_degree, - shape_info.n_q_points_1d); + Eval eval(AlignedVector >(), + shape_info.shape_gradients_eo, + shape_info.shape_hessians_eo, + shape_info.fe_degree, + shape_info.n_q_points_1d); // These avoid compiler warnings; they are only used in sensible context // but compilers typically cannot detect when we access something like @@ -783,33 +792,31 @@ namespace internal for (unsigned int c=0; c(values_dofs[c], gradients_quad[c][0]); - if (dim >= 2) + if (dim > 1) eval.template gradients<1,true,false>(values_dofs[c], gradients_quad[c][d1]); - if (dim >= 3) + if (dim > 2) eval.template gradients<2,true,false>(values_dofs[c], gradients_quad[c][d2]); } - if (evaluate_lapl == true) + if (evaluate_hessians == true) { eval.template hessians<0,true,false> (values_quad[c], hessians_quad[c][0]); if (dim > 1) { - eval.template gradients<0,true,false> (values_dofs[c], hessians_quad[c][d2]); - eval.template gradients<1,true,false> (hessians_quad[c][d2], hessians_quad[c][d3]); + // re-use grad_x already in gradients + eval.template gradients<1,true,false> (gradients_quad[c][0], hessians_quad[c][d3]); eval.template hessians<1,true,false> (values_dofs[c], hessians_quad[c][d1]); } if (dim > 2) { - // note that grad x is already in hessians_quad[c][d2] - eval.template gradients<2,true,false> (hessians_quad[c][d2], hessians_quad[c][d4]); - - eval.template gradients<1,true,false> (values_dofs[c], hessians_quad[c][d2]); - eval.template gradients<2,true,false> (hessians_quad[c][d2], hessians_quad[c][d5]); + // re-use grad_x already in gradients + eval.template gradients<2,true,false> (gradients_quad[c][0], hessians_quad[c][d4]); + eval.template gradients<2,true,false> (gradients_quad[c][d1], hessians_quad[c][d5]); eval.template hessians<2,true,false> (values_dofs[c], hessians_quad[c][d2]); } } @@ -819,22 +826,22 @@ namespace internal template inline void - FEEvaluationImplSpectral + FEEvaluationImplCollocation ::integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, VectorizedArray *values_dofs[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], VectorizedArray *, - const bool integrate_val, - const bool integrate_grad) + const bool integrate_values, + const bool integrate_gradients) { typedef EvaluatorTensorProduct > Eval; - Eval eval (shape_info.shape_val_evenodd, - shape_info.shape_gra_evenodd, - shape_info.shape_hes_evenodd, - shape_info.fe_degree, - shape_info.n_q_points_1d); + Eval eval(AlignedVector >(), + shape_info.shape_gradients_eo, + shape_info.shape_hessians_eo, + shape_info.fe_degree, + shape_info.n_q_points_1d); // These avoid compiler warnings; they are only used in sensible context // but compilers typically cannot detect when we access something like @@ -844,18 +851,18 @@ namespace internal for (unsigned int c=0; c(gradients_quad[c][0], values_dofs[c]); else eval.template gradients<0,false,false>(gradients_quad[c][0], values_dofs[c]); - if (dim >= 2) + if (dim > 1) eval.template gradients<1,false,true>(gradients_quad[c][d1], values_dofs[c]); - if (dim >= 3) + if (dim > 2) eval.template gradients<2,false,true>(gradients_quad[c][d2], values_dofs[c]); } } diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 906fe53911..3eb419d967 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -2741,6 +2741,8 @@ namespace internal res = vector_access (const_cast(vec), index); } + // variant where VectorType::value_type is the same as Number -> can call + // gather template void process_dof_gather (const unsigned int *indices, VectorType &vec, @@ -2750,6 +2752,8 @@ namespace internal res.gather(vec.begin(), indices); } + // variant where VectorType::value_type is not the same as Number -> must + // manually load the data template void process_dof_gather (const unsigned int *indices, VectorType &vec, @@ -2807,6 +2811,8 @@ namespace internal vector_access (vec, index) += res; } + // variant where VectorType::value_type is the same as Number -> can call + // scatter template void process_dof_gather (const unsigned int *indices, VectorType &vec, @@ -2814,7 +2820,6 @@ namespace internal internal::bool2type) const { // TODO: enable scatter path when indices are fixed - //#if DEAL_II_COMPILER_VECTORIZATION_LEVEL < 3 #if 1 for (unsigned int v=0; v::n_array_elements; ++v) @@ -2828,6 +2833,8 @@ namespace internal #endif } + // variant where VectorType::value_type is not the same as Number -> must + // manually append all data template void process_dof_gather (const unsigned int *indices, VectorType &vec, @@ -5408,9 +5415,9 @@ FEEvaluation else { if (fe_degree+1 == n_q_points_1d && - this->data->element_type == internal::MatrixFreeFunctions::tensor_gausslobatto) + this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation) { - internal::FEEvaluationImplSpectral + internal::FEEvaluationImplCollocation ::evaluate(*this->data, &this->values_dofs[0], this->values_quad, this->gradients_quad, this->hessians_quad, this->scratch_data, evaluate_val, evaluate_grad, evaluate_lapl); @@ -5418,7 +5425,7 @@ FEEvaluation else if (fe_degree+1 == n_q_points_1d && this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric) { - internal::FEEvaluationImplTransformSpectral + internal::FEEvaluationImplTransformToCollocation ::evaluate(*this->data, &this->values_dofs[0], this->values_quad, this->gradients_quad, this->hessians_quad, this->scratch_data, evaluate_val, evaluate_grad, evaluate_lapl); @@ -5520,9 +5527,9 @@ FEEvaluation else { if (fe_degree+1 == n_q_points_1d && - this->data->element_type == internal::MatrixFreeFunctions::tensor_gausslobatto) + this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation) { - internal::FEEvaluationImplSpectral + internal::FEEvaluationImplCollocation ::integrate(*this->data, &this->values_dofs[0], this->values_quad, this->gradients_quad, this->scratch_data, integrate_val, integrate_grad); @@ -5530,7 +5537,7 @@ FEEvaluation else if (fe_degree+1 == n_q_points_1d && this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric) { - internal::FEEvaluationImplTransformSpectral + internal::FEEvaluationImplTransformToCollocation ::integrate(*this->data, &this->values_dofs[0], this->values_quad, this->gradients_quad, this->scratch_data, integrate_val, integrate_grad); diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index 95b981f43b..eafeda2f72 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -40,11 +40,42 @@ namespace internal */ enum ElementType { - tensor_general, - tensor_symmetric, - truncated_tensor, - tensor_symmetric_plus_dg0, - tensor_gausslobatto + /** + * Tensor product shape function where the shape value evaluation in the + * quadrature point corresponds to the identity operation and no + * interpolation needs to be performed (collocation approach, also + * called spectral evaluation). This is for example the case for an + * element with nodes in the Gauss-Lobatto support points and + * integration in the Gauss-Lobatto quadrature points of the same order. + */ + tensor_symmetric_collocation = 0, + /** + * Symmetric tensor product shape functions fulfilling a Hermite + * identity with values and first derivatives zero at the element end + * points in 1D. + */ + tensor_symmetric_hermite = 1, + /** + * Usual tensor product shape functions whose shape values and + * quadrature points are symmetric about the midpoint of the unit + * interval 0.5 + */ + tensor_symmetric = 2, + /** + * Tensor product shape functions without further particular properties + */ + tensor_general = 3, + /** + * Polynomials of complete degree rather than tensor degree which can be + * described by a truncated tensor product + */ + truncated_tensor = 4, + /** + * Tensor product shape functions that are symmetric about the midpoint + * of the unit interval 0.5 that additionally add a constant shape + * function according to FE_Q_DG0. + */ + tensor_symmetric_plus_dg0 = 5 }; /** @@ -128,33 +159,35 @@ namespace internal * even-odd scheme where the symmetries in shape_values are used for * faster evaluation. */ - AlignedVector > shape_val_evenodd; + AlignedVector > shape_values_eo; /** * Stores the shape gradients in a different format, namely the so- * called even-odd scheme where the symmetries in shape_gradients are * used for faster evaluation. */ - AlignedVector > shape_gra_evenodd; + AlignedVector > shape_gradients_eo; /** * Stores the shape second derivatives in a different format, namely the * so-called even-odd scheme where the symmetries in shape_hessians are * used for faster evaluation. */ - AlignedVector > shape_hes_evenodd; + AlignedVector > shape_hessians_eo; /** - * Stores the shape gradients of the spectral element space - * FE_DGQ<1>(Quadrature<1>) for faster evaluation in even-odd format. + * Stores the shape gradients of the shape function space associated to + * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). For + * faster evaluation only the even-odd format is necessary. */ - AlignedVector > shape_grad_spectral_eo; + AlignedVector > shape_gradients_collocation_eo; /** - * Stores the shape hessians of the spectral element space - * FE_DGQ<1>(Quadrature<1>) for faster evaluation in even-odd format. + * Stores the shape hessians of the shape function space associated to + * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). For + * faster evaluation only the even-odd format is necessary. */ - AlignedVector > shape_hessian_spectral_eo; + AlignedVector > shape_hessians_collocation_eo; /** * Stores the indices from cell DoFs to face DoFs. The rows go through @@ -234,16 +267,17 @@ namespace internal /** * Check whether we have symmetries in the shape values. In that case, - * also fill the shape_???_evenodd fields. + * also fill the shape_???_eo fields. */ bool check_1d_shapes_symmetric(const unsigned int n_q_points_1d); /** * Check whether symmetric 1D basis functions are such that the shape - * values form a diagonal matrix, which allows to use specialized - * algorithms that save some operations. + * values form a diagonal matrix, i.e., the nodal points are collocated + * with the quadrature points. This allows for specialized algorithms + * that save some operations in the evaluation. */ - bool check_1d_shapes_gausslobatto(); + bool check_1d_shapes_collocation(); }; @@ -257,6 +291,7 @@ namespace internal const FiniteElement &fe_in, const unsigned int base_element_number) : + element_type(tensor_general), fe_degree (0), n_q_points_1d (0), n_q_points (0), @@ -267,9 +302,8 @@ namespace internal reinit (quad, fe_in, base_element_number); } - - } // end of namespace MatrixFreeFunctions + } // end of namespace internal DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 9d7b32818b..45918788d9 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -212,40 +212,42 @@ namespace internal if (element_type == tensor_general && check_1d_shapes_symmetric(n_q_points_1d)) { - if (check_1d_shapes_gausslobatto()) - element_type = tensor_gausslobatto; + if (check_1d_shapes_collocation()) + element_type = tensor_symmetric_collocation; else { element_type = tensor_symmetric; - // get spectral evaluation points + // get gradient and Hessian transformation matrix for the + // polynomial space associated with the quadrature rule + // (collocation space) if (fe_degree+1 == n_q_points_1d) { const unsigned int stride = fe_degree/2+1; - shape_grad_spectral_eo.resize((fe_degree+1)*stride); - shape_hessian_spectral_eo.resize((fe_degree+1)*stride); + shape_gradients_collocation_eo.resize((fe_degree+1)*stride); + shape_hessians_collocation_eo.resize((fe_degree+1)*stride); FE_DGQArbitraryNodes<1> fe(quad.get_points()); for (unsigned int i=0; i<(fe_degree+1)/2; ++i) for (unsigned int q=0; q bool - ShapeInfo::check_1d_shapes_gausslobatto() + ShapeInfo::check_1d_shapes_collocation() { if (dofs_per_cell != n_q_points) return false; @@ -441,11 +443,11 @@ namespace internal memory += MemoryConsumption::memory_consumption(shape_values); memory += MemoryConsumption::memory_consumption(shape_gradients); memory += MemoryConsumption::memory_consumption(shape_hessians); - memory += MemoryConsumption::memory_consumption(shape_val_evenodd); - memory += MemoryConsumption::memory_consumption(shape_gra_evenodd); - memory += MemoryConsumption::memory_consumption(shape_hes_evenodd); - memory += MemoryConsumption::memory_consumption(shape_grad_spectral_eo); - memory += MemoryConsumption::memory_consumption(shape_hessian_spectral_eo); + memory += MemoryConsumption::memory_consumption(shape_values_eo); + memory += MemoryConsumption::memory_consumption(shape_gradients_eo); + memory += MemoryConsumption::memory_consumption(shape_hessians_eo); + memory += MemoryConsumption::memory_consumption(shape_gradients_collocation_eo); + memory += MemoryConsumption::memory_consumption(shape_hessians_collocation_eo); memory += face_indices.memory_consumption(); for (unsigned int i=0; i<2; ++i) {