From: Martin Kronbichler Date: Thu, 28 May 2009 09:57:02 +0000 (+0000) Subject: Avoid some unnecessary work in FE_Q setup. Give compiler some hint in FEPolyTensor. X-Git-Tag: v8.0.0~7652 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9308945206e801cab53d1d42679da3e90cd1fc01;p=dealii.git Avoid some unnecessary work in FE_Q setup. Give compiler some hint in FEPolyTensor. git-svn-id: https://svn.dealii.org/trunk@18885 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/deal.II/source/fe/fe_poly_tensor.cc index a6cc1f226b..9abae0427a 100644 --- a/deal.II/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc @@ -330,8 +330,12 @@ FE_PolyTensor::get_data ( data->shape_values[i][k] = values[i]; else for (unsigned int i=0; idofs_per_cell; ++i) - for (unsigned int j=0; jdofs_per_cell; ++j) - data->shape_values[i][k] += inverse_node_matrix(j,i) * values[j]; + { + Tensor<1,dim> add_values; + for (unsigned int j=0; jdofs_per_cell; ++j) + add_values += inverse_node_matrix(j,i) * values[j]; + data->shape_values[i][k] = add_values; + } } if (flags & update_gradients) @@ -341,8 +345,12 @@ FE_PolyTensor::get_data ( data->shape_grads[i][k] = grads[i]; else for (unsigned int i=0; idofs_per_cell; ++i) - for (unsigned int j=0; jdofs_per_cell; ++j) - data->shape_grads[i][k] += inverse_node_matrix(j,i) * grads[j]; + { + Tensor<2,dim> add_grads; + for (unsigned int j=0; jdofs_per_cell; ++j) + add_grads += inverse_node_matrix(j,i) * grads[j]; + data->shape_grads[i][k] = add_grads; + } } } diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index e3fa057a36..8b27f0e2f0 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -1530,6 +1530,60 @@ FE_Q::initialize_embedding () const std::vector &index_map= this->poly_space.get_numbering(); + const double zero_threshold = 2e-13*this->degree*this->degree*dim; + + // precompute subcell interpolation + // matrix + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++j) + { + const Point p_subcell + = FE_Q_Helper::generate_unit_point (index_map[j], this->dofs_per_cell, + dealii::internal::int2type()); + const double + subcell_value = this->poly_space.compute_value (i, p_subcell); + + // cut off values that are + // too small. note that we + // have here Lagrange + // interpolation functions, + // so they should be zero + // at almost all points, + // and one at the others, + // at least on the + // subcells. so set them to + // their exact values + // + // the actual cut-off value + // is somewhat fuzzy, but + // it works for + // 2e-13*degree^2*dim (see + // above), which is kind of + // reasonable given that we + // compute the values of + // the polynomials via an + // degree-step recursion + // and then multiply the + // 1d-values. this gives us + // a linear growth in + // degree*dim, times a + // small constant. + if (std::fabs(subcell_value) < zero_threshold) + subcell_interpolation(j, i) = 0.; + else if (std::fabs(subcell_value-1) < zero_threshold) + subcell_interpolation(j, i) = 1.; + else + // we have put our + // evaluation + // points onto the + // interpolation + // points, so we + // should either + // get zeros or + // ones! + Assert (false, ExcInternalError()); + } + for (unsigned int ref=0; ref::isotropic_refinement; ++ref) for (unsigned int child=0; child::n_children(RefinementCase(ref+1)); ++child) { @@ -1543,66 +1597,43 @@ FE_Q::initialize_embedding () = FE_Q_Helper::generate_unit_point (index_map[j], this->dofs_per_cell, dealii::internal::int2type()); const Point p_cell = - GeometryInfo::child_to_cell_coordinates (p_subcell, child, RefinementCase(ref+1)); + GeometryInfo::child_to_cell_coordinates (p_subcell, child, + RefinementCase(ref+1)); for (unsigned int i=0; idofs_per_cell; ++i) { const double - cell_value = this->poly_space.compute_value (i, p_cell), - subcell_value = this->poly_space.compute_value (i, p_subcell); - - // cut off values that - // are too small. note - // that we have here - // Lagrange - // interpolation - // functions, so they - // should be zero at - // almost all points, - // and one at the - // others, at least on - // the subcells. so set - // them to their exact - // values + cell_value = this->poly_space.compute_value (i, p_cell); + + // cut off values that are + // too small. note that we + // have here Lagrange + // interpolation functions, + // so they should be zero + // at almost all points, + // and one at the others, + // at least on the + // subcells. so set them to + // their exact values // - // the actual cut-off - // value is somewhat - // fuzzy, but it works - // for - // 1e-14*degree*dim, - // which is kind of - // reasonable given - // that we compute the - // values of the - // polynomials via an - // degree-step - // recursion and then - // multiply the - // 1d-values. this - // gives us a linear - // growth in + // the actual cut-off value + // is somewhat fuzzy, but + // it works for + // 2e-13*degree^2*dim (see + // above), which is kind of + // reasonable given that we + // compute the values of + // the polynomials via an + // degree-step recursion + // and then multiply the + // 1d-values. this gives us + // a linear growth in // degree*dim, times a // small constant. - if (std::fabs(cell_value) < 2e-13*this->degree*this->degree*dim) + if (std::fabs(cell_value) < zero_threshold) cell_interpolation(j, i) = 0.; else cell_interpolation(j, i) = cell_value; - - if (std::fabs(subcell_value) < 2e-13*this->degree*this->degree*dim) - subcell_interpolation(j, i) = 0.; - else - if (std::fabs(subcell_value-1) < 2e-13*this->degree*this->degree*dim) - subcell_interpolation(j, i) = 1.; - else - // we have put our - // evaluation - // points onto the - // interpolation - // points, so we - // should either - // get zeros or - // ones! - Assert (false, ExcInternalError()); } } @@ -1631,7 +1662,7 @@ FE_Q::initialize_embedding () // here for (unsigned int i=0; idofs_per_cell; ++i) for (unsigned int j=0; jdofs_per_cell; ++j) - if (std::fabs(this->prolongation[ref][child](i,j)) < 2e-13*this->degree*dim) + if (std::fabs(this->prolongation[ref][child](i,j)) < zero_threshold) this->prolongation[ref][child](i,j) = 0.; // and make sure that the row @@ -1643,7 +1674,7 @@ FE_Q::initialize_embedding () double sum = 0; for (unsigned int col=0; coldofs_per_cell; ++col) sum += this->prolongation[ref][child](row,col); - Assert (std::fabs(sum-1.) < 2e-13*this->degree*this->degree*dim, + Assert (std::fabs(sum-1.) < zero_threshold, ExcInternalError()); } }