From: Martin Kronbichler Date: Tue, 28 Mar 2017 13:24:45 +0000 (+0200) Subject: Put new code into separate function. X-Git-Tag: v9.0.0-rc1~1768^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1ac01aa6d4fcadca67fd9a60c90776ca10751023;p=dealii.git Put new code into separate function. Do not use function pointers. --- diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 893bd2ac44..337b5a9618 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -77,17 +77,23 @@ namespace internal - // This struct performs the evaluation of function values, gradients and - // Hessians for tensor-product finite elements. The operation is used for - // both the symmetric and non-symmetric case, which use different apply - // functions 'values', 'gradients' in the individual coordinate - // directions. The apply functions for values are provided through one of - // the template classes EvaluatorTensorProduct which in turn are selected - // from the MatrixFreeFunctions::ElementType template argument. - // - // There is a specialization made for Gauss-Lobatto elements further down - // where the 'values' operation is identity, which allows us to write - // shorter code. + /** + * This struct performs the evaluation of function values, gradients and + * Hessians for tensor-product finite elements. The operation is used for + * both the symmetric and non-symmetric case, which use different apply + * functions 'values', 'gradients' in the individual coordinate + * directions. The apply functions for values are provided through one of + * 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. + * + * @author Katharina Kormann, Martin Kronbichler, 2012, 2014, 2017 + */ template struct FEEvaluationImpl @@ -207,46 +213,6 @@ namespace internal const unsigned int d4 = dim>2?4:0; const unsigned int d5 = dim>2?5:0; - // check if we can go through the spectral evaluation option which is - // faster than the standard one - if (fe_degree+1 == n_q_points_1d && - (type == MatrixFreeFunctions::tensor_symmetric || - type == MatrixFreeFunctions::tensor_general) && - evaluate_lapl == false) - { - Eval eval_grad(shape_info.shape_values, - variant == evaluate_evenodd ? shape_info.shape_grad_spectral_eo : - shape_info.shape_grad_spectral, - shape_info.shape_hessians, - shape_info.fe_degree, - shape_info.n_q_points_1d); - for (unsigned int c=0; c(values_dofs[c], values_quad[c]); - else if (dim == 2) - { - eval.template values<0,true,false>(values_dofs[c], gradients_quad[c][0]); - eval.template values<1,true,false>(gradients_quad[c][0], values_quad[c]); - } - else if (dim == 3) - { - eval.template values<0,true,false>(values_dofs[c], values_quad[c]); - eval.template values<1,true,false>(values_quad[c], gradients_quad[c][0]); - eval.template values<2,true,false>(gradients_quad[c][0], values_quad[c]); - } - if (evaluate_grad == true) - { - eval_grad.template gradients<0,true,false>(values_quad[c], gradients_quad[c][0]); - if (dim >= 2) - eval_grad.template gradients<1,true,false>(values_quad[c], gradients_quad[c][d1]); - if (dim >= 3) - eval_grad.template gradients<2,true,false>(values_quad[c], gradients_quad[c][d2]); - } - } - return; - } - switch (dim) { case 1: @@ -445,48 +411,6 @@ namespace internal const unsigned int d1 = dim>1?1:0; const unsigned int d2 = dim>2?2:0; - // check if we can go through the spectral evaluation option which is - // faster than the standard one - if (fe_degree+1 == n_q_points_1d && - (type == MatrixFreeFunctions::tensor_symmetric || - type == MatrixFreeFunctions::tensor_general)) - { - Eval eval_grad(shape_info.shape_values, - variant == evaluate_evenodd ? shape_info.shape_grad_spectral_eo : - shape_info.shape_grad_spectral, - shape_info.shape_hessians, - shape_info.fe_degree, - shape_info.n_q_points_1d); - for (unsigned int c=0; c(gradients_quad[c][0], values_quad[c]); - else - eval_grad.template gradients<0,false,false>(gradients_quad[c][0], values_quad[c]); - if (dim >= 2) - eval_grad.template gradients<1,false,true>(gradients_quad[c][d1], values_quad[c]); - if (dim >= 3) - eval_grad.template gradients<2,false,true>(gradients_quad[c][d2], values_quad[c]); - } - if (dim == 1) - eval.template values<0,false,false>(values_quad[c], values_dofs[c]); - else if (dim == 2) - { - eval.template values<0,false,false>(values_quad[c], gradients_quad[c][0]); - eval.template values<1,false,false>(gradients_quad[c][0], values_dofs[c]); - } - else if (dim == 3) - { - eval.template values<0,false,false>(values_quad[c], gradients_quad[c][0]); - eval.template values<1,false,false>(gradients_quad[c][0], values_quad[c]); - eval.template values<2,false,false>(values_quad[c], values_dofs[c]); - } - } - return; - } - switch (dim) { case 1: @@ -609,12 +533,20 @@ namespace internal } } - // This a specialization for "spectral" elements like Gauss-Lobatto elements - // where the 'values' operation is identity, which allows us to write - // shorter code. - template - struct FEEvaluationImpl + + + /** + * This struct performs the evaluation of function values, gradients and + * Hessians for tensor-product finite elements. This a specialization for + * symmetric basis functions with the same number of quadrature points as + * degrees of freedom. In that case, we can first transform to a spectral + * basis and then perform the evaluation of the first and second derivatives + * in spectral space, using the identity operation for the shape values. + * + * @author Katharina Kormann, Martin Kronbichler, 2017 + */ + template + struct FEEvaluationImplTransformSpectral { static void evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, @@ -637,148 +569,200 @@ namespace internal const bool integrate_grad); }; - template + template inline void - FEEvaluationImpl - ::evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, - VectorizedArray *values_dofs[], - VectorizedArray *values_quad[], - 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) + FEEvaluationImplTransformSpectral + ::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 , + const bool evaluate_grad, + const bool evaluate_lapl) { 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_val (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(shape_info.shape_values, + shape_info.shape_grad_spectral_eo, + shape_info.shape_hessian_spectral_eo, + shape_info.fe_degree, + shape_info.n_q_points_1d); - // These avoid compiler errors; they are only used in sensible context but + // These avoid compiler warnings; they are only used in sensible context but // compilers typically cannot detect when we access something like // gradients_quad[2] only for dim==3. const unsigned int d1 = dim>1?1:0; - const unsigned int d2 = dim>2?2:0; - const unsigned int d3 = dim>2?3:0; - const unsigned int d4 = dim>2?4:0; - const unsigned int d5 = dim>2?5:0; + const unsigned int d2 = dim>2?2:d1; + const unsigned int d3 = d1+d2; + const unsigned int d4 = dim>2?4:d3; + const unsigned int d5 = dim>2?5:d4; - switch (dim) + for (unsigned int c=0; c(values_dofs[c], values_quad[c]); + else if (dim == 2) { - if (evaluate_grad == true) - eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]); - if (evaluate_lapl == true) - eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]); + eval_val.template values<0,true,false>(values_dofs[c], gradients_quad[c][0]); + eval_val.template values<1,true,false>(gradients_quad[c][0], values_quad[c]); } - break; - - case 2: - if (evaluate_val == true) + else if (dim == 3) { - std::memcpy (values_quad[0], values_dofs[0], - Eval::dofs_per_cell * n_components * - sizeof (values_dofs[0][0])); + eval_val.template values<0,true,false>(values_dofs[c], values_quad[c]); + eval_val.template values<1,true,false>(values_quad[c], gradients_quad[c][0]); + eval_val.template values<2,true,false>(gradients_quad[c][0], values_quad[c]); } + + // apply derivatives in spectral space if (evaluate_grad == true) - for (unsigned int comp=0; comp (values_dofs[comp], - gradients_quad[comp][0]); - // grad y - eval.template gradients<1,true,false> (values_dofs[comp], - gradients_quad[comp][d1]); - } + { + eval.template gradients<0,true,false>(values_quad[c], gradients_quad[c][0]); + if (dim >= 2) + eval.template gradients<1,true,false>(values_quad[c], gradients_quad[c][d1]); + if (dim >= 3) + eval.template gradients<2,true,false>(values_quad[c], gradients_quad[c][d2]); + } if (evaluate_lapl == true) - for (unsigned int comp=0; comp (values_dofs[comp], - hessians_quad[comp][0]); - // hess y - eval.template hessians<1,true,false> (values_dofs[comp], - hessians_quad[comp][d1]); - - // grad x grad y - eval.template gradients<0,true,false> (values_dofs[comp], scratch_data); - eval.template gradients<1,true,false> (scratch_data, hessians_quad[comp][d1+d1]); - } - break; + { + 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]); + 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]); - case 3: - if (evaluate_val == true) + 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]); + eval.template hessians<2,true,false> (values_quad[c], hessians_quad[c][d2]); + } + } + } + } + + template + inline + void + FEEvaluationImplTransformSpectral + ::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) + { + typedef EvaluatorTensorProduct > Eval; + Eval eval_val (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(shape_info.shape_values, + shape_info.shape_grad_spectral_eo, + shape_info.shape_hessian_spectral_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 + // gradients_quad[2] only for dim==3. + const unsigned int d1 = dim>1?1:0; + const unsigned int d2 = dim>2?2:0; + + 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) + eval.template gradients<1,false,true>(gradients_quad[c][d1], values_quad[c]); + if (dim >= 3) + eval.template gradients<2,false,true>(gradients_quad[c][d2], values_quad[c]); + } + + // transform back to the usual space + if (dim == 1) + eval_val.template values<0,false,false>(values_quad[c], values_dofs[c]); + else if (dim == 2) + { + eval_val.template values<0,false,false>(values_quad[c], gradients_quad[c][0]); + eval_val.template values<1,false,false>(gradients_quad[c][0], values_dofs[c]); + } + else if (dim == 3) + { + eval_val.template values<0,false,false>(values_quad[c], gradients_quad[c][0]); + eval_val.template values<1,false,false>(gradients_quad[c][0], values_quad[c]); + eval_val.template values<2,false,false>(values_quad[c], values_dofs[c]); } - if (evaluate_grad == true) - for (unsigned int comp=0; comp (values_dofs[comp], - gradients_quad[comp][0]); - // grad y - eval.template gradients<1,true,false> (values_dofs[comp], - gradients_quad[comp][d1]); - // grad y - eval.template gradients<2,true,false> (values_dofs[comp], - gradients_quad[comp][d2]); - } - if (evaluate_lapl == true) - for (unsigned int comp=0; comp (values_dofs[comp], - hessians_quad[comp][0]); - // grad y - eval.template hessians<1,true,false> (values_dofs[comp], - hessians_quad[comp][d1]); - // grad y - eval.template hessians<2,true,false> (values_dofs[comp], - hessians_quad[comp][d2]); - - VectorizedArray *temp1 = scratch_data; - // grad xy - eval.template gradients<0,true,false> (values_dofs[comp], temp1); - eval.template gradients<1,true,false> (temp1, hessians_quad[comp][d3]); - // grad xz - eval.template gradients<2,true,false> (temp1, hessians_quad[comp][d4]); - // grad yz - eval.template gradients<1,true,false> (values_dofs[comp], temp1); - eval.template gradients<2,true,false> (temp1, hessians_quad[comp][d5]); - } - break; - default: - AssertThrow(false, ExcNotImplemented()); } } - template + + + /** + * 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. + * + * @author Katharina Kormann, 2012 + */ + template + struct FEEvaluationImplSpectral + { + static + void evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs[], + VectorizedArray *values_quad[], + 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); + + static + void integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + VectorizedArray *scratch_data, + const bool integrate_val, + const bool integrate_grad); + }; + + template inline void - FEEvaluationImpl - ::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) + FEEvaluationImplSpectral + ::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) { typedef EvaluatorTensorProduct > Eval; @@ -788,77 +772,92 @@ namespace internal shape_info.fe_degree, shape_info.n_q_points_1d); - // These avoid compiler errors; they are only used in sensible context but - // compilers typically cannot detect when we access something like + // These avoid compiler warnings; they are only used in sensible context + // but compilers typically cannot detect when we access something like // gradients_quad[2] only for dim==3. const unsigned int d1 = dim>1?1:0; - const unsigned int d2 = dim>2?2:0; + const unsigned int d2 = dim>2?2:d1; + const unsigned int d3 = d1+d2; + const unsigned int d4 = dim>2?4:d3; + const unsigned int d5 = dim>2?5:d4; - if (integrate_val == true) - std::memcpy (values_dofs[0], values_quad[0], - Eval::dofs_per_cell * n_components * - sizeof (values_dofs[0][0])); - switch (dim) + for (unsigned int c=0; c(values_dofs[c], gradients_quad[c][0]); + if (dim >= 2) + eval.template gradients<1,true,false>(values_dofs[c], gradients_quad[c][d1]); + if (dim >= 3) + eval.template gradients<2,true,false>(values_dofs[c], gradients_quad[c][d2]); + } + if (evaluate_lapl == true) + { + eval.template hessians<0,true,false> (values_quad[c], hessians_quad[c][0]); + if (dim > 1) { - if (integrate_val == 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]); + 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]); + 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]); + eval.template hessians<2,true,false> (values_dofs[c], hessians_quad[c][d2]); } } + } + } - break; - case 2: - if (integrate_grad == true) - for (unsigned int comp=0; comp (gradients_quad[comp][0], - values_dofs[comp]); - else - eval.template gradients<0, false, false> (gradients_quad[comp][0], - values_dofs[comp]); - - // grad y - eval.template gradients<1, false, true> (gradients_quad[comp][d1], - values_dofs[comp]); - } - break; + template + inline + void + FEEvaluationImplSpectral + ::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) + { + 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); - case 3: - if (integrate_grad == true) - for (unsigned int comp=0; comp (gradients_quad[comp][0], - values_dofs[comp]); - else - eval.template gradients<0, false, false> (gradients_quad[comp][0], - values_dofs[comp]); - - // grad y - eval.template gradients<1, false, true> (gradients_quad[comp][d1], - values_dofs[comp]); - - // grad z - eval.template gradients<2, false, true> (gradients_quad[comp][d2], - values_dofs[comp]); - } - break; + // These avoid compiler warnings; they are only used in sensible context + // but compilers typically cannot detect when we access something like + // gradients_quad[2] only for dim==3. + const unsigned int d1 = dim>1?1:0; + const unsigned int d2 = dim>2?2:0; - default: - AssertThrow(false, ExcNotImplemented()); + 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) + eval.template gradients<1,false,true>(gradients_quad[c][d1], values_dofs[c]); + if (dim >= 3) + 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 a8e8bef3e8..143002fc7c 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -2086,30 +2086,6 @@ private: */ void check_template_arguments(const unsigned int fe_no, const unsigned int first_selected_component); - - /** - * Function pointer for the evaluate function - */ - void (*evaluate_funct) (const internal::MatrixFreeFunctions::ShapeInfo &, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - 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); - - /** - * Function pointer for the integrate function - */ - void (*integrate_funct)(const internal::MatrixFreeFunctions::ShapeInfo &, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, - const bool evaluate_val, - const bool evaluate_grad); }; @@ -5214,94 +5190,6 @@ FEEvaluation ::check_template_arguments(const unsigned int fe_no, const unsigned int first_selected_component) { - const unsigned int fe_degree_templ = fe_degree != -1 ? fe_degree : 0; - const unsigned int n_q_points_1d_templ = n_q_points_1d > 0 ? n_q_points_1d : 1; - if (fe_degree == -1) - { - if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0) - { - evaluate_funct = internal::FEEvaluationImpl::evaluate; - integrate_funct = internal::FEEvaluationImpl::integrate; - } - else if (this->data->element_type == internal::MatrixFreeFunctions::truncated_tensor) - { - evaluate_funct = internal::FEEvaluationImpl::evaluate; - integrate_funct = internal::FEEvaluationImpl::integrate; - } - else - { - evaluate_funct = internal::FEEvaluationImpl::evaluate; - integrate_funct = internal::FEEvaluationImpl::integrate; - } - } - else - switch (this->data->element_type) - { - case internal::MatrixFreeFunctions::tensor_symmetric: - evaluate_funct = - internal::FEEvaluationImpl::evaluate; - integrate_funct = - internal::FEEvaluationImpl::integrate; - break; - - case internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0: - evaluate_funct = - internal::FEEvaluationImpl::evaluate; - integrate_funct = - internal::FEEvaluationImpl::integrate; - break; - - case internal::MatrixFreeFunctions::tensor_general: - evaluate_funct = - internal::FEEvaluationImpl::evaluate; - integrate_funct = - internal::FEEvaluationImpl::integrate; - break; - - case internal::MatrixFreeFunctions::tensor_gausslobatto: - evaluate_funct = - internal::FEEvaluationImpl::evaluate; - integrate_funct = - internal::FEEvaluationImpl::integrate; - break; - - case internal::MatrixFreeFunctions::truncated_tensor: - evaluate_funct = - internal::FEEvaluationImpl::evaluate; - integrate_funct = - internal::FEEvaluationImpl::integrate; - break; - - default: - AssertThrow(false, ExcNotImplemented()); - } - (void)fe_no; (void)first_selected_component; @@ -5484,11 +5372,89 @@ FEEvaluation Assert(this->matrix_info != 0 || this->mapped_geometry->is_initialized(), ExcNotInitialized()); - // Select algorithm matching the element type at run time (the function - // pointer is easy to predict, so negligible in cost) - evaluate_funct (*this->data, &this->values_dofs[0], this->values_quad, - this->gradients_quad, this->hessians_quad, this->scratch_data, - evaluate_val, evaluate_grad, evaluate_lapl); + // Select algorithm matching the element type at run time + const unsigned int fe_degree_templ = fe_degree != -1 ? fe_degree : 0; + const unsigned int n_q_points_1d_templ = n_q_points_1d > 0 ? n_q_points_1d : 1; + if (fe_degree == -1) + { + if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0) + { + internal::FEEvaluationImpl + ::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); + } + else if (this->data->element_type == internal::MatrixFreeFunctions::truncated_tensor) + { + internal::FEEvaluationImpl + ::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); + } + else + { + internal::FEEvaluationImpl + ::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); + } + } + else + { + if (fe_degree+1 == n_q_points_1d && + this->data->element_type == internal::MatrixFreeFunctions::tensor_gausslobatto) + { + internal::FEEvaluationImplSpectral + ::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); + } + else if (fe_degree+1 == n_q_points_1d && + this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric) + { + internal::FEEvaluationImplTransformSpectral + ::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); + } + else if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric) + { + internal::FEEvaluationImpl + ::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); + } + else if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0) + { + internal::FEEvaluationImpl + ::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); + } + else if (this->data->element_type == internal::MatrixFreeFunctions::tensor_general) + { + internal::FEEvaluationImpl + ::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); + } + else if (this->data->element_type == internal::MatrixFreeFunctions::truncated_tensor) + { + internal::FEEvaluationImpl + ::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); + } + else + AssertThrow(false, ExcNotImplemented()); + } #ifdef DEBUG if (evaluate_val == true) @@ -5518,11 +5484,89 @@ FEEvaluation Assert(this->matrix_info != 0 || this->mapped_geometry->is_initialized(), ExcNotInitialized()); - // Select algorithm matching the element type at run time (the function - // pointer is easy to predict, so negligible in cost) - integrate_funct (*this->data, this->values_dofs, this->values_quad, - this->gradients_quad, this->scratch_data, - integrate_val, integrate_grad); + // Select algorithm matching the element type at run time + const unsigned int fe_degree_templ = fe_degree != -1 ? fe_degree : 0; + const unsigned int n_q_points_1d_templ = n_q_points_1d > 0 ? n_q_points_1d : 1; + if (fe_degree == -1) + { + if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0) + { + internal::FEEvaluationImpl + ::integrate(*this->data, &this->values_dofs[0], this->values_quad, + this->gradients_quad, this->scratch_data, + integrate_val, integrate_grad); + } + else if (this->data->element_type == internal::MatrixFreeFunctions::truncated_tensor) + { + internal::FEEvaluationImpl + ::integrate(*this->data, &this->values_dofs[0], this->values_quad, + this->gradients_quad, this->scratch_data, + integrate_val, integrate_grad); + } + else + { + internal::FEEvaluationImpl + ::integrate(*this->data, &this->values_dofs[0], this->values_quad, + this->gradients_quad, this->scratch_data, + integrate_val, integrate_grad); + } + } + else + { + if (fe_degree+1 == n_q_points_1d && + this->data->element_type == internal::MatrixFreeFunctions::tensor_gausslobatto) + { + internal::FEEvaluationImplSpectral + ::integrate(*this->data, &this->values_dofs[0], this->values_quad, + this->gradients_quad, this->scratch_data, + integrate_val, integrate_grad); + } + else if (fe_degree+1 == n_q_points_1d && + this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric) + { + internal::FEEvaluationImplTransformSpectral + ::integrate(*this->data, &this->values_dofs[0], this->values_quad, + this->gradients_quad, this->scratch_data, + integrate_val, integrate_grad); + } + else if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric) + { + internal::FEEvaluationImpl + ::integrate(*this->data, &this->values_dofs[0], this->values_quad, + this->gradients_quad, this->scratch_data, + integrate_val, integrate_grad); + } + else if (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0) + { + internal::FEEvaluationImpl + ::integrate(*this->data, &this->values_dofs[0], this->values_quad, + this->gradients_quad, this->scratch_data, + integrate_val, integrate_grad); + } + else if (this->data->element_type == internal::MatrixFreeFunctions::tensor_general) + { + internal::FEEvaluationImpl + ::integrate(*this->data, &this->values_dofs[0], this->values_quad, + this->gradients_quad, this->scratch_data, + integrate_val, integrate_grad); + } + else if (this->data->element_type == internal::MatrixFreeFunctions::truncated_tensor) + { + internal::FEEvaluationImpl + ::integrate(*this->data, &this->values_dofs[0], this->values_quad, + this->gradients_quad, this->scratch_data, + integrate_val, integrate_grad); + } + else + AssertThrow(false, ExcNotImplemented()); + } #ifdef DEBUG this->dof_values_initialized = true; diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index 2cda877895..95b981f43b 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -151,10 +151,10 @@ namespace internal AlignedVector > shape_grad_spectral_eo; /** - * Stores the shape gradients of the spectral element space - * FE_DGQ<1>(Quadrature<1>) for faster evaluation in standard format. + * Stores the shape hessians of the spectral element space + * FE_DGQ<1>(Quadrature<1>) for faster evaluation in even-odd format. */ - AlignedVector > shape_grad_spectral; + AlignedVector > shape_hessian_spectral_eo; /** * Stores the indices from cell DoFs to face DoFs. The rows go through diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index eb6a2e96b6..9d7b32818b 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -209,23 +209,47 @@ namespace internal this->face_gradient[1][i] = fe->shape_grad(my_i,q_point)[0]; } - // get spectral evaluation points - if (fe_degree+1 == n_q_points_1d) - { - FE_DGQArbitraryNodes<1> fe(quad.get_points()); - shape_grad_spectral.resize(n_dofs_1d*n_dofs_1d); - for (unsigned int i=0; i fe(quad.get_points()); + for (unsigned int i=0; i<(fe_degree+1)/2; ++i) + for (unsigned int q=0; q