From: kronbichler Date: Mon, 11 Jul 2011 12:14:39 +0000 (+0000) Subject: Fix bug in FE_Q with arbitrary nodes: the hp line/quad identities need to be computed... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=adf7b5e5c56124098ec6debfe408889a77a28ce6;p=dealii-svn.git Fix bug in FE_Q with arbitrary nodes: the hp line/quad identities need to be computed from support points, not equidistant points. Also, the embedding matrices for arbitrary nodes do not need to be computed with FE_Tools. They can be generated in the usual case if we use the support points directly and not equidistant points. git-svn-id: https://svn.dealii.org/trunk@23939 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/source/fe/fe_q.cc b/deal.II/source/fe/fe_q.cc index 3a5088109b..a316bbe31b 100644 --- a/deal.II/source/fe/fe_q.cc +++ b/deal.II/source/fe/fe_q.cc @@ -69,107 +69,6 @@ namespace FE_Q_Helper out[in[i]]=i; return out; } - - - // given an integer N, compute its - // integer square root (if it - // exists, otherwise give up) -#if defined(DEAL_II_ANON_NAMESPACE_BUG) && defined(DEAL_II_ANON_NAMESPACE_LINKAGE_BUG) - static -#endif - inline - unsigned int int_sqrt (const unsigned int N) - { - for (unsigned int i=0; i<=N; ++i) - if (i*i == N) - return i; - Assert (false, ExcInternalError()); - return numbers::invalid_unsigned_int; - } - - - // given an integer N, compute its - // integer cube root (if it - // exists, otherwise give up) -#if defined(DEAL_II_ANON_NAMESPACE_BUG) && defined(DEAL_II_ANON_NAMESPACE_LINKAGE_BUG) - static -#endif - inline - unsigned int int_cuberoot (const unsigned int N) - { - for (unsigned int i=0; i<=N; ++i) - if (i*i*i == N) - return i; - Assert (false, ExcInternalError()); - return numbers::invalid_unsigned_int; - } - - - // given N, generate i=0...N-1 - // equidistant points in the - // interior of the interval [0,1] -#if defined(DEAL_II_ANON_NAMESPACE_BUG) && defined(DEAL_II_ANON_NAMESPACE_LINKAGE_BUG) - static -#endif - inline - Point<1> - generate_unit_point (const unsigned int i, - const unsigned int N, - const dealii::internal::int2type<1> ) - { - Assert (i(i*h); - } - - - // given N, generate i=0...N-1 - // equidistant points in the domain - // [0,1]^2 -#if defined(DEAL_II_ANON_NAMESPACE_BUG) && defined(DEAL_II_ANON_NAMESPACE_LINKAGE_BUG) - static -#endif - inline - Point<2> - generate_unit_point (const unsigned int i, - const unsigned int N, - const dealii::internal::int2type<2> ) - { - Assert (i=4, ExcInternalError()); - - const unsigned int N1d = int_sqrt(N); - const double h = 1./(N1d-1); - - return Point<2> (i%N1d * h, - i/N1d * h); - } - - - - // given N, generate i=0...N-1 - // equidistant points in the domain - // [0,1]^3 -#if defined(DEAL_II_ANON_NAMESPACE_BUG) && defined(DEAL_II_ANON_NAMESPACE_LINKAGE_BUG) - static -#endif - inline - Point<3> - generate_unit_point (const unsigned int i, - const unsigned int N, - const dealii::internal::int2type<3> ) - { - Assert (i=8, ExcInternalError()); - - const unsigned int N1d = int_cuberoot(N); - const double h = 1./(N1d-1); - - return Point<3> (i%N1d * h, - (i/N1d)%N1d * h, - i/(N1d*N1d) * h); - } - } } @@ -625,7 +524,7 @@ FE_Q::FE_Q (const Quadrature<1> &points) // Fill prolongation matrices with // embedding operators - FETools::compute_embedding_matrices (*this, this->prolongation); + initialize_embedding (); // Fill restriction matrices initialize_restriction(); @@ -1134,26 +1033,37 @@ hp_line_dof_identities (const FiniteElement &fe_other) const // FE_Nothing if (const FE_Q *fe_q_other = dynamic_cast*>(&fe_other)) { - // dofs are located along - // lines, so two dofs are - // identical if they are - // located at identical - // positions. note that for - // elements of orders p and q, - // nodes are located on edges - // at locations (i+1)/p and - // (j+1)/q, so we need to find - // those combinations i,j for - // which (i+1)/p == (j+1)/q, - // i.e. (i+1)*q == (j+1)*p + // dofs are located along lines, so two + // dofs are identical if they are + // located at identical positions. if + // we had only equidistant points, we + // could simple check for similarity + // like (i+1)*q == (j+1)*p, but we + // might have other support points + // (e.g. Gauss-Lobatto + // points). Therefore, read the points + // in unit_support_points for the first + // coordinate direction. We take the + // lexicographic ordering of the points + // in the first direction (i.e., + // x-direction), which we access + // between index 1 and p-1 (index 0 and + // p are vertex dofs). const unsigned int p = this->degree; const unsigned int q = fe_q_other->degree; std::vector > identities; + const std::vector &index_map_inverse= + this->poly_space.get_numbering_inverse(); + const std::vector &index_map_inverse_other= + fe_q_other->poly_space.get_numbering_inverse(); + for (unsigned int i=0; iunit_support_points[index_map_inverse[i+1]][0]- + fe_q_other->unit_support_points[index_map_inverse_other[j+1]][0]) + < 1e-14) identities.push_back (std::make_pair(i,j)); return identities; @@ -1187,26 +1097,35 @@ hp_quad_dof_identities (const FiniteElement &fe_other) cons if (const FE_Q *fe_q_other = dynamic_cast*>(&fe_other)) { // this works exactly like the line - // case above, except that now we - // have to have two indices i1, i2 - // and j1, j2 to characterize the - // dofs on the face of each of the - // finite elements. since they are - // ordered in lexicographic order, - // the rest is rather - // straightforward + // case above, except that now we have + // to have two indices i1, i2 and j1, + // j2 to characterize the dofs on the + // face of each of the finite + // elements. since they are ordered + // lexicographically along the first + // line and we have a tensor product, + // the rest is rather straightforward const unsigned int p = this->degree; const unsigned int q = fe_q_other->degree; std::vector > identities; + const std::vector &index_map_inverse= + this->poly_space.get_numbering_inverse(); + const std::vector &index_map_inverse_other= + fe_q_other->poly_space.get_numbering_inverse(); + for (unsigned int i1=0; i1unit_support_points[index_map_inverse[i1+1]][0]- + fe_q_other->unit_support_points[index_map_inverse_other[j1+1]][0]) + < 1e-14) && - ((i2+1)*q == (j2+1)*p)) + (std::fabs(this->unit_support_points[index_map_inverse[i2+1]][0]- + fe_q_other->unit_support_points[index_map_inverse_other[j2+1]][0]) + < 1e-14)) identities.push_back (std::make_pair(i1*(p-1)+i2, j1*(q-1)+j2)); @@ -1564,9 +1483,6 @@ FE_Q::initialize_embedding () // value eps is used a threshold to // decide when certain evaluations of the // Lagrange polynomials are zero or one. - const std::vector &index_map= - this->poly_space.get_numbering(); - const double eps = 1e-13*this->degree*this->degree*this->degree*this->degree*dim; unsigned n_ones = 0; @@ -1586,9 +1502,7 @@ FE_Q::initialize_embedding () 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 Point p_subcell = this->unit_support_points[j]; const double subcell_value = this->poly_space.compute_value (i, p_subcell); @@ -1629,9 +1543,7 @@ FE_Q::initialize_embedding () // generate a point on the // child cell and evaluate the // shape functions there - const Point p_subcell - = FE_Q_Helper::generate_unit_point (index_map[j], this->dofs_per_cell, - dealii::internal::int2type()); + const Point p_subcell = this->unit_support_points[j]; const Point p_cell = GeometryInfo::child_to_cell_coordinates (p_subcell, child, RefinementCase(ref+1));