From: kronbichler Date: Tue, 26 Jul 2011 13:05:02 +0000 (+0000) Subject: Improved initialization of embedding and restriction: need to evaluate a tensor produ... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d6942cd2eaceef89daf7add6f7b2d7ac468ff4f3;p=dealii-svn.git Improved initialization of embedding and restriction: need to evaluate a tensor product polynomial a lot of times on tensor product points. Hence, do the polynomial evaluation in one dimension at a time. This makes code slightly more complicated to read, but it evaluates much faster, especially for high orders. Now, FE_Q can be initialized for order 15 in 3D in less than a second. Before it took hundreds of seconds. git-svn-id: https://svn.dealii.org/trunk@23961 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/source/fe/fe_q.cc b/deal.II/source/fe/fe_q.cc index a316bbe31b..d5328724d8 100644 --- a/deal.II/source/fe/fe_q.cc +++ b/deal.II/source/fe/fe_q.cc @@ -20,7 +20,6 @@ #include #include -#include #include DEAL_II_NAMESPACE_OPEN @@ -69,6 +68,83 @@ namespace FE_Q_Helper out[in[i]]=i; return out; } + + + + // in initialize_embedding() and + // initialize_restriction(), want to undo + // tensorization on inner loops for + // performance reasons. this clears a + // dim-array + template +#ifdef DEAL_II_ANON_NAMESPACE_BUG + static +#endif + inline + void + zero_indices (unsigned int indices[dim]) + { + for (unsigned int d=0; d +#ifdef DEAL_II_ANON_NAMESPACE_BUG + static +#endif + inline + void + increment_indices (unsigned int indices[dim], + const unsigned int dofs1d) + { + ++indices[0]; + for (int d=0; d +#ifdef DEAL_II_ANON_NAMESPACE_BUG + static +#endif + inline + std::vector > + generate_poly_space1d (const std::vector > &unit_support_points, + const std::vector &index_map_inverse, + const unsigned int dofs1d) + { + AssertDimension (Utilities::fixed_power (dofs1d), + unit_support_points.size()); + std::vector > points1d (dofs1d); + for (unsigned int i=0; i(unit_support_points[j](0)); + for (unsigned int d=1; d::FE_Q (const unsigned int degree) initialize_unit_support_points (); initialize_unit_face_support_points (); - // compute constraint, embedding - // and restriction matrices + // reinit constraints initialize_constraints (); + + // Reinit the vectors of restriction and + // prolongation matrices to the right sizes + // and compute the matrices this->reinit_restriction_and_prolongation_matrices(); - initialize_embedding (); - initialize_restriction (); + initialize_embedding(); + initialize_restriction(); initialize_quad_dof_index_permutation(); } @@ -513,20 +592,14 @@ FE_Q::FE_Q (const Quadrature<1> &points) initialize_unit_support_points (points); initialize_unit_face_support_points (points); - // compute constraint, embedding - // and restriction matrices + // reinit constraints Implementation::initialize_constraints (points, *this); - // Reinit the vectors of - // restriction and prolongation - // matrices to the right sizes + // Reinit the vectors of restriction and + // prolongation matrices to the right sizes + // and compute the matrices this->reinit_restriction_and_prolongation_matrices(); - - // Fill prolongation matrices with - // embedding operators - initialize_embedding (); - - // Fill restriction matrices + initialize_embedding(); initialize_restriction(); initialize_quad_dof_index_permutation(); @@ -958,7 +1031,7 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe for (unsigned int i=0; idofs_per_face; ++i) sum += interpolation_matrix(j,i); - Assert (std::fabs(sum-1) < 2e-13*this->degree*this->degree*dim, + Assert (std::fabs(sum-1) < 2e-13*this->degree*dim, ExcInternalError()); } } @@ -1483,54 +1556,54 @@ FE_Q::initialize_embedding () // value eps is used a threshold to // decide when certain evaluations of the // Lagrange polynomials are zero or one. - const double eps = 1e-13*this->degree*this->degree*this->degree*this->degree*dim; - - unsigned n_ones = 0; - // precompute subcell interpolation - // information, which will give us a - // vector of permutations. it actually is - // a matrix (the inverse of which we'd - // need to multiply the celL - // interpolation matrix with), but since - // we use Lagrangian basis functions - // here, we know that each basis function - // will just one at one node and zero on - // all the others. this makes this - // process much cheaper. - std::vector subcell_permutations (this->dofs_per_cell, - deal_II_numbers::invalid_unsigned_int); - for (unsigned int i=0; idofs_per_cell; ++i) - for (unsigned int j=0; jdofs_per_cell; ++j) - { - const Point p_subcell = this->unit_support_points[j]; - const double - subcell_value = this->poly_space.compute_value (i, p_subcell); + const double eps = 1e-15*this->degree*dim; - if (std::fabs(subcell_value-1) < eps) - { - subcell_permutations[i] = j; - // in debug mode, still want to check - // whether we're not getting any strange - // results with more than one 1 per row. -#ifndef DEBUG - break; -#else - n_ones++; -#endif - } - else - Assert (std::fabs(subcell_value) < eps, +#ifdef DEBUG + // in DEBUG mode, check that the evaluation of + // support points in the current numbering + // gives the identity operation + for (unsigned int i=0; idofs_per_cell; ++i) + { + Assert (std::fabs (1.-this->poly_space.compute_value + (i, this->unit_support_points[i])) < eps, + ExcInternalError()); + for (unsigned int j=0; jdofs_per_cell; ++j) + if (j!=i) + Assert (std::fabs (this->poly_space.compute_value + (i, this->unit_support_points[j])) < eps, ExcInternalError()); + } +#endif + + // to efficiently evaluate the polynomial at + // the subcell, make use of the tensor product + // structure of this element and only evaluate + // 1D information from the polynomial. This + // makes the cost of this function almost + // negligible also for high order elements + const unsigned int dofs1d = this->degree+1; + std::vector > + subcell_evaluations (dim, Table<2,double>(dofs1d, dofs1d)); + const std::vector &index_map_inverse = + this->poly_space.get_numbering_inverse(); + + // recreate 1D polynomials + std::vector > poly_space1d = + FE_Q_Helper::generate_poly_space1d (this->unit_support_points, + index_map_inverse, dofs1d); + + // helper value: step size how to walk through + // diagonal and how many points we have left + // apart from the first dimension + unsigned int step_size_diag = 0; + { + unsigned int factor = 1; + for (unsigned int d=0; ddofs_per_cell, - ExcDimensionMismatch(n_ones, this->dofs_per_cell)); - for (unsigned int i=0; idofs_per_cell; ++i) - Assert (subcell_permutations[i] < this->dofs_per_cell, - ExcInternalError()); + } // next evaluate the functions // for the different refinement @@ -1538,20 +1611,20 @@ FE_Q::initialize_embedding () for (unsigned int ref=0; ref::isotropic_refinement; ++ref) for (unsigned int child=0; child::n_children(RefinementCase(ref+1)); ++child) { - for (unsigned int j=0; jdofs_per_cell; ++j) + // go through the points in diagonal to + // capture variation in all directions + // simultaneously + for (unsigned int j=0; j p_subcell = this->unit_support_points[j]; + const unsigned int diag_comp = index_map_inverse[j*step_size_diag]; + const Point p_subcell = this->unit_support_points[diag_comp]; const Point p_cell = 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); + for (unsigned int i=0; i::initialize_embedding () // // the actual cut-off value is somewhat // fuzzy, but it works for - // 2e-13*degree^2*dim (see above), which + // 2e-13*degree*dim (see above), which // is kind of reasonable given that we // compute the values of the polynomials // via an degree-step recursion and then @@ -1579,25 +1652,63 @@ FE_Q::initialize_embedding () // all we need to do is to switch the // rows we write the data into. moreover, // cut off very small values here - if (std::fabs(cell_value) < eps) - this->prolongation[ref][child](subcell_permutations[j],i) = 0; - else - this->prolongation[ref][child](subcell_permutations[j],i) = - cell_value; + if (std::fabs(cell_value) < eps) + subcell_evaluations[d](j,i) = 0; + else + subcell_evaluations[d](j,i) = cell_value; + } + } + + // now expand from 1D info. block innermost + // dimension (x_0) in order to avoid difficult + // checks at innermost loop + unsigned int j_indices[dim]; + FE_Q_Helper::zero_indices (j_indices); + for (unsigned int j=0; jdofs_per_cell; j+=dofs1d) + { + unsigned int i_indices[dim]; + FE_Q_Helper::zero_indices (i_indices); + for (unsigned int i=0; idofs_per_cell; i+=dofs1d) + { + double val_extra_dim = 1.; + for (unsigned int d=1; dprolongation[ref][child](j,i) = + // this->poly_space.compute_value (i, p_cell); + for (unsigned int jj=0; jjprolongation[ref][child](j_ind,index_map_inverse[i+ii]) + = val_extra_dim * subcell_evaluations[0](jj,ii); + } + + // update indices that denote the tensor + // product position. a bit fuzzy and therefore + // not done for innermost x_0 direction + FE_Q_Helper::increment_indices (i_indices, dofs1d); } + Assert (i_indices[dim-1] == 1, ExcInternalError()); + FE_Q_Helper::increment_indices (j_indices, dofs1d); } // and make sure that the row sum is // 1. this must be so since for this // element, the shape functions add up to // one +#ifdef DEBUG for (unsigned int row=0; rowdofs_per_cell; ++row) { double sum = 0; for (unsigned int col=0; coldofs_per_cell; ++col) sum += this->prolongation[ref][child](row,col); - Assert (std::fabs(sum-1.) < this->degree*eps, ExcInternalError()); + Assert (std::fabs(sum-1.) < eps, ExcInternalError()); } +#endif } } @@ -1630,10 +1741,7 @@ FE_Q::initialize_restriction () // function there. rather than // doing it for the shape functions // on the mother cell, we take the - // interpolation points there are - // also search which shape function - // corresponds to it (too lazy to - // do this mapping by hand) + // interpolation points there // // note that the interpolation // point of a shape function can be @@ -1653,33 +1761,22 @@ FE_Q::initialize_restriction () // (compute on a later child), so // we don't have to care about this - const double eps = 1e-13*this->degree*this->degree*this->degree*this->degree*dim; - const std::vector &index_map_inverse= + const double eps = 1e-15*this->degree*dim; + const std::vector &index_map_inverse = this->poly_space.get_numbering_inverse(); + + // recreate 1D polynomials for faster + // evaluation of polynomial + const unsigned int dofs1d = this->degree+1; + std::vector > poly_space1d = + FE_Q_Helper::generate_poly_space1d (this->unit_support_points, + index_map_inverse, dofs1d); + std::vector > evaluations1d (dofs1d); + for (unsigned int i=0; idofs_per_cell; ++i) { - const Point p_cell = this->unit_support_points[index_map_inverse[i]]; - unsigned int mother_dof = 0; - for (; mother_dofdofs_per_cell; ++mother_dof) - { - const double val - = this->poly_space.compute_value(mother_dof, p_cell); - if (std::fabs (val-1.) < eps) - // ok, this is the right - // dof - break; - else - // make sure that all - // other shape functions - // are zero there - Assert (std::fabs(val) < eps, ExcInternalError()); - } - // check also the shape - // functions after that - for (unsigned int j=mother_dof+1; jdofs_per_cell; ++j) - Assert (std::fabs (this->poly_space.compute_value(j, p_cell)) - < eps, - ExcInternalError()); + unsigned int mother_dof = index_map_inverse[i]; + const Point p_cell = this->unit_support_points[mother_dof]; // then find the children on // which the interpolation @@ -1687,18 +1784,33 @@ FE_Q::initialize_restriction () for (unsigned int ref=RefinementCase::cut_x; ref<=RefinementCase::isotropic_refinement; ++ref) for (unsigned int child=0; child::n_children(RefinementCase(ref)); ++child) { - // first initialize this - // column of the matrix - for (unsigned int j=0; jdofs_per_cell; ++j) - this->restriction[ref-1][child](mother_dof, j) = 0.; - - // then check whether this + // check whether this // interpolation point is // inside this child cell const Point p_subcell = GeometryInfo::cell_to_child_coordinates (p_cell, child, RefinementCase(ref)); if (GeometryInfo::is_inside_unit_cell (p_subcell)) { + // same logic as in initialize_embedding to + // evaluate the polynomial faster than from + // the tensor product: since we evaluate all + // polynomials, it is much faster to just + // compute the 1D values for all polynomials + // before and then get the dim-data. + for (unsigned int j=0; j (j_indices); + double sum_check = 0; + for (unsigned int j = 0; jdofs_per_cell; j += dofs1d) + { + double val_extra_dim = 1.; + for (unsigned int d=1; d::initialize_restriction () // get more than one nonzero // value per row. Still, the // values should sum up to 1 - double sum_check = 0; - for (unsigned int child_dof = 0; child_dofdofs_per_cell; - ++child_dof) - { - const double val - = this->poly_space.compute_value(child_dof, p_subcell); - if (std::fabs (val-1.) < eps) - this-> restriction[ref-1][child](mother_dof,child_dof)=1.; - else if(std::fabs(val) > eps) - this->restriction[ref-1][child](mother_dof,child_dof)= val; - - sum_check += val; + const double val + = val_extra_dim * evaluations1d[jj][0]; + const unsigned int child_dof = + index_map_inverse[j+jj]; + if (std::fabs (val-1.) < eps) + this->restriction[ref-1][child](mother_dof,child_dof)=1.; + else if(std::fabs(val) > eps) + this->restriction[ref-1][child](mother_dof,child_dof)=val; + sum_check += val; + } + FE_Q_Helper::increment_indices (j_indices, dofs1d); } - Assert (std::fabs(sum_check-1) < this->degree*eps, + Assert (std::fabs(sum_check-1) < eps, ExcInternalError()); } }