From: heister Date: Wed, 16 Jan 2013 16:26:13 +0000 (+0000) Subject: new versions of fe_q_dg0 that support embedding/restriction (thanks Daniel Arndt) X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=586a0fe87b71ecb55745e1f55990b00f922b6861;p=dealii-svn.git new versions of fe_q_dg0 that support embedding/restriction (thanks Daniel Arndt) git-svn-id: https://svn.dealii.org/trunk@28075 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/fe/fe_q_dg0.h b/deal.II/include/deal.II/fe/fe_q_dg0.h index 6c26a94714..cfca992e01 100644 --- a/deal.II/include/deal.II/fe/fe_q_dg0.h +++ b/deal.II/include/deal.II/fe/fe_q_dg0.h @@ -33,9 +33,12 @@ DEAL_II_NAMESPACE_OPEN * provides all values and derivatives of the shape functions. In case a * quadrature rule is given, the constructor creates a * TensorProductPolynomialsConst object that includes the tensor product of @p - * Lagrange polynomials with the * support points from @p points and a locally + * Lagrange polynomials with the support points from @p points and a locally * constant function. * + * Furthermore the constructor fills the @p interface_constrains, the @p + * prolongation (embedding) and the @p restriction matrices. + * *

Numbering of the degrees of freedom (DoFs)

* * The original ordering of the shape functions represented by the @@ -535,6 +538,20 @@ private: */ void initialize_constraints (); + /** + * Initialize the embedding + * matrices. Called from the + * constructor. + */ + void initialize_embedding (); + + /** + * Initialize the restriction + * matrices. Called from the + * constructor. + */ + void initialize_restriction (); + /** * Initialize the * @p unit_support_points field diff --git a/deal.II/source/fe/fe_q_dg0.cc b/deal.II/source/fe/fe_q_dg0.cc index ce21feb269..f6bae5220c 100644 --- a/deal.II/source/fe/fe_q_dg0.cc +++ b/deal.II/source/fe/fe_q_dg0.cc @@ -57,6 +57,68 @@ namespace FE_Q_DG0_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 + inline + void + zero_indices (unsigned int indices[dim]) + { + for (unsigned int d=0; d + inline + void + increment_indices (unsigned int indices[dim], + const unsigned int dofs1d) + { + ++indices[0]; + for (int d=0; d + 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)+1, + 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_DG0 (const unsigned int degree) // reinit constraints initialize_constraints (); - //TODO not working - // Reinit the vectors of restriction and - // prolongation matrices to the right sizes - // and compute the matrices - /*this->reinit_restriction_and_prolongation_matrices(); - if (dim == spacedim) - { - FETools::compute_embedding_matrices (*this, this->prolongation); - FETools::compute_projection_matrices (*this, this->restriction); - } - else - { - FE_Q_DG0 tmp (degree); - this->prolongation = tmp.prolongation; - this->restriction = tmp.restriction; - }*/ - + // 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_quad_dof_index_permutation(); } @@ -522,12 +573,13 @@ FE_Q_DG0::FE_Q_DG0 (const Quadrature<1> &points) // reinit constraints Implementation::initialize_constraints (points, *this); - // Reinit the vectors of restriction and - // prolongation matrices to the right sizes - // and compute the matrices + // 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_quad_dof_index_permutation(); + initialize_quad_dof_index_permutation(); } @@ -1617,6 +1669,297 @@ FE_Q_DG0::initialize_constraints () Implementation::initialize_constraints (points, *this); } + +template +void +FE_Q_DG0::initialize_embedding () +{ + // compute the interpolation matrices in much the same way as we do for + // the constraints. it's actually simpler here, since we don't have this + // weird renumbering stuff going on. The trick is again that the + // interpolation matrix is formed by a permutation of the indices of the + // cell matrix. The value eps is used a threshold to decide when certain + // evaluations of the Lagrange polynomials are zero or one. + const double eps = 1e-15*this->degree*dim; + + #ifdef DEBUG + // in DEBUG mode, check that the evaluation of support points (except for + // the discontinuous node) in the current numbering gives the identity + // operation + for (unsigned int i=0; idofs_per_cell-1; ++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-1; ++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, the last entry is not used + std::vector > poly_space1d = + FE_Q_DG0_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; d::isotropic_refinement; ++ref) + for (unsigned int child=0; + child::n_children(RefinementCase(ref+1)); + ++child) + { + // 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[diag_comp]; + const Point p_cell = + GeometryInfo::child_to_cell_coordinates + (p_subcell, child, RefinementCase(ref+1)); + for (unsigned int i=0; i (j_indices); + for (unsigned int j=0; jdofs_per_cell-1; j+=dofs1d) + { + unsigned int i_indices[dim]; + FE_Q_DG0_Helper::zero_indices (i_indices); + for (unsigned int i=0; idofs_per_cell-1; 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_DG0_Helper::increment_indices (i_indices, dofs1d); + } + Assert (i_indices[dim-1] == 1, ExcInternalError()); + FE_Q_DG0_Helper::increment_indices (j_indices, dofs1d); + } + + // the discontinuous node is simply mapped on the discontinuous node + // on the child cell + for (unsigned int j=0; jdofs_per_cell-1; j++) + { + this->prolongation[ref][child](j,this->dofs_per_cell-1) = 0.; + this->prolongation[ref][child](this->dofs_per_cell-1,j) = 0.; + } + this->prolongation[ref][child](this->dofs_per_cell-1, + this->dofs_per_cell-1) = 1.; + + // make sure that the row sum is one. this must be so since for this + // element, the continuous shape functions add up to one and the + // discontinuous node is mapped to the discontinuous node on the + // child cell + #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.) < eps, ExcInternalError()); + } + #endif + } +} + + + +template +void +FE_Q_DG0::initialize_restriction () +{ + // for these Lagrange interpolation polynomials, construction of the + // restriction matrices is relatively simple. the reason is that the + // interpolation points on the mother cell are (except for the case with + // arbitrary nonequidistant nodes and for the discontinuous node) always also + // interpolation points for some shape function on one or the other child, + // because we have chosen equidistant Lagrange interpolation points for the + // polynomials + // + // so the only thing we have to find out is: for each shape function on + // the mother cell, which is the child cell (possibly more than one) on + // which it is located, and which is the corresponding shape function there. + // rather than doing it for the shape functions on the mother cell, we take + // the interpolation points there + // + // note that the interpolation point of a shape function can be on the + // boundary between subcells. in that case, restriction from children to + // mother may produce two or more entries for a dof on the mother cell. + // however, this doesn't hurt: since the element is continuous (apart from + // the discontinuous node), the contribution from each child should yield + // the same result, and since the element (apart from the discontinuous + // node) is non-additive we just overwrite one value (compute one one child) + // by the same value (compute on a later child), so we don't have to care + // about this. + // for the discontinuous node we just take the mean of all the child cells' + // contributions. + + 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, the last + // entry is not used + const unsigned int dofs1d = this->degree+1; + std::vector > poly_space1d = + FE_Q_DG0_Helper::generate_poly_space1d (this->unit_support_points, + index_map_inverse, dofs1d); + std::vector > evaluations1d (dofs1d); + + for (unsigned int i=0; idofs_per_cell-1; ++i) + { + 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 point is located + for (unsigned int ref=RefinementCase::cut_x; + ref<=RefinementCase::isotropic_refinement; ++ref) + for (unsigned int child=0; + child::n_children(RefinementCase(ref)); + ++child) + { + // 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-1; j += dofs1d) + { + double val_extra_dim = 1.; + for (unsigned int d=1; drestriction[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_DG0_Helper::increment_indices (j_indices, dofs1d); + } + Assert (std::fabs(sum_check-1) < eps, + ExcInternalError()); + } + } + } + // for the discontinuous node we take the mean of all the child cells' + // contributions. so there is also a diagonal entry equal to the inverse of + // n_children. + for (unsigned int ref=RefinementCase::cut_x; + ref<=RefinementCase::isotropic_refinement; ++ref) + for (unsigned int child=0; + child::n_children(RefinementCase(ref)); + ++child) + { + for (unsigned int i=0; idofs_per_cell-1; ++i) + { + this->restriction[ref-1][child](i,this->dofs_per_cell-1)=0.; + this->restriction[ref-1][child](this->dofs_per_cell-1,i)=0.; + } + this->restriction[ref-1][child] + (this->dofs_per_cell-1,this->dofs_per_cell-1) = + 1./(double)GeometryInfo::n_children + (RefinementCase(ref)); + } + +} + + + //--------------------------------------------------------------------------- // Data field initialization //---------------------------------------------------------------------------