* @p FiniteElementData.
*/
static std::vector<unsigned int> 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:
- *
- * 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.
- *
- * 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.
- *
- * 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.
- *
- * 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.
- *
- * This function is called from
- * the constructor.
- */
- static
- std::vector<unsigned int>
- lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe_data,
- const unsigned int degree);
/**
* This is an analogon to the
- * previous function, but working
- * on faces. Called from the
+ * FETools::lexicographic_to_hierarchic_numbering
+ * function, but working on
+ * faces. Called from the
* constructor.
*/
static
//---------------------------------------------------------------------------
#include <fe/fe_q.h>
+#include <fe/fe_tools.h>
#ifdef HAVE_STD_STRINGSTREAM
# include <sstream>
get_dpo_vector(degree),1, degree).dofs_per_cell, std::vector<bool>(1,true))),
face_index_map(FE_Q_Helper::invert_numbering(face_lexicographic_to_hierarchic_numbering (degree)))
{
- this->poly_space.set_numbering(FE_Q_Helper::invert_numbering(
- lexicographic_to_hierarchic_numbering (*this, degree)));
-
+ std::vector<unsigned int> renumber (this->dofs_per_cell);
+ FETools::lexicographic_to_hierarchic_numbering (*this, renumber);
+ this->poly_space.set_numbering(FE_Q_Helper::invert_numbering(renumber));
+
// finally fill in support points
// on cell and face
initialize_unit_support_points ();
}
-
-template <int dim>
-std::vector<unsigned int>
-FE_Q<dim>::lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe_data,
- const unsigned int degree)
-{
- std::vector<unsigned int> renumber (fe_data.dofs_per_cell);
-
- const unsigned int n = degree+1;
-
- if (degree == 0)
- {
- 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;
- };
-
- if (degree > 0)
- {
- Assert (fe_data.dofs_per_vertex == 1, ExcInternalError());
- for (unsigned int i=0; i<GeometryInfo<dim>::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, degree };
- index = values[i];
- break;
- };
-
- case 2:
- {
- const unsigned int values[GeometryInfo<2>::vertices_per_cell]
- = { 0, degree, n*degree+degree, n*degree };
- index = values[i];
- break;
- };
-
- case 3:
- {
- const unsigned int values[GeometryInfo<3>::vertices_per_cell]
- = { 0, degree,
- n*n*degree + degree, n*n*degree,
- n*degree, n*degree+degree,
- n*n*degree + n*degree+degree, n*n*degree + n*degree};
- index = values[i];
- break;
- };
-
- default:
- Assert(false, ExcNotImplemented());
- }
-
- Assert (index<renumber.size(), ExcInternalError());
- renumber[index] = i;
- }
- };
-
- // for degree 2 and higher: Lines,
- // quads, hexes etc also carry
- // degrees of freedom
- if (degree > 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 (int i=0; i<static_cast<signed int>(GeometryInfo<dim>::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 = degree;
- break;
- // y=1 x=z=0
- case 202:
- case 304: case 307:
- tensorstart = n*degree;
- break;
- // x=z=1 y=0
- case 310:
- tensorstart = n*n*degree+degree;
- break;
- // z=1 x=y=0
- case 302: case 311:
- tensorstart = n*n*degree;
- break;
- // x=y=1 z=0
- case 305:
- tensorstart = n*degree+degree;
- break;
- // y=z=1 x=0
- case 306:
- tensorstart = n*n*n-n;
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
-
- for (unsigned int jx = 1; jx<degree ;++jx)
- {
- unsigned int tensorindex = tensorstart + jx * incr;
- Assert (tensorindex<renumber.size(), ExcInternalError());
- renumber[tensorindex] = index++;
- }
- }
-
- for (int i=0; i<static_cast<signed int>(GeometryInfo<dim>::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)
- {
- case 0:
- tensorstart = 0; incx = 1;
- if (dim==2)
- incy = n;
- else
- incy = n*n;
- break;
- case 1:
- tensorstart = n*degree; incx = 1; incy = n*n;
- break;
- case 2:
- tensorstart = 0; incx = 1; incy = n;
- break;
- case 3:
- tensorstart = degree; incx = n; incy = n*n;
- break;
- case 4:
- tensorstart = n*n*degree; incx = 1; incy = n;
- break;
- case 5:
- tensorstart = 0; incx = n; incy = n*n;
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
-
- for (unsigned int jy = 1; jy<degree; jy++)
- for (unsigned int jx = 1; jx<degree ;++jx)
- {
- unsigned int tensorindex = tensorstart
- + jx * incx + jy * incy;
- Assert (tensorindex<renumber.size(), ExcInternalError());
- renumber[tensorindex] = index++;
- }
- }
-
- if (GeometryInfo<dim>::hexes_per_cell > 0)
- for (int i=0; i<static_cast<signed int>(GeometryInfo<dim>::hexes_per_cell); ++i)
- {
- unsigned int index = fe_data.first_hex_index;
-
- for (unsigned int jz = 1; jz<degree; jz++)
- for (unsigned int jy = 1; jy<degree; jy++)
- for (unsigned int jx = 1; jx<degree; jx++)
- {
- const unsigned int tensorindex = jx + jy*n + jz*n*n;
- Assert (tensorindex<renumber.size(), ExcInternalError());
- renumber[tensorindex]=index++;
- }
- }
- }
-
- return renumber;
-}
-
-
-
template <int dim>
std::vector<unsigned int>
FE_Q<dim>::face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
{
- const FiniteElementData<dim-1> fe_data(FE_Q<dim-1>::get_dpo_vector(degree),1);
- return FE_Q<dim-1>::lexicographic_to_hierarchic_numbering (fe_data, degree);
+ const FiniteElementData<dim-1> face_data(FE_Q<dim-1>::get_dpo_vector(degree),1);
+ std::vector<unsigned int> face_renumber (face_data.dofs_per_cell);
+ FETools::lexicographic_to_hierarchic_numbering (face_data, face_renumber);
+ return face_renumber;
}