From 5e4c1296f0515e6bad6c68908b9583bddd3d2758 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Tue, 9 Dec 2008 07:51:24 +0000 Subject: [PATCH] The previous implementation of FEValues::get_function_values/gradients/hessians() was rather inefficient. The submitted version is more than twice as fast for get_function_values and 20% faster for get_function_hessians. Changed the order of the loops, read FEValues data using pointers. git-svn-id: https://svn.dealii.org/trunk@17903 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/fe_values.cc | 545 +++++++++++++++---------- 1 file changed, 335 insertions(+), 210 deletions(-) diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index af0ae8f103..fcfa232d29 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -546,10 +546,13 @@ void FEValuesBase::get_function_values ( // elements, so no need to check // for non-primitivity of shape // functions - for (unsigned int point=0; pointshape_value(shape_func, point)); + for (unsigned int shape_func=0; shape_funcshape_values(shape_func, 0); + for (unsigned int point=0; point::get_function_values ( // elements, so no need to check // for non-primitivity of shape // functions - for (unsigned int point=0; pointshape_value(shape_func, point)); + for (unsigned int shape_func=0; shape_funcshape_values(shape_func, 0); + for (unsigned int point=0; point::get_function_values ( // not. if it is, then set its only // non-zero component, otherwise // loop over components - for (unsigned int point=0; pointis_primitive(shape_func)) - values[point](fe->system_to_component_index(shape_func).first) - += (dof_values(shape_func) * shape_value(shape_func, point)); + { + const double *shape_value_ptr = &this->shape_values(shape_func, 0); + const unsigned int comp = fe->system_to_component_index(shape_func).first; + for (unsigned int point=0; pointget_nonzero_components(shape_func)[c] == false) + continue; + + const unsigned int + row = (this->shape_function_to_row_table[shape_func] + + + std::count (fe->get_nonzero_components(shape_func).begin(), + fe->get_nonzero_components(shape_func).begin()+c, + true)); + + const double *shape_value_ptr = &this->shape_values(row, 0); + + for (unsigned int point=0; point::get_function_values ( // non-zero component, otherwise // loop over components for (unsigned int mc = 0; mc < component_multiple; ++mc) - for (unsigned int point=0; pointis_primitive(shape_func)) - values[point](fe->system_to_component_index(shape_func).first - +mc * n_components) - += (fe_function(indices[shape_func+mc*dofs_per_cell]) - * shape_value(shape_func, point)); + { + const double *shape_value_ptr = &this->shape_values(shape_func, 0); + const unsigned int comp = fe->system_to_component_index(shape_func).first + + mc * n_components; + for (unsigned int point=0; pointget_nonzero_components(shape_func)[c] == false) + continue; + + const unsigned int + row = (this->shape_function_to_row_table[shape_func] + + + std::count (fe->get_nonzero_components(shape_func).begin(), + fe->get_nonzero_components(shape_func).begin()+c, + true)); + + const double *shape_value_ptr = &this->shape_values(row, 0); + const unsigned int comp = c + mc * n_components; + + for (unsigned int point=0; point::get_function_values ( // not. if it is, then set its only // non-zero component, otherwise // loop over components - if (quadrature_points_fastest) - for (unsigned int mc = 0; mc < component_multiple; ++mc) - for (unsigned int point=0; pointis_primitive(shape_func)) - values[fe->system_to_component_index(shape_func).first - +mc * n_components][point] - += (fe_function(indices[shape_func+mc*dofs_per_cell]) - * shape_value(shape_func, point)); - else - for (unsigned int c=0; cis_primitive(shape_func)) - values[point][fe->system_to_component_index(shape_func).first - +mc * n_components] - += (fe_function(indices[shape_func+mc*dofs_per_cell]) - * shape_value(shape_func, point)); - else - for (unsigned int c=0; cis_primitive(shape_func)) + { + const double *shape_value_ptr = &this->shape_values(shape_func, 0); + const unsigned int comp = fe->system_to_component_index(shape_func).first + + mc * n_components; + if (quadrature_points_fastest) + for (unsigned int point=0; pointget_nonzero_components(shape_func)[c] == false) + continue; + + const unsigned int + row = (this->shape_function_to_row_table[shape_func] + + + std::count (fe->get_nonzero_components(shape_func).begin(), + fe->get_nonzero_components(shape_func).begin()+c, + true)); + + const double *shape_value_ptr = &this->shape_values(row, 0); + const unsigned int comp = c + mc * n_components; + + if (quadrature_points_fastest) + for (unsigned int point=0; point::get_function_gradients ( // elements, so no need to check // for non-primitivity of shape // functions - for (unsigned int point=0; point tmp = this->shape_grad(shape_func,point); - tmp *= dof_values(shape_func); - gradients[point] += tmp; - }; + for (unsigned int shape_func=0; shape_func *shape_gradient_ptr = &this->shape_gradients[shape_func][0]; + for (unsigned int point=0; point void FEValuesBase::get_function_gradients ( const InputVector& fe_function, const VectorSlice >& indices, - std::vector > &values) const + std::vector > &gradients) const { Assert (this->update_flags & update_gradients, ExcAccessToUninitializedField()); // This function fills a single @@ -869,11 +929,11 @@ void FEValuesBase::get_function_gradients ( ExcDimensionMismatch(indices.size(), dofs_per_cell)); // This vector has one entry for // each quadrature point - Assert (values.size() == n_quadrature_points, - ExcDimensionMismatch(values.size(), n_quadrature_points)); + Assert (gradients.size() == n_quadrature_points, + ExcDimensionMismatch(gradients.size(), n_quadrature_points)); // initialize with zero - std::fill_n (values.begin(), n_quadrature_points, Tensor<1,spacedim>()); + std::fill_n (gradients.begin(), n_quadrature_points, Tensor<1,spacedim>()); // add up contributions of trial // functions. note that here we @@ -881,10 +941,13 @@ void FEValuesBase::get_function_gradients ( // elements, so no need to check // for non-primitivity of shape // functions - for (unsigned int point=0; pointshape_grad(shape_func, point)); + for (unsigned int shape_func=0; shape_func *shape_gradient_ptr = &this->shape_gradients[shape_func][0]; + for (unsigned int point=0; point::get_function_gradients ( // not. if it is, then set its only // non-zero component, otherwise // loop over components - for (unsigned int point=0; pointis_primitive (shape_func)) + for (unsigned int shape_func=0; shape_funcis_primitive(shape_func)) { - Tensor<1,spacedim> tmp = this->shape_grad(shape_func,point); - tmp *= dof_values(shape_func); - gradients[point][fe->system_to_component_index(shape_func).first] - += tmp; + const Tensor<1,spacedim> *shape_gradient_ptr + = &this->shape_gradients[shape_func][0]; + const unsigned int comp = fe->system_to_component_index(shape_func).first; + for (unsigned int point=0; point tmp = this->shape_grad_component(shape_func,point,c); - tmp *= dof_values(shape_func); - gradients[point][c] += tmp; + if (fe->get_nonzero_components(shape_func)[c] == false) + continue; + + const unsigned int + row = (this->shape_function_to_row_table[shape_func] + + + std::count (fe->get_nonzero_components(shape_func).begin(), + fe->get_nonzero_components(shape_func).begin()+c, + true)); + + const Tensor<1,spacedim> *shape_gradient_ptr + = &this->shape_gradients[row][0]; + + for (unsigned int point=0; point void FEValuesBase::get_function_gradients ( const InputVector& fe_function, const VectorSlice >& indices, - std::vector > >& values, + std::vector > >& gradients, bool quadrature_points_fastest) const { const unsigned int n_components = fe->n_components(); @@ -975,19 +1054,19 @@ void FEValuesBase::get_function_gradients ( // initialized to the correct sizes if (quadrature_points_fastest) { - Assert (values.size() == result_components, - ExcDimensionMismatch(values.size(), result_components)); - for (unsigned i=0;i::get_function_gradients ( Assert (this->update_flags & update_values, ExcAccessToUninitializedField()); // initialize with zero - for (unsigned i=0;i()); + for (unsigned i=0;i()); // add up contributions of trial // functions. now check whether the @@ -1008,32 +1087,50 @@ void FEValuesBase::get_function_gradients ( // not. if it is, then set its only // non-zero component, otherwise // loop over components - if (quadrature_points_fastest) - for (unsigned int mc = 0; mc < component_multiple; ++mc) - for (unsigned int point=0; pointis_primitive(shape_func)) - values[fe->system_to_component_index(shape_func).first - +mc * n_components][point] - += fe_function(indices[shape_func+mc*dofs_per_cell]) - * shape_grad(shape_func, point); - else - for (unsigned int c=0; cis_primitive(shape_func)) - values[point][fe->system_to_component_index(shape_func).first - +mc * n_components] - += fe_function(indices[shape_func+mc*dofs_per_cell]) - * shape_grad(shape_func, point); + for (unsigned int mc = 0; mc < component_multiple; ++mc) + for (unsigned int shape_func=0; shape_funcis_primitive(shape_func)) + { + const Tensor<1,spacedim> *shape_gradient_ptr + = &this->shape_gradients[shape_func][0]; + const unsigned int comp = fe->system_to_component_index(shape_func).first + + mc * n_components; + + if (quadrature_points_fastest) + for (unsigned int point=0; pointget_nonzero_components(shape_func)[c] == false) + continue; + + const unsigned int + row = (this->shape_function_to_row_table[shape_func] + + + std::count (fe->get_nonzero_components(shape_func).begin(), + fe->get_nonzero_components(shape_func).begin()+c, + true)); + + const Tensor<1,spacedim> *shape_gradient_ptr + = &this->shape_gradients[row][0]; + const unsigned int comp = c + mc * n_components; + + if (quadrature_points_fastest) + for (unsigned int point=0; point void FEValuesBase:: get_function_hessians (const InputVector &fe_function, - std::vector > &hessians) const + std::vector > &hessians) const { Assert (fe->n_components() == 1, ExcDimensionMismatch(fe->n_components(), 1)); @@ -1070,13 +1167,13 @@ get_function_hessians (const InputVector &fe_function, // elements, so no need to check // for non-primitivity of shape // functions - for (unsigned int point=0; point tmp = this->shape_hessian(shape_func,point); - tmp *= dof_values(shape_func); - hessians[point] += tmp; - }; + for (unsigned int shape_func=0; shape_func *shape_hessians_ptr = &this->shape_hessians[shape_func][0]; + for (unsigned int point=0; point void FEValuesBase::get_function_hessians ( const InputVector& fe_function, const VectorSlice >& indices, - std::vector > &values) const + std::vector > &hessians) const { Assert (this->update_flags & update_second_derivatives, ExcAccessToUninitializedField()); // This function fills a single @@ -1098,11 +1195,11 @@ void FEValuesBase::get_function_hessians ( ExcDimensionMismatch(indices.size(), dofs_per_cell)); // This vector has one entry for // each quadrature point - Assert (values.size() == n_quadrature_points, - ExcDimensionMismatch(values.size(), n_quadrature_points)); + Assert (hessians.size() == n_quadrature_points, + ExcDimensionMismatch(hessians.size(), n_quadrature_points)); // initialize with zero - std::fill_n (values.begin(), n_quadrature_points, Tensor<2,spacedim>()); + std::fill_n (hessians.begin(), n_quadrature_points, Tensor<2,spacedim>()); // add up contributions of trial // functions. note that here we @@ -1110,10 +1207,13 @@ void FEValuesBase::get_function_hessians ( // elements, so no need to check // for non-primitivity of shape // functions - for (unsigned int point=0; pointshape_hessian(shape_func, point)); + for (unsigned int shape_func=0; shape_func *shape_hessians_ptr = &this->shape_hessians[shape_func][0]; + for (unsigned int point=0; point void FEValuesBase:: get_function_hessians (const InputVector &fe_function, - std::vector > > &second_derivs, - bool quadrature_points_fastest) const + std::vector > > &hessians, + bool quadrature_points_fastest) const { - Assert (n_quadrature_points == second_derivs.size(), - ExcDimensionMismatch(second_derivs.size(), n_quadrature_points)); + Assert (n_quadrature_points == hessians.size(), + ExcDimensionMismatch(hessians.size(), n_quadrature_points)); const unsigned int n_components = fe->n_components(); - for (unsigned i=0;iupdate_flags & update_hessians, ExcAccessToUninitializedField()); Assert (present_cell.get() != 0, @@ -1148,55 +1248,62 @@ get_function_hessians (const InputVector &fe_function, present_cell->get_interpolated_dof_values(fe_function, dof_values); // initialize with zero - for (unsigned i=0;i()); + for (unsigned i=0;i()); // add up contributions of trial // functions - if (quadrature_points_fastest) - for (unsigned int point=0; pointis_primitive(shape_func)) - { - Tensor<2,spacedim> tmp(shape_hessian(shape_func,point)); - tmp *= dof_values(shape_func); - second_derivs[fe->system_to_component_index(shape_func).first][point] - += tmp; - } - else - for (unsigned int c=0; c tmp = this->shape_hessian_component(shape_func,point,c); - tmp *= dof_values(shape_func); - second_derivs[c][point] += tmp; - } - else - for (unsigned int point=0; pointis_primitive(shape_func)) + for (unsigned int shape_func=0; shape_funcis_primitive(shape_func)) + { + const Tensor<2,spacedim> *shape_hessian_ptr + = &this->shape_hessians[shape_func][0]; + const unsigned int comp = fe->system_to_component_index(shape_func).first; + + if (quadrature_points_fastest) + for (unsigned int point=0; point tmp(shape_hessian(shape_func,point)); - tmp *= dof_values(shape_func); - second_derivs[point][fe->system_to_component_index(shape_func).first] - += tmp; + if (fe->get_nonzero_components(shape_func)[c] == false) + continue; + + const unsigned int + row = (this->shape_function_to_row_table[shape_func] + + + std::count (fe->get_nonzero_components(shape_func).begin(), + fe->get_nonzero_components(shape_func).begin()+c, + true)); + + const Tensor<2,spacedim> *shape_hessian_ptr + = &this->shape_hessians[row][0]; + + if (quadrature_points_fastest) + for (unsigned int point=0; point tmp = this->shape_hessian_component(shape_func,point,c); - tmp *= dof_values(shape_func); - second_derivs[point][c] += tmp; - } + } } template template -void FEValuesBase::get_function_hessians ( +void FEValuesBase::get_function_hessians ( const InputVector& fe_function, const VectorSlice >& indices, - std::vector > >& values, + std::vector > >& hessians, bool quadrature_points_fastest) const { Assert (this->update_flags & update_second_derivatives, ExcAccessToUninitializedField()); @@ -1221,19 +1328,19 @@ void FEValuesBase::get_function_hessians ( // initialized to the correct sizes if (quadrature_points_fastest) { - Assert (values.size() == result_components, - ExcDimensionMismatch(values.size(), result_components)); - for (unsigned i=0;i::get_function_hessians ( const unsigned int component_multiple = result_components / n_components; // initialize with zero - for (unsigned i=0;i()); + for (unsigned i=0;i()); // add up contributions of trial // functions. now check whether the @@ -1252,32 +1359,50 @@ void FEValuesBase::get_function_hessians ( // not. if it is, then set its only // non-zero component, otherwise // loop over components - if (quadrature_points_fastest) - for (unsigned int mc = 0; mc < component_multiple; ++mc) - for (unsigned int point=0; pointis_primitive(shape_func)) - values[fe->system_to_component_index(shape_func).first - +mc * n_components][point] - += fe_function(indices[shape_func+mc*dofs_per_cell]) - * shape_hessian(shape_func, point); - else - for (unsigned int c=0; cis_primitive(shape_func)) - values[point][fe->system_to_component_index(shape_func).first - +mc * n_components] - += fe_function(indices[shape_func+mc*dofs_per_cell]) - * shape_hessian(shape_func, point); + for (unsigned int mc = 0; mc < component_multiple; ++mc) + for (unsigned int shape_func=0; shape_funcis_primitive(shape_func)) + { + const Tensor<2,spacedim> *shape_hessian_ptr + = &this->shape_hessians[shape_func][0]; + const unsigned int comp = fe->system_to_component_index(shape_func).first + + mc * n_components; + + if (quadrature_points_fastest) + for (unsigned int point=0; pointget_nonzero_components(shape_func)[c] == false) + continue; + + const unsigned int + row = (this->shape_function_to_row_table[shape_func] + + + std::count (fe->get_nonzero_components(shape_func).begin(), + fe->get_nonzero_components(shape_func).begin()+c, + true)); + + const Tensor<2,spacedim> *shape_hessian_ptr + = &this->shape_hessians[row][0]; + const unsigned int comp = c + mc * n_components; + + if (quadrature_points_fastest) + for (unsigned int point=0; point