From 368199ac285be876c80e583128676277a4c8dc3f Mon Sep 17 00:00:00 2001 From: hartmann Date: Fri, 5 Aug 2005 08:47:47 +0000 Subject: [PATCH] Remove lexicographic_to_hierarchic_numbering. Implement hierarchic_to_fe_q_hierarchical_numbering analogous to FETools::hierarchic_to_lexicographic_numbering. Rename face_lexicographic_to_hierarchic_numbering to face_fe_q_hierarchical_to_hierarchic_numbering. git-svn-id: https://svn.dealii.org/trunk@11227 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/fe/fe_q_hierarchical.h | 120 +++--- .../deal.II/source/fe/fe_q_hierarchical.cc | 402 ++++++++---------- 2 files changed, 222 insertions(+), 300 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_q_hierarchical.h b/deal.II/deal.II/include/fe/fe_q_hierarchical.h index 040bb30fe5..62dea1042d 100644 --- a/deal.II/deal.II/include/fe/fe_q_hierarchical.h +++ b/deal.II/deal.II/include/fe/fe_q_hierarchical.h @@ -228,7 +228,7 @@ template class MappingQ; * Note the reverse ordering of degrees of freedom on the left and upper * line. * - * @author Brian Carnes, 2002, Ralf Hartmann 2004 + * @author Brian Carnes, 2002, Ralf Hartmann 2004, 2005 */ template class FE_Q_Hierarchical : public FE_Poly,dim> @@ -315,75 +315,63 @@ class FE_Q_Hierarchical : public FE_Poly,dim> static std::vector get_dpo_vector(const unsigned int degree); /** - * Map tensor product data to - * shape function numbering. This - * function is actually an alike - * replica of the respective - * function in the FETools - * class, but is kept for three - * reasons: + * The numbering of the degrees + * of freedom in continous finite + * elements is hierarchic, + * i.e. in such a way that we + * first number the vertex dofs, + * in the order of the vertices + * as defined by the + * triangulation, then the line + * dofs in the order and + * respecting the direction of + * the lines, then the dofs on + * quads, etc. * - * 1. It only operates on a - * FiniteElementData - * structure. This is ok in the - * present context, since we can - * control which types of - * arguments it is called with - * because this is a private - * function. However, the - * publicly visible function in - * the FETools class needs - * to make sure that the - * FiniteElementData object - * it works on actually - * represents a continuous finite - * element, which we found too - * difficult if we do not pass an - * object of type FE_Q() - * directly. + * The dofs associated with 1d + * hierarchical polynomials are + * ordered with the vertices + * first ($phi_0(x)=1-x$ and + * $phi_1(x)=x$) and then the + * line dofs (the higher degree + * polynomials). The 2d and 3d + * hierarchical polynomials + * originate from the 1d + * hierarchical polynomials by + * tensor product. In the + * following, the resulting + * numbering of dofs will be + * denoted by `fe_q_hierarchical + * numbering`. * - * 2. If we would call the - * publicly available version of - * this function instead of this - * one, we would have to pass a - * finite element - * object. However, since the - * construction of an entire - * finite element object can be - * costly, we rather chose to - * retain this function. + * This function constructs a + * table which fe_q_hierarchical + * index each degree of freedom + * in the hierarchic numbering + * would have. * - * 3. Third reason is that we - * want to call this function for - * faces as well, by just calling - * this function for the finite - * element of one dimension - * less. If we would call the - * global function instead, this - * would require us to construct - * a second finite element object - * of one dimension less, just to - * call this function. Since that - * function does not make use of - * hanging nodes constraints, - * interpolation and restriction - * matrices, etc, this would have - * been a waste. Furthermore, it - * would have posed problems with - * template instantiations. + * This function is anologous to + * the + * FETools::hierarchical_to_lexicographic_numbering() + * function. However, in contrast + * to the fe_q_hierarchical + * numbering defined above, the + * lexicographic numbering + * originates from the tensor + * products of consecutive + * numbered dofs (like for + * LagrangeEquidistant). * - * To sum up, the existence of - * this function is a compromise - * between simplicity and proper - * library design, where we have - * chosen to weigh the simplicity - * aspect a little more than - * proper design. + * It is assumed that the size of + * the output argument already + * matches the correct size, + * which is equal to the number + * of degrees of freedom in the + * finite element. */ static - std::vector - lexicographic_to_hierarchic_numbering (const FiniteElementData &fe_data, - const unsigned int degree); + std::vector hierarchic_to_fe_q_hierarchical_numbering ( + const FiniteElementData &fe); /** * This is an analogon to the @@ -392,7 +380,7 @@ class FE_Q_Hierarchical : public FE_Poly,dim> */ static std::vector - face_lexicographic_to_hierarchic_numbering (const unsigned int degree); + face_fe_q_hierarchical_to_hierarchic_numbering (const unsigned int degree); /** * Initialize two auxiliary @@ -463,6 +451,6 @@ void FE_Q_Hierarchical<1>::initialize_unit_face_support_points (); template <> std::vector -FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int); +FE_Q_Hierarchical<1>::face_fe_q_hierarchical_to_hierarchic_numbering (const unsigned int); #endif diff --git a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc index bb6c9145e6..6bc2343da6 100644 --- a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc +++ b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc @@ -46,10 +46,10 @@ FE_Q_Hierarchical::FE_Q_Hierarchical (const unsigned int degree) get_dpo_vector(degree),1, degree).dofs_per_cell, false), std::vector >(FiniteElementData( get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector(1,true))), - face_renumber(face_lexicographic_to_hierarchic_numbering (degree)) + face_renumber(face_fe_q_hierarchical_to_hierarchic_numbering (degree)) { - this->poly_space.set_numbering(invert_numbering( - lexicographic_to_hierarchic_numbering (*this, degree))); + this->poly_space.set_numbering( + hierarchic_to_fe_q_hierarchical_numbering(*this)); // The matrix @p{dofs_cell} contains the // values of the linear functionals of @@ -339,7 +339,6 @@ initialize_embedding_and_restriction (const std::vector > &do const std::vector &renumber= this->poly_space.get_numbering(); - for (unsigned int c=0; c::children_per_cell; ++c) { this->prolongation[c].reinit (this->dofs_per_cell, this->dofs_per_cell); @@ -575,251 +574,186 @@ FE_Q_Hierarchical::get_dpo_vector(const unsigned int deg) template std::vector FE_Q_Hierarchical:: -lexicographic_to_hierarchic_numbering (const FiniteElementData &fe_data, - const unsigned int degree) +hierarchic_to_fe_q_hierarchical_numbering (const FiniteElementData &fe) { - std::vector renumber (fe_data.dofs_per_cell, 0); - - const unsigned int n = degree+1; + Assert (fe.n_components() == 1, ExcInternalError()); + std::vector h2l(fe.dofs_per_cell); + // polynomial degree + const unsigned int degree = fe.dofs_per_line+1; + // number of grid points in each + // direction + const unsigned int n = degree+1; - if (degree == 0) + // 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) { - Assert ( (fe_data.dofs_per_vertex == 0) && - ((fe_data.dofs_per_line == 0) || (dim == 1)) && - ((fe_data.dofs_per_quad == 0) || (dim == 2)) && - ((fe_data.dofs_per_hex == 0) || (dim == 3)), - ExcInternalError() ); - renumber[0] = 0; - }; + case 1: + { + for (unsigned int i=0; i 0) - { - Assert (fe_data.dofs_per_vertex == 1, ExcInternalError()); - for (unsigned int i=0; i::vertices_per_cell; ++i) - { - unsigned int index = 0; - // Find indices of vertices. - // Unfortunately, somebody - // switched the upper corner - // points of a quad. The same - // person decided to find a very - // creative numbering of the - // vertices of a hexahedron. - // Therefore, this looks quite - // sophisticated. - // - // NB: This same person - // claims to have had good - // reasons then, but seems to - // have forgotten about - // them. At least, the - // numbering was discussed - // with the complaining - // person back then when all - // began :-) - switch (dim) - { - case 1: - { - const unsigned int values[GeometryInfo<1>::vertices_per_cell] - = { 0, 1 }; - index = values[i]; - break; - }; - - case 2: - { - const unsigned int values[GeometryInfo<2>::vertices_per_cell] - = { 0, 1, n + 1, n }; - index = values[i]; - break; - }; - - case 3: - { - const unsigned int values[GeometryInfo<3>::vertices_per_cell] - = { 0, 1, - n * n + 1, n * n, - n, n + 1, - n * n + n + 1, n * n + n}; - index = values[i]; - break; - }; + break; + } + + case 2: + { + // Example: degree=3 + // + // hierarchical numbering: + // 3 8 9 2 + // 11 14 15 7 + // 10 12 13 6 + // 0 4 5 1 + // + // fe_q_hierarchical numbering: + // 4 6 7 5 + // 12 14 15 13 + // 8 10 11 9 + // 0 2 3 1 + unsigned int next_index = 0; + // first the four vertices + h2l[next_index++] = 0; + h2l[next_index++] = 1; + h2l[next_index++] = n+1; + h2l[next_index++] = n; + // bottom line + for (unsigned int i=0; i 1) - { - Assert (fe_data.dofs_per_line == degree-1, ExcInternalError()); - Assert ((fe_data.dofs_per_quad == (degree-1)*(degree-1)) || - (dim < 2), ExcInternalError()); - Assert ((fe_data.dofs_per_hex == (degree-1)*(degree-1)*(degree-1)) || - (dim < 3), ExcInternalError()); - - for (unsigned int i=0; i::lines_per_cell; ++i) - { - unsigned int index = fe_data.first_line_index - + i*fe_data.dofs_per_line; - unsigned int incr = 0; - unsigned int tensorstart = 0; - // This again looks quite - // strange because of the odd - // numbering scheme. - switch (i+100*dim) - { - // lines in x-direction - case 100: - case 200: case 202: - case 300: case 302: case 304: case 306: - incr = 1; - break; - // lines in y-direction - case 201: case 203: - case 308: case 309: case 310: case 311: - incr = n; - break; - // lines in z-direction - case 301: case 303: case 305: case 307: - incr = n * n; - break; - default: - Assert(false, ExcNotImplemented()); - } - switch (i+100*dim) - { - // x=y=z=0 - case 100: - case 200: case 203: - case 300: case 303: case 308: - tensorstart = 0; - break; - // x=1 y=z=0 - case 201: - case 301: case 309: - tensorstart = 1; - break; - // y=1 x=z=0 - case 202: - case 304: case 307: - tensorstart = n; - break; - // x=z=1 y=0 - case 310: - tensorstart = n * n + 1; - break; - // z=1 x=y=0 - case 302: case 311: - tensorstart = n * n; - break; - // x=y=1 z=0 - case 305: - tensorstart = n + 1; - break; - // y=z=1 x=0 - case 306: - tensorstart = n * n + n; - break; - default: - Assert(false, ExcNotImplemented()); - } - - for (unsigned int jx = 2; jx<=degree ;++jx) - { - unsigned int tensorindex = tensorstart + jx * incr; - Assert (tensorindex < renumber.size(), ExcInternalError()); - renumber[tensorindex] = index++; - } - } - - for (int i=0; i(GeometryInfo::quads_per_cell); ++i) - { - unsigned int index = fe_data.first_quad_index+i*fe_data.dofs_per_quad; - unsigned int tensorstart = 0; - unsigned int incx = 0; - unsigned int incy = 0; - switch (i) - { - // z=0 (dim==2), y=0 (dim==3) - case 0: - tensorstart = 0; incx = 1; - if (dim==2) - incy = n; - else - incy = n * n; - break; - // y=1 - case 1: - tensorstart = n; incx = 1; incy = n * n; - break; - // z=0 - case 2: - tensorstart = 0; incx = 1; incy = n; - break; - // x=1 - case 3: - tensorstart = 1; incx = n; incy = n * n; - break; - // z=1 - case 4: - tensorstart = n * n; incx = 1; incy = n; - break; - // x=0 - case 5: - tensorstart = 0; incx = n; incy = n * n; - break; - default: - Assert(false, ExcNotImplemented()); - } - - for (unsigned int jy = 2; jy<=degree; jy++) - for (unsigned int jx = 2; jx<=degree ;++jx) - { - unsigned int tensorindex = tensorstart - + jx * incx + jy * incy; - Assert (tensorindex < renumber.size(), ExcInternalError()); - renumber[tensorindex] = index++; - } - } - - if (GeometryInfo::hexes_per_cell > 0) - for (int i=0; i(GeometryInfo::hexes_per_cell); ++i) - { - unsigned int index = fe_data.first_hex_index; - - for (unsigned int jz = 2; jz<=degree; jz++) - for (unsigned int jy = 2; jy<=degree; jy++) - for (unsigned int jx = 2; jx<=degree; jx++) - { - const unsigned int tensorindex = jx + jy * n + jz * n * n; - Assert (tensorindex < renumber.size(), ExcInternalError()); - renumber[tensorindex]=index++; - } - } + default: + Assert (false, ExcNotImplemented()); } - - return renumber; + return h2l; } - template std::vector FE_Q_Hierarchical:: -face_lexicographic_to_hierarchic_numbering (const unsigned int degree) +face_fe_q_hierarchical_to_hierarchic_numbering (const unsigned int degree) { FiniteElementData fe_data(FE_Q_Hierarchical::get_dpo_vector(degree),1); - return FE_Q_Hierarchical::lexicographic_to_hierarchic_numbering (fe_data, - degree); + return invert_numbering(FE_Q_Hierarchical:: + hierarchic_to_fe_q_hierarchical_numbering (fe_data)); } @@ -827,7 +761,7 @@ face_lexicographic_to_hierarchic_numbering (const unsigned int degree) template <> std::vector -FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int) +FE_Q_Hierarchical<1>::face_fe_q_hierarchical_to_hierarchic_numbering (const unsigned int) { return std::vector (); } -- 2.39.5