]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
To avoid code duplication, use FETools::lexicographic_to_hierarchic_numbering functio...
authorhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 13 Jun 2005 11:49:17 +0000 (11:49 +0000)
committerhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 13 Jun 2005 11:49:17 +0000 (11:49 +0000)
git-svn-id: https://svn.dealii.org/trunk@10858 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_q.h
deal.II/deal.II/source/fe/fe_q.cc

index a18e4cabf48621f7798fc40824adc718860991ac..8ca83af7810c0cad38ee6641a1a655eecef92003 100644 (file)
@@ -328,85 +328,12 @@ class FE_Q : public FE_Poly<TensorProductPolynomials<dim>,dim>
                                      * @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
index 22e879c1ba5aee515a1deea38450b8bfea1c2be0..e02fc15f0edfd475114ff02c6151d794f71cddfd 100644 (file)
@@ -12,6 +12,7 @@
 //---------------------------------------------------------------------------
 
 #include <fe/fe_q.h>
+#include <fe/fe_tools.h>
 
 #ifdef HAVE_STD_STRINGSTREAM
 #  include <sstream>
@@ -431,9 +432,10 @@ FE_Q<dim>::FE_Q (const unsigned int degree)
                    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 ();
@@ -665,245 +667,14 @@ FE_Q<dim>::get_dpo_vector(const unsigned int deg)
 }
 
 
-
-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;
 }
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.