From: kronbichler Date: Wed, 15 May 2013 08:21:03 +0000 (+0000) Subject: Simplify get_interpolation_matrix X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=363b21c7fc3131914054e671cfd9aba302047927;p=dealii-svn.git Simplify get_interpolation_matrix git-svn-id: https://svn.dealii.org/trunk@29516 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/source/fe/fe_q_base.cc b/deal.II/source/fe/fe_q_base.cc index 331131daf3..60c7fd251c 100644 --- a/deal.II/source/fe/fe_q_base.cc +++ b/deal.II/source/fe/fe_q_base.cc @@ -449,11 +449,8 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, FullMatrix &interpolation_matrix) const { // this is only implemented, if the source FE is also a Q element - typedef FE_Q_Base FEQ; - typedef FiniteElement FEL; - - AssertThrow (dynamic_cast(&x_source_fe) != 0, - typename FEL::ExcInterpolationNotImplemented()); + AssertThrow ((dynamic_cast *>(&x_source_fe) != 0), + (typename FiniteElement::ExcInterpolationNotImplemented())); Assert (interpolation_matrix.m() == this->dofs_per_cell, ExcDimensionMismatch (interpolation_matrix.m(), @@ -466,42 +463,41 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, const FE_Q_Base &source_fe = dynamic_cast&>(x_source_fe); - // compute the interpolation matrices in much the same way as we do for the - // embedding matrices from mother to child. - FullMatrix cell_interpolation (this->dofs_per_cell, - this->dofs_per_cell); - FullMatrix source_interpolation (this->dofs_per_cell, - source_fe.dofs_per_cell); - // only evaluate Q dofs const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); const unsigned int source_q_dofs_per_cell = Utilities::fixed_power(source_fe.degree+1); + + // evaluation is simply done by evaluating the other FE's basis functions on + // the unit support points (FE_Q has the property that the cell + // interpolation matrix is a unit matrix, so no need to evaluate it and + // invert it) for (unsigned int j=0; j p = this->unit_support_points[j]; - for (unsigned int i=0; ipoly_space.compute_value (i, p); + + // FE_Q element evaluates to 1 in unit support point and to zero in all + // other points by construction + Assert(std::abs(this->poly_space.compute_value (j, p)-1.)<1e-13, + ExcInternalError()); for (unsigned int i=0; idofs_per_cell) { AssertDimension(source_q_dofs_per_cell+1, source_fe.dofs_per_cell); - cell_interpolation(q_dofs_per_cell, q_dofs_per_cell) = 1.; - source_interpolation(q_dofs_per_cell, source_q_dofs_per_cell) = 1.; + for (unsigned int i=0; idegree*dim; - // cut off very small values + const double eps = 2e-13*this->degree*dim; for (unsigned int i=0; idofs_per_cell; ++i) for (unsigned int j=0; j &x_source_fe, Assert (dim > 1, ExcImpossibleInDim(1)); // this is only implemented, if the source FE is also a Q element - typedef FE_Q_Base FEQ; - typedef FiniteElement FEL; - AssertThrow (dynamic_cast(&x_source_fe) != 0, - typename FEL::ExcInterpolationNotImplemented()); + AssertThrow ((dynamic_cast *>(&x_source_fe) != 0), + (typename FiniteElement:: + ExcInterpolationNotImplemented())); Assert (interpolation_matrix.n() == this->dofs_per_face, ExcDimensionMismatch (interpolation_matrix.n(), @@ -552,8 +547,8 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, // in that case might lead to problems in the hp procedures, which use this // method. Assert (this->dofs_per_face <= source_fe.dofs_per_face, - typename FEL:: - ExcInterpolationNotImplemented ()); + (typename FiniteElement:: + ExcInterpolationNotImplemented ())); // generate a quadrature with the unit support points. This is later based // as a basis for the QProjector, which returns the support points on the diff --git a/deal.II/source/fe/fe_tools.cc b/deal.II/source/fe/fe_tools.cc index 0589b9a3b4..b29fd46f7d 100644 --- a/deal.II/source/fe/fe_tools.cc +++ b/deal.II/source/fe/fe_tools.cc @@ -2630,22 +2630,15 @@ namespace FETools const unsigned int dofs_per_cell = fe.dofs_per_cell; // polynomial degree const unsigned int degree = fe.dofs_per_line+1; - // number of grid points in each - // direction + // number of grid points in each direction const unsigned int n = degree+1; - // the following lines of code are - // somewhat odd, due to the way the - // hierarchic numbering is - // organized. if someone would - // really want to understand these - // lines, you better draw some - // pictures where you indicate the - // indices and orders of vertices, - // lines, etc, along with the - // numbers of the degrees of - // freedom in hierarchical and - // lexicographical order + // the following lines of code are somewhat odd, due to the way the + // hierarchic numbering is organized. if someone would really want to + // understand these lines, you better draw some pictures where you + // indicate the indices and orders of vertices, lines, etc, along with the + // numbers of the degrees of freedom in hierarchical and lexicographical + // order switch (dim) { case 1: